# Magnetic ordering and dynamics in monolayers and bilayers of chromium
trihalides: atomistic simulations approach

S. Stagraczyński  Faculty of Physics, Adam Mickiewicz University in Poznań,
ul. Uniwersytetu Poznańskiego 2, 61-614 Poznań, Poland  P. Baláž  FZU –
Institute of Physics of the Czech Academy of Sciences, Na Slovance 1999/2, 182
21 Prague 8, Czech Republic  M. Jafari  Faculty of Physics, Adam Mickiewicz
University in Poznań, ul. Uniwersytetu Poznańskiego 2, 61-614 Poznań, Poland
J. Barnaś  Faculty of Physics, Adam Mickiewicz University in Poznań, ul.
Uniwersytetu Poznańskiego 2, 61-614 Poznań, Poland  A. Dyrdał

###### Abstract

We analyze magnetic properties of monolayers and bilayers of chromium
trihalides, CrI3, in two different stacking configurations: AA and
rhombohedral ones. Our main focus is on the corresponding Curie temperatures,
hysteresis curves, equilibrium spin structures, and spin wave excitations. To
obtain all these magnetic characteristic, we employ the atomistic spin
dynamics and Monte Carlo simulation techniques. The model Hamiltonian includes
isotropic exchange coupling, magnetic anisotropy, and Dzyaloshinskii-Moriya
interaction. Though the latter is relatively weak in CrI3, we consider a more
general case assuming also an enhancement of Dzyaloshinskii-Moriya interaction
in the corresponding Janus structures and by external electric fields. An
important issue of the analysis is the correlation between hysteresis curves
and spin configurations in the system, as well as formation of the skyrmion
textures.

## Introduction

Van der Waals magnetic (vdW) materials are currently of great interest due to
expected applications in two-dimensional (2D) nanoelectronics and spintronics.
These materials are built from weakly coupled 2D monolayers, and therefore
they can be relatively easily obtained in the multilayered form [1], down to a
single monolayer. The discovery of ferromagnetism in monolayers of CrI3 [2]
and Cr2Ge2Te6 [3] initiated an enormous interest in magnetic 2D vdW crystals
[4, 5, 6]. There are several groups of well-known magnetic vdW materials,
including chromium trihalides CrX3 (X= I, Cl, Br, I) [2, 7], metal tribromides
MBr3 (M=Mn, Cu, Fe, V), chromium-based ternary tellurides Cr2XTe6 (X = Ge, P)
[3], transition metal dichalcogenides MnX2 and VX2 (X = S, Se, Te) in two
different phases denoted as 2H and 1T, [8, 9, 10], and also others. The number
of known 2D magnetic van der Waals materials is growing rapidly, and they have
various structural, electronic, magnetic, and transport properties. There is
currently extensive search for new materials, with new properties, and new
perspectives for applications, as these materials seem to be ideal for
building ultra-thin (of atomic scale) nanoelectronic devices (like spin valves
or nonvolatile memory elements) for, e.g., quantum computing, spintronics,
magnonics, opto-electronics, and others[4, 11]. However, magnetic properties
of most of the known 2D van der Waals systems still require further
investigations. The magnetic ground state of 2D vdW crystals depends, among
others, on their crystallographic phase, stacking geometry, substrate, and on
possible twisting of adjacent monolayers. In many cases, the magnetic ground
state survives only at low temperatures, while at higher temperatures the
systems become paramagnetic. Therefore, much efforts are focused on searching
for materials with critical temperatures above the room temperature.

Another group of 2D materials are layered molecular magnets. Indeed, a number
of molecular magnets have been synthesized in the monolayer and bilayer forms
[12, 13, 14]. Both 2D vdW materials and 2D molecular magnets offer a unique
opportunity to test various theoretical models of magnetic interactions and
magnetic physical phenomena in two dimensions. As for the interactions, this
includes exchange symmetric coupling, Dzyaloshinskii-Moriya Interaction (DMI),
magnetic anisotropy, higher-order interactions, like two-ion magnetic
anisotropy (or Kitaeev term), and also others [7, 15]. Many of these materials
have the ground state spin orientation perpendicular to the monolayers or
bilayers. In some cases, the magnetic anisotropy enforces in-plane spin
orientation. Moreover, these materials also allow to observe various
topological phases and phase transitions, like strictly 2D skyrmions or
topological Berezinsky-Kosterlitz-Thouless phase transitions [16, 12].

Among various groups of vdW materials, there are materials which are promising
for applications, like for instance transition metal dichalcogenides,
especially Vanadium dichalcogenides VX2 (X=S, Se, Te) [8, 9, 10], and chromium
trihalides, e.g., chromium iodine and chromium chlorine, CrI3 and CrCl3 [2, 7,
15]. In this paper we will focus on the latter group, specifically on CrI3.
Monolayers of this material are ferromagnetic, with the magnetic moments
oriented perpendicularly to the layer plane [2, 17]. This material is very
interesting for fundamental research [18, 19]. However, the corresponding
Curie temperature is around 40 K [2], which is too low for most of practical
applications. In the bilayer form, the two monolayers are coupled
antiferromagnetically. However magnetic properties can be externally tuned by
electric field [20, 21] or strains [22, 23, 24, 25, 26, 27, 28], and
interlayer coupling can be even changed from antiferromagnetic to
ferromagnetic one [25]. Moreover, magnetic properties of CrI3 in a multilayer
form depend on the stacking geometry [29, 30, 31]. This material seems to be
ideal for testing various theoretical models of spin configuration, including
topological skyrmion textures [32, 33, 34, 35], as well as excited states
(spin waves) which display certain topological features [7, 36, 37, 38, 39,
40]. It was pointed out [35], that DMI in CrI3 is too small to generate
skyrmion textures. However, it was also predicted that DMI in Janus monolayers
of chromium trihalides Cr(I,X)3 can be significantly larger than that in CrI3,
and can generate magnetic skyrmion textures [33]. Therefore, we assume a more
general case, which admits larger values of DMI (e.g. in Janus structures).

In this paper we analyze magnetic properties of CrI3, and to do this we employ
simulation methods like atomistic spin dynamic (ASD) technique and Monte Carlo
simulations [41, 42, 43]. From this we determine such magnetic properties like
Curie temperature, hysteresis curves, (quasi)equilibrium spin configurations,
or excited states (spin waves). We limit considerations to monolayers and
bilayers of CrI3. As some results on the monolayers are already available in
the literature, much less work has been done on the bilayer structures of
CrI3. Atomic structure of the bilayers of CrI3 is shown in Fig.1 for two
different stacking of the monolayers. This figure shows top and side views.
Atomic structure of a single monolayer is effectively shown by the top part of
Fig.1(a).

Figure 1: Atomic structure of the bilayers of CrI3 in the AA stacking of
layers without offset (a), and in the rhombohedral stacking, where the offset
is half of the unit cell vectors $\vec{a}/2$, $\vec{b}/2$ (b). The
$z$-component of unit vectors is kept the same for both cases. The Cr atoms of
the top (bottom) monolayer are indicated in blue (green). The top panels
present the top view, while the bottom panels show the side view. Note, the
top part of (a) shows effectively atomic structure of a monolayer. The I atoms
are indicated in the top parts in pink, while for clarity reasons they are not
indicated in the bottom parts. Intralayer and interlayer exchange couplings
are also marked in the top (bottom) parts.

