Title

Double-Q and quadruple-Q instabilities at low-symmetry ordering wave vectors under tetragonal symmetry

Author(s)Hayami, Satoru

Citation

Physical Review B, 110(1), 014404
https://doi.org/10.1103/PhysRevB.110.014404

Issue Date

2024-07-01

Doc URL

http://hdl.handle.net/2115/92971

Rights©2024 American Physical SocietyTypearticleFile Information

PhysRevB.110.014404.pdf

Instructions for use

Hokkaido University Collection of Scholarly and Academic Papers : HUSCAPPHYSICAL REVIEW B 110, 014404 (2024)

Double-Q and quadruple-Q instabilities at low-symmetryordering wave vectors under tetragonal symmetry
Satoru Hayami*

Graduate School of Science, Hokkaido University, Sapporo 060-0810, Japan(Received 18 February 2024; revised 19 May 2024; accepted 20 June 2024; published 1 July 2024)
Multiple-Q states are expressed as a superposition of spin density waves at multiple ordering wave vectors,
which results in unconventional complicated spin textures, such as skyrmions, hedgehogs, and vortexes. We
investigate the multiple-Q instability by focusing on the low-symmetry ordering wave vectors in momentum
space. By systematically performing the simulated annealing for effective spin models with various orderingwave vectors on a two-dimensional square lattice, we classify the magnetic phase diagram into four types
according to the position of the ordering wave vectors. Three out of four cases lead to a plethora of isotropic
multiple-Q instabilities yielding collinear, coplanar, and noncoplanar double-Q and quadruple-Q magnetic
phases, while the remaining case leads to an anisotropic double-Q instability when the multiple-spin interaction
is introduced. Our results indicate that exotic multiple-Q phases distinct from the skyrmion crystal phase are
expected when the ordering wave vectors lie on the low-symmetry positions in the Brillouin zone.
DOI: 10.1103/PhysRevB.110.014404

I. INTRODUCTION

Multiple-Q states characterized by a superposition of multiple spin density waves have attracted much attention as asource of complicated magnetic structures [1]. Depending on
the type of constituent spin density waves, a variety of spin
textures have been realized [2], such as magnetic skyrmion
crystals [3–11], hedgehog crystals [12–17], meron crystals[18–22], and tetra-axial vortex crystals [23]. For example,
a two-dimensional superposition of the three spiral waves
on a triangular lattice leads to a skyrmion crystal, while athree-dimensional superposition on a cubic lattice leads toa hedgehog crystal. These spin textures often induce topologically nontrivial physical phenomena like the topologicalHall effect [24–31], which would be promising for future
spintronics applications [32–34].
In a classical spin system with a fixed spin length at eachsite, only the isotropic exchange interactions in the form of
Si · S j , where Si stands for the classical spin at site i, are notenough to stabilize the above multiple-Q states in the groundstate. In the case of the skyrmion crystal, additional factors,
such as the Dzyaloshinskii-Moriya interaction [35,36], thermal fluctuations [37], multiple-spin interaction [38], dipolarinteraction [39], and magnetic anisotropy [40–42], are needed
to stabilize the ground-state skyrmion crystal. This is because
the spin configurations as a result of the multiple-Q superposition lead to a nonzero intensity at higher-harmonic wave
vectors, which gives rise to the energy loss compared witha single-Q (1Q) spiral state under the constraint on the spin
length.
Meanwhile, some multiple-Q states do not possess theintensity at higher-harmonic wave vectors depending on the*Contact author: hayami@phys.sci.hokudai.ac.jp2469-9950/2024/110(1)/014404(10)position of the ordering wave vectors in the Brillouin zone.
One of the typical examples is a triple-Q state on a triangularlattice, whose ordering wave vectors lie on the special points
in the Brillouin zone [43–46]. Similar situations also occur
for other lattice structures, when the ordering wave vectors
lie on the special points in the Brillouin zone, as found in
the pyrochlore structure [47], simple cubic lattice [48,49],
checkerboard structure [50], and square lattice [51]. In thesecases, the energies between the 1Q and multiple-Q states
are degenerate with each other within the bilinear isotropic
exchange interaction; the degeneracy is lifted by considering infinitesimally small multiple-spin interactions [46,52]
and thermal fluctuations [53]. Furthermore, another example
has recently been clarified for a system with specific lowsymmetry ordering wave vectors on a square lattice, where
double-Q (2Q) and quadruple-Q (4Q) states without the intensity at higher-harmonic wave vectors are realized [54,55].
Although such a multiple-Q instability has been shown in a
specific situation, there have been no systematic investigations
so far, especially for the case with low-symmetry orderingwave vectors in the Brillouin zone. It is highly desired toclarify the relationship between the positions of the orderingwave vectors and the resultant multiple-Q instability, which
would help to explore materials with complicated magnetic
structures in both theory and experiments.
In this paper, we systematically investigate the multiple-Q
instability on a square lattice by focusing on the lowsymmetry ordering wave vectors that cause magnetic frustration [56,57]. The analysis is based on the numerical simulated
annealing for a simple spin model with the Heisenberg-typeexchange interaction, multiple-spin interaction, and easy-axis
two-spin interaction, where the multiple-spin interaction is
supposed to be small. As a result, we classify the magnetic
phase diagrams into four cases depending on the positionsof the ordering wave vectors. We find that the ordering wave

014404-1©2024 American Physical Society

SATORU HAYAMIPHYSICAL REVIEW B 110, 014404 (2024)

vectors Q = (Qx , Qy ) satisfying the condition (i) Qx = π /2
or Qy = π /2 or (ii) Qx + Qy = π result in isotropic 2Q and
4Q instabilities at zero field, while the anisotropic 2Q instability is caused in the case of other low-symmetry orderingwave vectors; all cases do not require the DzyaloshinskiiMoriya interaction [57]. We also investigate the field-induced
multiple-Q states in all cases. Our results show that a further
intriguing multiple-Q state is expected when ordering wave
vectors lie on the low-symmetry positions in the Brillouinzone.
The rest of this paper is organized as follows. In Sec. II,
the spin model on the square lattice is introduced, and thenumerical method is outlined. Then we discuss the multipleQ instability in the spin model with different low-symmetryordering wave vectors in Sec. III. We describe the four cases
according to the positions of the ordering wave vectors one byone. Section IV is devoted to a summary of this paper.
II. MODEL AND METHOD

