Introduction

The skyrmion is a topological configuration of a continuous field. Although it was originally proposed to explain hadrons in the particle theory1,2, it has turned out to be realized in various forms in condensed matter physics3. One possible realization was discovered in magnets, in the form of skyrmion-like magnetic textures4,5,6,7,8. The magnetic skyrmions often exist in their stable form, the so-called skyrmion crystal (SkX), which is a periodic array of the particle-like magnetic skyrmions. Importantly, such a SkX is approximately expressed as a superposition of multiple helical spin density waves, and hence, it can be regarded as an interference pattern by the multiple helices. It has attracted enormous attention since the swirling magnetic texture generates an emergent magnetic field through the Berry phase mechanism and results in peculiar transport phenomena, such as the topological Hall effect3,9,10.

Similar to an isolated skyrmion, the SkX is characterized by three quantities: skyrmion number, vorticity, and helicity3. However, as the SkX is regarded as an interference pattern, it has another degree of freedom, which has been overlooked in the previous researches, the phases of the constituent waves. This is exemplified for three scalar waves in Fig. 1a, b, where a phase shift in one of the three waves leads to a different interference pattern with different symmetry. Such a phase degree of freedom exists in all the interference phenomena, except for linearly independent waves in continuous space for which a phase shift is equivalent to a spatial translation. The SkX appears not in a continuous field but for spins on a discrete lattice, which leads to a further variety of the interference patterns by the discretization, even for the linearly independent waves. A shift of the relative phases changes not only magnetic textures but also emergent magnetic fields, and hence, transport properties, but such an interesting possibility has not been elucidated thus far.

Fig. 1: Phase shift and interference patterns.
figure 1

a, b Superpositions of three density waves with the wave vectors Q1, Q2, and Q3. b is the figure generated from a with the phase shift of π/2 in Q1, which breaks sixfold rotational symmetry in a. c, d superpositions of three spirals waves: the nsk = 1 skyrmion crystal (SkX1) (c) and the meron-antimeron crystal (MAX) (d). e, f Superpositions of three sinusoidal waves: the nsk = 2 skyrmion crystal (SkX2) (e) and the tetra-axial vortex crystal (TVX) (f). d, f are generated from (c) and (e), respectively, with the phase shift of Θ = π/2, and both of them break sixfold rotational symmetry similar to b. In c–f, the upper and lower planes show the spin \({{{{{{{{\bf{S}}}}}}}}}_{i}=({S}_{i}^{x},{S}_{i}^{y},{S}_{i}^{z})\) [the color scale indicates \({S}_{i}^{z}\), and the arrows indicate \(({S}_{i}^{x},{S}_{i}^{y})\)] and the scalar chirality χR, respectively. The spin textures in the magnetic unit cell are shown in the insets of c–f.

Full size image

In this study, we theoretically unveil the effect of phase shifts in the SkX and propose how to control the phase degree of freedom. Considering an itinerant electron model on a triangular lattice, we show that the SkX turns into a tetra-axial vortex crystal (TVX) or a meron-antimeron crystal (MAX) by a phase shift of π/2. The phase-shifted states have distinct properties from the SkX: The SkX exhibits a net scalar chirality leading to the topological Hall effect, while the TVX and MAX exhibit a staggered one that does not lead to the topological Hall effect, but induces nonreciprocal transport phenomena that do not require the spin-orbit coupling. We find that such a phase shift can be caused by several different mechanisms, such as exchange interactions between the localized spins, thermal fluctuations, and long-range chirality interactions. Our results open another route to a further variety of magnetic textures which have been overlooked in skyrmion-hosting materials.

Results