##  1 Model and Method

###  1.1 Structure and spin Hamiltonian

Magnetic properties of chromium trihalides remarkably depend on the number of
monolayers in the sample, and also on the geometry and offset of their
stacking. This difference is especially significant for small numbers of
layers in the system. The most stable stacking geometries are the AA and
Rhombohedral ones, see Fig.1. The following analysis is limited to both
monolayers and bilayers of CrI3. Magnetic moments in CrI3 are predominantly
localized on the chromium atoms, while contributions from iodine atoms are
rather small and may be omitted. However, the iodine atoms contribute to
exchange interactions between the magnetic Cr atoms (e.g. via superexchange),
and thus have a significant impact on magnetic properties of CrI3. The Cr
atoms in an individual monolayer form a hexagonal (honeycomb) lattice, and
positions of these atoms have been taken from the unit cell of bulk CrI3.

Accordingly, magnetic properties of CrI3 will be described by an effective
spin Hamiltonian attributed to the corresponding hexagonal lattice of Cr
atoms, which we split into two terms,

$\mathcal{H}=\mathcal{H}_{0}+\mathcal{H}_{\rm DM},$ (1)

where $\mathcal{H}_{0}$ includes the exchange interactions, magnetic
anisotropies, and the Zeemann energy in an external magnetic field
$\mathbf{B}$. The second term, $\mathcal{H}_{\rm DM}$, takes into account
antisymmetric exchange, known also as Dzyaloshinskii-Moriya Interaction (DMI).

Assuming isotropic exchange interactions, we write the first term in the form

$\displaystyle\mathcal{H}_{0}=-\sum\limits_{ij}J_{ij}{\mathbf{S}}_{i}\cdot{\mathbf{S}}_{j}-k\sum\limits_{i}\left(S^{z}_{i}\right)^{2}-\mu_{S}\sum\limits_{i}{\mathbf{B}}\cdot{\mathbf{S}}_{i},$
(2)

where the sum over $i$ and $j$ means the summation over lattice sites. Here,
the positive exchange parameters $J_{ij}$ correspond to ferromagnetic (FM)
coupling, whereas negative ones correspond to antiferromagnetic (AFM)
coupling. In the following we take nonzero exchange parameters between
nearest-neighbours, second-nearest and third-nearest neighbors. In a more
general case the exchange parameters can be anisotropic. The second term
describes the magnetic anisotropy. Here, positive constant $k$ corresponds to
the easy-axis anisotropy along the $z$-direction (normal to the layer), while
negative $k$ corresponds to the easy-($xy$)plane anisotropy. Furthermore,
$\mu_{s}$ is the atomic magnetic moment of Cr atoms. Here, and in the
following, if not stated otherwise, the vector $\mathbf{S}$ is a unit vector
along the orientation of a local spin moment [41]. Accordingly, one needs to
adjust the exchange parameters, which usually are related to spin Hamiltonian,
where $\mathbf{S}$ is the spin operator.

The second term in Eq.(1) includes the Dzyaloshinskii-Moriya interactions (or
alternatively antisymmetric exchange). Existence of nonzero components of this
interaction depends on symmetry of the system. In van der Waals materials one
can easily tune externally the symmetry and thus also the components of DMI,
e.g., by an external mechanical strain or by electric field normal to the
layers. In the case under consideration, the DMI Hamiltonian can be written in
the form [15]

$\mathcal{H}_{\rm DM}=-\sum\limits_{\left\langle
i,j\right\rangle}{\mathbf{d}}^{NN}_{ij}\cdot\left(\mathbf{S}_{i}\times\mathbf{S}_{j}\right)-\sum\limits_{\left\langle\left\langle
i,j\right\rangle\right\rangle}\mathbf{d}^{NNN}_{ij}\cdot\left(\mathbf{S}_{i}\times\mathbf{S}_{j}\right),$
(3)

where the first term includes interaction between intralayer nearest-
neighbours (NN), while the second term between next nearest-neighbours (NNN).
Here, ${\mathbf{d}}^{NN}_{ij}$ and ${\mathbf{d}}^{NNN}_{ij}$ are the
corresponding Dzyaloshinskii-Moriya vectors [15]. This vector for NN has the
following form:
$\mathbf{d}^{NN}_{ij}=d^{NN}_{xy}\hat{\mathrm{z}}\times\hat{\mathrm{l}}_{ij}+\tau_{ij}^{NN}d^{NN}_{z}\hat{\mathrm{z}}$,
where $\hat{\mathrm{z}}$ is a unit vector along the axis $z$ (normal to the
layer), $\hat{\mathrm{l}}_{ij}$ is the unit vector oriented from the site $i$
to site $j$, while $d^{NN}_{xy}$ and $d^{NN}_{z}$ are parameters and
$\tau_{ij}^{NN}=1(-1)$ for the AB (BA) link. In turn, in the second term
$\mathbf{d}^{NNN}_{ij}=\tau^{NNN}_{ij}\left(d_{xy}^{NNN}\hat{\mathrm{l}}_{ij}+\nu_{ij}d_{z}^{NNN}\hat{\mathrm{z}}\right)$,
with $\tau^{NNN}_{ij}=+1$ for AA-link and $-1$ for BB-link, while
$\nu^{NNN}_{ij}=\pm 1$ and alternates its sign on the consecutive NNN links of
a chosen site [15].

The DMI between nearest-neighbors in pristine CrI3 vanishes, while the DMI
between the next-nearest-neighbors is rather small,
$d_{z}^{NNN}=-8.8\mu\mathrm{eV}$ and $d_{xy}^{NNN}=73.0\mu\mathrm{eV}$ [15].
However, DMI can be tuned by an external electric field normal to the layer,
so the parameters $d_{xy}$ and $d_{z}$ can be nonzero and relatively large
[15]. It can be also tuned by external strain. Moreover, the DMI can be
significantly enhanced in the corresponding Janus monolayers of chromium
trihalides Cr(I,X)3 [33]. All this makes the DMI parameters between NNs and
NNNs tunable in a quite considerable range. Finally,we note, that in the
Hamiltonian assumed here we neglect such terms, like anisotropic exchange and
higher-order anisotropy terms (e.g., two-ion anisotropy or Kitaev term).

###  1.2 Numerical procedures

Our numerical calculations are based on the atomistic spin dynamics (ASD)
method implemented into the Vampire code [41]. This method allows for taking
into account all details of the system geometry, and also for keeping insight
into behaviour of the spin moments during the simulation process. The ASD
simulations are carried out by solving the coupled system of atomistic
stochastic (at finite temperature) Landau–Lifshitz-Gilbert (LLG) equations,

$\frac{\partial\mathbf{S}_{i}}{\partial
t}=-\frac{\gamma}{1+\alpha^{2}}\left[\mathbf{S}_{i}\times\mathbf{B}^{\mathrm{eff}}_{i}+\alpha\mathbf{S}_{i}\times(\mathbf{S}_{i}\times\mathbf{B}^{\mathrm{eff}}_{i})\right],$
(4)