To investigate the multiple-Q instability at low-symmetryordering wave vectors in a systematic way, we consider an
effective spin model with a momentum-resolved interaction
on a two-dimensional square lattice under the space group
P4/mmm [58], which is given by
H = −Jz
SQν · S−Qν + I z SQ
Sz
ν −Qν
ν


2
K 
z
z
+
SQν · S−Qν + I z SQS−H
Siz . (1)
ν −Qν
N ν
i
where Si = (Six , Siy , Siz ) is the localized spin at site i; thespin length is fixed to be |Si | = 1. The Qν component ofyzx
the spin is given by SQν = (SQ
, SQ
, SQ
), which is obtainedννν
via the Fourier transformation to Si ; ν is the index of theordering wave vectors Qν . We take the lattice constant of
the square lattice as unity. The first term represents the bilinear exchange interaction with the coupling constant J. We
suppose that the first term originates from the RudermanKittel-Kasuya-Yosida (RKKY) interaction [59–61]; this termis microscopically obtained in the lowest-order perturbative expansions in terms of the Kondo coupling in the
Kondo lattice model. We further introduce the easy-axis-type
anisotropic form factor in the interaction, I z > 0, which arises
from the relativistic spin-orbit coupling; in the case of the Gdand Eu-based compounds without orbital angular momenta,
such anisotropy originates from the anisotropic coupling
to itinerant magnets, which is effectively renormalized into
the anisotropic spin interaction. On the other hand, we neglect other magnetic anisotropy, such as the bond-dependent
magnetic anisotropy, for simplicity [62]. The second term represents the biquadratic exchange interaction with the coupling
constant K, which is obtained as the higher-order contribution of the RKKY interaction in the perturbative expansion
in the Kondo lattice model [38,46,52], where N stands forthe number of total spins in the system. We deal with K asa perturbation term since K appears as a higher-order termthan J in the expansion. In the following, we set J = 1 as theenergy unit of the model and K = 0 or 0.02. The third termO
FIG. 1. Grouping of multiple-Q instability at low-symmetry ordering wave vectors under tetragonal symmetry in the first Brillouinzone. The green hexagons, black circles, red squares, and blue triangle stand for Cases I–IV, respectively.

in Eq. (1) stands for the Zeeman coupling under an external
magnetic field along the z (out-of-plane) direction.
In the model in Eq. (1), we only consider the interaction
at wave vectors that mainly contribute to the ground-stateenergy, which enables us to examine the multiple-Q instability
efficiently [38]. We choose the ordering wave vectors locatedat the low-symmetry positions in the first Brillouin zone.
Since we consider the tetragonal P4/mmm symmetry, it isenough to consider the region for qy < qx < π and 0 < qy < π(denoted as region I) as the independent wave vectors, as
shown in Fig. 1. We set Q1 to be the wave vector in region
I. From the symmetry, there are seven symmetry-related wave
vectors to Q1 in the first Brillouin zone, which are definedby Q2 = R(π /2)Q1 , Q3 = M(x)Q1 , Q4 = R(π /2)Q3 , Q5 =−Q1 , Q6 = −Q2 , Q7 = −Q3 , and Q8 = −Q4 , where R(π /2)
denotes the rotational operation by π /2 around the z axis andM(x) denotes the mirror operation with respect to the xz plane.
As for Q1 , we consider 21 sets by supposing the system size with N = 162 , which are denoted as the symbolsof the green hexagons, black circles, red squares, and blue
triangle in Fig. 1; the different symbols mean the different
multiple-Q instabilities, as detailed in Sec. III. From the microscopic viewpoint, the interactions at the low-symmetryordering wave vectors become dominant when the nesting of
the Fermi surface in itinerant electron systems occurs at thecorresponding wave vectors, which have been recently foundin EuNiGe3 [63–65]. We ignore the effect of the interactionsat other wave vectors since their contributions to the totalenergy are negligible in determining the ground-state spin
configuration. A similar phenomenological approach has beenused for other models, where several mechanisms causing themultiple-Q instability have been found [39,66–70].
We calculate the phase diagram while varying I z and H
for 21 sets of ordering wave vectors based on the simulated

014404-2

DOUBLE-Q AND QUADRUPLE-Q INSTABILITIES AT …

PHYSICAL REVIEW B 110, 014404 (2024)