Let us start by classifying noncoplanar spin textures according to the type of constituent waves and the relative phases. First, we consider a superposition of spiral spin textures represented by \({{{{{{{{\bf{S}}}}}}}}}_{i}^{{{{{{{{\rm{spiral}}}}}}}}}=\mathop{\sum }\nolimits_{\nu = 1}^{3}\left(\sin {{{{{{{{\mathcal{Q}}}}}}}}}_{\nu }\cos {\phi }_{\nu },\sin {{{{{{{{\mathcal{Q}}}}}}}}}_{\nu }\sin {\phi }_{\nu },-\cos {{{{{{{{\mathcal{Q}}}}}}}}}_{\nu }\right)\), where \({{{{{{{{\mathcal{Q}}}}}}}}}_{\nu }={{{{{{{{\bf{Q}}}}}}}}}_{\nu }\cdot {{{{{{{{\bf{r}}}}}}}}}_{i}+{\theta }_{\nu }\) and \({\phi }_{\nu }=\frac{2}{3}\pi (\nu -1)\); Qν and θν are the wave vector and the phase of the νth spiral, respectively, and ri is the position vector for site i. In the following analyses, as an archetype, we consider a two-dimensional triangular lattice system with threefold rotationally symmetric wave vectors with spiral pitch Q: Q1 = (Q, 0), \({{{{{{{{\bf{Q}}}}}}}}}_{2}=(-Q/2,\sqrt{3}Q/2)\), and \({{{{{{{{\bf{Q}}}}}}}}}_{3}=(-Q/2,-\sqrt{3}Q/2)\) satisfying ∑νQν = 0. The spin texture \({{{{{{{{\bf{S}}}}}}}}}_{i}^{{{{{{{{\rm{spiral}}}}}}}}}\) is modulated by shifting the phases θν. We demonstrate the situation for θ1 = θ2 = θ3 by changing the total phase Θ = ∑νθν. Figure 1c, d displays the spin textures with Θ = 0 and π/2. In each figure, the upper and lower planes show the textures of spin and scalar spin chirality, respectively; the latter is defined by χR = Si ⋅ (Sj × Sk), where R represents the position vector at the center of a triangle with sites i, j, k in the counterclockwise order. The case with Θ = 0 in Fig. 1c is a periodic array of skyrmions with the skyrmion number of one (nsk = 1), which we call the SkX1. It retains sixfold rotational symmetry in both spin and chirality, and has a nonzero net chirality leading to the topological Hall effect. Meanwhile, the case with Θ = π/2 in Fig. 1d is a staggered arrangement of merons and antimerions (half skyrmions and antiskyrmions11), which we call the MAX. In this state, the rotational symmetry is reduced to threefold. Moreover, the meron and antimeron carry the skyrmion number of  +1/2 and −1/2, respectively, and hence, the total skyrmion number is zero in the MAX; accordingly, the net value of χR, \({\chi }^{\rm total}=\frac{1}{N}{\sum }_{{{{{{{{\bf{R}}}}}}}}}{\chi }_{{{{{{{{\bf{R}}}}}}}}}\) where N is the number of lattice sites, also vanishes and the MAX does not show the topological Hall effect. We note that the MAX was proposed as a candidate for the unidentified magnetic state next to the SkX1 found in a triangular-lattice magnet Gd2PdSi311.

Next, we consider a superposition of sinusoidal waves12,13,14, which is represented by \({{{{{{{{\bf{S}}}}}}}}}_{i}^{\sin }=\left(\cos {{{{{{{{\mathcal{Q}}}}}}}}}_{1},\cos {{{{{{{{\mathcal{Q}}}}}}}}}_{2},\cos {{{{{{{{\mathcal{Q}}}}}}}}}_{3}\right)\). Similar to \({{{{{{{{\bf{S}}}}}}}}}_{i}^{{{{{{{{\rm{spiral}}}}}}}}}\), different Θ gives different spin and chirality textures, as shown in Fig. 1e, f (the spin frame is rotated for better visibility). The spin texture with Θ = 0 in Fig. 1e is the other SkX called the SkX2, in which each skyrmion has the skyrmion number of two (nsk = 2). In this state, while the spin texture has threefold rotational symmetry, the chirality χR is sixfold and the net value χtotal is nonzero, similar to the SkX1 in Fig. 1c. The phase shift by π/2 lowers the symmetry from sixfold to threefold, as shown in Fig. 1f; χR has a staggered configuration with no net scalar chirality, similar to the MAX in Fig. 1d. In this state, the spin texture is given by a periodic array of four types of vortices; the vortex axes, which are defined by the vorticity for xy, yz, and zx components of spins point to four corners of the tetrahedron (see Supplementary Information). Hence, we call the Θ = π/2 state the TVX.

