# Skyrmion crystal formation and temperature - magnetic field phase diagram of
the frustrated tirangular-lattice Heisenberg magnet with easy-axis masugnetic
anisotropy

Hikaru Kawamura<EMAIL_ADDRESS>Molecular Photoscience
Research Center, Kobe University, Kobe, Hyogo 657-8501, Japan

###### Abstract

The nature of the skyrmion-crystal (SkX) formation and various multiple-$q$
phases encompassing the SkX phase are investigated by extensive Monte Carlo
simulations on the frustrated $J_{1}$-$J_{3}$ triangular-lattice Heisenberg
model with the weak easy-axis magnetic anisotropy. Phase diagram in the
temperature $T$ vs. magnetic-field $H$ plane are constructed, leading to a
rich variety of multiple-$q$ phases. The anisotropy stabilizes the SkX state
down to $T=0$ at intermediate fields, while in the lower-field range the SkX
state becomes only metastable, and new multiple-$q$ states with a broken
$C_{3}$ symmetry are instead stabilized. Implications to experiments are
discussed.

###### pacs:

Much attention has recently been paid to various types of topologically
protected nanoscale spin textures in magnets, e.g., vortex, skyrmion and
hedgehog, from both fundamental interest in topology-related physics and
possible applications to spintronics. Skymion, a swirling noncoplanar spin
texture characterized by an integer topological charge whose constituent spin
directions wrap a sphere in spin space, has got special attention. In
magnetically ordered state, skyrmion is often stabilized as a periodic array
called the skyrmion crystal (SkX). At an earlier stage, the SkX state was
discussed for non-centrosymmetric magnets as induced by the antisymmetric
Dzyaloshinskii-Moriya (DM) interaction Muhlbauer ; Neubauer ; Munzer ; Yu2010
; Yu2011 . In 2012, it was theoretically proposed that the “symmetric” SkX is
also possible in certain class of frustrated centrosymmetric magnets without
the DM interaction, where the size of constituent skyrmions can be varied
continuously from very small to infinitely large (corresponding to the
continuum limit) by tuning the extent of frustration OkuboChungKawamura . An
interesting characteristic of frustration-induced symmetric skyrmion is that,
due to the underlying chiral degeneracy, both skyrmion and antiskyrmion of
mutually opposite signs of topological charge, or the scalar chirality, are
equally possible, leading to the unique and rich electromagnetic responses
OkuboChungKawamura .

In Ref.OkuboChungKawamura , the SkX was identified in a simplified model,
i.e., the frustrated $J_{1}$-$J_{3}$ ($J_{1}$-$J_{2}$) isotropic Heisenberg
model on the triangular lattice as a triple-$q$ state stabilized by magnetic
fields and thermal fluctuations. Subsequent experiment successfully observed
the SkX for centrosymmetric triangular-lattice metallic magnet, e.g.,
Gd2PdSi3, accompanied by the pronounced topological Hall effect Kurumaji .
Recent Monte Carlo (MC) simulation indicated that the SkX could also be
stabilized in the standard RKKY system with only the bilinear interaction
modelling weak-coupling metals, where the oscillating nature of the RKKY
interaction bears frustration Mitsumoto2021 ; Mitsumoto2022 .

Of course, real material possesses various perturbative interactions not taken
into account in a simplified model OkuboChungKawamura , e.g., the three-
dimensionality (interplane coupling), the magnetic anisotropy, quantum
fluctuations, etc. In particular, experiment has indicated that the SkX can be
stabilized even at zero temperature ($T=0$), where the effect of certain
perturbative interactions, e.g., the magnetic anisotropy, was argued to play a
role Kurumaji ; LeonovMostovoy . Possible mechanism leading to the $T=0$ SkX
state was theoretically discussed in the literature, including the biquadratic
interaction arising from the higher-order perturbation beyond the second-order
(strong-coupling effect in itinerant metals) OzawaHayamiMotome2017 ;
HayamiOzawaMotome2017 ; WangBatista , quantum spin fluctuations Rosch , etc.
Among them, the magnetic anisotropy prevails in real magnets, both classical
and quantum, and generally exists even at the spin-bilinear order.

On the basis of the ground-state phase diagram of the frustrated
$J_{1}$-$J_{2}$ triangular Heisenberg model obtained by the simulated
annealing, it was theoretically suggested that the easy-axis magnetic
anisotropy stabilized the SkX state even at $T=0$ LeonovMostovoy . While the
effect of magnetic anisotropy on the SkX formation was examined further by
various authors LinHayami2016 ; HayamiLinBatista ; HayamiMotome2019 ; Wang2020
; Wang2021 ; HayamiXY ; Hayami2022 ; Rosch , most of them concentrated on the
$T=0$ properties, with few studies on the temperature ($T$) vs. magnetic-field
($H$) phase diagram (see LinHayami2016 ; HayamiLinBatista ; HayamiXY ,
however). Even concerning with the $T=0$ properties, the proposed magnetic-
anisotropy stabilization of the SkX state might deserve further careful
examination, since the numerical method employed, e.g., the simulated
annealing, might capture the metastable SkX state, while such a metastable,
not truly stable SkX state was indeed reported under certain annealing
conditions even experimentally Oike ; Karube .