where $\gamma$ is the gyromagnetic ratio, $\alpha$ stands for the Gilbert
damping parameter, and $\mathbf{B}^{\mathrm{eff}}_{i}$ is the effective
magnetic field acting on the spin at site $i$. This effective field
$\mathbf{B}^{\mathrm{eff}}_{i}$ consists of two terms: one is deterministic
and is related to the spin Hamiltonian, while the second one is a stochastic
term, that is determined by thermal fluctuations,

$\mathbf{B}^{\mathrm{eff}}_{i}=\frac{-1}{M_{s}}\frac{\partial\mathscr{H}}{\partial\mathbf{S}_{i}}+\mathbf{\xi}^{\rm(th)}_{i}.$
(5)

In the low energy regime, the stochastic term, $\mathbf{\xi}^{\rm(th)}_{i}$ is
usually modelled by the Gaussian thermal noise.

To verify the results on spin dynamics obtained by the atomistic LLG
equations, we also used other methods, like the method based on the spin-spin
correlation function and the Monte Carlo simulation technique. We find a very
good agreement between the results obtained by the atomistic LLG equations and
those obtained by the Monte Carlo simulations. As the former technique is more
appropriate for dynamical properties, the latter one is more convenient for
static properties. However, the static properties are also accessible from
spin dynamics, but require much higher computational efforts. In addition, to
study dynamical properties we also use the atomistic approach implemented
within the Uppsala ASD code [44].

An important step in the simulations process, especially when simulating
hysteresis loops, critical temperatures, or magnon dispersion relations, is
the preparation of the relevant initial spin configuration. To do this we use
the field-cooling procedure. We start this procedure always with the random
spin state at a high temperature (above the Curie temperature). After
performing one million of Monte Carlo equilibration steps (or time steps for
LLG equation), the temperature becomes decreased nearly adiabatically towards
the minimum temperature (i.e. temperature, at which the required simulations
are to be performed, e.g. 15 K in most of our simulations of the hysteresis
loops).

In the following we will apply the above mentioned methods to determine the
Curie (Néel) temperature, as well as hysteresis loops at a finite (nonzero)
temperature. The simulations will be performed for monolayers, as well as for
bilayers in two different stacking geometries, i.e., in the AA and rhmbohedral
configurations. In addition we will also consider dynamical properties of the
systems under consideration, especially spin wave excitations.

##  2 Results

###  2.1 Critical temperatures

The reported experimental data indicate on a ferromagnetic ground state in
CrI3 monolayer [2]. Moreover, FM ordering is also expected in multilayers of
chromium trihalides, (i.e. in corresponding bulk material). However, the
bilayer of CrI3 was reported to display antiferromagnetic interlayer exchange
coupling [2], so such bilayers should also reveal not only features typical of
ferromagnets, but also some properties typical of antiferromagnets. The
interlayer coupling can be tuned (including also sign change), e.g., by an
external gate voltage [20, 21] or strain[23, 24, 25, 26, 27, 28]. Therefore,
one of our main objectives in the case of bilayers is to investigate the
impact of the layer stacking on the critical temperature. The analyzed
geometries (AA and rhombohedral) are visualized in Fig. 1. Second, in
antiferromagnetically coupled bilayers one may expect certain features related
to the breakdown of ferromagnetic order within the monolayers and of
antiferromagnetic order related to the interlayer coupling (Curie and Néel
temperatures.

a)b) Figure 2: (a) Simulated (dots) magnetization vs. temperature for the
monolayer and ferromagnetic bilayers in the AA and Rh-stacking. The Curie
temperatures have been determined from fitting to the formula
$(1-T/T_{C})^{\beta}$ (solid lines). From this fitting follows: monolayer:
$T_{C}=33.93$K, $\beta=0.32$; FM-coupled bilayer in Rhombohedral stacking:
$T_{C}=37.64$K, $\beta=0.31$; FM-coupled bilayer in the AA stacking:
$T_{C}=39.64$K, $\beta=0.32$. (b) Simulated mean specific heat vs. temperature
for the antiferromagnetic bilayers in the AA and Rh-stacking.. We added here
also the data for the ferromagnetic monolayer, to show that the Curie
temperatures determined from simulations of the specific heat and from
simulations of the magnetization are equal. The parameters assumed in
simulations: The exchange couplings (after Ref. 7) $J_{1}=4.020$meV,
$J_{2}=0.643$meV, $J_{3}=0.016$meV), while for the interlayer exchange in
bilayers we assumed $J_{z}=J_{2}$ (a) and $J_{z}=-J_{2}$ (b). The nearest-
neigbour DMI, $d^{NN}_{xy}=0.201$meV and $d^{NN}_{z}=0$, the magnetic
anisotropy constant $k=0.109$ meV, and the chromium magnetic moment
$\mu_{Cr}=3.36\mu_{B}$. The parameters of the NNN DMI as in the text.

The simulated data (points) for both monolayer and ferromagnetic bilayer
(ferromagnetic interlayer coupling) in the AA and Rh stacking geometries are
shown in Fig. 2(a). We included there the exchange parameters between nearest
neighbours, $J_{1}$, second-nearest neighbours, $J_{2}$, third-nearest
neighbours, $J_{3}$, and interlayer exchange coupling between nearest
neighbours, $J_{z}$, as described in the figure caption. The other parameters
are also described there. To determine the Curie temperatures, the
corresponding data have been fitted to the formula $(1-T/T_{C})^{\beta}$, and
from this fitting we found $T_{C}\approx 33.93$ K for the monolayer,
$T_{C}\approx 37.64$ K for the Rh-stacked bilayer, and $T_{C}\approx 39.64$ K
for the AA-stacked bilayer. The Curie temperature for the monolayer is lower
than that for the bilayers, and the bilayer in the AA-stacking has Curie
temperature higher than that for the Rh-stacking. For the assumed parameters,
the DMI has no significant impact on the magnetization vs. temperature profile
and on Curie temperature.

The critical temperatures (both Curie and Néel) also can be determined from
magnetic susceptibility or specific heat. In the case of ferromagnetic
monolayer and ferromagnetically coupled bilayers, both these methods give
results which are consistent with those obtained from the above magnetization
analysis. Therefore, to determine the critical temperature in bilayers with
antiferromagnetic interlayer coupling, we used the Vampire package [43, 42] to
simulate the specific heat as a function of temperature, $C_{v}=[\left\langle
U^{2}\right\rangle-\left\langle U\right\rangle^{2}]/k_{B}T^{2}$, where $U$ is
the internal energy. The specific heat has a cusp at the critical temperature,
and this property is used to determine the Curie temperature, $T_{c}$, in
ferromagnetic systems, and the Néel temperature in antiferromagnets.

The numerical results on the specific heat for bilayers with antiferromagnetic
coupling in the AA and Rh stackings and also for the ferromagnetic monolayer
are shown in Fig. 2(b). The critical temperatures for the bilayers with
antiferromagnetic interlayer coupling are clearly seen and are roughly the
same as Curie temperatures for the corresponding bilayers with ferromagnetic
interlayer coupling. This is because absolute value of interlayer coupling is
remarkably smaller than the intralayer exchange coupling constant $J_{1}$, and
absolute magnitude of the interlayer coupling is assumed the same in bilayers
with ferro- and antiferromagnetic interlayer coupling. The curve for
ferromagnetic monolayer is presented in Fig. 2(b) just to show that the
methods used in simulations of the results in Fig. 2(a) and Fig. 2(b) give the
same critical temperatures.