Thus, in both cases, the phase shift changes not only the spin texture but also the symmetry and topology. In particular, the net value of the scalar chirality χtotal is sensitively dependent on Θ; the two types of SkXs at Θ = 0 have nonzero values and cause the topological Hall effect, while the MAX and TVX at Θ = π/2 have no net value and do not show the topological Hall effect. Interestingly, however, the breaking of sixfold rotational symmetry in χR in the MAX and TVX leads to Fermi surface deformations as discussed below, which can induce direction-dependent nonreciprocal transport phenomena without the spin-orbit coupling.

The optimal values of the phases θν will be determined by multiple factors, such as lattice geometry and interactions between the spins. In the previous studies, the SkXs with Θ = 0 are stabilized, e.g., by the Dzyaloshinskii−Moriya (DM)6,15, four-spin16,17,18,19, frustrated20,21,22, and spin-charge interactions14,23,24,25 on various lattices. The key question addressed here is what is the relevant parameter to cause a phase shift that leads to switching of magnetic, topological, and transport properties. In the following, we unveil three different mechanisms for such a phase shift, by taking an archetypal model for itinerant magnets hosting SkXs, the Kondo lattice model on a triangular lattice where both the SkX1 and SkX2 appear in the ground state (see “Methods”).

We first demonstrate that a phase shift can be caused by introducing the exchange interactions between the localized spins described by \({{{{{{{{\mathcal{H}}}}}}}}}^{{{{{{{{\rm{loc}}}}}}}}}={\sum }_{ij}{J}_{ij}{{{{{{{{\bf{S}}}}}}}}}_{i}\cdot {{{{{{{{\bf{S}}}}}}}}}_{j}\) to the original Kondo lattice Hamiltonian \({{{{{{{\mathcal{H}}}}}}}}\). Considering the first-, second-, and third-neighbor interactions, J1, J2, and J3 for Jij, respectively, we perform a variational calculation to determine the ground-state phase diagram of the Hamiltonian \({{{{{{{\mathcal{H}}}}}}}}+{{{{{{{{\mathcal{H}}}}}}}}}^{{{{{{{{\rm{loc}}}}}}}}}\) at zero field (see “Methods”). Figure 2a, b shows the results on the J1−J2 and J1−J3 planes. While the SkX2 is stable in the wide range of parameters, we find two topological phase transitions: One is to the TVX while increasing J1 and decreasing J2 (increasing J3) and the other is to the SkX1 while decreasing J1 and decreasing J2 (increasing J3), as shown in Fig. 2a (Fig. 2b). The former transition from the SkX2 to the TVX is accompanied by the phase shift with Θ = π/2.

Fig. 2: Variational ground states of the Kondo lattice model with the exchange interactions between the localized spins up to third neighbors, J1, J2, and J3.
figure 2

Variational phase diagrams on the J1−J2 (a) and J1−J3 (b) planes. c Spin (left) and chirality (right) structure factors for the SkX2 (top), TVX (middle), and SkX1 (bottom). The circles represent the peak positions at Q1, Q2, and Q3, and the arrows indicate the subdominant (dominant) peak positions in the spin (chirality).

Full size image

The phase transitions among these three states are understood from the higher harmonics in the spin structure. We show the spin structure factors in momentum (q) space for the SkX2, TVX, and SkX1 in the left panels of Fig. 2c (see “Methods”). Although the dominant peaks appear at Q1, Q2, and Q3 commonly in the three phases, subdominant peaks are found at different q among them: \({{{{{{{\bf{q}}}}}}}}={{{{{{{{\bf{Q}}}}}}}}}_{\nu }-{{{{{{{{\bf{Q}}}}}}}}}_{\nu ^{\prime} }\) (\(\nu \,\ne\, \nu ^{\prime}\)) in the SkX2, q = 3Qν in the TVX, and q = 0 in the SkX1. Thus, considering that the Fourier transform of \({{{{{{{{\mathcal{H}}}}}}}}}^{{{{{{{{\rm{loc}}}}}}}}}\) is written as ∑qJqSq ⋅ S−q, the (anti)ferromagnetic interactions giving J0 < 0 (\({J}_{3{{{{{{{{\bf{Q}}}}}}}}}_{\nu }} < \, 0\)) tend to prefer the SkX1 (TVX).