Under these circumstances, we study by extensive MC simulations the SkX
formation and the $T$-$H$ phase diagram of the frustrated $J_{1}$-$J_{3}$
Heisenberg model on the triangular lattice with the easy-axis magnetic
anisotropy, an anisotropic extention of the isotropic model of
Ref.OkuboChungKawamura . We wish to clarify how the $T$-$H$ phase diagram of
the isotropic model changes by the magnetic anisotropy, paying special
attention to the questions of whether the SkX state is truly stabilized at
$T=0$, whether some new phases appear induced by the anisotropy, and if any,
the nature of these phases. We then find that the SkX state is stabilized in
intermediate fields at $T=0$, while its stability range is considerably
reduced compared with that obtained by the simulated annealing, and in the
region where the SkX phase turns out to be only metastable two new anisotropy-
induced multiple-$q$ phases with a broken $C_{3}$ symmetry emerge as stable
phases.

We consider the $J_{1}$-$J_{3}$ classical Heisenberg model on the two-
dimensional triangular lattice with the easy-axis uniaxial anisotropy. The
Hamiltonian is given by

$\displaystyle\mathcal{H}$
$\displaystyle=-J_{1}\sum_{\Braket{i,j}_{1}}(S_{ix}S_{jx}+S_{iy}S_{jy}+\gamma
S_{iz}S_{jz})$ $\displaystyle-$ $\displaystyle J_{3}$
$\displaystyle\sum_{\Braket{i,j}_{3}}(S_{ix}S_{jx}+S_{iy}S_{jy}+\gamma
S_{iz}S_{jz})-H\sum_{i}S_{iz},$ (1)

where $J_{1}>0$ is the ferromagnetic nearest-neighbor coupling, $J_{3}<0$ the
antiferromagnetic third-neighbor coupling,
${\bm{S}}_{i}=(S_{ix},S_{iy},S_{iz})$ a three-component unit vector at site
$i$, magnetic field is applied along the easy axis with $H$ the magnetic-field
intensity, and $\gamma$ the uniaxial exchange anisotropy parameter. We assume
a rather weak easy-axis anisotropy and set $\gamma=1.1$, i.e., 10% anisotropy.
Following Ref.OkuboChungKawamura , we set $J_{1}/J_{3}=-1/3$, and $J_{1}$, $T$
and $H$ are given in units of $|J_{3}|$, hereafter.

MC simulation based on the standard heat-bath method combined with the over-
relaxation method is performed. In addition, fully equilibrated temperature-
exchange simulations are also made in the higher-$T$ range. The lattice is a
$L\times L$ triangular lattice with $L=144,180,216,288$ with periodic boundary
conditions. Unit MC step consists of one heat-bath and $L$ over-relaxation
sweeps. Typically, each run contains $2\times 10^{5}$ MC steps per spin (MCS)
at each temperature, the first half discarded for thermalization.

To reach a given $(T^{*},H^{*})$ state, together with the field-cooling (FC)
run, i.e., the gradual cooling simulated-annealing run at fixed $H^{*}$,
various other computation protocols are tried by combining $H$\- and
$T$-sweeps in search for the stable state. Since, at sufficiently low
$T=T_{0}$, a truly stable state should have the lowest energy among several
metastable states generated by different protocols, it can be determined by
comparing their energies. One standard protocol might be the zero-field
cooling (ZFC) run to ($T_{0},H^{*}$), i.e., gradual cooling in zero field (or
weak fields of $H\lesssim 1.5$) to a low $T=T_{0}$ (we set here $T_{0}=0.1$,
0.05) followed by the gradual increase of $H$ to $H^{*}$ at $T=T_{0}$. Such
ZFC runs are repeated 10$\sim$20 times in search for the lower-energy state by
changing the random numbers and the way of $H$ application. If the ZFC
protocol yields a stable state with the lowest energy at ($T_{0},H^{*}$),
gradual warming run from that state is also performed to higher $T=T^{*}$ at
fixed $H^{*}$ (sometimes further cooling run also made). Consistency is then
checked by confirming the obtained state to be compatible with that obtained
by the $T$-exchange simulations at moderately high $T$.

---
Figure 1:  The phase diagram of the model with $J_{1}/J_{3}=-1/3$ and
$\gamma=1.1$ in the temperature vs. magnetic-field plane. The notations
($mq,nq$), D and U are explained in the text.

The $T$-$H$ phase diagram obtained in this way is shown in Fig.1. It contains
ten distinct ordered phases. Although it might look rather complicated, all
the phases appearing in the isotropic phase diagram OkuboChungKawamura  also
appear here. The triple-$q$ SkX state, described as ($3q,3q$) in Fig.1, is
stabilized down to $T=0$ in the intermediate field range, where $m$ and $n$
(=1,2,3) in ($mq,nq$) represent the number of the dominant (quasi-)Bragg peaks
except for ${\bf q}=0$ in the transverse ($S_{xy}$) and longitudinal ($S_{z}$)
spin structure factors, $S_{\perp}({\bm{q}})$ and $S_{\parallel}({\bm{q}})$
mqnq . The SkX state is essentially of the same type as the one of the
isotropic model, as demonstrated in the real-space spin and scalar-chirality
configurations of Figs. 2(c, d). The SkX state is characterized by the nonzero
total scalar chirality $\chi_{tot}>0$, leading to the topological Hall effect.
The $T$-dependence of the specific heat and $\chi_{tot}$ at $H=3.5$ are shown
in Figs. 2(a, b), the definitions of $\chi_{tot}$ and $S({\bf q})$ being given
in Appendix A. In the higher-$T$ region, the so-called $Z$ phase, the random
domain state of the SkX and the anti-SkX OkuboChungKawamura , also appears as
the collinear triple-$q$ (D,$3q$) state (D means disordered, i.e., the absence
of sharp (quasi-)Bragg peak in $S({\bm{q}})$): See Fig. 11.

