# Melting and Freezing of a Skyrmion Lattice

Dmitry A. Garanin, Jorge F. Soriano, and Eugene M. Chudnovsky  Physics
Department, Herbert H. Lehman College and Graduate School, The City University
of New York, 250 Bedford Park Boulevard West, Bronx, New York 10468-1589, USA

(September 3, 2024)

###### Abstract

We report comprehensive Monte-Carlo studies of the melting of skyrmion
lattices in systems of small, medium, and large sizes with the number of
skyrmions ranging from $10^{3}$ to over $10^{5}$. Large systems exhibit
hysteresis similar to that observed in real experiments on the melting of
skyrmion lattices. For sufficiently small systems which achieve thermal
equilibrium, a fully reversible sharp solid-liquid transition on temperature
with no intermediate hexatic phase is observed. A similar behavior is found on
changing the magnetic field that provides the control of pressure in the
skyrmion lattice. We find that on heating the melting transition occurs via a
formation of grains with different orientations of hexagonal axes. On cooling,
the fluctuating grains coalesce into larger clusters until a uniform
orientation of hexagonal axes is slowly established. The observed scenario is
caused by collective effects involving defects and is more complex than a
simple picture of a transition driven by the unbinding and annihilation of
dislocation and disclination pairs.

##  I Introduction

Recently, several groups have reported theoretical and experimental studies of
the melting of skyrmion lattices (SkL) in ferromagnetic films [1, 2, 3, 4, 5,
6, 7, 8]. They have been driven by the interest to the fundamental problem of
crystallinity and melting transition in two dimensions (2D), which generated
controversies for almost a century of theoretical and experimental research,
see, e.g., Refs. [9, 10] for review. This paper contains a detailed study of
the melting of the skyrmion lattice. The brief history of the subject is
outlined below.

The study of 2D crystals began with the work of Peierls [11] who had shown
that in 2D the average square deviation of atoms $\langle{\bf u}^{2}\rangle$
from equilibrium positions in a crystal diverges with the size of the system
$L$ as $\ln L$ due to thermal fluctuations. Thirty years later, Mermin and
Wagner [12], and independently Hohenberg, [13] argued that the long-range
order cannot exist in 2D not only for crystals but also for magnets (no net
magnetization) and superfluids (no Bose condensation). This, however,
contradicted the numerical work on hard disks [14] performed in the early
1960s, which showed a melting transition between solid and liquid states.
Also, contrary to the theory, the superfluidity in 2D helium films was
demonstrated experimentally in the late 1960s [15]. It was subsequently
realized that while fluctuations indeed grew in 2D as $\ln L$, they did not
prevent the system from developing rigidity at low temperatures, which for a
2D crystal meant finite elastic moduli.

In 1971 Berezinskii demonstrated theoretically that the destruction of the
long-range order in thin magnetic films and 2D crystals should occur through
the creation of topological defects [16]. One year later, Kosterlitz and
Thouless (KT) [17] arrived at the same conclusion for 2D solids and
superfluids. The KT results for the temperature dependence of the superfluid
density in a helium film were confirmed by Bishop and Reppy in the late 1970s
[18]. Soon, however, Halperin and Nelson [19, 20], and independently Young
[21], realized that as far as the solids were concerned the KT theory was
missing the angular terms in the energy of interaction between dislocations.
By including such terms in the Hamiltonian, they showed that the melting of a
2D crystal occurred into an anisotropic liquid with the exponential decay of
translational correlations but only the algebraic decay of orientational
correlations, which they called a hexatic liquid. The latter, according to the
KTHNY theory, melted into the isotropic liquid due to the creation of
disclinations via a second transition at a higher temperature.

In contrast with the KTHNY theory, however, early simulations that used
molecular dynamics, see, e.g., Ref. [22], emphasized the importance of short-
wave phonons that can drive a conventional first-order melting in 2D as they
do in 3D via the Lindemann melting criterion [23]. (It was recently suggested
[24] that in 2D the Lindemann criterion essentially coincides with the KTHNY
condition for spontaneous creation of dislocations). In addition, in the early
1980s, Chui [25] and independently Saito [26], considering the 2D melting
problem in terms of interacting dislocations, argued that the KTHNY scenario
of melting via the spontaneous creation of independent dislocation pairs by
thermal agitation was self-consistent only at a sufficiently large energy of
the dislocation core, compared to their interaction energy. As the core energy
is lowered, dislocations are created in clusters. Close to the transition
temperature they aggregate into the grain boundaries that break the 2D crystal
into differently oriented grains. Such a possibility was first acknowledged in
the earlier studies of the melting of a 2D electron crystal [27]. The energy
of the dislocation core is difficult to measure directly and it requires a
computation at the level of individual atoms. The lack of that knowledge
creates the biggest uncertainty about the nature of the melting transition in
a 2D system.

Numerous experimental studies performed on 2D crystals in the 1980s and 1990s,
see e.g. Ref. [28], showed a first-order melting transition, as in 3D, with no
sign of the hexatic phase. The existence of the latter at low temperatures due
to disorder was predicted [29, 30] and observed in vortex lattices of type-II
superconductors [31], but convincing evidence of the thermally induced KTHNY
two-step melting emerged from the study of lattices of colloidal particles
[32, 33] only in the 2000s. Subsequently, it was established by numerous
studies that the melting scenario was not universal. Besides the energy of the
dislocation core, it depends on the interaction potential, the shape of the
interacting particles, and the symmetry of the crystal lattice, see, e.g.,
Refs. [34, 35, 36, 37] and references therein.

The interest in skyrmions in magnetic films, besides the beauty of the physics
[38] and imagery associated with them, has been driven by the prospect of
developing topologically protected information technology [39, 40, 41, 42,
43]. In materials lacking inversion symmetry, individual skyrmions are
stabilized by Dzyaloshinskii-Moriya interaction (DMI) [44, 45, 46, 47, 48,
49]. Other mechanisms of stabilization of skyrmions include frustrated
exchange interactions [50, 51], magnetic anisotropy [52, 53], disorder [54],
and geometrical confinement [55].

Theoretical research on 2D lattices of magnetic vortices in systems with
Dzyaloshinskii-Moriya (DMI) interaction was pioneered by Bogdanov and Hubert
[45]. They obtained the magnetic phase diagram separating the conventional
laminar domain structure from the magnetic state with a periodic arrangement
of vortices. The first evidence of regular skyrmion lattices was reported in
FeCoSi films by Yu et al. [56] with the help of the real-space Lorentz
transmission electron microscopy. It was followed by an avalanche of papers
demonstrating transitions between uniformly magnetized states, laminar
domains, and skyrmion crystals on the magnetic field and temperature, see,
e.g., Refs. [57, 58] and references therein.

Figure 1: Triangular skyrmion lattice in the system with free boundary
condition at $T=0$. Spin components are color-coded, $s_{z}=-1$ green,
$s_{z}=1$ orange. White arrows show the in-plane spin components $s_{x}$ and
$s_{y}$.

In recent years it was realized that the skyrmion lattice (see Fig. 1)
provides another rich system for the study of 2D melting. An advantage here is
that the melting transition can be achieved on both temperature and the
magnetic field. Conflicting results have been obtained so far by theorists and
experimentalists. Nishikawa et al. [2] used Monte Carlo simulations for the
lattice-spins model to show a direct melting into an isotropic 2D liquid, with
no evidence of the intermediate hexatic phase. Monte-Carlo studies [7] of
lattices of up to $10^{5}$ skyrmions considered as point particles with
repulsive interaction obtained from a microscopic spin model [59] rendered a
similar conclusion. However, measurements of Huang et al. [4] performed on a
nano-slab of Cu2OSeO3 by using cryo-Lorentz transmission electron microscopy
reported a two-step KTHNY melting transition. Such a transition was also
observed in the simulations of skyrmion lattices in a GaV4S8 spinel by Baláž
et al. [5].