It is noteworthy that the different superpositions of the constituent waves also give rise to the different peak positions in the q-resolved scalar chirality shown in the right panels of Fig. 2c (see “Methods”). As mentioned above, both the SkX2 and SkX1 have the dominant peak at q = 0 reflecting nonzero χtotal, while the TVX has the dominant peaks at 2Qν with equal intensities and no weight at q = 0.

The second mechanism to cause the phase shift is thermal fluctuations. We study the finite-temperature behavior of the Kondo lattice Hamiltonian \({{{{{{{\mathcal{H}}}}}}}}\) by performing the Langevin dynamics simulations with the kernel polynomial method (see “Methods”)26,27,28. Figure 3a shows the results at zero magnetic field H = 0 where the ground state is the SkX214. We find two phase transitions at T1 ≃ 0.0055 and T2 ≃ 0.009. The transition at T1 is characterized by the onset of χ0, suggesting that the SkX2 remains stable up to T1. Note that the true magnetic long-range order is limited to zero temperature in the present two-dimensional system due to the Mermin−Wagner theorem29; the state for 0 < T < T1 is a chiral spin liquid having the SkX2 spin texture with a finite correlation length (but much longer than the system size). The real-space spin configuration at T = 0 is shown in Fig. 3c, which corresponds to that in Fig. 1e. Meanwhile, the transition at T2 appears to be signaled by the onset of the higher harmonics \({\chi }_{2{{{{{{{{\bf{Q}}}}}}}}}_{\nu }}\) as plotted in Fig. 3a. A snapshot of the real-space spin configuration at T = 0.006 is shown in Fig. 3d, which well reproduces the spin texture for the TVX in Fig. 1f. From these results, we conclude that the low- and intermediate-temperature phases are chiral spin liquids with the SkX2 and TVX spin textures, respectively, and the transition at T1 is associated with the phase shift of π/2 between them.

Fig. 3: Finite-temperature behaviors of the Kondo lattice model.
figure 3

a,b Temperature dependences of the scalar chirality, \({\chi }_{0}^{2}\) and \({\chi }_{2{{{{{{{{\bf{Q}}}}}}}}}_{\nu }}^{2}\) at zero magnetic field H = 0 (a) and H = 0.004 (b). The dashed (dash-dotted) line represents the transition temperature T1 (T2) in a and \({T}_{1}^{\prime}\) (\({T}_{2}^{\prime}\)) in b. c–e Real-space spin configurations in the SkX2 at T = 0 and H = 0 (c), the TVX at T = 0.006 and H = 0 (d), and the SkX1 at T = 0 and H = 0.004 (e). The color scale indicates \({S}_{i}^{z}\), while the arrows indicate \(({S}_{i}^{x},{S}_{i}^{y})\). The insets display the Fermi surfaces in each state; the sixfold rotational symmetry is broken in all the three states. See the main text.

Full size image