---
Figure 2:  The temperature and size dependence of (a) the specific heat, and
of (b) the total scalar chirality, at $H=3.5$. The data of the FC runs and of
the temperature-exchange runs ($T_{{\rm ex}}$) are shown. The typical real-
space configurations of (c) the transverse components of the spin, and of (d)
the scalar chirality, in the SkX state at $H=3.5$ and $T=0.0001$ are shown for
a common part of the $L=180$ lattice.

The single-$q$ spiral states also appear both in zero (or sufficiently weak)
field and in higher fields. The first type, the ($1q,1q$) state, is a vertical
spiral (VS) induced by the easy-axis anisotropy, exhibiting a $90^{\circ}$
rotation from the conical spiral (CS) stabilized in the isotropic model
OkuboChungKawamura . In zero field, this VS state in the $T\rightarrow 0$
limit is vertically coplanar, but it becomes weakly noncoplanar in nonzero
fields, the latter corresponding to the “M state” of Ref.LeonovMostovoy . By
contrast, the high-field single-$q$ state, the ($1q$, U) state (U means
uniform, i.e., only the ${\bm{q}}={\bm{0}}$ peak in $S_{\parallel}(\bm{q})$),
is essentially the same CS as that of the isotropic model OkuboChungKawamura .
In the high-field region, there also appears the double-$q$ state, the
($2q,1q$) state, essentially the same as that of the isotropic model
OkuboChungKawamura . Further details of these states are given in Appendices B
and C.

Between the SkX phase at intermediate $H$ and the VS phase at low enough $H$,
there appear two new phases absent in the isotropic model, i.e., the ($1q,2q$)
and the ($2q,2q$) phases, which persist even in the $T\rightarrow 0$ limit.
The associated $S_{\perp}(\bm{q})$ and $S_{\parallel}(\bm{q})$ are
respectively given in Figs. 3(a, b) and 3(e, f). Note that “$1q$” (“$2q$”)
here means $S({\bm{q}})$ possesses three (quasi-)Bragg peaks, but the $C_{3}$
symmetry is broken resulting in one (two) pair of higher-intensity peaks mqnq
. In both phases, $\chi_{tot}$ vanishes, meaning the absence of the
topological Hall effect.

An important caveat might be in order here: If one makes FC simulated-
annealing runs in the field range $2.0\lesssim H\lesssim 3.6$ to low $T$, one
ends up with the triple-$q$ SkX state even if one makes a very slow cooling.
By contrast, if one makes a ZFC run to ($T_{0},H^{*}$) with $H^{*}$ in the
relevant range, one generally ends up with the $C_{3}$-symmetry broken state.
The energy ($e$) comparison indicates that, for the field $1.35\lesssim
H\lesssim 2.6$ the ($2q,2q$) state reached by the ZFC run is stable ($e$ is
lower than $e$ of the SkX state reached by the FC run by $\sim$0.54% and by
$\sim$0.27% at $H=2$ and 2.5, respectively, well beyond the typical error bar
of order 0.001%); for $2.6\lesssim H\lesssim 3.2$ the ($1q,2q$) state is
stable ($e$ is lower than $e$ of the SkX state by $\sim$0.12% at $H=3$); but
for $3.2\lesssim H\lesssim 4.1$ the ($3q,3q$) SkX state is stable ($e$ is
lower than $e$ of the ($2q,2q$)/($1q,2q$) states by $\sim$0.08% and by
$\sim$0.12% at $H=3.5$ and 4, respectively). In fact, the energy difference
between the ($2q,2q$) and ($1q,2q$) states is rather small of order of the
error bar, although we have observed a clear phase transition between these
two states with varying $T$: See Figs. 4(b) and 6(b).

Our observation then indicates that the easy-axis anisotropy stabilizes the
triple-$q$ SkX state even at $T=0$ at intermediate $H$ ($3.2\lesssim H\lesssim
4.1$), where its stability range is considerably reduced compared with that
obtained by the simulated annealing, and the truly stable state in the
lower-$H$ region turns out to be the ($1q,2q$) state for $2.6\lesssim
H\lesssim 3.2$, and the ($2q,2q$) states for $1.35\lesssim H\lesssim 2.6$. To
the author’s knowledge, these two states, the ($1q,2q$) and ($2q,2q$) states,
are new unnoticed so far. The observed strong hysteretic effect might give the
reason why these states were not reported in the $T=0$ phase diagram
constructed by, e.g., the simulated annealing LeonovMostovoy .

---
Figure 3:  The spin structure factors (a) $S_{\perp}({\bm{q}})$, and (b)
$S_{\parallel}({\bm{q}})$, and the projected spin (c) ($S_{x},S_{y}$), and (d)
($S_{x},S_{z}$) components, in the ($1q,2q$) phase at $H=3$ and $T=0.0001$.
(e)- (h): The corresponding plots in the ($2q,2q$) phase at $H=1.5$ and
$T=0.0001$. The lattice size is $L=180$.

Let us further look into the nature of these new phases. Although the
($1q,2q$) state is a vertical coplanar state as shown in Figs. 3(c, d) qandS ,
it is not a simple VS ($1q,1q$) state. In fact, as can be seen from the
vector-chirality ($\kappa_{x},\kappa_{y}$) projection shown in Fig. 4(c),
where ${\bm{\kappa}}$ is the vector chirality defined on each upward triangle
by ${\bm{\kappa}}=(2/3\sqrt{3})\sum_{<ij>}{\bm{S}}_{i}\times{\bm{S}}_{j}$ (the
summation taken over three clockwise bonds on each triangle) changing sign
under the spatial inversion, ${\bm{\kappa}}$ mostly exhibits parallel
alignment in the direction perpendicular to the coplanar spin plane, but some
${\bm{\kappa}}$ exhibits antiparallel alignment in the opposite direction,
indicating that the spins rotate mostly in a certain (say, clockwise)
direction, but occasionally rotate in an opposite (say, anti-clockwise)
direction. Closer inspection reveals that such a counter-rotation occurs when
the spins stay in the vicinity of the $H$-direction to gain the Zeeman energy.