Nonequilibrium effects hinder comparison with the analytical theory of the 2D
melting that assumes thermal equilibrium. For example, thermal studies of
skyrmion lattices in a stack of Ta(5)/Co20Fe60B20(0.9)/Ta(0.08)/MgO(2)/Ta(5)
by Zavorka et al. [3], while providing good data on the orientational
correlations, did not produce consistent results on the nature of the phase
transition due to long equilibration times. McGray et al. [6] reported
hysteresis in the temperature dependence of the orientational order in the
skyrmion lattices of Fe3GeTe2 similar to that in the field-cooled/zero-field-
cooled (FC/ZFC) magnetization studies of the arrays of single-domain magnetic
particles in a solid matrix [60]. According to Meisenheimer et al. [8],
disorder and pinning of skyrmions inhibited simple scenarios of the melting
transition in a polar van der Waals magnet (Fe0.5Co0.5)5GeTe2.

Such a broad range of melting scenarios in the skyrmion lattices can be
explained by several factors. One factor can be the difference in the
dislocation core energy for different materials. The slow equilibration in
both experiments and numerical work is another factor. We try to overcome the
latter problem in our Monte Carlo simulations while also elucidating
nonequilibrium effects that can be seen in real experiments. Ambrose and
Stamps [1] mentioned the distortion of the skyrmion profile on approaching the
melting transition, which they observed in the Monte Carlo studies that used
parameters of FeCoSi and FeGe. It was confirmed that such effect does exist at
the melting transition where it it responsible for the kink in the temperature
dependence of the magnetization [64]. Since the number of magnetic systems in
which skyrmion lattices have been observed experimentally increases rapidly,
in our studies we stuck to the most common range of parameters without
specifying a concrete material. The goal is to elucidate features of the
skyrmion matter that may be seen in the melting and freezing of other 2D
systems.

Apart from elucidating the thermodynamics of the melting/freezing transition,
we investigate dislocations which spontaneously emerge when a 2D skyrmion
solid is heated. The picture that we observe in Monte Carlo simulations and
illustrate with images in the paper and with movies in the supplemental
material, clearly points toward the Chui-Saito melting scenario [25, 26] in
which the dislocations emerge in clusters on approaching the melting
transition. At the transition, they aggregate into grain boundaries that break
the 2D solid with a quasi-long-range orientational order into a system of
fluctuating differently oriented grains. The transition is well defined by the
abrupt change in the orientational order parameter. However, the
orientationally ordered grains of fluctuating shape and size are also seen in
the liquid phase. We study their properties with an eye on real experiments
that can image the structure of the skyrmion matter near the melting/freezing
transition.

The paper is organized as follows. The model of the skyrmion lattice, the
order parameter, and the relevant thermodynamics are explained in Sec. II. The
numerical method is described in Sec. III. Thermodynamic quantities, including
the vicinity of the melting/freezing transition, are computed in Sec. IV. The
characteristics of the defects that spontaneously emerge on rising temperature
are computed and illustrated by images in Sec. V. The nonequilibrium effects
are studied in Sec. VI. Our conclusions and suggestions for experiments are
presented in Sec. VII.

##  II The model

###  II.1 Skyrmions in the lattice-spins model

We consider skyrmions in the model of ferromagnetically coupled three-
component classical spin vectors ${\bf s}_{i}$ of length $1$ on a square
lattice with the energy given by

$\displaystyle\mathcal{H}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{ij}J_{ij}\mathbf{s}_{i}\cdot\mathbf{s}_{j}-H\sum_{i}s_{iz}$
(1)
$\displaystyle-A\sum_{i}\left[(\mathbf{s}_{i}\times\mathbf{s}_{i+\delta_{x}})_{x}+(\mathbf{s}_{i}\times\mathbf{s}_{i+\delta_{y}})_{y}\right].$

The nearest neighbors exchange interaction of strength $J>0$ favors
ferromagnetic ordering and incorporates the actual length of the spin. The
stabilizing field $H$ in the second (Zeeman) term is applied in the negative
z-direction, $H<0$. The third term in Eq. (1) is the Dzyaloshinskii-Moriya
interaction (DMI) of the Bloch type of strength $A$, with $\delta_{x}$ and
$\delta_{y}$ standing for to the nearest lattice site in the positive $x$ or
$y$ direction. In this configuration, the dominant direction of spins is down
and that in the skyrmions is up. The DMI of the Néel type is described by the
term
$(\mathbf{s}_{i}\times\mathbf{s}_{i+\delta_{x}})_{y}-(\mathbf{s}_{i}\times\mathbf{s}_{i+\delta_{y}})_{x}$
and leads to similar results.

Two-dimensional configurations of the spin field are characterized by the
topological charge [38]

$Q=\int\frac{dxdy}{4\pi}\>{\bf s}\cdot\frac{\partial{\bf s}}{\partial
x}\times\frac{\partial{\bf s}}{\partial y}$ (2)

that takes quantized values $Q=0,\pm 1,\pm 2,...$. For the pure-exchange
model, the analytical solution for topological configurations with a given
value of $Q$ was found by Belavin and Polyakov [38] (see also Ref. [59] for a
more transparent approach). Due to the scale invariance of the exchange
interaction in 2D, the energy of a skyrmion or antiskyrmion ($Q=\pm 1$) with
respect to the uniform state, $\Delta E_{\mathrm{BP}}=4\pi J$, does not depend
on their size $\lambda$. In a continuous spin-field model, the conservation of
the topological charge prevents skyrmions from decaying. Finite lattice
spacing $a$ breaks this invariance by adding a term of the order
$-\left(a/\lambda\right)^{2}$ to the energy, which leads to the skyrmion
collapse [61]. This result was generalized for topological structures with any
$Q$ in Ref. [62].

In the presence of the DMI, only skyrmions, but not antiskyrmions, exist. To
find a single-skyrmion spin configuration numerically, one can start with any
state with $Q=1$ and perform energy minimization. The energy of a single
skyrmion with respect to the uniform ferromagnetic state becomes negative for
$A$ large enough and $|H|$ small enough, making the uniform state unstable to
the creation of skyrmions. However, skyrmions compete with the laminar
domains. A stable skyrmion solution exists in the field interval
$H_{s}\leq|H|\leq H_{c}$. Below $H_{s}\simeq 0.55A^{2}/J$ skyrmions become
unstable against converting into a laminar domain structure [7]. Above
$H_{c}\simeq 0.97A^{3/2}/J^{1/2}$ skyrmions collapse [65]. The ratio of these
fields is $H_{c}/H_{s}\simeq 1.76\left(J/A\right)^{1/2}$. In this paper, we
use the value $A/J=0.2$ which is typical for materials with nanoscale
skyrmions, see, e.g. [49, 63]. In this case, $H_{c}/J=0.034$ and
$H_{s}/J=0.022$. Figure 1 shows a Bloch-skyrmion lattice in our model, Eq.(1),
obtained by the numerical energy minimization at $T=0$ for $A/J=0.2$ and
$H/J=-0.025$.

###  II.2 Skyrmions as point particles

There is a short-range repulsion between skyrmions due to DMI [59]

$U(r)\simeq F\exp\left(-\frac{r}{\delta_{H}}\right),\qquad F\equiv
60J\left(\frac{A^{2}}{JH}\right)^{2},$ (3)

where $r$ is the distance between the skyrmions’ centers and
$\delta_{H}=a\sqrt{J/|H|}$ is the magnetic length. In sufficiently dense
skyrmion lattices that we study here, this interaction dominates, whereas
dipole-dipole repulsion of skyrmions can be neglected. This result allows
considering a skyrmion lattice as a system of point particles with a repulsive
interaction, Eq. (3), that is computationally easier than using the full model
of lattice spins, Eq. (1).

The number of skyrmions $N_{S}$ in the system at $T=0$ can be found by the
minimization of the total energy. For a triangular lattice of skyrmions of the
period $a_{S}$ the unit-cell area is $a_{S}^{2}\sqrt{3}/2$ and $N_{S}$ can be
expressed as