When the interlayer exchange coupling is antiferromagnetic, the total magnetic
moment of the bilayer vanishes in the ground state. As already mentioned
above, the absolute magnitude of the interlayer coupling is remarkably smaller
than the intralayer exchange coupling $J_{1}$. Thus, with increasing
temperature, the antiferromagnetic moment (Néel vector) vanishes first and we
have a system of two ferromagnetic layers (with average moments either aligned
or antialigned). When the temperature grows further the magnetic moment in
each monolayer (and thus in the whole system) drops to zero at the Curie
temperature as determined above. Thus, on the simulated magnetization vs
temperature curves one could expect two cusps, one at a lower temperature that
corresponds to the breakdown of the antiferromagnetic alignment of the two
monolayers, and another one at higher Curie temperature. However, in the
simulation data on the magnetization vs temperature, there is no signature of
the breakdown of antiferromagnetic configuration. Apparently, it would be
resolved in a multilayer system [45].

Figure 3: Simulated normalized values of $S^{z}$ for both monolayers in the
AFM-coupled Rh-stacked bilayer, presented as a function of temperature. These
normalized values of $S^{z}$ for the top and bottom layers are shown as the
yellow and blue dots, respectively. The solid red and blue lines indicate the
average of $S^{z}$ for both layers. The grey area between these lines
indicates the temperature range where reorientation of the magnetic ordering
between the layers does not happen in the absence of external magnetic field.
Reorientation begins at the border between the grey and green areas, and the
temperature at which this takes place (around 34 K) we assign as the Néel
temperature $T_{N}$. When an external magnetic field normal to the bilayer is
applied, the reorientations start at lower temperature, as shown in the inset,
and also marked schematically on the grey area.

In order to get some information on the temperature at which the
antiferromagnetic alignment of the two ferromagnetic monolayers is
destabilized by thermal fluctuations, we analyzed average values of the total
spin $z$-component of each monolayer separately as a function of temperature.
Results of the simulations are shown in Fig. 3, where the simulated normalised
values of $S^{z}$ for the two monolayers are presented with blue and yellow
dots, while the average values are shown by the solid blue and red lines. With
increasing temperature, fluctuations around the average values become larger,
and at a certain temperature (recognized as the relevant Néel temperature) the
average spin moments of the two monolayers start to change sign. The Néel
temperature was identified as a temperature at which spin moments of the
monolayers begin to fluctuate between the up and down orientations, see Fig.
3. The corresponding Néel temperature is lower than the Curie temperature by
about 3.5 K, as one might expect.

###  2.2 Hysteresis loops and spin textures

Hysteresis curves are important characteristics as they include information on
magnetization variation with increasing external magnetic field and thus also
on possible field-induced phase transitions. We have determined the hysteresis
loops from the Monte-Carlo (MC) simulations as well as from the atomistic
Landau-Lifshitz-Gilbert (LLG) equation. The starting point for each simulation
was a random state at temperature $T=80\mathrm{K}$ and cooling field
$H_{\mathrm{fc}}=1\mathrm{T}$ applied perpendicularly to the sample plane. The
field-cooling procedure was used to cool the system’s temperature down to
$T=15\mathrm{K}$, and then the simulations of a hysteresis loop was initiated.

####  2.2.1 Monolayers

Figure 4:  Hysteresis loops for monolayers at $T=15$K, easy-axis (normal to
the layer) anisotropy parameter $k=0.109$meV, and with the NN DMI parameters:
(a) $d_{xy}=0$ and $d_{z}/J_{1}$ as indicated; (b) $d_{z}=0$ and
$d_{xy}/J_{1}$ indicated; and (c) ($d_{xy}/J_{1}$, $d_{z}/J_{1}$) as
indicated. The other parameters as in Fig. 2. The presented results are
average over 50 realizations.

Figure 4 shows the hysteresis curves for fixed values of DMI between NNNs (as
described above in the text), and different values of the components $d_{z}$
and $d_{xy}$ of the DM vector between NNs. When $d_{xy}=0$ , then increasing
the component $d_{z}$ leads mainly to reduced coercivity of the hysteresis
curves, while their shape remains similar to those for ferromagnetic layers,
see Figure 4(a). A more interesting impact of DMI on the hysteresis curves
appears for the $d_{xy}$ component of the DMI vector. The corresponding
hysteresis curves for monolayers are shown in Figure 4(b) for different values
of $d_{xy}$ and $d_{z}=0$. For small values of $d_{xy}/J_{1}$, the hysteresis
curves are qualitatively similar to those for zero $d_{xy}$, but the
corresponding coercivity becomes reduced with increasing $d_{xy}$. When
$d_{xy}$ grows further, some new features appear in the curves, see e.g. the
curves for $d_{xy}/J_{1}=0.3$ and $d_{xy}/J_{1}=0.35$, which clearly
distinguish these loops from those typical for ferromagnetic systems. The
interplay of both components of NN DMI vector is presented in Figure 4(c),
where both components are nonzero and comparable. For the assumed parameters
the corresponding hysteresis loops contain clearly features of both $d_{z}$
and $d_{xy}$ components of the DMI vector. From the hysteresis curves shown
above one may conclude that DMI modifies magnetic texture of the system and
this modification leads to the observed features in the hysteresis loops. In
the curves shown here, the most pronounced features are for
$d_{xy}=0.35J_{1}$. Therefore we analyse now in more details this specific
situation ($d_{z}=0$ and $d_{xy}=0.35J_{1}$), and for clarity, we also omit
here the NNN DMI.

Figure 5: (a) Hysteresis loop for the monolayer with $d_{xy}=0.35J_{1}$,
$d_{z}=0$, $T=15$K, and vanishing NNN DMI. The solid black line is the average
over 50 realizations, while the shaded area represents the standard deviation.
The circle points represent realizations that are closest to the average
(solid line), for which the corresponding magnetic textures are shown in part
(b) of this figure. (b) The magnetic textures corresponding to the points A,
B, C, D, E and F, marked in the hysteresis loop (a). The results of MC
simulations are shown in the top row. Then, the final set of MC simulations
was used as initial step for simulations based on the atomistic LLG equation,
and the middle row shows the magnetic textures after $1$ns of spin-dynamics
using the LLG Heun method with damping factor $\lambda=0.1$. Then the results
from the middle row were used in the next simulation for cooling the
temperature from $15$K down to $0$K. The spin textures upon cooling down to 0K
are shown in the bottom row. Upon cooling, the magnetic patterns (stripes,
skyrmions) have regular edges. (c) Magnified single skyrmion obtained upon
cooling the magnetic pattern corresponding to the point F ($B=0.68\mathrm{T}$)
in (b). Here, the spin structure of the skyrmion is clearly resolved. This
spin structure also confirms that the created skyrmions are of Néel-type, as
presented by the vertical crossection through the skyrmion center, shown on
the right side.