annealing, which minimizes the energy of the model in Eq. (1)at a low temperature. We consider the two-dimensional squarelattice consisting of N = 162 spins under periodic boundary
conditions. We confirmed that the following results are notchanged by the larger system size with N = 322 . Startingfrom a random spin configuration at the temperature T0 =1.5, we gradually reduce the temperature with a rate Tn+1 =0.999999Tn in each Monte Carlo sweep up to the final temperature T = 0.0001, where Tn is the nth-step temperature. At the
final temperature, we perform 105 –106 Monte Carlo sweeps
for measurements.
For the obtained spin configuration, we calculate the spin
structure factor, which is given by1  η ηS S exp[iq · (ri − r j )],(2)
Ssηη (q) =N i, j i jfor η = x, y, z; ri is the position vector at site i and q is
the wave vector in the first Brillouin zone. The total spin
structure factor is given by Ss (q) = Ssxx (q) + Ssyy (q) + Sszz (q).
It is noted that Ssηη (Q1 ) = Ssηη (Q5 ), Ssηη (Q2 ) = Ssηη (Q6 ),Ssηη (Q3 ) = Ssηη (Q7 ), and Ssηη (Q4 ) = Ssηη (Q8 ). We also calculate the Qν component of the magnetic moments from the spin
structure factor, which is given by ηη
Ss (Qν )
η
.
(3)
mQν =
N
η
In addition, we define the in-plane component of mQ
as
ν
y 2
⊥ 2
x 2
)
=
(m
)
+
(m
)
.
The
net
magnetization
along
the
(mQ
Qν
Qν
ν
field direction is given by
1  z
S.
(4)
Mz =
N i i
III. RESULTS

We show the multiple-Q instabilities in the model with
21 sets of low-symmetry ordering wave vectors. By systematically performing the simulated annealing, we find thatmultiple-Q instabilities are classified into four cases, which
are denoted by the different symbols in Fig. 1. We describe
them one by one in the following Secs. III A–III D.

(a) 1.0(c) 1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.0
0.0

0.5

1.0

1.5

0.0
0.0

2.0

(b) 0.5

0.5

1.0

1.5

2.0

0.5

1.0

1.5

2.0

(d)
0.4

0.4
0.3
0.3
0.2

0.2

0.1

0.1
0.0
0.0

0.5

1.0

1.5

2.0

0.0
0.0

FIG. 2. H dependence of (a) and (c) the magnetization M z andη 2) for ν = 1, 2 and(b) and (d) the squared magnetic moments (mQνη =⊥, z in the model belonging to Case I with Q1 = (5π /8, π /8)
and Q2 = (−π /8, 5π /8) at K = 0. The model parameter I z is set
as (a) and (b) I z = 0 and (c) and (d) I z = 0.1. The vertical lines
represent the phase boundaries between different magnetic phases.
η 2
η 2) = 0 and (mQ) = 0.
It is noted that (mQ34

the spin rotational symmetry under I z = H = 0. When theeffect of H is introduced, the spiral plane is fixed on the xy
plane to gain the Zeeman energy. In other words, the spin
⊥ 2
) = 0state is characterized by the 1Q conical state with (mQ1z 2
and (mQ1 ) = 0, whose real-space spin configuration is shownin Fig. 3(a). This state remains stable up to the saturation
field H = 2, as shown in Figs. 2(a) and 2(b). Then this statecontinuously turns into the fully polarized state for H  2.
By considering the effect of I z , the magnetization curve
exhibits the jump at H ∼ 0.6, as shown in the case of I z = 0.1
in Fig. 2(c), which means that an additional magnetic phaseappears in the presence of I z . In the low-field region, theA. Case I

We consider the case where the ordering wave vectors lie
at low-symmetry positions in the Brillouin zone except forQ1x = π /2, Q1y = π /2, and Q1x + Q1y = π , as shown by thegreen hexagons in Fig. 1. For the system size with N = 162 ,there are 12 independent ordering wave vectors. Among them,
we show the results for Q1 = (5π /8, π /8), although different
choices of the ordering wave vector for Q1 do not affect thequalitative result.
First, we discuss the result in the absence of I z and K.
Figures 2(a) and 2(b) show the H dependence of the fieldinduced magnetization M z and the squared magnetic moments
η 2(mQ) for ν = 1, 2 and η =⊥, z, respectively. In Fig. 2(b),νη 2) for betwe show the results by appropriately sorting (mQν
ter readability. At zero field (H = 0), the 1Q spiral state isstabilized. The spin configuration is given by Si = (cos Q1 ·ri , sin Q1 · ri , 0), where the spiral plane is arbitrary owing to(a) 1Q conical(b) 1Q vertical spiral

1

0

-1
FIG. 3. Real-space spin configurations of (a) the 1Q conical stateat I z = 0 and H = 0.35 and (b) the 1Q vertical spiral state at I z = 0.1
and H = 0.15. The arrows represent the in-plane spin components,
while the color represents the out-of-plane spin component.

014404-3

SATORU HAYAMIPHYSICAL REVIEW B 110, 014404 (2024)(a) 1.0(b) 0.5

0.8

0.4

0.6

0.3

0.4

0.2

0.2

0.1

(a)

(b)

1

0

0.0
0.0

0.5

1.0

1.5

0.0
0.0

2.0

-1
0.5

1.0

1.5

(c)

2.0

η 2
FIG. 4. H dependence of (a) M z and (b) (mQ) for ν = 1, 2 andνη =⊥, z at K = 0.02 and I z = 0. The vertical lines represent thephase boundaries between different magnetic phases. It is noted thatη 2
η 2
) = 0 and (mQ) = 0.
(mQ3
4

η 2
behavior of (mQ) is different from that at I z = 0; nonzeroνz 2(mQν ) appears, which indicates that the spiral plane lies on
the xz or yz plane corresponding to the 1Q vertical spiralstate. Indeed, one finds such a tendency in the real-space spin
configuration in Fig. 3(b). The 1Q vertical spiral state shows
the phase transition into the 1Q conical state with increasing
H, as shown in Fig. 2(d). Thus, no multiple-Q instability
occurs in the absence of K; this is because the energy of the
1Q state is lower than the multiple-Q state in Case I.
Next, we introduce the small biquadratic interaction K =0.02 at I z = 0. As shown in Fig. 4(a), the magnetization looks
continuous against H, which is like that in Fig. 2(a). Meanz 2⊥ 2
while, (mQ) becomes nonzero in addition to (mQ) with12different magnitudes, as shown in Fig. 4(b), which indicatesthe emergence of the anisotropic 2Q state. The spin configuration is characterized by a superposition of the conical spiralwave in the Q1 component and the sinusoidal wave with thez-spin modulation in the Q2 component. Such a feature can beseen in the real-space spin configuration in Fig. 5(a) and thespin structure factor in Fig. 5(b), where the Q2 componentof magnetic moments is dominant and the Q1 component
is subdominant. It is noted that small intensities at higherharmonic wave vectors, such as Q2 + 2Q1 , exist owing to thesuperposition of the 2Q spin density waves at Q1 and Q2 .
A similar anisotropic 2Q state has been clarified in the case

(d)

8

0

4

0

12

0

0

6

0

0

FIG. 6. Real-space spin configurations of the 2Q state at K = 0
and I z = 0 for (a) H = 0.3 and (b) H = 0.35 in the model belonging
to Case II with Q1 = (π /2, π /8). The arrows represent the in-plane
spin components, while the color represents the out-of-plane spincomponent. (c) and (d) Square root of the spin structure factor corresponding to (a) and (b).

where the ordering wave vectors lie in the high-symmetric100 and 110 lines [38,71].
B. Case II

0

In this section, we consider the case where either the x or ycomponent of the ordering wave vectors takes π /2. There are
six possibilities for N = 162 , as denoted by the black circles
in Fig. 1. We specifically set Q1 = (π /2, π /8), although we
confirmed that similar results are obtained for other choices of
Q1 , as detailed below.
In contrast to Case I, we find that a 2Q state has thesame energy as the 1Q conical state at I z = 0 and K = 0,
irrespective of H. We show the real-space spin configurations
and spin structure factors for nonzero H in Figs. 6(a)–6(d),which are obtained from the simulated annealing as one of
the lowest-energy states. As shown in Figs. 6(c) and 6(d), this2Q state is characterized by a superposition of Q1 and Q3[Fig. 6(d)] or that of Q2 and Q4 [Fig. 6(c)]. Thus, constituent
ordering wave vectors are connected by the mirror symmetry
rather than the fourfold rotational symmetry. We find the spin
ansatz in this 2Q state, which is given by⎡ ⎤1 − (M z )2 c1 cos Q1 − 1 − c12 sin Q3⎢⎥⎢⎥
Si = ⎢ 1 − (M z )2  − c sin Q + 1 − c2 cos Q ⎥, (5)
1
1
3 ⎦
⎣
1
Mz

FIG. 5. (a) Real-space spin configurations of the anisotropic 2Q
state at K = 0.02, I z = 0, and H = 0.1. The arrows represent thein-plane spin components, while the color represents the out-of-plane
spin component. (b) Square root of the spin structure factor corresponding to (a).

where Qν = Qν · ri , M z represents the magnetization, and c1
represents the numerical coefficient satisfying 0  c1  1. A
similar spin ansatz is also obtained for a 2Q state with nonzero
Q2 and Q4 components.
The expressions in Eq. (5) exhibit two characteristic points.
The first is that the spin configuration in Eq. (5) connects to

(a)

(b)

1

0

-1

12

0

6

0

014404-4

DOUBLE-Q AND QUADRUPLE-Q INSTABILITIES AT …

PHYSICAL REVIEW B 110, 014404 (2024)(b) 2Q II(a) 2Q I

fully polarized state

1

0

2Q II

-1
(c)

(d)

8

2Q I

8

0

4

FIG. 7. Magnetic phase diagram in the plane of I z and H in themodel belonging to Case II with Q1 = (π /2, π /8) at K = 0.02.
0that in the 1Q conical state when setting c1 = 0 or 1. Thus,
this state is regarded as a superposition of two conical spiralwaves, i.e., a 2Q conical state. The second is that the energy ofthis state is unchanged for 0  c1  1; the 1Q and 2Q conicalstates are energetically degenerate. This is also understoodfrom the fact that the obtained 2Q states do not have theintensities at higher-harmonic wave vectors like Q1 + Q3 andQ2 + Q4 in the spin structure factor, as shown in Figs. 6(c)and 6(d). Indeed, from the expression in Eq. (5), one finds thatthe interference term arising from the normalization condition|Si | = 1 is summarized as sin(Q1 + Q3 ), and it identicallyvanishes. Such a degeneracy is not lifted by I z . To lift the degeneracy between the 1Q and 2Q spin states, the introductionof the multiple-spin interaction K is required.
By introducing K, the energy of the 2Q state becomes
lower than that of the 1Q state since K tends to enhance
the instability toward the multiple-Q states compared withthe 1Q state [38,46,52]. For small K = 0.02, we construct themagnetic phase diagram against I z and H in Fig. 7. There are
two phases below the saturation magnetic field at H = 2: One
is the 2Q I in the low-field region, and the other is the 2Q
II state in the high-field region. In contrast with Case I, no
single-Q state appears in the phase diagram.
In the low-field region, the 2Q I state appears, whose
stability region becomes wider for larger I z . The real-spacespin configuration in this state is shown in Fig. 8(a), where
both xy and z spins exhibit the 2Q modulation. Such a 2Q
feature is also found in the spin structure factor in Fig. 8(c),
z 2
z 2
⊥ 2
⊥ 2
where the relations of (mQ) = (mQ) and (mQ
) = (mQ)1
3
1
3
are satisfied. In contrast with the case at K = 0 in Figs. 6(c)and 6(d), there are intensities at the high-harmonic wave
vectors of Qν ; subdominant peak structures are found at 3Q1
and 3Q3 , which are attributed to the fact that the additional
out-of-plane modulations of both Q1 and Q3 have different
intensities from the in-plane modulations owing to I z and H,
which results in the appearance of the elliptic spiral plane tohave the 3Qν component.
When H increases, the 2Q I state turns into the 2Q IIstate with a jump of the magnetization like the behavior in
Fig. 2(c). The spin configuration in this√ state corresponds tothat in Eq. (5) with setting c1 = 1/ 2; there is no z-spin

0

0

4

0

0

FIG. 8. Real-space spin configurations of (a) the 2Q I state for
H = 0.05 and (b) the 2Q II state for H = 1 at K = 0.02 and I z = 0.2.
The arrows represent the in-plane spin components, while the color
represents the out-of-plane spin component. (c) and (d) Square rootof the spin structure factor corresponding to (a) and (b).

modulation, in contrast with the 2Q I state. The real-spacespin configuration and the spin structure factor are shownin Figs. 8(b) and 8(d), respectively. This state continuously
changes into the fully polarized state.
It is noted that the emergence of the 2Q state in Case II isqualitatively different from that in Case I since, in the absenceof K, the energy of the 2Q state is degenerate with that of the
1Q state in the former, while their energies are different in thelatter. This difference is attributed from the fact that the latter2Q spin structure possesses the intensity at high-harmonic
wave vectors in the spin structure factor, while the former
does not. In this sense, the instability toward the 2Q state isprominent for Case II. Indeed, only the 2Q state appears in
the phase diagram even for small K, as shown in Fig. 7. On
the other hand, the 1Q state remains stable for finite H in CaseI, as shown in Fig. 4(b).
The 2Q I and 2Q II states at Q1 = (π /2, π /8) also appear
for other choices of the ordering wave vectors once either
Q1x or Q1y takes π /2. We show the snapshots of the realspace spin configurations of the 2Q I and 2Q II states in thecases of Q1 = (π /2, π /4) in Fig. 9(a), Q1 = (π /2, 3π /8) in
Fig. 9(b), Q1 = (5π /8, π /2) in Fig. 9(c), Q1 = (3π /4, π /2)
in Fig. 9(d), and Q1 = (7π /8, π /2) in Fig. 9(e), where the left
(right) panel shows the spin configuration of the 2Q I (2Q II)
state. The spin structure factors in these states are qualitatively
like those in Figs. 8(c) and 8(d). Thus, the instability towardthe 2Q spin configurations is common when the ordering wave
vectors include the π /2 modulation.
C. Case III

In this section, we discuss the situation where the ordering wave vectors satisfy Q1x + Q1y = π except for Q1 =(3π /4, π /4), as denoted by the red squares in Fig. 1; we

014404-5

SATORU HAYAMI(a)

PHYSICAL REVIEW B 110, 014404 (2024)
(a)

1

(b)

1

0

0

-1
(c)

-1
(b)

(d)

8

6

1
0

4

0

3

0
0

-1
(c)

1

0

-1
(d)

1

0

-1
(e)

1

0

-1

FIG. 9. Real-space spin configurations of (left panel) the 2Q
I state for H = 0.1 and (right panel) the 2Q II state for H = 1
at K = 0.02 and I z = 0.2, which are obtained in the models withdifferent positions of the ordering wave vectors. The Q1 vector isgiven by (a) Q1 = (π /2, π /4), (b) Q1 = (π /2, 3π /8), (c) Q1 =(5π /8, π /2), (d) Q1 = (3π /4, π /2), and (e) Q1 = (7π /8, π /2).
The arrows represent the in-plane spin components, while the color
represents the out-of-plane spin component.

0

0

0

FIG. 10. Real-space spin configurations of (a) the 2Q state for
H = 0.1 and (b) the 4Q state for H = 0.15 at K = 0 and I z = 0 in themodel belonging to Case III with Q1 = (5π /8, 3π /8). The arrows
represent the in-plane spin components, while the color represents
the out-of-plane spin component. (c) and (d) Square root of the spin
structure factor corresponding to (a) and (b).

denote this case as Case III. We specifically set Q1 =(5π /8, 3π /8) without loss of generality.
For K = 0 and I z = 0, we find that the 1Q, 2Q, and 4Q
states are energetically degenerate. The 2Q state is represented
by a superposition of the conical spiral waves at Q1 and Q4 ,whose spin ansatz is given by⎡ ⎤1 − (M z )2 c1 cos Q1 − 1 − c12 sin Q4⎢⎥⎢⎥
Si = ⎢ 1 − (M z )2  − c sin Q + 1 − c2 cos Q ⎥. (6)
1
1
4 ⎦
⎣
1
Mz
Although the expression is like that in Eq. (5), the constituent
ordering wave vectors are different from each other; the ordering wave vectors in the present 2Q state are connectedby the vertical mirror plane on the [110] line, while those in
the 2Q state in Eq. (5) are connected by the vertical mirror
plane on the [100] line. Meanwhile, this 2Q state does notpossess the intensity at high-harmonic wave vectors in the spin
structure factor like the situation in Sec. III B. Like Case II, Si
is normalized as |Si | = 1 for all i without the normalization
constant owing to the relation sin(Q1 + Q4 ) = 0. The realspace spin configuration and spin structure factor are shownin Figs. 10(a) and 10(c), respectively.
The spin ansatz of the 4Q state is given by⎤⎡A(cos Q1 − sin Q2 − sin Q3 − cos Q4 )⎥⎢Si = ⎣A(sin Q1 + cos Q2 + cos Q3 + sin Q4 )⎦, (7)Mzwhere A = 1 − (M z )2 /2. In contrast with the 2Q state in
Eq. (6), the intensities at Q1 –Q4 are the same as each other,which indicates the isotropic 4Q state; the isotropic intensity

014404-6

DOUBLE-Q AND QUADRUPLE-Q INSTABILITIES AT …

PHYSICAL REVIEW B 110, 014404 (2024)(b) 4Q II(a) 4Q I

1

0

-1
(c)

(d)

6

8

0

FIG. 11. Magnetic phase diagram in the plane of I z and H in themodel belonging to Case III with Q1 = (5π /8, 3π /8) at K = 0.02.

of the spin structure factor is found in Fig. 10(d). It is notedthat there is no intensity at higher-harmonic wave vectors,
which is understood from the expression in Eq. (7); thefactor of cos(Q2 − Q3 ) − cos(Q1 + Q4 ) + sin(Q1 − Q2 ) +sin(Q1 − Q3 ) + sin(Q2 + Q4 ) + sin(Q3 + Q4 ) appearing
from the normalization condition of the spin length isexpressed
asc2 sin π xi sin π yi + c3 (sin π xi + sin π yi )for ri = (xi , yi ), which vanishes for all i (c2 and c3 arecoefficients). Since the modulations of Q1 –Q4 occur in thexy plane, this state is regarded as a 4Q conical state, whose
spin configuration is shown in Fig. 10(b). Like Case II, theintroduction of K lifts the degeneracy among the 1Q and
multiple-Q states.
Next, we turn on I z and K. Figure 11 shows the lowtemperature phase diagram with changing I z and H for small
K = 0.02. By introducing K, the 4Q state becomes the groundstate in the whole region in the phase diagram except for
H = 2, where the fully polarized state appears. There are two
types of the 4Q states, which are denoted as the 4Q I statestabilized in the low-field region and the 4Q II state stabilized
in the high-field region.
The 4Q I state is characterized by a superposition of spin
density waves at Q1 –Q4 with equal intensity, as shown by
the spin structure factor in Fig. 12(c). The real-space spin
configuration is shown in Fig. 12(a), where a complicatednoncoplanar spin texture happens, although there is no net
scalar spin chirality in contrast with the skyrmion crystal.
There are intensities at high-harmonic wave vectors, such as−Q1 + 2Q4 , owing to a superposition of multiple spin density waves in both xy- and z-spin components, as shown in
Fig. 12(c). By increasing H, the 4Q I state changes into the
4Q II state with a jump of the magnetization. The spin configuration of the 4Q II state corresponds to that in Eq. (7), where
the real-space snapshot is presented in Fig. 12(b), and the spin
structure factor in momentum space is shown in Fig. 12(d).
This phase continuously changes into the fully polarized state
at H = 2.
The qualitatively same phase diagram is obtained for another ordering wave vector, i.e., Q1 = (7π /8, π /8). We show
the real-space spin configurations of the 4Q I and 4Q II

3

0

0

0

4

0

0

FIG. 12. Real-space spin configurations of (a) the 4Q I state for
H = 0.1 and (b) the 4Q II state for H = 1 at K = 0.02 and I z = 0.2.
The arrows represent the in-plane spin components, while the color
represents the out-of-plane spin component. (c) and (d) Square rootof the spin structure factor corresponding to (a) and (b).

states in the left and right panels of Fig. 13, respectively.
Compared with the spin configuration in Fig. 12(a), the spin
configuration of the 4Q I state in the left panel of Fig. 13
seems to be different, although the spin structure factor shows
similar profiles to each other. In this way, the instability toward the 4Q state is expected when the ordering wave vectors
belong to Case III.
D. Case IV

Finally, let us discuss the multiple-Q instability in CaseIV, where the ordering wave vector lies at Q1 = (3π /4, π /4),
as denoted by the blue triangle in Fig. 1. Since the magnetic
phase diagram at this ordering wave vector has already beenshown in previous literature [55], we briefly mention the difference from other cases.

1

0

-1
FIG. 13. Real-space spin configurations of (left panel) the 4Q
I state for H = 0.1 and (right panel) the 4Q II state for H = 1 at
K = 0.02 and I z = 0.2, which are obtained for different orderingwave vectors Q1 = (7π /8, π /8). The arrows represent the in-plane
spin components, while the color represents the out-of-plane spincomponent.

014404-7

SATORU HAYAMIPHYSICAL REVIEW B 110, 014404 (2024)

In Case IV, a checkerboard bubble lattice with the collinear
4Q spin configuration appears in the low-field region, whilethe 4Q conical state, whose spin configuration is like that in
Eq. (7), appears in the high-field region when introducing
a small biquadratic interaction K [55]. In other words, theinstability toward the 4Q state occurs at Q1 = (3π /4, π /4)
like Case III. On the other hand, the spin configurationsin the low-field region are different between Cases III and
IV; the noncoplanar spin configuration emerges in Case III[Fig. 12(a)], while the collinear spin configuration emerges inCase IV.

IV. SUMMARY

To summarize, we have investigated the multiple-Q instability on the two-dimensional square lattice by focusing
on the situation where the ordering wave vectors lie on the
low-symmetry position in the Brillouin zone. By performing
numerical simulations based on the simulated annealing forthe spin model with various ordering wave vectors, we findthat the multiple-Q instabilities are classified into four cases.
In Case I (Sec. III A), the ground state corresponds to thesingle-Q state in the bilinear spin model. The anisotropic
double-Q state is induced by considering the effect of thebiquadratic interaction (K). In Case II (Sec. III B), the singleQ and double-Q states are energetically degenerate in the[1] Y. A. Izyumov, Sov. Phys. Usp. 27, 845 (1984).
[2] B. Göbel, I. Mertig, and O. A. Tretiakov, Phys. Rep. 895, 1(2021).
[3] T. H. R. Skyrme, Nucl. Phys. 31, 556 (1962).
[4] A. N. Bogdanov and D. A. Yablonskii, Sov. Phys. JETP 68,
101 (1989).
[5] A. Bogdanov and A. Hubert, Phys. Status Solidi B 186, 527(1994).
[6] U. K. Rößler, A. N. Bogdanov, and C. Pfleiderer, Nature(London) 442, 797 (2006).
[7] S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A.
Neubauer, R. Georgii, and P. Böni, Science 323, 915 (2009).
[8] X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y.
Matsui, N. Nagaosa, and Y. Tokura, Nature (London) 465, 901(2010).
[9] S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A.
Kubetzka, R. Wiesendanger, G. Bihlmayer, and S. Blügel,Nat. Phys. 7, 713 (2011).
[10] N. Nagaosa and Y. Tokura, Nat. Nanotechnol. 8, 899 (2013).
[11] Y. Tokura and N. Kanazawa, Chem. Rev. 121, 2857 (2021).
[12] N. Kanazawa, Y. Nii, X.-X. Zhang, A. Mishchenko, G. De
Filippis, F. Kagawa, Y. Iwasa, N. Nagaosa, and Y. Tokura,Nat. Commun. 7, 11622 (2016).
[13] N. Kanazawa, J. S. White, H. M. Rønnow, C. D. Dewhurst,
D. Morikawa, K. Shibata, T. Arima, F. Kagawa, A. Tsukazaki,
Y. Kozuka et al., Phys. Rev. B 96, 220414(R) (2017).
[14] Y. Fujishiro, N. Kanazawa, T. Nakajima, X. Z. Yu, K. Ohishi,
Y. Kawamura, K. Kakurai, T. Arima, H. Mitamura, A. Miyake
et al., Nat. Commun. 10, 1059 (2019).

bilinear spin model. The isotropic double-Q state is chosen
for small K. In Case III (Sec. III C), the single-Q, double-Q,
and quadruple-Q states are energetically degenerate in thebilinear spin model, and the quadruple-Q state is chosen as
the ground state when K is considered. In Case IV (Sec. III D),
the situation is like that in Case III, although the quadrupleQ spin configurations in the low-field region are different
from each other. Our present results indicate the possibilityof realizing further exotic nontopological multiple-Q states,
as observed in iron-based magnets [72–77], CeRh2 Si2 [78],and Na3 Co2 SbO6 [79] and predicted in the theoretical models [80,81], which will stimulate experimental explorationin future studies. The candidate materials in this paper are
the europium and gadolinium compounds with large total
angular momentum, such as EuPtSi [82–91], EuAl4 [92–96],
EuNiGe3 [63,64,97,98], GdRu2 Si2 [99–104], and GdRu2 Ge2[105], since they often host a variety of multiple-Q states.
ACKNOWLEDGMENTS

This paper was supported by JSPS KAKENHI Grants
No. JP21H01037, No. JP22H04468, No. JP22H00101, No.
JP22H01183, No. JP23H04869, No. JP23K20827, and No.
JP23K03288, and by JST PRESTO (Grant No. JPMJPR20L8)
and JST CREST (Grant No. JPMJCR23O4). Parts of the numerical calculations were performed with the supercomputing
systems in ISSP, the University of Tokyo.

[15] S. Ishiwata, T. Nakajima, J.-H. Kim, D. S. Inosov, N.
Kanazawa, J. S. White, J. L. Gavilano, R. Georgii, K. M.
Seemann, G. Brandl et al., Phys. Rev. B 101, 134406 (2020).
[16] K. Shimizu, S. Okumura, Y. Kato, and Y. Motome, Phys. Rev.
B 103, 054427 (2021).
[17] K. Aoyama and H. Kawamura, Phys. Rev. B 103, 014406(2021).
[18] S.-Z. Lin, A. Saxena, and C. D. Batista, Phys. Rev. B 91,
224407 (2015).
[19] X. Z. Yu, W. Koshibae, Y. Tokunaga, K. Shibata, Y. Taguchi,
N. Nagaosa, and Y. Tokura, Nature (London) 564, 95 (2018).
[20] S. Hayami and Y. Motome, Phys. Rev. Lett. 121, 137202(2018).
[21] S. Hayami and R. Yambe, Phys. Rev. B 104, 094425 (2021).
[22] K. Chen, Q. Luo, Z. Zhou, S. He, B. Xi, C. Jia, H.-G. Luo, andJ. Zhao, New J. Phys. 25, 023006 (2023).
[23] S. Hayami, T. Okubo, and Y. Motome, Nat. Commun. 12, 6927
(2021).
[24] K. Ohgushi, S. Murakami, and N. Nagaosa, Phys. Rev. B 62,
R6065 (2000).
[25] R. Shindou and N. Nagaosa, Phys. Rev. Lett. 87, 116801(2001).
[26] G. Tatara and H. Kawamura, J. Phys. Soc. Jpn. 71, 2613
(2002).
[27] A. Neubauer, C. Pfleiderer, B. Binz, A. Rosch, R. Ritz, P. G.
Niklowitz, and P. Böni, Phys. Rev. Lett. 102, 186602 (2009).
[28] N. Kanazawa, Y. Onose, T. Arima, D. Okuyama, K. Ohoyama,
S. Wakimoto, K. Kakurai, S. Ishiwata, and Y. Tokura, Phys.
Rev. Lett. 106, 156603 (2011).

014404-8

DOUBLE-Q AND QUADRUPLE-Q INSTABILITIES AT …

PHYSICAL REVIEW B 110, 014404 (2024)[29] K. Hamamoto, M. Ezawa, and N. Nagaosa, Phys. Rev. B 92,
115417 (2015).
[30] T. Kurumaji, T. Nakajima, M. Hirschberger, A. Kikkawa, Y.
Yamasaki, H. Sagayama, H. Nakao, Y. Taguchi, T.-H. Arima,
and Y. Tokura, Science 365, 914 (2019).
[31] B. Göbel, A. Mook, J. Henk, and I. Mertig, Phys. Rev. B 95,
094413 (2017).
[32] N. Romming, C. Hanneken, M. Menzel, J. E. Bickel, B.
Wolter, K. von Bergmann, A. Kubetzka, and R. Wiesendanger,
Science 341, 636 (2013).
[33] A. Fert, V. Cros, and J. Sampaio, Nat. Nanotechnol. 8, 152(2013).
[34] A. Fert, N. Reyren, and V. Cros, Nat. Rev. Mater. 2, 17031(2017).
[35] I. Dzyaloshinsky, J. Phys. Chem. Solids 4, 241 (1958).
[36] T. Moriya, Phys. Rev. 120, 91 (1960).
[37] T. Okubo, S. Chung, and H. Kawamura, Phys. Rev. Lett. 108,017206 (2012).
[38] S. Hayami, R. Ozawa, and Y. Motome, Phys. Rev. B 95,
224424 (2017).
[39] O. I. Utesov, Phys. Rev. B 103, 064414 (2021).
[40] A. O. Leonov and M. Mostovoy, Nat. Commun. 6, 8275
(2015).
[41] D. Amoroso, P. Barone, and S. Picozzi, Nat. Commun. 11,5784 (2020).
[42] R. Yambe and S. Hayami, Sci. Rep. 11, 11184 (2021).
[43] I. Martin and C. D. Batista, Phys. Rev. Lett. 101, 156402(2008).
[44] Y. Akagi and Y. Motome, J. Phys. Soc. Jpn. 79, 083711 (2010).
[45] S. Kumar and J. van den Brink, Phys. Rev. Lett. 105, 216405(2010).
[46] S. Hayami and Y. Motome, Phys. Rev. B 90, 060402(R)(2014).
[47] G.-W. Chern, Phys. Rev. Lett. 105, 226403 (2010).
[48] J. L. Alonso, J. A. Capitán, L. A. Fernández, F. Guinea, and V.
Martín-Mayor, Phys. Rev. B 64, 054408 (2001).
[49] S. Hayami, T. Misawa, Y. Yamaji, and Y. Motome, Phys. Rev.
B 89, 085124 (2014).
[50] J. W. F. Venderbos, M. Daghofer, J. van den Brink, and S.
Kumar, Phys. Rev. Lett. 109, 166405 (2012).
[51] D. F. Agterberg and S. Yunoki, Phys. Rev. B 62, 13816 (2000).
[52] Y. Akagi, M. Udagawa, and Y. Motome, Phys. Rev. Lett. 108,096401 (2012).
[53] G.-W. Chern and C. D. Batista, Phys. Rev. Lett. 109, 156801(2012).
[54] S. Hayami, R. Ozawa, and Y. Motome, Phys. Rev. B 94,
024424 (2016).
[55] S. Hayami, Phys. Rev. B 108, 094415 (2023).
[56] J. Wang, Supercond. Sci. Technol. 35, 125006 (2022).
[57] M. Georgiou, I. Rousochatzakis, and J. J. Betouras, Phys. Rev.
B 108, 024412 (2023).
[58] S. Hayami and Y. Motome, J. Phys.: Condens. Matter 33,
443001 (2021).
[59] M. A. Ruderman and C. Kittel, Phys. Rev. 96, 99 (1954).
[60] T. Kasuya, Prog. Theor. Phys. 16, 45 (1956).
[61] K. Yosida, Phys. Rev. 106, 893 (1957).
[62] R. Yambe and S. Hayami, Phys. Rev. B 106, 174437 (2022).
[63] D. Singh, Y. Fujishiro, S. Hayami, S. H. Moody, T. Nomoto,
P. R. Baral, V. Ukleev, R. Cubitt, N.-J. Steinke, D. J. Gawryluk
et al., Nat. Commun. 14, 8050 (2023).

[64] T. Matsumura, K. Kurauchi, M. Tsukagoshi, N. Higa, H.
Nakao, M. Kakihana, M. Hedo, T. Nakama, and Y. Ōnuki,
J. Phys. Soc. Jpn. 93, 074705 (2024).
[65] S. Hayami, Phys. Rev. B 109, 054422 (2024).
[66] S. Hayami and Y. Motome, Phys. Rev. B 103, 024439 (2021).
[67] Z. Wang, Y. Su, S.-Z. Lin, and C. D. Batista, Phys. Rev. B 103,
104408 (2021).
[68] S. Hayami, Phys. Rev. B 105, 174437 (2022).
[69] S. Hayami, J. Phys. Soc. Jpn. 91, 023705 (2022).
[70] S. Hayami and Y. Kato, J. Magn. Magn. Mater. 571, 170547(2023).
[71] R. Ozawa, S. Hayami, K. Barros, G.-W. Chern, Y. Motome,
and C. D. Batista, J. Phys. Soc. Jpn. 85, 103703 (2016).
[72] J. Lorenzana, G. Seibold, C. Ortix, and M. Grilli, Phys. Rev.
Lett. 101, 186402 (2008).
[73] S. Ducatman, N. B. Perkins, and A. Chubukov, Phys. Rev.
Lett. 109, 157206 (2012).
[74] M. N. Gastiasoro and B. M. Andersen, Phys. Rev. B 92,
140506(R) (2015).
[75] J. Allred, K. Taddei, D. Bugaris, M. Krogstad, S. Lapidus, D.
Chung, H. Claus, M. Kanatzidis, D. Brown, J. Kang et al.,Nat. Phys. 12, 493 (2016).
[76] M. H. Christensen, D. D. Scherer, P. Kotetes, and B. M.
Andersen, Phys. Rev. B 96, 014523 (2017).
[77] M. H. Christensen, B. M. Andersen, and P. Kotetes, Phys. Rev.
X 8, 041022 (2018).
[78] H. Saito, F. Kon, H. Hidaka, H. Amitsuka, C. Kwanghee, M.
Hagihala, T. Kamiyama, S. Itoh, and T. Nakajima, Phys. Rev.
B 108, 094440 (2023).
[79] Y. Gu, X. Li, Y. Chen, K. Iida, A. Nakao, K. Munakata, V. O.
Garlea, Y. Li, G. Deng, I. A. Zaliznyak et al., Phys. Rev. B
109, L060410 (2024).
[80] T. Shimokawa and H. Kawamura, Phys. Rev. Lett. 123,057202 (2019).
[81] P. Balla, Y. Iqbal, and K. Penc, Phys. Rev. Res. 2, 043278(2020).
[82] M. Kakihana, D. Aoki, A. Nakamura, F. Honda, M.
Nakashima, Y. Amako, S. Nakamura, T. Sakakibara, M.
Hedo, T. Nakama et al., J. Phys. Soc. Jpn. 87, 023701(2018).
[83] K. Kaneko, M. D. Frontzek, M. Matsuda, A. Nakao, K.
Munakata, T. Ohhara, M. Kakihana, Y. Haga, M. Hedo, T.
Nakama et al., J. Phys. Soc. Jpn. 88, 013702 (2019).
[84] C. Tabata, T. Matsumura, H. Nakao, S. Michimura, M.
Kakihana, T. Inami, K. Kaneko, M. Hedo, T. Nakama, and Y.
Ōnuki, J. Phys. Soc. Jpn. 88, 093704 (2019).
[85] M. Kakihana, D. Aoki, A. Nakamura, F. Honda, M.
Nakashima, Y. Amako, T. Takeuchi, H. Harima, M. Hedo, T.
Nakama et al., J. Phys. Soc. Jpn. 88, 094705 (2019).
[86] A. K. Mishra and V. Ganesan, Phys. Rev. B 100, 125113(2019).
[87] T. Takeuchi, M. Kakihana, M. Hedo, T. Nakama, and Y. Ōnuki,
J. Phys. Soc. Jpn. 88, 053703 (2019).
[88] T. Sakakibara, S. Nakamura, S. Kittaka, M. Kakihana, M.
Hedo, T. Nakama, and Y. Ōnuki, J. Phys. Soc. Jpn. 88, 093701(2019).
[89] T. Takeuchi, M. Kakihana, M. Hedo, T. Nakama, and Y. Ōnuki,
J. Phys. Soc. Jpn. 89, 093703 (2020).
[90] S. Hayami and R. Yambe, J. Phys. Soc. Jpn. 90, 073705(2021).

014404-9

SATORU HAYAMIPHYSICAL REVIEW B 110, 014404 (2024)[91] T. Sakakibara, S. Nakamura, S. Kittaka, M. Kakihana, M.
Hedo, T. Nakama, and Y. Onuki, J. Phys. Soc. Jpn. 90, 064701(2021).
[92] R. Takagi, N. Matsuyama, V. Ukleev, L. Yu, J. S. White, S.
Francoual, J. R. L. Mardegan, S. Hayami, H. Saito, K. Kaneko
et al., Nat. Commun. 13, 1472 (2022).
[93] X. Y. Zhu, H. Zhang, D. J. Gawryluk, Z. X. Zhen, B. C. Yu,
S. L. Ju, W. Xie, D. M. Jiang, W. J. Cheng, Y. Xu et al., Phys.
Rev. B 105, 014423 (2022).
[94] M. Gen, R. Takagi, Y. Watanabe, S. Kitou, H. Sagayama,
N. Matsuyama, Y. Kohama, A. Ikeda, Y. Ōnuki, T. Kurumaji
et al., Phys. Rev. B 107, L020410 (2023).
[95] S. Hayami, J. Phys.: Mater. 6, 014006 (2023).
[96] R. Yang, C.-C. Le, P. Zhu, Z.-W. Wang, T. Shang, Y.-M.
Dai, J.-P. Hu, and M. Dressel, Phys. Rev. B 109, L041113(2024).
[97] R. J. Goetsch, V. K. Anand, and D. C. Johnston, Phys. Rev. B
87, 064406 (2013).
[98] X. Fabrèges, A. Gukasov, P. Bonville, A. Maurya, A.
Thamizhavel, and S. K. Dhar, Phys. Rev. B 93, 214414(2016).

[99] N. D. Khanh, T. Nakajima, X. Yu, S. Gao, K. Shibata, M.
Hirschberger, Y. Yamasaki, H. Sagayama, H. Nakao, L. Peng
et al., Nat. Nanotechnol. 15, 444 (2020).
[100] Y. Yasui, C. J. Butler, N. D. Khanh, S. Hayami, T. Nomoto, T.
Hanaguri, Y. Motome, R. Arita, T. h. Arima, Y. Tokura et al.,Nat. Commun. 11, 5925 (2020).
[101] N. D. Khanh, T. Nakajima, S. Hayami, S. Gao, Y. Yamasaki,
H. Sagayama, H. Nakao, R. Takagi, Y. Motome, Y. Tokura
et al., Adv. Sci. 9, 2105452 (2022).
[102] N. Matsuyama, T. Nomura, S. Imajo, T. Nomoto, R. Arita, K.
Sudo, M. Kimata, N. D. Khanh, R. Takagi, Y. Tokura et al.,Phys. Rev. B 107, 104421 (2023).
[103] G. D. A. Wood, D. D. Khalyavin, D. A. Mayoh, J. Bouaziz,
A. E. Hall, S. J. R. Holt, F. Orlandi, P. Manuel, S. Blügel, J. B.
Staunton et al., Phys. Rev. B 107, L180402 (2023).
[104] S. V. Eremeev, D. Glazkova, G. Poelchen, A. Kraiker, K. Ali,
A. V. Tarasov, S. Schulz, K. Kliemt, E. V. Chulkov, V. S.
Stolyarov et al., Nanoscale Adv. 5, 6678 (2023).
[105] H. Yoshimochi, R. Takagi, J. Ju, N. Khanh, H. Saito,
H. Sagayama, H. Nakao, S. Itoh, Y. Tokura, T. Arima et al.,Nat. Phys. 20, 1001 (2024).

014404-10