$N_{S}=\frac{2}{\sqrt{3}}\frac{S}{a_{S}^{2}},$ (4)

where $S$ is the total area of the system. Each skyrmion in the lattice
interacts with its six nearest neighbors, whereas interactions with further
neighbors can be neglected. The energy per skyrmion in the perfect lattice is

$E_{0}=\varDelta E+3F\exp\left(-\frac{a_{S}}{\delta_{H}}\right),$ (5)

where $\Delta E<0$ is the skyrmion’s core energy computed from the lattice-
spins model [7]. The equilibrium state at $T=0$ can be obtained by the
minimization of the total energy, $\mathcal{E}=N_{S}E_{0}$ with respect to
$a_{S}$, as described in Ref. [7]. For $A/J=0.2$ and $H/J=-0.025$ one has
$\Delta E/J=-4.23$, $\delta_{H}=6.32a$, and $F/J=154$ that results in
$a_{S}=38.5a$. The interaction energy between nearest neighbors
$U_{0}=F\exp\left(-a_{S}/\delta_{H}\right)=0.349J$ defines the scale for the
skyrmion-lattice melting temperature $T_{m}/J\simeq 0.12$.

In the temperature range below as well as above melting the number of
skyrmions in the system does not change because both annihilating a skyrmion
or adding a skyrmion to the system requires overcoming large energy barriers
[7]. In the processes under investigation $\Delta E$ plays no role. One can
consider metastable states in which $N_{S}$ deviates from the value found by
the energy minimization.

The quantity describing the orientation of the hexagon of the nearest
neighbors of any skyrmion $i$ in a 2D lattice is the local hexagonality

$\Psi_{i}=\frac{1}{6}\sum_{j}\exp(6i\theta_{ij}).$ (6)

Here the summation is carried out over the six nearest neighbors $j$ of the
skyrmion $i$ (without the distance cut-off), $\theta_{ij}$ is the angle
between the $ij$ bond and any fixed direction in the lattice. If $\theta$ is
counted from the direction of the $x$-axis (that is our choice), and two of
the sides of the perfect hexagon coincide with the $x$-axis (as in the figures
above), all terms in the sum are equal to $1$, and thus $\Psi_{i}=1$. We call
this horizontal orientation of hexagons and label it by the red color.
Rotation of the hexagons by 30∘ from the horizontal orientation results in the
vertical orientation for which $\Psi_{i}=-1$. For any other orientation of a
perfect hexagon, $\Psi_{i}$ is a complex number of modulus 1. The angle
$\phi_{i}$ by which the hexagon $i$ is rotated from its initial horizontal
orientation is related to the phase angle $\Theta_{i}$ in
$\Psi_{i}=\left|\Psi_{i}\right|e^{i\Theta_{i}}$ as $\phi_{i}=\Theta_{i}/6$.

At finite temperatures, orientations of the bonds fluctuate and the condition
$|\Psi_{i}|=1$ no longer holds. The quality of hexagons can be described by
global hexagonality value defined as

$V_{6}\equiv\sqrt{\frac{1}{N_{S}}\sum_{i}|\Psi_{i}|^{2}}.$ (7)

At temperatures well above the melting transition, when even the short-range
order is completely destroyed, the orientations of the bonds and the angles
$\theta_{ij}$ become random. In this limit $V_{6}=\sqrt{1/6}$ as each particle
has six nearest neighbors on average. The common orientation of hexagons in
the lattice is described by the complex order parameter which is defined as
local hexagonality averaged over the system:

$\Psi=\frac{1}{N_{S}}\sum_{i}\Psi_{i}.$ (8)

We compute the averages of these quantities over the Monte Carlo process
obtaining

$V_{6}\equiv\left\langle\sqrt{\frac{1}{N_{S}}\sum_{i}|\Psi_{i}|^{2}}\right\rangle,\qquad
O_{6}=\left\langle\Psi\right\rangle.$ (9)

In addition, we define

$G_{6}=\sqrt{\left\langle\left|\Psi\right|^{2}\right\rangle}.$ (10)

For the hexagonality value one has $V_{6}^{2}>0.5$ in the solid phase and
$V_{6}^{2}<0.5$ in the liquid phase. In the solid phase
$G_{6}\cong\left|O_{6}\right|$. In the liquid phase $O_{6}$ reduces to zero by
a long thermalization, whereas $G_{6}$ remains finite and scales with
$1/\sqrt{N_{S}}$. One can formally introduce the susceptibility

$\chi=\frac{N_{S}}{3T}\left(\left\langle\left|\Psi\right|^{2}\right\rangle-\left|\left\langle\Psi\right\rangle\right|^{2}\right)$
(11)

that contains both $O_{6}$ and $G_{6}$ and peaks at the melting/freezing
point.

Similarly, one can compute the heat capacity via energy fluctuations:

$C=\frac{d\left\langle E\right\rangle}{dT}=\frac{N_{S}}{T}\left(\left\langle
E^{2}\right\rangle-\left\langle E\right\rangle^{2}\right)$ (12)

with the energy per skyrmion given by

$E=\Delta E+\frac{F}{2N_{S}}\sum_{ij}e^{-r_{ij}/\delta_{H}},$ (13)

where $\Delta E$ does not contribute.

The quantity $G_{6}$ also provides the information about orientational
correlations in the system:

$G_{6}^{2}=\left\langle\frac{1}{N_{S}^{2}}\sum_{ij}\Psi_{i}\Psi_{j}^{*}\right\rangle=\frac{V_{6}^{2}}{N_{S}}\sum_{i-j}C_{6,i-j},$
(14)

where

$C_{6,ij}=\frac{1}{V_{6}^{2}}\left\langle\Psi_{i}\Psi_{j}^{*}\right\rangle.$
(15)

is the correlation function of the hexagonalities. Here the factor $V_{6}^{2}$
is introduced to account for distortions of hexagons at nonzero temperatures
and enforce $C_{6,ii}=1$. If correlations are long, one can use continuous
approximation

$G_{6}^{2}=\frac{V_{6}^{2}}{N_{S}}\int\frac{d^{2}r}{a_{S}^{2}\sqrt{3}/2}C_{6}(r),$
(16)

see Eq. (4). For $C_{6}(r)=\exp\left[-\left(r/R_{6}\right)^{p}\right]$ one
obtains

$G_{6}^{2}=\pi
K_{p}\frac{2V_{6}^{2}R_{6}^{2}}{N_{S}\sqrt{3}a_{S}^{2}}\Longrightarrow\frac{R_{6}}{a_{S}}=\frac{G_{6}}{V_{6}}\sqrt{\frac{\sqrt{3}N_{S}}{2\pi
K_{p}}},$ (17)

where $K_{1}=2$ and $K_{2}=1$. Theory of critical phenomena considers
correlation functions of the type $C(r)\propto r^{2-d-\eta}e^{-r/R}$, where
$d$ is the dimensionality of the space. With this Ansatz, one obtains
$G_{6}^{2}\propto R_{6}^{2-\eta}$ and further the scaling relation
$\gamma=\nu(2-\eta)$. The critical index $\eta$ could be extracted from the
behavior of the correlation function at criticality. However, in our case,
there is a very steep temperature dependence at the melting/freezing point
consistent with only a zero value of the critical index $\beta$ for the
orientational order parameter. Thus the state of the system at criticality is
difficult to define and we are not using the critical Ansatz for the
orientational CF.

##  III Numerical method

Throughout the paper, we consider the system with periodic boundary conditions
(pbc) and sizes $L_{x}\cong L_{y}$ compatible with the lattice without
distortions (see Ref. [7] for details). We studied the small system having the
linear size 30 in the SkL units $a_{S}$ and containing $N_{S}=1080$ skyrmions,
the medium system with the linear size 100 and 1$1600$ skyrmions, and the
large system with the linear size 300 and 104400 skyrmions.