By contrast, the ($2q,2q$) state is a noncoplanar state as shown in Figs. 3(g,
h). As can be seen from the ($\kappa_{x},\kappa_{z}$) projection of Fig. 4(d),
${\bm{\kappa}}$ is dominated by the horizontal (perpendicular to the field)
component, which suggests that the associated noncoplanar spin configuration
is basically “vertical”. If one compares this with the corresponding plot for
another double-$q$ state, i.e., the ($2q,1q$) state in high fields,
${\bm{\kappa}}$ in the latter rather lacks in the horizontal component (see
Fig. 4(f)), consistently with the meron-like “conical” character of its spin
state (see Fig. 4(e)). Hence, the two double-$q$ states, the ($2q,2q$) and
($2q,1q$) states, are different kinds of states, i.e., “vertical” vs.
“conical”.

---
Figure 4:  The temperature dependence of (a) the specific heat for various
sizes, and of (b) the intensities of the three (quasi-) Bragg peaks of
$S_{\perp}({\bm{q}})$ and $S_{\parallel}({\bm{q}})$ for $L=288$, at $H=2$. The
data are taken by the gradual warming run from the $T=T_{0}$ state as
explained in the text, while in (a) the data of the temperature-exchange run
($T_{{\rm ex}}$) are also added. The projected vector-chirality (c)
($\kappa_{x},\kappa_{y}$) components in the ($1q,2q$) phase at $H=3$ and
$T=0.0001$, and (d) ($\kappa_{x},\kappa_{z}$) components in the ($2q,2q$)
phase at $H=1.5$ and $T=0.0001$. The projected (e) spin ($S_{x},S_{z}$), and
(f) vector-chirality ($\kappa_{x},\kappa_{z}$) components, in the ($2q,1q$)
phase at $H=4.5$ and $T=0.0001$.