The appearance of the TVX at finite temperature is explained by an effective chirality interaction as follows. At the mean-field level, the entropic contributions are in general given in the form of n th-order magnetic interactions as \(T{\sum }_{{{{{{{{{\bf{q}}}}}}}}}_{1}\cdots {{{{{{{{\bf{q}}}}}}}}}_{n}}({{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{q}}}}}}}}}_{1}}\cdot {{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{q}}}}}}}}}_{2}})\cdots ({{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{q}}}}}}}}}_{n-1}}\cdot {{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{q}}}}}}}}}_{n}})\delta ({{{{{{{{\bf{q}}}}}}}}}_{1}+\cdots +{{{{{{{{\bf{q}}}}}}}}}_{n})\)30,31. Among them, the lowest-order contribution to the phase shift appears in the sixth order. By considering \({{{{{{{{\bf{S}}}}}}}}}_{i}^{\sin }\), the dominant entropic contribution is given as \(T{{{{{{{\rm{Re}}}}}}}}[({{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{1}}\cdot {{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{1}})({{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{2}}\cdot {{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{2}})({{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{3}}\cdot {{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{3}})]=T{{{{{{{\rm{Re}}}}}}}}[{\{{{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{1}}\cdot ({{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{2}}\times {{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{3}})\}}^{2}]\propto T\cos 2{{\Theta }}\). Thus, the sixth-order entropic term, which has the form of chirality-chirality interactions, depends on Θ and tends to stabilize the TVX at finite temperature.

Next, let us discuss the case of SkX1. Figure 3b represents the results under a magnetic field, where the ground state is the SkX1. Similar to the zero-field case in Fig. 3a, two phase transitions occur at \({T}_{1}^{\prime}\simeq 0.0045\) and \({T}_{2}^{\prime}\simeq 0.0085\). The low-temperature state below \({T}_{1}^{\prime}\) is the SkX1 (with quasi-long-range order in the xy components), whose spin configuration is shown in Fig. 3e14. Recalling the phase shift of π/2 from the SkX2 to the TVX in the zero-field case, one may expect that the SkX1 changes into the MAX by raising temperature, but we find that the intermediate phase for \({T}_{1}^{\prime} \, < \, T \, < \, {T}_{2}^{\prime}\) is a different superposition of three sinusoidal waves which has χtotal = 0. This is because the Zeeman energy gain in the MAX is not sufficient to overcome the obtained state. We note that the threefold rotational symmetry is broken in the intermediate phase; hence, we call this phase the anisotropic 3Q state (see Supplementary Information). This phase transition from the SkX1 to the anisotropic 3Q state is also accounted for by the sixth-order entropic contribution, similar to that from the SkX2 to the TVX at zero field, as will be shown below.

We display the Fermi surfaces in the SkX2, TVX, and SkX1 in the inset of Fig. 3c–e, respectively. The Fermi surface in the TVX in Fig. 3d is threefold rotationally symmetric, meaning the breaking of the sixfold rotational symmetry of the triangular lattice, as expected from the above discussion. This leads to a nonreciprocal transport in itinerant electrons. Notably, there appear threefold rotationally symmetric Fermi surfaces even in the SkX2 in Fig. 3c (very weakly broken) and SkX1 in Fig. 3e. This is not due to the shift of Θ but by phase-locking at (θ1, θ2, θ3) = (π/3, −π/3, 0) (any permutation is allowed) so that the skyrmion cores avoid the lattice sites (the values of θν depend on Q; see Supplementary Information). Thus, the results indicate that the individual phase θν, as well as the total phase Θ, are relevant in the actual discrete lattice systems.

The third mechanism is higher-order spin interactions, inferred from the above entropic mechanism. In general, the kinetic motion of itinerant electrons induces effective spin interactions, which can be explicitly derived by perturbation expansion in terms of the spin-charge coupling in the Kondo lattice model. The lowest-order contribution is a bilinear interaction called the Ruderman−Kittel−Kasuya−Yosida interaction32,33,34, and the next fourth-order biquadratic interaction was shown to be relevant to stabilize the SkXs35. We here consider a higher-order six-spin contribution given by \(\tilde{L}[{\{{{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{1}}\cdot ({{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{2}}\times {{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{3}})\}}^{2}+{{{{{{{\rm{H.c.}}}}}}}}]\) (see “Methods” and Supplementary Information). It is worthy to note that this has a similar form to the above six-spin entropic term, and it is the lowest-order contribution whose energy depends on Θ under ∑νQν = 0 in the perturbation expansion. This chirality interaction is different from those discussed in the previous studies that stabilize noncoplanar spin states but appear to be irrelevant to the phase shift36,37.

To clarify the effect of the chirality interaction, we investigate the ground-state phase diagram of the effective spin Hamiltonian \({{{{{{{{\mathcal{H}}}}}}}}}^{{{{{{{{\rm{eff}}}}}}}}}\) by variational calculations and the simulated annealing (see “Methods”). We find that the SkX2 gives the lowest energy for 0 ≤ L < 1, while the TVX does for larger L, as shown in Fig. 4a. The optimal θν are obtained as (θ1, θ2, θ3) = (π/3, −π/3, 0) for the former and (π/3, π/6, 0) for the latter (any permutation is allowed), both of which are consistent with those in the Kondo lattice model above. Thus, the chirality interaction prefers the Θ = π/2 states to the Θ = 0 ones, namely, it brings about the phase shift in SkXs.

Fig. 4: Ground states of the effective spin model.
figure 4

a L dependence of the variational energy per site, E, at J = 1 and K = 0.4. The solid line shows the result by the simulated annealing down to temperature T = 10−4 for N = 962. b Phase diagram on the L−J1 plane obtained by the simulated annealing.

Full size image

Furthermore, we find that the value of L necessary for the phase shift can be largely reduced by combining the first mechanism by the exchange interactions for the localized spins. Figure 4b shows the ground-state phase diagram of the Hamiltonian \({{{{{{{{\mathcal{H}}}}}}}}}^{{{{{{{{\rm{eff}}}}}}}}}+{{{{{{{{\mathcal{H}}}}}}}}}^{{{{{{{{\rm{loc}}}}}}}}}\) with J1 only obtained by the simulated annealing. While increasing J1, the phase boundary between the SkX2 and TVX is shifted to a smaller L rapidly. We note that the SkX2 turns into the SkX1 while decreasing J1, whose tendency is similar to the results in Fig. 2a, b. In addition, we obtain the phase transition from the SkX1 to the anisotropic 3Q state while increasing L, similar to the result in Fig. 3b, which also indicates that L plays a similar role to the temperature.

Discussion

We have theoretically guided a new direction of exploring further exotic magnetic states beyond the SkXs by using the phase degree of freedom among the constituent spin density waves. The phase shift turns the SkXs into other states represented by the TVX, which are characterized by staggered emergent magnetoelectric fields and breaking of the lattice rotational symmetry. We unveiled three microscopic mechanisms which drive such a phase shift in itinerant magnets: the exchange interactions between the localized spins, thermal fluctuations, and the long-range chirality interactions.

Our results indicate that the skyrmion-based physics arising from nonzero net chirality can be switched on and off by changing the relative phases among the constituent waves. Furthermore, the lowering of the rotational symmetry by the phase shift induces nonreciprocal transport even in centrosymmetric systems without the spin-orbit coupling. These features open a new direction of emergent electromagnetism by the phase shift engineering. This would be realized in centrosymmetric skyrmion-hosting materials where the multiple-spin interactions rooted in the spin-charge coupling might play an important role11,38,39. While it is not straightforward to identify the phase shift by diffraction techniques such as the neutron scattering and the resonant x-ray scattering, our findings suggest that the angle-resolved photoemission spectroscopy and transport measurements will give good probes to detect the phase shift and new phases like the TVX.

Furthermore, the concept of the phase shift is not limited to the field of skyrmionics but ubiquitously useful for a variety of topological spin crystals, such as vortex crystals and hedgehog crystals, since they are characterized by the multiple-Q spin density waves. Our results suggest that the overlooked phase degree of freedom can induce further interesting topological phase transitions, unconventional electronic structures, topological properties, and conductive phenomena, which will stimulate future exploration of functional spintronics materials in both experiment and theory.

Methods

Kondo lattice model

We consider the Kondo lattice model on a triangular lattice, whose Hamiltonian is given by

$${{{{{{{\mathcal{H}}}}}}}}=-\mathop{\sum}\limits_{i,j,\sigma }{t}_{ij}{c}_{i\sigma }^{{{{\dagger}}} }{c}_{j\sigma }+{J}_{{{{{{{{\rm{K}}}}}}}}}\mathop{\sum}\limits_{i}{{{{{{{{\bf{s}}}}}}}}}_{i}\cdot {{{{{{{{\bf{S}}}}}}}}}_{i}-H\mathop{\sum}\limits_{i}{S}_{i}^{z}.$$
(1)

The first term represents the kinetic energy of itinerant electrons, where \({c}_{i\sigma }^{{{{\dagger}}} }\) (ciσ) is a creation (annihilation) operator of an itinerant electron at site i and spin σ. The second term represents the exchange coupling between itinerant electron spins \({{{{{{{{\bf{s}}}}}}}}}_{i}=(1/2){\sum }_{\sigma ,\sigma ^{\prime} }{c}_{i\sigma }^{{{{\dagger}}} }{{{{{{{{\boldsymbol{\sigma }}}}}}}}}_{\sigma \sigma ^{\prime} }{c}_{i\sigma ^{\prime} }\) [σ = (σx, σy, σz) is the vector of Pauli matrices] and classical localized spins Si with ∣Si∣ = 1. The third term represents the Zeeman coupling to an external magnetic field H. In the calculations, we take the model parameters common to those in ref. 14: the nearest- and third-neighbor hoppings, t1 = 1 and t3 = −0.85, respectively, JK = 1, and the chemical potential μ =−3.5, which gives the characteristic wave vectors at Qν (ν = 1, 2, 3) with Q = π/3 in the main text (we take the lattice constant unity). In this parameter set, the ground state at zero field becomes the SkX2 in Fig. 1e, and that in a field becomes the SkX1 in Fig. 1c14. Note that the SkX2 is stable for other values of Qν while changing the hopping parameters and the electron filling14; our arguments on the phase shift in the main text are not limited to Q = π/3 but can be applied to such other cases. To identify the spin and chirality structure, we compute the spin structure factor

$${S}_{s}({{{{{{{\bf{q}}}}}}}})=\frac{1}{N}\mathop{\sum}\limits_{\alpha =x,y,z}\mathop{\sum}\limits_{j,l}{S}_{j}^{\alpha }{S}_{l}^{\alpha }{{{{{{\rm{e}}}}}}}^{{{{{{\rm{i}}}}}}{{{{{{{\bf{q}}}}}}}}\cdot ({{{{{{{{\bf{r}}}}}}}}}_{j}-{{{{{{{{\bf{r}}}}}}}}}_{l})},$$
(2)

and the chirality structure factor

$${S}_{\chi }({{{{{{{\bf{q}}}}}}}})=\frac{1}{N}\mathop{\sum}\limits_{\mu }\mathop{\sum}\limits_{{{{{{{{\bf{R}}}}}}}},{{{{{{{\bf{R}}}}}}}}^{\prime} \in \mu }{\chi }_{{{{{{{{\bf{R}}}}}}}}}{\chi }_{{{{{{{{\bf{R}}}}}}}}^{\prime} }{{{{{{\rm{e}}}}}}}^{{{{{{\rm{i}}}}}}{{{{{{{\bf{q}}}}}}}}\cdot ({{{{{{{\bf{R}}}}}}}}-{{{{{{{\bf{R}}}}}}}}^{\prime} )},$$
(3)

respectively, where μ = (u, d) represent upward and downward triangles, respectively.

Variational calculation for the Kondo lattice model

In the variational calculations in Fig. 2a, b, we compare the energy of the following spin textures as the variational states: the triple spiral and sinusoidal crystals in Fig. 1c–f while varying three θν under the constraint ∣Si∣ = 1 at each site, the single-Q spiral state characterized by \({{{{{{{{\bf{S}}}}}}}}}_{i}=(\cos {{{{{{{\bf{q}}}}}}}}\cdot {{{{{{{{\bf{r}}}}}}}}}_{i},\sin {{{{{{{\bf{q}}}}}}}}\cdot {{{{{{{{\bf{r}}}}}}}}}_{i},0)\) where q = (0, 0), (Q1, 0), (4π/3, 0) denoted as 1Q K and \((0,2\pi /\sqrt{3})\) denoted as 1Q M, and the conical (1Q C) state characterized by \({{{{{{{{\bf{S}}}}}}}}}_{i}=(1/{N}_{m})(\cos {{{{{{{{\bf{Q}}}}}}}}}_{1}\cdot {{{{{{{{\bf{r}}}}}}}}}_{i},\sin {{{{{{{{\bf{Q}}}}}}}}}_{1}\cdot {{{{{{{{\bf{r}}}}}}}}}_{i},{a}_{z})\) where Nm is the normalization factor and az is the variational parameter. We assume Q = π/3 in the variational calculations, since it was shown that the spin states with Q = π/3 give the lowest grand potential by performing the unbiased Langevin dynamics simulations with the kernel polynomial method when the exchange interactions between the localized spins are zero14. The phase diagram is obtained for the system size N = 962 under the periodic boundary condition.

Finite-temperature calculation for the Kondo lattice model

We adopt the Langevin dynamics simulation with the kernel polynomial method27 to study the finite-temperature properties of the Kondo lattice model in Fig. 3. In the kernel polynomial method, we expand the density of states by up to 2000th order of the Chebyshev polynomials with 162 random vectors. In the Langevin dynamics, we use a projected Heun scheme for 1000−5000 steps with the time interval Δτ = 2. The simulations are done for N = 602, 722, 962, and 1202 sites, and the thermal averages are taken for 100−800 samplings after the thermalization. In the main text, we show the results for N = 962 and 1202, as those for N ≥ 722 are convergent within the error bars. The data at different temperatures are obtained independently starting from different random states. In the simulations, we compute the q components of the scalar chirality \({\chi }_{{{{{{{{\bf{q}}}}}}}}}=\sqrt{{S}_{\chi }({{{{{{{\bf{q}}}}}}}})/N}\) at q = 0, 2Q1, 2Q2, and 2Q3, as plotted in Fig. 3a.

Effective spin model

An effective spin model, which is derived from the Kondo lattice model in equation (1), is given by35

$${{{{{{{{\mathcal{H}}}}}}}}}^{{{{{{{{\rm{eff}}}}}}}}}=2\mathop{\sum }\limits_{\nu = 1}^{3}\left[-J{{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{\nu }}\cdot {{{{{{{{\bf{S}}}}}}}}}_{-{{{{{{{{\bf{Q}}}}}}}}}_{\nu }}+\tilde{K}{({{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{\nu }}\cdot {{{{{{{{\bf{S}}}}}}}}}_{-{{{{{{{{\bf{Q}}}}}}}}}_{\nu }})}^{2}\right]+\tilde{L}\left[{\left\{{{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{1}}\cdot ({{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{2}}\times {{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{3}})\right\}}^{2}+{{{{{{{\rm{H.c.}}}}}}}}\right],$$
(4)

where \({{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{\nu }}=(1/\sqrt{N}){\sum }_{i}{{{{{{{{\bf{S}}}}}}}}}_{i}{e}^{i{{{{{{{{\bf{Q}}}}}}}}}_{\nu }\cdot {{{{{{{{\bf{r}}}}}}}}}_{i}}\). The first two terms describe bilinear and biquadratic interactions, which are derived by second- and fourth-order perturbation expansions in terms of the spin-charge coupling, respectively35; J > 0 and \(\tilde{K}=K/N \; > \; 0\), and N denotes the number of sites. The J and K terms provide a minimal effective model for the Kondo lattice model, stabilizing the SkX2 at zero field (Fig. 3c) and the SkX1 at finite fields (Fig. 3e)35. Meanwhile, the third term with \(\tilde{L}=L/{N}^{2}\) represents an interaction between the scalar spin chirality composed of \({{{{{{{{\bf{S}}}}}}}}}_{{{{{{{{{\bf{Q}}}}}}}}}_{\nu }}\) (see Supplementary Information).

Variational calculation and simulated annealing for the effective spin model

In the variational calculations in Fig. 4a, we compare the energy of the four states in Fig. 1c–f while varying three θν under the constraint ∣Si∣ = 1 at each site. The results are for J = 1 and K = 0.4 while changing L for the system with N = 962 sites. In the simulated annealing in Fig. 4a, b, our simulations are carried out with the standard Metropolis local updates in real space, by reducing the temperature successively, from T = 0.1−1.0 to ≃ 10−4 with the cooling rate of 0.99995−0.99999. The final temperature is typically T = 10−4. Each phase is identified by its spin and chirality configurations. We also start the simulations from the spin structures obtained at low temperatures to determine the phase boundaries between different magnetic states.