To equilibrate the system at a particular temperature and a magnetic field, we
used the classical Metropolis Monte Carlo (MC) method. Particles were
successively moved by the random vector $\Delta\mathbf{r}$ of the length
$\left|\Delta\mathbf{r}\right|=\mathrm{dispmax\times rand}$, where rand is a
random real number in the range $(0,1)$ and dispmax is the maximal
displacement chosen automatically to maximize the effective displacement
$\mathrm{dispeff}=0.5\mathrm{dispmax}\times\mathrm{acceptancerate}$. It turns
out that near and above melting the maximum of dispeff on dispmax is very flat
and there is a tendency for dispmax to become comparable with $a_{S}$. Because
of that, we introduced a limitation
$\mathrm{dispmax}<\mathrm{dispmaxmax}\equiv 0.3a_{S}$. After each trial move,
the energy change $\Delta E_{\mathrm{trial}}$ was computed and the move was
accepted if $\exp\left(-\Delta E_{\mathrm{trial}}/T\right)>\mathrm{rand}$.

To properly deal with the interacting neighbors, we introduced the interaction
cutoff $r_{\mathrm{cutoff}}=0.95\sqrt{3}a_{S}$ that is just shy of the next-
nearest distance in the triangular lattice, $\sqrt{3}a_{S}$, and subdivided
the system into square bins of the size
$1.05r_{\mathrm{cutoff}}+\mathrm{dispmaxmax}$. For any skyrmion, interacting
neighbors are looked up in the same bin and in the eight surrounding bins.
This drastically reduces the procedure of finding the interacting neighbors
that becomes not straightforward in the liquid phase where the particles are
moving by large distances and the interacting neighbors change. The fast
binning procedure was repeated after each Monte Carlo step (MCS), that is,
after updating all particles in the system.

Figure 2: Parallelized relay computation – a four-core example with lowering
temperature. Figure 3: The number of MCS to thermalize the system strongly
increases in the region of melting/freezing.

For computations in a broad region across the phase transition, it is very
important to employ a Monte Carlo routine with an automated stopping
criterion, as the number of MCS’ necessary for equilibrating the system and
measuring the physical quantities changes by orders of magnitude. Routines of
this kind were used in some preceding work, e.g., in Ref. [7], but here we
formulated the definitive algorithm. First, one chooses the control quantity
(CQ) which is nonzero in both phases. One suitable quantity is the energy but
here we use $V_{6}$ defined by Eq. (7). The MC simulation is performed in
blocks of MCBlock Monte Carlo steps, typical values being
$\mathrm{MCBlock}=50-100$. After the completion of each block, the control
quantity and other physical quantities are computed. The next computational
parameter is $\mathrm{MCMass}<1$ which defines the volume of the computed data
used for data processing and finally for the measurement. We form the lists of
the most recent data of the length $L=\mathrm{Round}[\mathrm{MCMass}\times
N_{\mathrm{val}}]$, where $N_{\mathrm{val}}=\mathrm{MCS}/\mathrm{MCBlock}$ is
the number of computed values of the physical quantities. That is, we use the
fraction MCMass of all computed data, keeping the most recent results, while
older results are discarded. Both $N_{\mathrm{val}}$ and $L$ increase in the
course of the simulation. Then, for the control quantity we compute the moving
averages over the lists of length $L$ using the triangular averaging kernel
$K_{L,i}$

$\mathrm{QCmean}=\sum_{i=1}^{L}K_{L,i}\mathrm{QCList}_{i},\qquad\sum_{i=1}^{L}K_{L,i}=1,$
(18)

and form the list of the same length $L$ from these averages, CQmeanList. The
triangular averaging kernel yields a rather smooth list of data which is used
to quantify the evolution of the control quantity. As $K_{L,i}$ has a maximum
at $i=L/2$, the values from the middle of the averaging interval $L$ make the
largest contribution.

We define the slope of the CQ per one MCS as

$\mathrm{CQ}_{\mathrm{slope}}=\frac{\mathrm{First[CQmeanList]}-\mathrm{Last[CQmeanList]}}{L\times\mathrm{MCBlock}}$
(19)

Then we define the slope parameter

$\mathrm{SlopePar}=\mathrm{CQ}_{\mathrm{slope}}/\left(\mathrm{CQ}_{0}\times\mathrm{MCSens}\right).$
(20)

Here $\mathrm{\mathrm{CQ}_{0}}$ is the typical value of the control quantity,
$\mathrm{CQ}_{0}=1$ for $V_{6}$, and MCSens is the sensitivity of our Monte
Carlo simulation, a very small number. We form the list of the most recent
values of $\left|\mathrm{SlopePar}\right|$ of the same length $L$. The
simulation is stopped when the mean value of this list becomes smaller than
one – this is the most important new feature of the method.

Figure 4: The temperature dependences of the energy per skyrmion (Upper panel.
Inset shows the magnification of the transition region with a negligible
hysteresis) and the heat capacity per skyrmion obtained from the energy
fluctuations, Eq. (12) (Lower panel).

Figure 5: Orientational order parameter $O_{6}$ vs temperature (Upper panel.
Inset shows a jump with a negligible hysteresis) and the temperature
dependence of the orientational quantities on cooling (Lower panel).

One more improvement makes an accent on the critical region that is defined as
the region where $V_{6,\mathrm{floor}}<V_{6}<V_{6,\mathrm{ceiling}}$. We used
$V_{6,\mathrm{floor}}^{2}=0.38$ and $V_{6,\mathrm{ceiling}}^{2}=0.69$. These
bracket the region in which melting/freezing takes place, $V_{6}^{2}\simeq
0.5$. Whenever the computed value of $V_{6}$ falls into this interval,
SlopePar is multiplied by $\mathrm{MCSensFact}=10$. As a result, the stopping
criterion becomes more stringent in the critical region than otherwise. The
method always converges, even in the case of slow large-amplitude stationary
fluctuations because of the increasing $L$ in the definition of
$\mathrm{CQ}_{\mathrm{slope}}$, Eq. (19). After stopping, the latest of the
moving averages of the control quantity, Last[CQmeanList], is output. For the
other quantities, the mean values over the lists of $L$ most recent values are
output. The smaller the parameter MCSens, the longer simulation is done and
the more extensive is the averaging of the computed quantities. Fig. 3
illustrates the number of MCS needed to achieve thermalization for the system
of 11600 skyrmions with $\mathrm{MCSens}=10^{-7}$. The number of MCS increases
in the critical region and peaks at the melting/freezing point $T/J=0.117$.
Since the main computing time is spent near melting/freezing anyway, one can
make computations in a broad range of $T$ with the same $T$-interval. One does
not need to rarify the points away from the phase transition region because
these points cost a relatively short computing time.

The Monte Carlo algorithm presented above is superior to the version used in
Ref. [7]. Whereas the number of MCS done in Ref. [7] does not exceed one
million, here it goes up to several millions that provides a better
thermalization and averaging of stationary fluctuations near the transition.
The main focus here is on the system with $\sim 10^{4}$ skyrmions which is
sufficiently large to study thermodynamics. Here in computing temperature
dependences of the quantities, we used a finer $T$-grid than in Ref. [7].

Computations were performed with Wolfram Mathematica using compilation of the
core routines such as Monte Carlo and those computing the physical quantities.
Dependences on $T$ and $H$ were computed in parallelized cycles using the
number $N_{\mathrm{cores}}$ of the available processor cores and the so-called
relay method (see Fig. 2) in which the next $T$-point within the same thread
uses the previously obtained state as the initial condition. With the interval
between the neighboring $T$-points being $\delta T$, the interval between
those within the same thread is $\Delta T=N_{\mathrm{cores}}\delta T$ that is
sufficiently small if $\delta T$ is small and $N_{\mathrm{cores}}$ is not too
large. We mainly used 8 cores for these computations. The initial conditions
were perfect lattice or random placement of skyrmions.