In Figure 5(a) we show the average hysteresis curve (solid line) and the
standard deviation from the average (shaded areas) for $d_{xy}=0.35J_{1}$,
$d_{z}=0$ and 15 K. The corresponding spin patterns for the points A, B, C, D,
E and F, indicated on the hysteresis curve are shown in Figure 5(b). Here, the
results of MC simulations are shown in the top row. Then, the final set of MC
simulations was used as initial step for simulations based on the atomistic
LLG equation. The middle row shows the magnetic textures after $1$ns of spin-
dynamics using the LLG Heun method with damping factor $\lambda=0.1$. We have
changed the numerical integrator from Monte Carlo to LLG in order to check if
the spin texture was stable under both integrators. The results from the
middle row were subsequently used in the next simulation for the cooling
temperature from $15$K down to $0$K. As the patterns for $T=15$K have spin
textures irregular at the boundaries due to thermal fluctuations, the edges
become regular when cooling the system temperature down to 0 K, see bottom row
in Figure 5(b). As one can see, upon cooling the edges of various objects
become smooth and one can clearly distinguish there regular circular objects
and stripe structures. The circular objects are Néel-type skyrmions. To see
this in more details, in Figure 5(c) we zoomed out an individual circular
object, so now the spin structure associated with the skyrmion is clearly
visible. To show further details, on the right side of this figure we present
a crossection through the center of the skyrmion .

####  2.2.2 Bilayers

Hysteresis loops in bilayers depend of the nature of interlayer coupling, and
for ferromagnetic and antiferomagnetic couplings are usually remarkably
different. In addition, they also depend on the stacking geometry. The
hysteresis curves for the AA and Rhombohedral stackings in the absence of DM
interactions and for isotropic exchange interactions are shown in Fig.6(a) for
systems with ferromagnetic interlayer coupling and in Fig.6(b). for systems
with antiferromagnetic interlayer coupling.

Figure 6: Hysteresis loops for bilayers in the AA and Rh stacking with (a)
ferromagnetic interlayer coupling $J_{z}=J_{2}$, and (b) antiferromagnetic
interlayer coupling, $J_{z}=-J_{2}$. The out-of-plane anisotropy constant
$k=0.109$ meV, and temperature $T=15$K. NN and NNN terms of DMI are omitted
here.  Figure 7: (a) Hysteresis loop of a ferromagnetically coupled bilayer in
the Rh-stacking, calculated for $d_{xy}=0.35J_{1}$, $T=15$K, ferromagnetic
interlayer coupling $J_{z}=J_{2}$, and for constant NNN DMI as described in
the main text. The solid black line is the average over 50 realizations, while
the shaded area represents the standard deviation. The numbers 1 to 4 indicate
four exemplary hysteresis loops obtained in the numerical simulations. (b)
Spin patterns in FM-coupled Rh-stacked bilayers upon cooling, obtained for the
points indicated on the hysteresis loop in (a). Only spin pattern of the top
layer is shown here, as that of the bottom layer is approximately the same.
One can distinguish here stripes, skyrmions, and other patterns with regular
edges.

From this figure follows that the hysteresis curves for ferromagnetically
coupled bilayers are similar to typical hysteresis loops of ferromagnetic
layers, and are roughly independent of the stacking geometry. On the other
hand, the hysteresis curves for bilayers with AFM interlayer coupling are
significantly different from those for ferromagnetically coupled bilayers, and
also different for different stacking geometries. An important feature of the
hysteresis loops for antiferromagnetic interlayer coupling is a wide plateau
at zero magnetization, which corresponds to antialigned magnetizations of the
two monolayers. At a certain value of magnetic field there is a jump in
magnetization associated with transition to the spin-flop phase (note, the
magnetic field is normal to the layers). Then the system goes smoothly to the
saturated state. The magnetic field, at which the transition to the spin-flop
phase occurs with increasing magnetic field is different from the field at
which the system goes back to antialigned configuration when the magnetic
field decreases. This leads to minor loops clearly visible in Fig.6(b) on both
sides (positive and negative magnetic fields). These loops follow from the
magnetic anisotropy in the system, and also depend on the stacking geometry.
No such features appear in ferromagnetically coupled bilayers. The saturation
field also depends on the stacking configuration. This is because the
saturation field depends on the effective coupling between the two monolayers,
and this coupling depends, in turn, on the strengths of interlayer exchange
parameter, and on the number of interlayer nearest neighbours. Note, the
hysteresis curves for the bilayers in the AA and Rh stacking have been
calculated for the same interlayer exchange parameters. However, the numbers
of nearest neighbors in the Rh stacking is larger (but distance between them
is longer then in the AA stacking). Thus, for the assumed parameters the
saturation fields are different for both stacking geometries and is larger for
the Rh stacking.

The impact of both NN and NNN terms of DMI on the hysteresis curves in the Rh
stacking geometry and for ferromagnetic and antiferromagnetic interlayer
coupling is shown in Fig.7(a) and Fig.8(a), respectively, while the
corresponding spin textures are shown in Fig.7(b) and Fig.8(b). As already
discussed above, the DMI between NNNs is rather small in pristine CrI3. In
turn, the DMI between NN’s can be easily tuned by external gate and therefore
can be dominant (even if it vanishes in pristine CrI3), so in Fig.7(a) and
Fig.8(a) we assumed intrinsic values of NNN DMI (as given in the text) and
quite significant value of DMI between NNs, $d_{xy}=0.35J_{1}$ and $d_{z}=0$.
We focus on the impact of the component $d_{xy}$ of DMI as this component is
crucial in the skyrmion texture formation. As in monolayers, the DMI modifies
the hysteresis loops and also the corresponding spin patterns. These patterns
include areas of various shapes, with irregular edges due to thermal
fluctuations at finite temperatures. However, when the temperature is cooled
down to T=0K, these edges become regular and reveal among others stripe
structures and well defined skyrmion textures.

Figure 8: (a) Hysteresis loop for antiferromagnetically coupled Rh-stacked
bilayer, simulated for $d_{xy}=0.35J_{1}$, $d_{z}=0$, $T=15$K,
antiferromagnetic interlayer coupling $J_{z}=-J_{2}$, and NNN DMI as described
in the text. The solid black line is the average over 50 realizations, with
the shadow represents the standard deviation. (b) The spin patterns for three
values of the magnetic field, including the case with zero magnetic field.
These magnetic fields correspond to the three points indicated in the
hysteresis loop (a).

Consider first the bilayers with ferromagnetic interlayer coupling. The
corresponding spin patterns are shown in Fig.7(b) for a few points indicated
on the corresponding hysteresis curve, see Fig.7(a). In general, we find there
stripes, skyrmions, and also other structures. Only spin patterns of the top
monolayer are shown there as the spin patterns of the bottom monolayer are
practically the same due to ferromagnetic interlayer coupling. Certain, almost
not resolved difference may appear due to a shift of the two monolayers in the
Rh stacking. Accordingly, the skyrmions in both layers are coupled
ferromagnetically and have the same topological charge.

In bilayers with antiferromagnetic coupling, the corresponding spin patterns
are qualitatively similar to those in FM-coupled bilayers, and we find the
spin stripes, wavy stripes, and skyrmions, as shown in Fig.8(b). Though they
look qualitatively similar to those in FM-coupled bilayers, they are different
and this difference originates from different magnetic fields for which the
spin patterns have been simulated and also from opposite signs of the
interlayer coupling. IN the simulation results shown in Fig.8(b), we show the
spin patterns of both top and bottom layers. Note, for small magnetic fields
(H=0 and H=0.5 T), the spin patterns of both layers are practically opposite
due to antiferromagnetic interlayer coupling. The situation is different for
stronger magnetic field, H=2.28 T, which is strong enough to overcome te
interlayer exchange coupling. In that case the spin patterns of both layer
reveal individual skyrmion in the spin pattern of the top layer and also
individual skyrmion in the bottom layer. These two skyrmions, however, are not
coupled and appear in different positions.