In the higher-$T$ range, phases absent in the isotropic model also appear,
including the collinear single-$q$ (D,$1q$) and the collinear double-$q$
(D,$2q$) phases. In Fig.1, the collinear single-$q$ phase appears in two
distinct regions, i.e., at intermediate fields and at zero and weaker fields,
each represented by (D,$1q$) and (D,$1q$)’, which are not connected in the
phase diagram. Indeed, in the (D,$1q$) state, $S_{\parallel}({\bm{q}})$
possesses three pairs of (quasi-)Bragg peaks among which single pair
$\pm{\bm{q}}_{1}^{*}$ exceeds the other twos $\pm{\bm{q}}_{2}^{*}$ and
$\pm{\bm{q}}_{3}^{*}$ by factor of $2\sim 3$ in their intensities (refer to
Figs. 4(b), 6(b) and 10, while, in the (D,$1q$)’ state,
$S_{\parallel}({\bm{q}})$ possesses only one pair of (quasi-)Bragg peaks.

The richness of the phase diagram suggests that even a simple cut of the phase
diagram could yield many phases and phase transitions among them. We
demonstrate such richness by showing the $T$-dependence of physical quantities
at a representative field $H=2$. The data are taken by the gradual warming
runs from the $T=T_{0}$ state prepared by the ZFC run explained above. On
increasing $T$ from $T_{0}$, one encounters the ($2q,2q$), ($1q,2q$),
(D,$2q$), (D,$1q$), (D,$3q$) states before finally reaching the paramagnetic
state. We show in Fig.4 the $T$-dependences of (a) the specific heat, and of
(b) the intensities of the three relevant (quasi-)Bragg peaks
${\bm{q}}_{i}^{*}$ ($i=1,2,3$) of the spin structure factor which are ordered
according to their intensities comment . Similar data for $H=2.5$ are also
given in Fig. 6. As can be seen from the figures, the system indeed exhibits a
rich phase structure.

In view of the complicated appearance of the phase diagram, we now try to give
a rough and intuitive picture of the complicated phase diagram, together with
an intuitive reason why the easy-axis anisotropy energetically stabilizes the
SkX phase at intermediate fields. Depending on the relative strength of the
easy-axis anisotropy and the magnetic field, the phase diagram might be
divided into three regimes, i.e., [I] the high-field regime where the field
exceeds the anisotropy, [II] the low-field regime where the anisotropy exceeds
the field, and [III] the medium-filed region where both compete. In the region
[I] involving the double-$q$ ($2q,1q$) and the CS ($1q$,U) states, the spin
states tend to be “conical” induced by the field, while in the region [II]
involving the VS ($1q,1q$), the double-$q$ ($2q,2q$), and the single-$q$
($1q,2q$) states, the spin states tend to be “vertical” induced by the
anisotropy. Since the conical and vertical states compete with each other, the
states in between tend to be virtually “spherical”, setting the stage for
stabilization of the SkX state.

In summary, by means of extensive MC simulations on the frustrated
$J_{1}$-$J_{3}$ triangular-lattice Heisenberg model with the easy-axis
exchange anisotropy, we have constructed the $T$-$H$ phase diagram containing
a rich variety of multiple-$q$ phases. The easy-axis anisotropy stabilizes the
triple-$q$ SkX state down to $T=0$ at intermediate fields. As the field gets
weaker, the SkX state becomes only metastable, and new multiple-$q$ states
with a broken $C_{3}$ symmetry, the ($2q,2q$) and ($1q,2q$) states, are
instead stabilized. In the high-$T$ regime, in addition to the collinear
triple-$q$ phase ($Z$ phase), the collinear single-$q$ and double-$q$ states
absent in the isotropic model are stabilized by the easy-axis anisotropy.

Finally, we discuss experimental implications of the present result.
Concerning the stability of the SkX and the multiple-$q$ states encompassing
it, while the weak easy-axis magnetic anisotropy enhances the SkX formation
even at $T=0$, it often accompanies a strong hysteretic effect associated with
the $C_{3}$-breaking. Thus, in order to experimentally clarify the SkX-related
phase structure, one needs to examine carefully the possible dependence of the
state on the $T$-cooling/$H$-application protocols. Especially when different
final states are to be obtained by different protocols to a common ($T,H$),
one should determine which state is truly stable. Since the direct comparison
of the energies as we did in the present analysis would be difficult
experimentally, the long-time off-equilibrium measurements toward equilibrium
might eventually be required. Experimental distinction among ($1q,2q$),
($2q,2q$) and ($3q,3q$) from $S({\bm{q}})$ measurements might sometimes be not
easy due to the domain problem, whereas the absence of the topological Hall
effect in the former twos could be used as a signature to distinguish them
from the SkX state.

Of course, features of the phase diagram might well depend on the type and the
strength of the anisotropy, e.g., the $\gamma$ value, as well as on other
perturbative interactions not taken into account in the present model, e.g.,
the dipolar interaction, quantum fluctuations, higher-order exchange
interactions, etc. For example, since the energy difference between the two
new $C_{3}$-broken states, the ($2q,2q$) and ($1q,2q$) states, is rather
small, a small change in $\gamma$ and/or other perturbative effects might
affect their relative stability. While further theoretical and experimental
studies are desirable to fully clarify the effects of these perturbative
interactions, the present work might hopefully serve as a useful starting
reference.

The author is thankful to Prof. T. Sato, Prof. T. Kurumaji and Dr. K.
Mitsumoto for useful discussion. This study was supported by JSPS KAKENHI
Grants No.17H06137 and No.24K00572. We are thankful to ISSP, the University of
Tokyo, for providing us with CPU time.

##  Appendix A Definitions of physical quantities

In this section of appendix, we give definitions of several physical
quantities computed in our Monte Carlo (MC) simulations.

###  A.1 Specific heat

The specific heat is computed generally via the energy fluctuation. In the
vicinity of the first-order transition, to capture the latent-heat
contribution, we also compute it via the temperature ($T$) difference of the
energy per spin $\langle e\rangle$, i.e., $\Delta\langle e\rangle/\Delta T$,
where $\Delta T$ is taken to be 0.001.

###  A.2 Spin structure factors

The transverse and longitudinal spin structure factors, $S_{\perp}(\bm{q})$
and $S_{\parallel}(\bm{q})$, are defined by

$\displaystyle S_{\perp}(\bm{q})$ $\displaystyle=$
$\displaystyle\frac{1}{N}\left\langle\sum_{\mu=x,y}\left|\sum_{i=1}^{N}S_{i\mu}e^{-i\bm{q}\cdot\bm{r}_{i}}\right|^{2}\right\rangle,$
(2) $\displaystyle S_{\parallel}(\bm{q})$ $\displaystyle=$
$\displaystyle\frac{1}{N}\left\langle\left|\sum_{i=1}^{N}S_{iz}e^{-i\bm{q}\cdot\bm{r}_{i}}\right|^{2}\right\rangle,$
(3)

where $N$ is the number of spins, the summation over $i$ is taken over all
sites on the triangular lattice, while $\langle\cdots\rangle$ represents the
thermal average.

###  A.3 Total scalar chirality

The local scalar chirality is defined for the upward (downward) triangle by
$\chi_{\bigtriangleup(\bigtriangledown)}=\bm{S}_{i}\cdot\bm{S}_{j}\times\bm{S}_{k}~{}(i,j,k\in\bigtriangleup(\bigtriangledown))$.
The total scalar chirality is defined by

$\displaystyle\chi_{\rm tot}$
$\displaystyle=\frac{1}{2N}\left(\left\langle\left(\sum_{\bigtriangleup}\chi_{\bigtriangleup}+\sum_{\bigtriangledown}\chi_{\bigtriangledown}\right)^{2}\right\rangle\right)^{1/2},$
(4)

where the summation $\sum_{\bigtriangleup}$ ($\sum_{\bigtriangledown}$) runs
over all upward (downward) triangles on the triangular lattice.

###  A.4 Total vector chirality

We define the local vector chirality for each upward triangle by
${\bm{\kappa}}=(\kappa_{x},\kappa_{y},\kappa_{z})=\frac{2}{3\sqrt{3}}\sum_{\langle
i,j\rangle}({\bm{S}_{i}}\times{\bm{S}_{j}})$, where the summation is taken in
the clockwise direction over three bonds on an upward triangle. The transverse
and longitudinal components of the total vector chirality per plaquette,
$\kappa_{t}$ and $\kappa_{l}$, are then defined by

$\displaystyle\kappa_{t}$
$\displaystyle=\frac{1}{N}\left(\left\langle\left(\sum\kappa_{x}\right)^{2}+\left(\sum\kappa_{y}\right)^{2}\right\rangle\right)^{1/2},$
(5) $\displaystyle\kappa_{l}$
$\displaystyle=\frac{1}{N}\left(\left\langle\left(\sum\kappa_{z}\right)^{2}\right\rangle\right)^{1/2},$
(6)

where the summation is taken over all $N$ upward triangles on the triangular
lattice.

---
Figure 5:  The temperature and size dependence of (a) the specific heat, and
of (c) the total scalar chirality, in zero field $H=0$. The data of the FC
runs and of the temperature-exchange runs ($T_{{\rm ex}}$) are shown. (b) and
(d) are magnified views of the transition region depicted by the dashed-line
boxes in (a) and (c), respectively.

##  Appendix B The temperature dependence of physical quantities in magnetic
fields

In this section, we wish to show our MC data of the temperature ($T$)
dependence of several physical quantities in magnetic fields, which are not
shown in the main text.

###  B.1 $H=0$

We begin with the $H=0$ case. At $H=0$, the system exhibits, on decreasing
$T$, phase transitions from the paramagnetic to the collinear single-$q$
(D,$1q$)’ phase, and then, to the vertical spiral (VS) ($1q,1q$) phase. The
quantity which can be regarded as the order parameter of the VS order might be
the transverse component of the vector chirality, $\kappa_{t}$, defined by
Eq.(4).

We show in Figs. 5(a, c) the $T$-dependence of the specific heat and the
transverse and longitudinal vector chiralities, respectively. Double peaks of
the specific heat associated with the expected two transitions are observed.
The low-$T$ ordered phase is characterized by a nonzero $\kappa_{t}$ with a
vanishing $\kappa_{l}$ defined by Eq.(5), consistently with the VS ordering.
In the intermediate (D,$1q$)’ phase, $\kappa_{t}$ tends to decrease
systematically with the system size $L$, consistently with the expected
collinear ordering. Indeed, the spin structure factor of the intermediate
phase shown in Figs.10(e, f) below also supports the (D,$1q$)’ nature of the
intermediate phase.

---
Figure 6:  The temperature dependence of (a) the specific heat for various
sizes, and of (b) the intensities of the three (quasi-)Bragg peaks of
$S_{\perp}({\bm{q}})$ and $S_{\parallel}({\bm{q}})$ for $L=288$ which are
ordered according to their intensities. The field is $H=2.5$. The data are
taken by the gradual warming runs from the $T=T_{0}$ state prepared by the ZFC
run as explained in the main text. In (a), the data taken by the temperature-
exchange runs ($T_{{\rm ex}}$) are also shown.

###  B.2 $H=2.5$

To demonstrate the richness of the phase structure of the model, we have shown
in Figs. 4(a, b) of the main text the $T$-dependence of the specific heat and
the intensities of the (quasi-)Bragg peaks of the spin structure factor at a
representative field $H=2$, which exhibits a variety of multiple-$q$ phases,
i.e., the ($2q,2q$), ($1q,2q$), (D,$2q$), (D,$1q$), (D,$3q$) states on
increasing $T$ before finally reaching the paramagnetic state. In this
subsection, we show similar plots for a different field $H=2.5$. In Figs. 6(a,
b), we show respectively the $T$-dependences of the specific heat and of the
intensities of the three relevant (quasi-)Bragg peaks ${\bm{q}}_{i}^{*}$
($i=1,2,3$) of the transverse and longitudinal spin structure factors
$S_{\perp}({\bf q})$ and $S_{\parallel}({\bf q})$, which are ordered according
to their intensities. The data are taken by the gradual warming runs from the
$T=T_{0}$ state prepared by the ZFC run as explained in the main text. Similar
phase sequences as observed in Figs. 4(a,b) of the main text for $H=2$ are
also observed for $H=2.5$.

---
Figure 7:  The temperature dependence of physical properties at $H=4.5$. (a)
The specific heat for various sizes, and of (b) the intensities of the three
(quasi-)Bragg peaks of $S_{\perp}({\bm{q}})$ and $S_{\parallel}({\bm{q}})$ for
$L=288$ which are ordered according to their intensities. The data are taken
by the FC runs, while in (a) the data taken by the temperature-exchange runs
($T_{{\rm ex}}$) are also shown.

###  B.3 $H=4.5$

At a higher field $H=4.5$, the system exhibits on decreasing $T$ a phase
transition from the paramagnetic state to the double-$q$ ($2q,1q$) state. The
spin and the vector-chirality configurations at $H=4.5$ have been given in
Figs. 4(e, f) of the main text. We show in Figs. 7(a, b) the $T$-dependence of
the specific heat and the intensities of the three (quasi-)Bragg peaks of the
transverse and longitudinal spin structure factors $S_{\perp}({\bm{q}})$ and
$S_{\parallel}({\bm{q}})$, at $H=4.5$. The occurrence of a single transition
to the ($2q,1q$) ordered state can be seen from the figure.

---
Figure 8:  The temperature dependence of physical properties in the transition
region at $H=6$: (a) The specific heat, and (b) the transverse and
longitudinal vector chiralities, $\kappa_{t}$ and $\kappa_{l}$, for various
sizes. The data are taken by the FC runs.

###  B.4 $H=6$

At a still higher field $H=6$, the conical spiral (CS) ($1q$,U) state
intervenes the paramagnetic and the ($2q,1q$) states, as can be seen from the
$H$-$T$ phase diagram shown in Fig. 1 of the main text. To demonstrate this,
we show in Figs. 8(a, b) the $T$-dependence of the specific heat and of the
transverse and longitudinal vector chiralities, $\kappa_{t}$ and $\kappa_{l}$.
Note that the CS state is characterized by a nonzero $\kappa_{l}$ with a
vanishing $\kappa_{t}$, in contrast to the VS state characterized by a nonzero
$\kappa_{t}$ with a vanishing $\kappa_{l}$. The specific heat shown in Fig.
8(a) now exhibits double peaks, suggesting the appearance of an intermediate
phase. As can be seen from Fig. 8(b), the intermediate state is characterized
by a nonzero $\kappa_{l}$ with a vanishing $\kappa_{t}$, consistently with the
CS nature of the state. We note that the CS state intervening the
paramamagnetic and the ($2q,1q$) states have also been observed in the
isotropic model OkuboChungKawamura .

---
Figure 9:  The projected spin (a,d) ($S_{x},S_{y}$) and (b,e) ($S_{x},S_{z}$)
components, and (c,f) the vector-chirality ($\kappa_{x},\kappa_{y}$)
components in the VS ($1q,1q$) phase state at $T=0.0001$. The upper row
(a,b,c) corresponds to $H=0$, and the lower row (d,e,f) corresponds to $H=1$.
The lattice size is $L=180$.
---
Figure 10:  The transverse and longitudinal spin structure factors,
$S_{\perp}({\bm{q}})$ and $S_{\parallel}({\bm{q}})$, in [left] the (D,$2q$)
phase at $H=2$ and $T=0.35$, [middle] the (D,$1q$) phase at $H=2$ and
$T=0.44$, and [right] the (D,$1q$)’ phase at $H=0$ and $T=0.465$. The upper
row (a,c,e) represents $S_{\perp}({\bm{q}})$, while the lower row (b,d,f)
represents $S_{\parallel}({\bm{q}})$. The lattice size is $L=180$.

##  Appendix C  Properties of each ordered phase

In this section, we show some of the properties of each ordered phase not
given in the main text.

###  C.1 The ($1q,1q$) vertical-spiral phases

We begin with the VS state in zero and lower fields. In Fig. 9, we show the
projected plots of the spin (a) ($S_{x},S_{y}$) and (b) ($S_{x},S_{z}$), and
(c) the vector chirality ($\kappa_{x},\kappa_{y}$) configurations. As can be
seen from Figs. 9(a-c), the VS state at $H=0$ is a vertical coplanar state
with a definite rotation. At a finite field $H=1$, by contrast, the VS becomes
noncoplanar as can be seen from Figs. 9(d-f), as briefly mentioned in the main
text. Yet, the state does not exhibit a counter rotation exhibited by the
($1q,2q$) state as explained in the main text. Compare Figs. 9(c, f) with Fig.
4(c) of the main text.

###  C.2 The collinear (D,$1q$), (D,$2q$) phases

As mentioned in the main text, the present anisotropic model also exhibits at
higher temperatures the collinear phases not realized in the corresponding
isotropic model, including the two types of collinear single-$q$ phase, which
are called (D,$1q$) and (D,$1q$)’ phases each stabilized in intermediate- and
low(zero)-field regimes, and the collinear double-$q$ (D,$2q$) phase. In Fig.
10, we show the transverse and longitudinal spin structure factors,
$S_{\perp}({\bm{q}})$ and $S_{\parallel}({\bm{q}})$, for these three collinear
phases, i.e., (D,$2q$), (D,$1q$) and (D,$1q$)’ phases.

---
Figure 11:  The real-space scalar-chirality configuration and the spin
structure factors in the $Z$ (D,$3q$) state at $H=3.5$ and $T=0.39$. (a) The
color map of the typical real-space configurations of the scalar chirality,
and the transverse and longitudinal spin structure factors, (b)
$S_{\perp}({\bm{q}})$, and (c) $S_{\parallel}({\bm{q}})$. The lattice size is
$L=180$, while a part of the lattice is shown in (a).

###  C.3 The collinear (D,$3q$) phase

The $Z$ phase, i.e., the collinear triple-$q$ (D,$3q$) state, which also
exists in the isotropic model, appears in the anisotropic model, too. In fact,
the anisotropy enhances the stability of this phase considerably in the
$T$-$H$ phase diagram OkuboChungKawamura , while the fundamental character of
the phase remains the same as in the isotropic case. To demonstrate the
random-domain character consisting of the SkX and anti-SkX states, we show in
Fig. 11(a) a typical real-space configuration of the scalar chirality. In
Figs. 11(b,c), we show the corresponding spin structure factors,
$S_{\perp}({\bm{q}})$ and $S_{\parallel}({\bm{q}})$, respectively.

## References

  * (1) S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni, Skyrmion lattice in a chiral magnet, Science 323, 915 (2009).
  * (2) A. Neubauer, C. Pfleiderer, B. Binz, A. Rosch, R. Ritz, P. G. Niklowitz, and P. Böni, Topological Hall Effect in the A Phase of MnSi, Phys. Rev. Lett. 102, 186602 (2009).
  * (3) W. Münzer, A. Neubauer, T. Adams, S. Mühlbauer, C. Franz, F. Jonietz, R. Georgii, P. Böni, B. Pedersen, M. Schmidt, A. Rosch, and C. Pfleiderer, Skyrmion lattice in the doped semiconductor Fe1-xCoxSi, Phys. Rev. B 81, 041203(R) (2010).
  * (4) X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y. Matsui, N. Nagaosa, and Y. Tokura, Real-space observation of a two-dimensional skyrmion crystal, Nature (London) 465, 901 (2010).
  * (5) X. Z. Yu, N. Kanazawa, Y. Onose, K. Kimoto, W. Z. Zhang, S. Ishiwata, Y. Matsui, and Y. Tokura, Near room-temperature formation of a skyrmion crystal in thin-films of the helimagnet FeGe, Nat. Mater. 10, 106 (2011).
  * (6) T. Okubo, S. Chung, and H. Kawamura, Multiple-$q$ States and the Skyrmion Lattice of the Triangular-Lattice Heisenberg Antiferromagnet Under Magnetic Fields, Phys. Rev. Lett. 108, 017206 (2012).
  * (7) T. Kurumaji, T. Nakajima, M. Hirschberger, A. Kikkawa, Y. Yamasaki, H. Sagayama, H. Nakao, Y. Taguchi, T.-h. Arima, and Y. Tokura, Skyrmion lattice with a giant topological Hall effect in a frustrated triangular-lattice magnet, Science 365, 914 (2019).
  * (8) K. Mitsumoto and H. Kawamura, Replica symmetry breaking in the RKKY skyrmion-crystal system, Phys. Rev. B 104, 184432 (2021).
  * (9) K. Mitsumoto and H. Kawamura, Skyrmion crystal in the RKKY system on the two-dimensional triangular lattice, Phys. Rev. B 105, 094427 (2022).
  * (10) A. O. Leonov and M. Mostovoy, Multiply periodic states and isolated skyrmions in an anisotropic frustrated magnet, Nat. Commun. 6, 8275 (2015).
  * (11) R. Ozawa, S. Hayami, and Y. Motome, Zero-Field Skyrmions with a High Topological Number in Itinerant Magnets, Phys. Rev. Lett. 118, 147205 (2017).
  * (12) S. Hayami, R. Ozawa, and Y. Motome, Effective bilinear-biquadratic model for noncoplanar ordering in itinerant magnets, Phys. Rev. B. 95, 224424 (2017).
  * (13) Z. Wang and C. D. Batista, Skyrmion crystals in the triangular Kondo lattice model, SciPost Phys. 15, 161 (2023).
  * (14) V. Lohani, C. Hickey, J. Masell, and A. Rosch, Quantum Skyrmions in Frustrated Ferromagnets, Phys. Rev. X 9, 041063 (2019).
  * (15) S.-Z. Lin and S. Hayami, Ginzburg-Landau theory for skyrmions in inversion-symmetric magnets with competing interactions, Phys. Rev. B 93, 064430 (2016).
  * (16) S. Hayami, S.-Z. Lin, and C. D. Batista, Bubble and skyrmion crystals in frustrated magnets with easy-axis anisotropy, Phys. Rev. B 93, 184413 (2016).
  * (17) S. Hayami and Y. Motome, Effect of magnetic anisotropy on skyrmions with a high topological number in itinerant magnets, Phys. Rev. B 99, 094420 (2019).
  * (18) Z. Wang, Y. Su, S.-Z. Lin, and C. D. Batista, Skyrmion Crystal from RKKY Interaction Mediated by 2D Electron Gas, Phys. Rev. Lett. 124, 207201 (2020).
  * (19) Z. Wang, Y. Su, S.-Z. Lin, and C. D. Batista, Meron, skyrmion, and vortex crystals in centrosymmetric tetragonal magnets, Phys. Rev. B 103, 104408 (2021).
  * (20) S. Hayami, In-plane magnetic field-induced skyrmion crystal in frustrated magnets with easy-plane anisotropy, Phys. Rev. B 103, 224418 (2021).
  * (21) S. Hayami, Skyrmion crystals in centrosymmetric triangular magnets under hexagonal and trigonal single-ion anisotropy, J. Mag. Mag. Mater. 553, 169220 (2022).
  * (22) H. Oike, A. Kikkawa, N. Kanazawa, Y. Taguchi, M. Kawasaki, Y. Tokura, and F. Kagawa, Interplay between topological and thermodynamic stability in a metastable magnetic skyrmion lattice, Nat. Phys. 12, 62 (2016).
  * (23) K. Karube, J. S. White, N. Reynolds, J. L. Gavilano, H. Oike, A. Kikkawa, F. Kagawa, Y. Tokunaga, H. M. Rönnow, Y. Tokura, and Y. Taguchi, Robust metastable skyrmions and their triangular-square lattice structural transition in a high-temperature chiral magnet, Nat. Mater. 15, 1237 (2016).
  * (24) For example, if $S({\bm{q}})$ possesses three (quasi-)Bragg peaks at ${\bm{q}}_{1}$, ${\bm{q}}_{2}$, and ${\bm{q}}_{3}$ with the relation $S({\bm{q}}_{1})=S({\bm{q}}_{2})>S({\bm{q}}_{3})>0$, we call it $2q$.
  * (25) Note that the real-space ($x,y$) (or the ${\bm{q}}$-space ($q_{x},q_{y}$)) part and the spin-space ($S_{x},S_{y}$) part are totally uncoupled in the present model so that any spin configuration generated by aribitrary global spin rotation around the magnetic-field ($S_{z}$) axis is equally possible. In Figs. 3(c) and (d), just for illustation, we take a special choice that the coplanar spin place coincides with the ($S_{x},S_{z}$) plane. The same convention is employed also in Fig. 4(c).
  * (26) The ($2q,2q$) $\leftrightarrow$ ($1q,2q$) transition temperature is estimated to be $T_{c}\simeq 0.23$ from Fig.4(a) and (b). Even within the ($2q,2q$) phase, however, some step-like “structure” or “weak anomaly” exists in Fig.4(b) just below $T_{c}$, while the basic ($2q,2q$) character of $S({\bm{q}})$ remains the same. We regard such “structure” only as secondary one associated with, e.g., the possible incommensurability effect or the domain formation, although we cannot completely rule out the possibility that the another phase with the same ($2q,2q$) ordering pattern intervenes between the ($2q,2q$) and ($1q,2q$) phases just below $T_{c}$.