##  IV Thermodynamics of melting and freezing

Figure 6: Upper panel: Susceptibility associated with the orientational order
peaks at the melting/freezing transition. Lower panel: Critical dependence of
the inverse susceptibility above freezing.

Figure 7: Orientational correlation length $R_{6}$ computed using Eq. (17).

We have found that in our case SkL melting/freezing occurs according to the
grain scenario [25, 26] via only one melting/freezing transition between solid
and liquid phases. In Ref. [7] we concluded that this phase transition is
first order because the orientational order parameter $O_{6}$ jumped to zero
at the melting temperature. However, as we show here, this transition
possesses features of second-order transitions such as divergent heat capacity
and susceptibility. These are more reliable in the determination of the order
of the transition than a jump of the order parameter because in simulations it
is difficult to distinguish a jump from a very steep dependence.

Melting/freezing is a very slow process including slow relaxation and slow
large-amplitude stationary fluctuations. Large areas of the system are melting
and freezing again. Nevertheless, for the system of $10^{4}$ skyrmions
thermodynamic averages at equilibrium can be obtained by the Metropolis Monte
Carlo within $10^{6}-10^{7}$ Monte Carlo steps (MCS), see Fig. 3. With a full
thermalization, there is practically no hysteresis and the melting temperature
$T_{m}$ is practically equal to the freezing temperature $T_{f}$. Hysteresis
is observed in the cases of incomplete equilibration, for instance, if the
stopping criterion for thermalization is too liberal or if the temperature
changes at a too high rate.

Larger systems, such as $10^{5}$ skyrmions, tend to freeze into a polycrystal
state which then evolves very slowly toward the monocrystal state, see Sec.
VI. This process for large systems is even slower than melting/freezing.

Figure 8: Temperature dependence of the pressure in the skyrmion matter.

Figure 9: Freezing of the skyrmion lattice on decreasing the strength of the
magnetic field is a continuous transition with a small Mayer-Wood loop for the
pressure.

The energy per particle $E$ shows a singularity just above the transition, see
Fig. 4. The heat capacity $C$ per particle obtained from the energy
fluctuations, Eq. (12), and shown in the lower panel of Fig. 4 peaks in this
region, although the numerical data does not allow to reliably infer its
critical behavior. The singularity of $C$ is not a robust power law and
resembles a logarithmic divergence. In the low-temperature region the
analytical result is $C=k_{B}$ that is reproduced in the simulation.

The orientational order parameter $O_{6}$ exhibits a linear dependence at low
temperatures (that is typical for classical systems) and a jump or an
unresolved very steep dependence at melting/freezing, see Fig. 5. The average
hexagonality value $V_{6}$ decreases with the temperature and shows a steep
descent just above the melting/freezing point (lower panel of Fig. 5). The
orientational parameter $G_{6}$ is smooth above melting and decreases with the
temperature and the system size. In the skyrmion-liquid phase, $G_{6}\gg
O_{6}$, while in the solid phase $G_{6}\cong O_{6}$, thus the susceptibility
has a peak at melting/freezing, as is shown in the upper panel of Fig. 6. To
obtain this plot, we used $\mathrm{MCSence}=10^{-8}$ with
$\mathrm{MCSensFact}=10$. The inverse susceptibility computed for the systems
of $10^{3}$, $10^{4}$, and $10^{5}$ skyrmions and shown in the lower panel of
the same figure has an apparent critical behavior with the critical index
$\gamma\simeq 1.33$.

Fig. 7 shows the temperature dependence of the orientational correlation
length above melting, obtained from $G_{6}$, Eq. (10). The exponential Ansatz
for the orientational correlation function, see Eq. (17), yields a critical
divergence of the correlation length $R_{6}/a_{S}\simeq
0.42\left[(T-T_{m})/J\right]^{-\nu}$ with $\nu=2/5$.

Figure 10: Early and advanced stages of melting of the SkL slightly above the
melting temperature. The process goes via formation of grains with a different
orientations of hexagons (color coded). Grains are separated by lines of
dislocations shown in black.

The system of skyrmions is a system with a pure repulsion, thus it is always
under pressure that keeps the system together. One can calculate pressure as
the force acting perpendicularly across any line in the system per unit length
of this line. For a perfect lattice at $T=0$ taking into account only the six
nearest neighbors one obtains

$P=\frac{\sqrt{3}f}{a_{S}},\qquad
f=\left.\frac{dU(r)}{dr}\right|_{r=a_{S}}=\frac{F}{\delta_{H}}\exp\left(-\frac{a_{S}}{\delta_{H}}\right).$
(21)

where $f$ is the magnitude of the repulsion force between the two skyrmions.
The temperature dependence of the pressure shown in Fig. 8 is similar to that
of the energy, Fig. 4. In particular, melting causes a steep increase of the
pressure.

Usually pressure is computed for different concentrations of particles. For a
system of skyrmions, it is more convenient to control pressure by the applied
magnetic field as the repulsion energy depends on $H$ via the magnetic length
$\delta_{H}=\sqrt{J/|H|}$ and the exponential prefactor $F$. Increasing the
field strength reduces the interaction and causes melting of the skyrmion
lattice, as shown in Fig. 9. The dependence on $H$ is less steep than the
temperature dependence, and one can see that the phase transition on $H$ is
clearly continuous.

Increasing the strength of the magnetic field reduces the pressure in the
system. However, melting of the skyrmion lattice increases pressure that can
create a Mayer-Wood loop. This loop is difficult to see, though, as it is
masked by a strong regular part of the dependence $P(H)$ that is much more
pronounced than the $P(T)$ dependence. Still, one can see a small Mayer-Wood
loop in the lower panel of Fig. 9.

##  V Spontaneous emergence of lattice defects and formation of grains

Figure 11: Temperature dependence of the number of lattice defects across the
melting transition. Figure 12: Crystal fraction vs temperature for different
lengths of positional averaging.

Figure 13: Skyrmion matter after positional averaging at $T/J=0.12$ and 0.13.
Skyrmions belonging to the solid part are shown by black points and those
belonging to the liquid part are shown by gray points.

The process of melting of a skyrmion lattice goes via spontaneous formation of
grains with different orientation of hexagons, as suggested by Chui (see Fig.
10). Grains are separated from each other by the lines where the lattice is
broken and dislocations are clustered. Clustered dislocations can be
identified as skyrmions having the number of Delaunay nearest neighbors
different from six. Near the melting/freezing temperature, in the course of
the Monte Carlo simulation, grains appear and disappear, grain boundaries move
back and forth. This process is very long. Above the melting temperature, the
tendency of increasing the number of differently oriented grains dominates.
Below this temperature, to the contrary, grains tend to coalesce into bigger
ones, until there is a uniform orientation of hexagons in the system. We do
not observe the so-called “unbinding” of dislocations or disclinations – with
increasing the temperature lattice defects emerge bound in clusters.

Figure 14: Temperature dependence of the number of islands and grains (upper
panel) and of the average size of islands and grains (lower panel).

Fig. 11 shows the temperature dependence of the concentration of lattice
defects: dislocations, free dislocations, and free disclinations. Although
dislocations are defined via the Burgers vector, at nonzero temperatures, when
the lattice is deformed, it is more convenient to look at the Delaunay graph
showing the number of nearest neighbors for each particle. Particles having
five of seven nearest neighbors are disclinations. We count free disclinations
which are separated by a distance of least $2a_{S}$ from other defects. A pair
of 5- and 7-disclinations separated by the distance of order $a_{S}$ form a
dislocation. Again, we count free dislocations which are separated by a
distance of at least $2a_{S}$ from other defects. Also we count all
dislocations, including clustered ones. As can be seen, only a small fraction
of lattice defects are free while the majority of them are clustered.
Moreover, the number of free dislocations and disclinations decreases with
temperature above $T_{m}$. This happens because the total number of
dislocations increases and it becomes difficult to find a place not close to
any of them.