Finally, we also note here, that the hysteresis loops and the corresponding
spin patterns of the monolayers and bilayers, as shown above, are for
isotropic exchange interactions. These loops and spin patterns may become
remarkably modified when the exchange interaction has anisotropic
contributions, or more complex anisotropy terms (like two-ion anisotrop) occur
in the system. A typical example of such interactions is the so-called Kitaev
term.

###  2.3 Artificially created Skyrmion lattices

In the spin patterns shown above for monolayers, Figure 5, and for bilayers,
Fig.7(b) and Fig.8(b), only single skyrmion states have been found in the
corresponding simulation processes for the assumed parameters. Of course, for
other parameters, especially for stronger DMI, one can achieve stable multi-
skyrmion textures or skyrmion lattices [32, 33]. However, we show now that one
can also create arrays of (quasi)stable skyrmions by appropriate artificial
arrangement of initial spin states, by reversing manually the local
magnetization at the nodes of a certain lattice, see the top panel in Fig.9.
The parts (a) and (e) of this panel present random initial states, while the
others are arranged manually in the magnetic field 0.3T (b-d) and 0.748T
(f-h). The evolution of the initial states is analyzed by Monte Carlo
simulations, and the second and third rows of these columns show the spin
structure after 300k and 1M of Monte Carlo steps.

Figure 9: Spin patterns in FM-coupled Rh-stacked bilayers upon cooling for
high density of skyrmions. The left part (the columns (a-d)) is for $B$=0.3T,
while the right part (columns (e-h)) is for $B$=0.748T. The top row shows the
initial spin configurations: random in (a) and (e) and set artificially in the
others. The other two rows correspond to the spin structure of the bottom
layer after certain number of Monte Carlo steps, as indicated. The other
parameters: $d_{xy}=0.35J_{1}$, $J_{z}=J_{2}$, and NNN DMI as mentioned in the
text.

In case (a) we get at the and (after 1M steps) individual as well as coupled
skyrmions. In the case (b) the initially created square lattice evolves into a
square lattice of skyrmions with a size smaller than the size of the initially
created areas of reversed magnetization. In turn, in case (c) and (d) the
system after 1M of simulation steps evolves to smaller skyrmions arranged in a
defected hexagonal lattice. The origin of the defects might be twofold. First,
they appear naturaly, because of the square shape of the simulation cell,
which is incompatible with a hexagonal lattice. Second, defects might be
excited at increasing temperature. The columns (e) to (h) in Fig.9 correspond
to the same initial states as in columns (a) to (d), but simulations are in
the magnetic field of 0.748T. In all the four cases (e) to (h), after 1M of
Monte Carlo steps we get the arrays of skyrmions (of approximately equal
size). Interestingly, the square lattice in column (f) survives after 1M of
simulation steps. From the above results one may conclude that the simulated
skyrmion lattices may be stable or quasi-stable in time, or can evolve into a
disordered hexagonal lattice. This also shows that one can intentionally
create skyrmion lattices, which may be of some interest in skyrmionics.
However, their stabilization requires further investigations [46]

###  2.4 Dynamical properties

Now, we consider briefly the spin excitations (spin waves) in both monolayers
and bilayers of CrI3. Therefore, we assume zero magnetic filed to exclude
skyrmin texture formation when DMI is relatively large. First, we consider
spin waves in bilayers assuming intrinsic NNN DMI parameters as given in the
main text. As for the NN DMI we assume rather small its component $d_{xy}$,
i.e., $d_{xy}=0.05J_{1}$, which may be induced externally, e.g, by a gate
voltage or strain due to a substrate, while for the second component we assume
$d_{z}=0$. We focus on spin waves in bilayers in the AA and Rh stackings and
with FM and AFM interlayer exchange coupling. The corresponding dispersion
curves for the bilayers with FM and AFM interlayer exchange coupling are shown
in Fig.10. The dispersion curves have been calculated using ASD Vampire code
[42], and are shown there along the main paths in the Brillouin zone. Within
this approach, dynamics of individual spins is described in terms of Landau-
Lifshitz-Gilber equation. To determine the spin wave mode, the system is
externally excited and then Fourier analysis of the time evolution is used to
determine the spectrum. As in a monolayer there are two modes [7, 36, 37, 38,
39], one can expect four modes in bilayers due to splitting induced by the
exchange coupling between the two individual monolayers. However, the
interlayer coupling is rather small, so the splitting is not well resolved.
Indeed, in the case of FM interlayer coupling, the splitting of the modes are
not resolved in Figs10(a,b) for the AA and Rh stackings. The splitting is
clearly resolved in the case of AFM interlayer coupling in the AA stacking,
see Fig.10(c), while it is not well resolved in the Rh stacking, see Fig.10(d)
However, in the latter case there is a small gap in the spectrum, induced by
the NNN DMI. However, this gap is very small due to small values of NNN DMI
parameters. To check variation of this gap with NNN DMI we consider the spin
wave spectrum in the monolayer case, where the NNN DMI can be significantly
enhanced, e.g., in the corresponding chromium trihalide Janus structures [33].

Figure 10: Dispersion curves of spin waves in bilayers with NNN DMI included
(the corresponding parameters are given in the main text), and with NN DMI
$d_{xy}=0.05J_{1}$. Spin wave spectrum for the bilayer in the (a) AA stacking
with FM interlayer coupling $J_{z}=J_{2}$; (b) Rh stacking with FM interlayer
coupling $J_{z}=J_{2}$; (c) AA stacking with AFM interlayer coupling
$J_{z}=-J_{2}$; (d) Rh stacking with AFM interlayer coupling $J_{z}=-J_{2}$.
Other exchange parameters as given in the main text. Figure 11: Dispersion
curves of spin waves in a monolayer (a) without DMI, and (b) – (d) with NNN
DMI. The NNN DMI parameters used in the calculations: $d_{xy}^{\rm
NNN}=-0.146\times\eta\,{\rm meV}$, $d_{z}^{\rm NNN}=-0.0176\times\eta\ {\rm
meV}$, where (b) $\eta=1$, (c) $\eta=2$, and (d) $\eta=4$. The other
parameters: $J_{1}=1.224\,{\rm meV}$, $J_{2}=0.1986\,{\rm meV}$, and
$J_{3}=0.0027\,{\rm meV}$.

Now we analyse spin waves in monolayers [7, 36, 37, 38, 39], assuming
relatively large values of the NNN DMI parameters, and focus on the
topological features in the spectrum associated with NNN DMI. Though this term
in DMI is small in pristine CrI3, it can be remarkably enhanced, e.g., in the
chromium trihalide Janus structures [33], to produce remarkable effects in the
spin wave spectrum by creating a gap in the Dirac point K of the Brillouin
zone. The spin wave modes for monolayers are simulated within the ASD Uppsala
code [44] . As for exchange parameters, we note, that these parameters are not
well determined yet and there is some scatter in the available data [47, 15,
7]. We assume here the parameters consistent with those provided in Ref. [7].
The results of simulations are shown in Fig.10 for four different values of
the NNN DMI parameters. In the absence of NNN DMI, there is no gap in the
spectrum. However, with increasing DMI parameters the gap opens in the
spectrum at the Dirac point K of the Brillouin zone, and increases with
increasing parameters of NNN DMI.