As the temperature increases beyond the melting temperature, the regions
between the grains, where the lattice is broken, become wider. Here the
skyrmion matter become liquid, and the fraction of the liquid phase increases
with temperature. This can be quantified as follows. The temperature
dependence $V_{6}(T)$ shown in the lower panel of Fig. 5 suggests using the
criterion $V_{6}^{2}>0.5$ for solid and $V_{6}^{2}<0.5$ for liquid. Around the
melting/freezing transition, there are regions of both types, and by computing
their areas (or the numbers of skyrmions in them) one can define the
crystalline fraction. This quantity is shown in red in Fig. 12. One can see
that there is some crystalline fraction above melting as well as some liquid
fraction below melting. However, besides slow processes of grain dynamics,
there is a fast quasi-harmonic motion of the particles that creates lattice
deformations decreasing the hexagonality $V_{6}$. To make a more reliable
separation between liquid and solid parts of the system, one can perform a
short averaging of skyrmions’ positions over the Monte Carlo process – so
short that it averages out the quasiharmonic motion but does not affect the
slow grain dynamics. The corresponding results are shown in Fig. 12 for
several different values of the averaging Monte Carlo steps $N_{avr}$
performed after the system was thermalized. As the result of positional
averaging, the liquid part below $T_{m}$ disappears completely that confirms
the hypothesis about its origin as fully due to the quasi-harmonic motion. The
crystalline part above melting increases as well for the same reason. However,
too long averaging begins to reduce the crystalline part at higher
temperatures because here grain dynamics becomes faster. For revealing grains
in the melt, there is an optimum length of the averaging cycle
$N_{avr}=100-300$ MCS, as can be seen in Fig. 12.

One can enumerate the islands of solid within the liquid with the help of the
connected-component labeling algorithm. In this way, one can also find the
average island size (the average number of skyrmions in the island). This
algorithm can be generalized making it resolved with respect to the
orientation of hexagons. In the following we call regions of solid with an
unresolved hexagon orientation “island” and those with a resolved orientation
“grains”. One island can consist of several grains, so there are more grains
than islands. Some grains are separated from each other by dislocations-rich
grain boundaries. There are also continuous changes of grain orientations that
do not lead to creation of dislocations. To define grains in this situation,
we set intervals of $10^{\circ}$ for hexagon orientations. When the
orientation changes by more than 10 degrees, this is considered as another
grain. As the full interval of hexagon orientations is $60^{\circ}$, there are
only six different grain orientations that we take into account here.

Temperature dependences of the number of islands and grains, as well as that
of the average island and grain size, is shown in Fig. 14. Well above $T_{m}$
everything is liquid and the number of islands and grains goes to zero. With
lowering the temperature, grains are formed and their number increases.
However, with a further lowering temperature, islands begin to coalesce, and
their number decreases after reaching a maximum. In this region, the number of
grains still increases as the big coalesced islands consist of many grains
with different orientations. Finally, at the freezing point all grains acquire
the same orientation and become one grain, that is, a monocrystal. Below
freezing, there is one island and one grain. The size of islands and grains in
the lower panel of Fig. 14 monotonically increases with lowering temperature
starting from its minimal value of six particles in the island or grain. Below
$T_{m}$ the island size and the grain size are equal to the number of
skyrmions in the system.

It should be noted that the state with grains of solid within the melt is not
a mixed state because grains are volatile and there is no second temperature,
above the melting/freezing temperature, at which grains completely disappear.
In most of the temperature range, except for a close vicinity of $T_{m}$,
grains are very small to be considered as a thermodynamically stable phase.
Rather, grains are correlated regions above the ordering transition,
precursors of the emerging ordering. The average grain size gives an
estimation of the correlation length above the transition.

##  VI Nonequilibrium effects in the polycrystalline state

Figure 15: Nonequilibrium effects on heating and cooling of the skyrmion
lattice (with the initial random placement of skyrmions) have some similarity
to field-cooled and zero-field-cooled phenomena observed in arrays of magnetic
particles. The data was obtained on the system of 104400 skyrmions by stepwise
increasing and decreasing $T$ with $\delta T/J=0.001$ and the number of MCS at
each temperature 1000, 3000, and 10000. Upper panel: the orientational order
parameter is reduced because the state is polycrystalline. Lower panel: the
average grain size (the number of skyrmions in the grain) in the
polycrystalline state.

Figure 16: The polycrystalline state in the system of 104400 skyrmions at
$T=0$ with 10000 MCS per temperature point. The hexagon orientation is color-
coded as above. Skyrmions having 5 or 7 Delaunay nearest neighbors are shown
in black. Upper panel: small grains obtained by quenching from the initial
state with a random placement of skyrmions. Such states are then warmed up
that leads to coarsening of grains. Lower panel: Large grains obtained by
gradual lowering the temperature.

Large systems of skyrmions on lowering the temperature freeze into a
polycrystalline state. The pure system with periodic boundary conditions that
we are studying here then tends to evolve very slowly towards the
monocrystalline state via the motion of boundaries between differently
oriented grains. If there is only one flat boundary between the two large
grains, the system will stay in this state forever. Also even a weak pinning
or the effect of boundaries will prevent the system from reaching the
monocrystalline state. In particular, a square system with rigid boundaries is
a frustrated system because differently oriented rigid boundaries favor
different orientations of hexagons.

To elucidate the evolution of the polycrystalline state, we performed
stimulated annealing experiments on the large system with 104400 skyrmions by
stepwise increasing and decreasing the temperature with the $T$-step $\delta
T/J=0.001$ and the number of MCS at each temperature 1000, 3000, and 10000. In
the last case, the total number of MCS done in changing the temperature
between zero and $T/J=0.1$ is $10000\times 0.1/0.001=10^{6}$. Both for warming
and cooling, we start with a random placement of skyrmions. At $T=0$, the
system relaxes into a metastable state with small differently oriented grains
and thus a small value of the orientational order parameter $O_{6}$. The
process of grains coarsening is very slow at low temperatures. With increasing
$T$, the coarsening goes faster and $\left|O_{6}\right|$ increases up to the
point of melting before it finally goes down. This behavior has a similarity
to that observed in the zero-field-cooled experiments on magnetic systems.
However, there is no analog of the applied field that causes increasing of the
magnetization with increasing $T$, and
$\left|O_{6}\right|\propto\sqrt{\mathrm{GrainSize}/N_{S}}$ decreases with the
system’s size.

A more suitable quantity characterizing the polycrystalline state is the
average grain size, as defined in the preceding section. This quantity does
not depend on the system size for large systems. As the grain size can be
rather big, studying the polycrystalline state requires large systems. Fig. 15
shows the temperature dependences of $\left|O_{6}\right|$ and the grain size
in the system of 104400 skyrmions. With increasing the number of MCS done at
each $T$-point, both $\left|O_{6}\right|$ and the grain size increase,
although much more MCS is needed to reach the monocrystalline state with
$\left|O_{6}\right|=1$ in the limit $T\rightarrow 0$ and the grain size equal
to the system size.

Examples of the states of the system of $N_{S}=104400$ with 10000 MCS per
temperature point at $T=0$ are shown in Fig. 16. The upper panel shows small
grains obtained by quenching from the initial state with a random placement of
skyrmions. Such states are then warmed up that leads to the coarsening of
grains, see warming curves in Fig. 15. The lower panel shows large grains
obtained by gradual lowering the temperature, see cooling curves in Fig. 15.

##  VII Conclusions

Skyrmion lattices provide another field for testing the predictions of the
theory of 2D melting. They present an advantage over systems of particles due
to the possibility of continuously changing interaction between skyrmions by
changing the external magnetic field. Here we performed comprehensive Monte-
Carlo studies of the melting of 2D skyrmion lattice. The temperature and the
field dependence of the positional and orientational order parameters have
been computed in systems of various sizes ranging from $10^{3}$ to $10^{5}$
skyrmions, and images of skyrmion lattices near the melting transition have
been obtained. Skyrmions repel each other with a force that goes down
exponentially with the distance between the skyrmions. The decay length is
determined by the magnetic field. When the system is small enough, thermal
equilibrium is achieved, and we observe a sharp reversible solid-liquid
transition without an intermediate hexatic phase. Larger systems fail to
achieve thermal equilibrium and exhibit hysteresis on cycling the temperature.
For both equilibrium and non-equilibrium systems we find that the melting and
freezing transitions occur via the formation of grains with different
orientations of hexagonal axes. The departure from a simple picture of a
transition caused by the unbinding of dislocation and disclination pairs is
due to the clustering of the defects into the grain boundaries. We observe the
granular structure on both heating and cooling, with the only difference that
in the liquid phase it corresponds to the grains of a fluctuating shape while
in the solid phase it corresponds to a more stable polycrystal. On cooling the
fluctuating grains freeze into a polycrystal which behaves similarly to
ordinary crystals on annealing. Our findings show that the scenarios of
melting driven by the unbinding of free defects and melting driven by the
formation of grains are likely to be the two limiting cases of a more complex
scenario in which many-body interaction between the defects plays a crucial
role.

##  VIII Acknowledgments

This work was supported by Grant No. DE-FG02-93ER45487 funded by the U.S.
Department of Energy, Office of Science, and Grant No. FA9550-24-1-0090 funded
by the Air Force Office of Scientific Research.