##  3 Summary and conclusions

In this paper we have studied magnetic properties of chromium trihalides using
numerical simulations within the atomistic spin dynamics and Monte Carlo
methods. These simulation techniques have been used to determine critical
temperatures, hysteresis curves, and spin wave frequencies. These
characteristics have been analysed for monlayers as well as bilayers in two
(AA and rhombohedral) stackings and for ferromagnetic and antiferromagnetic
interlayer exchange coupling. The main focus was on the influence of
Dzyaloshinskii-Moriya interaction on the basic magnetic properties, especially
on the hysteresis curves, spin textures (stripe domains, skyrmion textures,
and others), and on the spin wave excitations, where Dzyaloshinskii-Moriya
coupling between next-nearest-neighbours leads to opening a gap in the
spectrum at the Dirac K point. Both, skyrmion formation and gap opening
require quite a relatively strong DMI between nearest-neighbours and next-
nearest-neighbours. Accordingly to observe these features one needs enhance
the DMI either by external electric field (or strain), or use chromium
trihalide Janus structures [33].

## References

  * [1] Wang, D. & Sanyal, B.  Systematic study of monolayer to trilayer cri3: Stacking sequence dependence of electronic structure and magnetism.  _The Journal of Physical Chemistry C_ 125, 18467–18473, DOI: 10.1021/acs.jpcc.1c04311 (2021).  https://doi.org/10.1021/acs.jpcc.1c04311.
  * [2] Huang, B. _et al._ Layer-dependent ferromagnetism in a van der waals crystal down to the monolayer limit.  _Nature_ 546, 270–273, DOI: 10.1038/nature22391 (2017).
  * [3] Gong, C. _et al._ Discovery of intrinsic ferromagnetism in two-dimensional van der waals crystals.  _Nature_ 546, 265–269, DOI: 10.1038/nature22060 (2017).
  * [4] Wang, Q. H. _et al._ The magnetic genome of two-dimensional van der waals materials.  _ACS nano_ 16, 6960–7079, DOI: 10.1021/acsnano.1c09150 (2022).
  * [5] Jain, A. _et al._ Commentary: The materials project: A materials genome approach to accelerating materials innovation.  _APL_ 1, 011002, DOI: 10.1063/1.4812323 (2013).
  * [6] Sarkar, S. & Kratzer, P.  Electronic correlation, magnetic structure and magnetotransport in few-layer cri3.  _Physical Review Materials_ 4, 104006, DOI: 10.1103/PhysRevMaterials.4.104006 (2020).
  * [7] Kartsev, A., Augustin, M., Evans, R. F. L., Novoselov, K. S. & Santos, E. J. G.  Biquadratic exchange interactions in two-dimensional magnets.  _npj Comput Mater_ 6, 150, DOI: 10.1038/s41524-020-00416-1 (2020).
  * [8] Chhowalla, M. _et al._ The chemistry of two-dimensional layered transition metal dichalcogenide nanosheets.  _Nature chemistry_ 5, 263–275, DOI: 10.1038/nchem.1589 (2013).
  * [9] Feng, J. _et al._ Metallic few-layered VS2 ultrathin nanosheets: high two-dimensional conductivity for in-plane supercapacitors.  _Journal of the American Chemical Society_ 133, 17832–17838, DOI: 10.1021/ja207176c (2011).
  * [10] Zhang, H., Liu, L.-M. & Lau, W.-M.  Dimension-dependent phase transition and magnetic properties of vs 2.  _Journal of Materials Chemistry A_ 1, 10821–10828, DOI: 10.1039/C3TA12098H (2013).
  * [11] Zhang, L., Zhou, J., Li, H., Shen, L. & Feng, Y. P.  Recent progress and challenges in magnetic tunnel junctions with 2d materials for spintronic applications.  _Applied Physics Reviews_ 8, 021308, DOI: 10.1063/5.0032538 (2021).
  * [12] Bałanda, M. _et al._ Magnetic ordering in the double-layered molecular magnet $\text{Cu}(\text{tetren})[\text{W}{(\text{CN})}_{8}]$: Single-crystal study.  _Phys. Rev. B_ 78, 174409, DOI: 10.1103/PhysRevB.78.174409 (2008).
  * [13] Kuzniak-Glanowska, E. _et al._ Engineering of the xy magnetic layered system with adeninium cations: Monocrystalline angle-resolved studies of nonlinear magnetic susceptibility.  _Inorganic Chemistry_ 60, 10186–10198, DOI: 10.1021/acs.inorgchem.1c00432 (2021).  PMID: 34232628, https://doi.org/10.1021/acs.inorgchem.1c00432.
  * [14] Konieczny, P., Pełka, R., Czernia, D. & Podgajny, R.  Rotating magnetocaloric effect in an anisotropic two-dimensional cuii[wv(cn)8]3– molecular magnet with topological phase transition: Experiment and theory.  _Inorganic Chemistry_ 56, 11971–11980, DOI: 10.1021/acs.inorgchem.7b01930 (2017).  PMID: 28915020, https://doi.org/10.1021/acs.inorgchem.7b01930.
  * [15] Jaeschke-Ubiergo, R., Morell, E. S. & Nunez, A. S.  Theory of magnetism in the van der waals magnet cri3.  _PRB_ 103, 174410, DOI: 10.1103/PhysRevB.103.174410 (2021).  2012.03099.
  * [16] Olsen, T.  Theory and simulations of critical temperatures in CrI3 and other 2d materials: Easy-axis magnetic order and easy-plane kosterlitz-thouless transitions.  _MRS Communications_ 9, 1142–1150, DOI: 10.1557/mrc.2019.117 (2019).
  * [17] Li, P. _et al._ Single-layer CrI3 grown by molecular beam epitaxy.  _Science Bulletin_ 65, 1064–1071, DOI: 10.1016/j.scib.2020.03.031 (2020).
  * [18] Lado, J. L. & Fernández-Rossier, J.  On the origin of magnetic anisotropy in two dimensional cri3.  _2D Materials_ 4, 035002, DOI: 10.1088/2053-1583/aa75ed (2017).
  * [19] Dupont, M. _et al._ Monolayer ${\mathrm{crcl}}_{3}$ as an ideal test bed for the universality classes of 2d magnetism.  _Phys. Rev. Lett._ 127, 037204, DOI: 10.1103/PhysRevLett.127.037204 (2021).
  * [20] Morell, E. S., León, A., Miwa, R. H. & Vargas, P.  Control of magnetism in bilayer cri3 by an external electric field.  _2D Materials_ 6, 025020, DOI: 10.1088/2053-1583/ab04fb (2019).
  * [21] Ghosh, S., Stojic, N. & Binggeli, N.  Structural and magnetic response of CrI3 monolayer to electric field.  _Physica B: Condensed Matter_ 570, 166–171, DOI: 10.1016/j.physb.2019.06.040 (2019).
  * [22] Ebrahimian, A., Dyrdał, A. & Qaiumzadeh, A.  Control of magnetic states and spin interactions in bilayer crcl3 with strain and electric fields: an ab initio study.  _Scientific Reports_ 13, 5336 (2023).
  * [23] León, A. M., González, J. W., Mejía-López, J., de Lima, F. C. & Morell, E. S.  Strain-induced phase transition in cri3 bilayers.  _2D Materials_ 7, 035008, DOI: 10.1088/2053-1583/ab8268 (2020).
  * [24] Xu, C., Wang, Q.-J., Xu, B. & Hu, J.  Effect of biaxial strain and hydrostatic pressure on the magnetic properties of bilayer CrI3.  _Frontiers of Physics_ 16, 53502, DOI: 10.1007/s11467-021-1073-x (2021).
  * [25] Zheng, F. _et al._ Tunable spin states in two-dimensional magnet CrI3.  _arXiv:1709.05472 [cond-mat]_ (2017).
  * [26] Ghosh, A. _et al._ Exotic magnetic and electronic properties of layered CrI3 single crystals under high pressure.  _Physical Review B_ 105, L081104, DOI: 10.1103/PhysRevB.105.L081104 (2022).  2108.00173.
  * [27] Wu, Z., Yu, J. & Yuan, S.  Strain-tunable magnetic and electronic properties of monolayer CrI ${}_{\textrm{3}}$.  _Physical Chemistry Chemical Physics_ 21, 7750–7755, DOI: 10.1039/C8CP07067A (2019).
  * [28] Li, T. _et al._ Pressure-controlled interlayer magnetism in atomically thin CrI3.  _Nature Materials_ 18, 1303–1308, DOI: 10.1038/s41563-019-0506-1 (2019).
  * [29] Guo, K. _et al._ Layer dependence of stacking order in nonencapsulated few-layer CrI3.  _Science China Materials_ 63, 413–420, DOI: 10.1007/s40843-019-1214-y (2020).
  * [30] Niu, B. _et al._ Coexistence of magnetic orders in two-dimensional magnet cri3.  _Nano Letters_ 20, 553–558, DOI: 10.1021/acs.nanolett.9b04282 (2020).
  * [31] Jiang, P. _et al._ Stacking tunable interlayer magnetism in bilayer CrI 3.  _Physical Review B_ 99, 144401, DOI: 10.1103/PhysRevB.99.144401 (2019).
  * [32] Behera, A. K., Chowdhury, S. & Das, S. R.  Magnetic skyrmions in atomic thin CrI3 monolayer.  _Applied Physics Letters_ 114, 232402 (2019).
  * [33] Xu, C., Feng, J., Xiang, H. & Bellaiche, L.  Magnetic skyrmion state in janus monolayers of chromium trihalides cr(i,x)3.  _Physical Review B_ 101, 060404, DOI: 10.1103/PhysRevB.101.060404 (2020).
  * [34] Liu, J., Shi, M., Mo, P. & Lu, J.  Electrical-field-induced magnetic Skyrmion ground state in a two-dimensional chromium tri-iodide ferromagnetic monolayer.  _AIP Advances_ 8, 055316, DOI: 10.1063/1.5030441 (2018).  https://pubs.aip.org/aip/adv/article-pdf/doi/10.1063/1.5030441/19733022/055316_1_online.pdf.
  * [35] Ghosh, S., Stojić, N. & Binggeli, N.  Comment on “Magnetic skyrmions in atomic thin CrI3 monolayer” [Appl. Phys. Lett. 114, 232402 (2019)].  _Applied Physics Letters_ 116, 086101, DOI: 10.1063/1.5128277 (2020).  https://pubs.aip.org/aip/apl/article-pdf/doi/10.1063/1.5128277/14531416/086101_1_online.pdf.
  * [36] Jin, W. _et al._ Raman fingerprint of two terahertz spin wave branches in a two-dimensional honeycomb ising ferromagnet.  _Nature Communications_ 9, 5122, DOI: 10.1038/s41467-018-07547-6 (2018).
  * [37] Aguilera, E., Jaeschke-Ubiergo, R., Vidal-Silva, N., Foa, L. & Núñez, A.  Topological magnonics in the two-dimensional van der waals magnet CrI3.  _Physical Review B_ 102, 024409, DOI: 10.1103/PhysRevB.102.024409 (2020).
  * [38] Zhang, X.-X. _et al._ Gate-tunable spin waves in antiferromagnetic atomic bilayers.  _Nature Materials_ 19, 838–842, DOI: 10.1038/s41563-020-0713-9 (2020).
  * [39] Hidalgo-Sacoto, R. _et al._ Magnon valley hall effect in CrI3-based vdW heterostructures.  _Physical Review B_ 101, 205425, DOI: 10.1103/PhysRevB.101.205425 (2020).
  * [40] Cenker, J. _et al._ Direct observation of two-dimensional magnons in atomically thin cri3.  _Nature Physics_ 17, 20–25 (2021).
  * [41] Evans, R. F. L. _et al._ Atomistic spin model simulations of magnetic nanomaterials.  _Journal of Physics: Condensed Matter_ 26, 103202, DOI: 10.1088/0953-8984/26/10/103202 (2014).
  * [42] Asselin, P. _et al._ Constrained monte carlo method and calculation of the temperature dependence of magnetic anisotropy.  _Phys. Rev. B_ 82, 054415, DOI: 10.1103/PhysRevB.82.054415 (2010).
  * [43] Evans, R. F. L., Atxitia, U. & Chantrell, R. W.  Quantitative simulation of temperature-dependent magnetization dynamics and equilibrium properties of elemental ferromagnets.  _Phys. Rev. B_ 91, 144425, DOI: 10.1103/PhysRevB.91.144425 (2015).
  * [44] Skubic, B., Hellsvik, J., Nordström, L. & Eriksson, O.  A method for atomistic spin dynamics simulations: implementation and examples.  _Journal of Physics: Condensed Matter_ 20, 315203, DOI: 10.1088/0953-8984/20/31/315203 (2008).
  * [45] Zhang, X. _et al._ Self-intercalation tunable interlayer exchange coupling in a synthetic van der waals antiferromagnet.  _Advanced Functional Materials_ 32, 2202977, DOI: https://doi.org/10.1002/adfm.202202977 (2022).  https://onlinelibrary.wiley.com/doi/pdf/10.1002/adfm.202202977.
  * [46] Wang, W., Sun, J. & Li, H.  Stabilization of skyrmions in two-dimensional systems with next-nearest-neighbor exchange interactions.  _Frontiers in Physics_ 10, DOI: 10.3389/fphy.2022.995902 (2022).
  * [47] Pizzochero, M., Yadav, R. & Yazyev, O. V.  Magnetic exchange interactions in monolayer cri3 from many-body wavefunction calculations.  _2D Materials_ 7, 035005, DOI: 10.1088/2053-1583/ab7cab (2020).

## Acknowledgements (not compulsory)

This work has been supported by the Norwegian Fi- nancial Mechanism 2014–2021
under the Polish Nor- wegian Research Project NCN GRIEG 2Dtronics no.
2019/34/H/ST3/00515.

## Author contributions statement

S.S. performed numerical simulations using Vampire code. P.B. performed
numerical simulations using ASD Uppsala code. M.J. supported numerical
calculations of S.S.. A.D. conceived of the idea. A.D and J.B supervised the
work and wrote the manuscript with input from all authors.