## References

  * [1] M. C. Ambrose and R. L. Stamps, Melting of hexagonal skyrmion states in chiral magnets, New Journal of Physics 15, 053003-(14) (2013).
  * [2] Y. Nishikawa, K. Hukushima, and W. Krauth, Solid-liquid transition of skyrmions in a two-dimensional chiral magnet, Physical Review B 99, 064435 (2019).
  * [3] J. Zázvorka, F. Dittrich, Y. Ge, N. Kerber, K. Raab, T. Winkler, K. Litzius, M. Veis, P. Virnau, and M. Kläui, Skyrmion lattice phases in thin film multilayers, Advanced Functional Materials 30, 2004037-(8) (2020).
  * [4] P. Huang, T. Schönenberger, M. Cantoni, L. Heinen, A. Magrez, A. Rosch, F. Carbone, and H. M. Rønnow, Melting of a skyrmion lattice to a skyrmion liquid via a hexatic phase, Nature Nanotechnology 15, 761-767 (2020).
  * [5] P. Baláž, M. Paściak, and J. Hlinka, Melting of Néel skyrmion lattice, Physical Review B 103, 174411-(8) (2021).
  * [6] A. R. C. McCray, Yue Li, R. Basnet, K. Pandey, J. Hu, D. P. Phelan, X. Ma, A. K. Petford-Long, and C. Phatak, Thermal hysteresis and ordering behavior of magnetic skyrmion lattices, Nanoletters 22, 7804-7810 (2022).
  * [7] D. A. Garanin and E. M. Chudnovsky, Polyhexatic and polycrystalline states of skyrmion lattices, Physical Review B 107, 014419-(17) (2023).
  * [8] P. Meisenheimer, H. Zhang, D. Raftrey, X. Chen, Y.-T. Shao, Y.-T. Chan, R. Yalisove, R. Chen. J. Yao, M. C. Scott, W. Wu, D. A. Muller, P. Fischer, R. J. Birgeneau, R. Ramesh, Ordering of room-temperature magnetic skyrmions in a polar van der Waals magnet, Nature Communications 14 3744-(9) (2023).
  * [9] J. M. Kosterlitz, Kosterlitz-Thouless physics: a review of key issues, Reports on Progress in Physics 79 026001-(60) (2016).
  * [10] V. N. Ryzhov, E. A. Gaidul, E. E. Tareeva, Y. D. Fomin, and E. N. Tsiok, Melting scenarios of two-dimensional systems: Possibilities of computer simulation, Journal of Theoretical and Experimental Physics (JETP) 137, 125-150 (2023).
  * [11] R. E. Peierls, Remarks on transition temperatures, Helvetica Physica Acta 7, 81-83 (1934).
  * [12] N. D. Mermin and H. Wagner H 1966 Absence of ferromagnetism or antiferromagnetism in one or two-dimensional isotropic Heisenberg models, Physical Review Letters 17 1133-1136 (1966).
  * [13] P. C. Hohenberg, Existence of long-range order in one and two dimensions Physical Review 158, 383-386 (1967).
  * [14] B. J. Adler and T. E. Wainwright, Phase transition in elastic disks, Physical Review 127, 369-361 (1962).
  * [15] R. S. Kagiwada, J. C. Fraser, I. Rudnick, and D. Bergman, Superflow in helium films, Physical Review Letters 22, 338-342 (1969).
  * [16] V. L. Berezinskii, Destruction of long range order in one-dimensional and two-dimensional systems having a continuous symmetry group. I. Classical systems, Journal of Experimental and Theoretical Physics (JETP) 32, 493-500 (1971).
  * [17] J. M. Kosterlitz and D. J. Thouless, Long range order and metastability in two-dimensional solids and superfluids, Journal of Physics C: Solid State Physics 5 L124-L126 (1972); Ordering, metastability and phase transitions in two-dimensional systems, Journal of Physics C: Solid State Physics 6 1181-1203 (1973).
  * [18] D. J. Bishop and J. D. Reppy, Study of the superfluid transition in two-dimensional 4He films Physical Review Letters 40 1727-1730 (1978); Physical Review B 22 5171-5185 (1980).
  * [19] B. I. Halperin and D. R. Nelson, Theory of two-dimensional melting, Physical Review Letters 41, 121-124 (1978).
  * [20] D. R. Nelson and B. I. Halperin, Dislocation-mediated melting in two dimensions, Physical Review B 19, 2457-2483 (1979).
  * [21] A. P. Young, Melting and the vector Coulomb gas in two dimensions, Physical Review B 19, 1855-1866 (1979).
  * [22] J. Q. Broughton, G. H. Gilmer, and J. D. Weeks, Molecular-dynamics study of melting in two dimensions: Inverse-twelfth-power interaction, Physical Review B 25, 4651-4669 (1982).
  * [23] F. Lindemann, About the Calculation of Molecular Own Frequencies, Zeitschrift für Physik 11, 14, 609-612 (1910).
  * [24] S. A. Khrapak, Lindenmann melting criterion in two dimensions, Physical Review Research 2, 012040(R)-(6) (2020).
  * [25] S. T. Chui, Grain-boundary theory of melting in two dimensions, Physical Review Letters 48, 933-935, (1982); Physical Review B 28, 178-194 (1983).
  * [26] Y. Saito, Monte Carlo studies of two-dimensional melting: Dislocation vector systems, Physical Review B 26, 6239-6253, (1982).
  * [27] D. S. Fisher, B. I. Halperin, and R. Morf, Defects in the two-dimensional electron solid and implications for melting, Physical Review B 20, 4692-4712 (1979).
  * [28] W. Janke and H. Kleinert, From first-order to two continuous melting transitions: Monte Carlo Study of a new 2D lattice-defect model, Physical Review Letters 61, 2344-2347 (1988); Erratum Physical Review Letters 62, 608 (1989).
  * [29] E.M. Chudnovsky, Structure of a solid film on an imperfect surface, Physical Review B 33, 245-250 (1986).
  * [30] E.M. Chudnovsky, Hexatic vortex glass in disordered superconductors, Physical Review B 40, (R)11355-11367 (1989).
  * [31] C. A. Murray, P. L. Gammel, D. J. Bishop, D. B. Mitzi, and A. Kapitulnik, Observation of a hexatic vortex glass in flux lattices of the high-$T_{c}$ superconductor Bi2.1Sr1.9Ca0.9Cu2O8+δ, Physical Review Letters 64, 2312-2315 (1990); D. G. Grier, C. A. Murray, C. A. Bolle, P. L. Gammel, D. J. Bishop, D. B. Mitzi, and A. Kapitulnik, Translational and bond-orientational order in the vortex lattice of the high-$T_{c}$ superconductor Bi2.1Sr1.9Ca0.9Cu2O8+δ, Physical Review Letters 66, 2270-2273 (1991).
  * [32] H. H. Günberg, P. Keim, K. Zahn, and G. Maret, Elastic behavior of a two-dimensional crystal near melting, Physical Review Letters 93 255703-(4) (2004).
  * [33] J. Zanghellini, P. Keim, and H. H. von Grünberg, The softening of two-dimensional colloidal crystals, Journal of Physics: Condensed Matter 17 S3579-S3586 (2005).
  * [34] J. Dietel and H. Kleinert, Triangular lattice model of two-dimensional defect melting, Physical Review B 73, 024113-(13) (2006).
  * [35] S. C. Kapfer and W. Krauth, Physical Review, Letters 114, 035702-(5) (2015).
  * [36] J. A. Anderson, J. Antonaglia, J. A. Millan, M. Engel, and S. Glotzer, Shape and symmetry determine two-dimensional melting transitions of hard regular polygons, Physical Review X 7, 021001-(14) (2017).
  * [37] E. N. Tsiok, Y. D. Fomin, E. A. Gaiduk, E. E. Tareyeva, V. N. Ryzhov, P. A. Libet, N. A. Dmitryuk, N. P. Kryuchkov, and S. O. Yurchenko, Journal of Chemical Physics 156, 114703-(11) (2022).
  * [38] A. A. Belavin and A. M. Polyakov, Metastable states of two-dimensional isotropic ferromagnets, Pis’ma Zh. Eksp. Teor. Fiz. 22 503–506 (1975) [JETP Lett. 22, 245–247 (1975)].
  * [39] N. Nagaosa and Y. Tokura, Topological properties and dynamics of magnetic skyrmions, Nature Nanotechnology 8, 899-911 (2013).
  * [40] X. Zhang, M. Ezawa, and Y. Zhou, Magnetic skyrmion logic gates: conversion, duplication and merging of skyrmions, Scientific Reports 5, 9400-(8) (2015).
  * [41] G. Finocchio, F. Büttner, R. Tomasello, M. Carpentieri, and M. Klaui, Magnetic skyrmions: from fundamental to applications, Journal of Physics D: Applied Physics. 49, 423001-(17) (2016).
  * [42] W. Jiang, G. Chen, K. Liu, J. Zang, S. G. E. te Velthuis, and A. Hoffmann, Skyrmions in magnetic multilayers, Physics Reports 704, 1-49 (2017).
  * [43] A. Fert, N. Reyren, and V. Cros, Magnetic skyrmions: advances in physics and potential applications, Nature Reviews Materials 2, 17031-(15) (2017).
  * [44] A. N. Bogdanov and D. A. Yablonskii, Thermodynamically stable “vortices” in magnetically ordered crystals. The mixed state of magnets. Soviet Physics JETP 68, 101-103 (1989).
  * [45] A. Bogdanov and A. Hubert, Thermodynamically stable magnetic vortex states in magnetic crystals, J. Magn. Magn. Mater. 138, 255-269 (1994).
  * [46] U. K. Rößler, N. Bogdanov, and C. Pfleiderer, Spontaneous skyrmion ground states in magnetic metals, Nature 442, 797-801 (2006).
  * [47] S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A. Kubetzka, R. Wiesendanger, G. Bihlmayer, and S. Blugel, Spontaneous atomic-scale magnetic skyrmion lattice in two dimensions, Nature Physics 7, 713-718 (2011).
  * [48] O. Boulle, J. Vogel, H. Yang, S. Pizzini, D. de Souza Chaves, A. Locatelli, T. O. Mentes, A. Sala, L. D. Buda-Prejbeanu, O. Klein, M. Belmeguenai, Y. Roussigné, A. Stahkevich, S. M. Chérif, L. Aballe, M. Foerster, M. Chshiev, S. Auffret, I. M. Miron, and G. Gaudin, Room-temperature chiral magnetic skyrmions in ultrathin magnetic nanostructures, Nature Nanotechnology 11, 449-454 (2016).
  * [49] A. O. Leonov, T. L. Monchesky, N. Romming, A. Kubetzka, A. N. Bogdanov, and R. Wiesendanger, The properties of isolated chiral skyrmions in thin magnetic films, New Journal of Physics 18, 065003-(16) (2016).
  * [50] A. O. Leonov and M. Mostovoy, Multiply periodic states and isolated skyrmions in an anisotropic frustrated magnet, Nature Communications 6, 8275-(8) (2015).
  * [51] X. Zhang, J. Xia, Y. Zhou, X. Liu, H. Zhang, and M. Ezawa, Skyrmion dynamics in a frustrated ferromagnetic film and current-induced helicity locking-unlocking transition, Nature Communications 8 1717-(10) (2017).
  * [52] B. A. Ivanov, A. Y. Merkulov, V. A. Stepanovich, C. E. Zaspel, Finite energy solitons in highly anisotropic two dimensional ferromagnets, Physical Review B 74, 224422-(17) (2006).
  * [53] S.-Z. Lin and S. Hayami, Ginzburg-Landau theory for skyrmions in inversion-symmetric magnets with competing interactions, Physical Review B 93, 064430-(16) (2016).
  * [54] E. M. Chudnovsky and D. A. Garanin, Skyrmion glass in a 2D Heisenberg ferromagnet with quenched disorder, New Journal of Physics 20 033006-(9) (2018).
  * [55] C. Moutafis, S. Komineas, and J. A. C. Bland, Dynamics and switching processes for magnetic bubbles in nanoelements, Physical Review B 79, 224429-(8) (2009).
  * [56] X. Yu, Y. Onose, N. Kanazawa, J. Park, J. Han, Y. Matsui, N. Nagaosa, and Y. Tokura, Real-space observation of a two-dimensional skyrmion crystal, Nature 465, 901-904 (2010).
  * [57] D. A. Garanin, E. M. Chudnovsky, S. Zhang, and X. X. Zhang, Thermal creation of skyrmions in ferromagnetic films with perpendicular anisotropy and Dzyaloshinskii-Moriya interaction, Journal of Magnetism and Magnetic Materials 493, 165724-(9).
  * [58] T. Dohi, R. M. Reeve, and M. Kläui, Thin film skyrmionics, Annual Review of Condensed Matter Physics 33, 73-95 (2022).
  * [59] D. Capic, D. A. Garanin, and E. M. Chudnovsky, Skyrmion-skyrmion interaction in a magnetic film, Journal of Physics: Condensed Matter 32, 415803-(14) (2020).
  * [60] E. M. Chudnovsky and J. Tejada, Lectures on Magnetsm, Rinton Press, Princeton, New Jersey (2006).
  * [61] L. Cai, E. M. Chudnovsky, and D. A. Garanin, Collapse of skyrmions in two-dimensional ferromagnets and antiferromagnets, Physical Review B 86, 024429-(4) (2012).
  * [62] D. Capic, D. A. Garanin, and E. M. Chudnovsky, Stabilty of biskyrmions in centrosymmetric magnetic films, Physical Review B 100, 014432 (2019).
  * [63] R. E. Camley and K. L. Livedey, Consequences of Dzyaloshinskii-Motiya interaction, Surface Science Reports 78, 100605-(33) (2023).
  * [64] D. A. Garanin and E. M. Chudnovsky, Solid-liquid transition in a skyrmion matter, arXiv:2403.00151.
  * [65] A. Derras-Chouk, E. M. Chudnovsky, and D. A. Garanin, Quantum collapse of a magnetic skyrmion, Physical Review B 98, 024423-(9) (2018).