# Magnetic domain wall dynamics under external electric field in bilayer CrI3

Sahar Izadi Vishkayi  School of Physics, Institute for Research in Fundamental
Sciences (IPM), Tehran 19395-5531, Iran  Reza Asgari<EMAIL_ADDRESS>School of
Physics, Institute for Research in Fundamental Sciences (IPM), Tehran
19395-5531, Iran  School of Physics, University of New South Wales,
Kensington, NSW 2052, Australia

###### Abstract

Motivated by manipulating the magnetic order of bilayer CrI3, we carry out
microscopic calculations to find the magnetic order and various magnetic
domains of the system in the presence of an electric field. Making use of
density functional simulations, a spin model Hamiltonian is introduced
consisting of isotropic exchange couplings, Dzyaloshinskii-Moriya (DM)
interaction, and on-site magnetic anisotropy. The spin dynamics of two well-
known states of bilayer CrI3, low temperature (LT) and high temperature (HT)
phases, are obtained by solving the Landau-Lifshitz-Gilbert equation. We show
that the magnetic texture is stacking-dependent in bilayer CrI3 and stable
magnetic domains can appear in the HT stack which are tunable by external
electric and magnetic fields. Therefore, we suggest that the HT phase
represents a promising candidate for data storage in the modern generation of
spintronic devices working on magnetic domain engineering.

##  I Introduction

Two-dimensional (2D) magnetic materials have been investigated by experimental
Huang  _et al._ (2017); Deng  _et al._ (2018); Gong  _et al._ (2017);
Gibertini  _et al._ (2019); Gong and Zhang (2019); Meseguer-Sánchez  _et al._
(2021) and theoretical Lado and Fernández-Rossier (2017); Morell  _et al._
(2019); Torelli and Olsen (2020); Frey  _et al._ (2019); Menichetti  _et al._
(2019) researchers to analyze their attractive properties from the initial
report of the successful synthesis of thin layer CrI3 Huang  _et al._ (2017).
The presence of a substantial magnetocrystalline anisotropy was discovered,
thus breaking the invariance under spin rotations, and the crystal magnetic
ordering was observed to be a layer-dependent phenomenon. Because of the
existence of metastable magnetic domains with ultra-thin domain walls (like
permanent nanoparticle magnets Hubert and Schäfer (2008)) which possess spin
dynamics Wahab  _et al._ (2021); Chen  _et al._ (2021); Akram  _et al._
(2021); Xu  _et al._ (2021), CrI3, in monolayer and multilayer forms, has
recently been described as a viable contender for the modern generation of
memory devices. The microscopic nature of these magnetic domains in the
presence of external fields is still poorly understood, offering a significant
challenge for theoretical materials science.

Two geometrically stable stacks of bilayer CrI3 structures are known as low-
temperature (LT) and high-temperature (HT) phases (in more than $210$ K)
Morell  _et al._ (2019); Jiang  _et al._ (2019); Kim  _et al._ (2019). The LT
phase is a ferromagnetic (FM) ground state in which the magnetic moments of
the Cr atoms are aligned in both layers, whereas the magnetic moments of the
Cr atoms in distinct layers of the HT phase are in an antiferromagnetic (AFM)
configuration Morell  _et al._ (2019); Lei  _et al._ (2021); Song  _et al._
(2019); Li  _et al._ (2019); Sivadas  _et al._ (2018). The competition between
orbitally dependent interlayer AFM super-super-exchange and interlayer FM
super-super-exchange causes the magnetic ordering. In the HT bilayer CrI3, a
vertical external electric field generates an AFM to an FM magnetic transition
Jiang  _et al._ (2018a); Morell  _et al._ (2019); Jiang  _et al._ (2018b).
Because of their sensitivity to the electric field, the HT phase of the
bilayer CrI3 can be more efficient in ultrafast memory devices, and it could
be the first device for the engineering of magnetic domains in the world of 2D
magnets. For a monolayer up to bulk structures, the Curie temperature of a
thin layer CrI3 ranges from $45$ K to $61$ K Huang  _et al._ (2017); Kashin
_et al._ (2020). The magnetic ordering temperature of a thin layer CrI3 is
raised by the number of layers, according to experiments Huang  _et al._
(2017). The Néel temperature of $45$ K is reported for the HT bilayer Huang
_et al._ (2018) and it is expected that the LT bilayer maintains a larger
Curie temperature due to strong FM ordering between the layers Wang and Sanyal
(2021).

Although the Ising Hamiltonian model is a classic model for explaining the
magnetic behavior of 2D CrI3, it cannot explain the experimentally observed
spin wave Chen  _et al._ (2018). To explain the spin canting physics, an
extended Heisenberg Hamiltonian containing Kitaev Lee  _et al._ (2020),
biquadratic anisotropy exchange term Wahab  _et al._ (2021), or
Dzyaloshinskii-Moria (DM) Vishkayi  _et al._ (2020); Jaeschke-Ubiergo  _et
al._ (2021) interaction terms was proposed. The spin wave in a CrI3 structure
Chen  _et al._ (2020) is described by a Heisenberg Hamiltonian with Kitaev or
DM interactions, and a small portion of the DM interactions in 2D CrI3 is
responsible for a topological gap Jaeschke-Ubiergo  _et al._ (2021).

The ability to alter and adjust magnetic domains of 2D bilayer magnets to make
them useful in technology and spintronics is a key question here. We focus on
the bilayer CrI3 in the presence of an external electric field Huang  _et al._
(2018); Jiang  _et al._ (2018a, b); Morell  _et al._ (2019); Lei  _et al._
(2021) to answer this specific topic.

In this study, we use a spin model Hamiltonian in the presence of external
fields, with isotropic exchange and DM interactions. We use density functional
theory (DFT) computations to derive a suitable spin model Hamiltonian for the
bilayer CrI3 that accurately reproduces the magnetic characteristics of the
system. The isotropic exchange coupling parameters of the bilayer CrI3 have
been investigated in the presence of an external electric field Lei  _et al._
(2021), but the Dzyaloshinskii-Moria (DM) interaction has yet to be
investigated. In contrast to what we found for a monolayer CrI3 Vishkayi  _et
al._ (2020), our calculations suggest that the DM interaction in the HT
bilayer CrI3 is substantially connected to the applied electric field. The
spin dynamics of the system are calculated using the Landau-Lifshitz-Gilbert
formalism Ellis  _et al._ (2015) with DFT-derived exchange parameters as
starting inputs. To obtain a temperature dependent magnetization, we also
employ the Monte Carlo Metropolis technique with $10^{6}$ time steps for a
$50\times 50$ nm2 layer.

The weak strength of the isotropic intralayer exchange coupling signifies the
magnetic ordering of the layers, which can be modified by an external
electrical field, according to our findings. The findings also reveal that an
electric field, like a magnetic field, alters the magnetic texture of the
system. Furthermore, by applying an electric field resulting from an increase
in the intralayer isotropic exchange coefficient, the temperature of the
magnetic order of the HT bilayer is raised. The electric field also affects
the magnetic texturing on the domain walls. The results reveal that the
magnetic textures in the studied systems can be manipulated using an external
electric field. To be more precise, in the LT bilayer, we show that the DM
interactions lead to the creation of a quasi-circular domain shape in the
layers, but its value is extremely smaller than the intralayer exchange
coupling to stabilize it in the system, therefore, the obtained magnetic
domain is unstable and disappearing quickly.

Figure 1: (Color online) The top view ((a) and (c)) and side view ((b) and
(d)) of the HT and LT bilayer CrI3 structures, respectively. The black
lozenges show the unit cell of bilayers. The blue (pink) circles represent
chromium (iodine) atoms. In (b) and (d), chromium atoms are labeled by 1-4 in
the unit cell. In (b) the bonding angles between the atoms of each layer are
shown by $\alpha$, $\beta$ and $\theta$. The crystal parameters are given in
Table I.

This paper is organized as follows. We commence with a description of our
theoretical formalism in Sec. II, followed by the details of the DFT
simulations, the spin model Hamiltonian and spin dynamics. Numerical results
for the spin-spin interaction parameters in the presence of an electric field,
the exchange coefficients in both the HT and LT phases and spin dynamics in
two cases are reported in Sec. III. We summarize our main findings in Sec. IV.

##  II Theory and Models

###  II.1 DFT calculations

We use Quantum Espresso package Giannozzi  _et al._ (2009) to simulate
density-functional based calculations. The generalized gradient exchange-
correlation within Perdew-Burke-Ernzerhof Perdew  _et al._ (1996) functional
is considered to expose the electronic and magnetic ground states of two
bilayer CrI3 structures. An energy cut-off of $50$ Ry for wave vector and an
$8\times 8\times 1$ grid mesh of $k$-points within the first Brillouin zone
are defined as the converged input parameters. To avoid repeats effects, we
consider a vacuum of $20$ Å along the $z$-direction. The van der Waals
Grimme-D2 Grimme (2006) correction is used to consider the interaction between
layers. DFT$+U$ calculations Cococcioni and De Gironcoli (2005) are also
needed to find correct magnetic ground states of CrI3 bilayers. At this point,
we consider $U=3$ eV as the on-site Hubbard parameter similar to Refs. Jiang
_et al._ (2019); Lei  _et al._ (2021); Leon  _et al._ (2020); Akram  _et al._
(2021). To obtain a reliable total ground-state energy, we maintain a high
degree of accuracy of $10^{-10}$ eV. The bilayers are relaxed until the
maximum force on each atom is $0.01$ eV/Å even under the electric field
conditions. It should be noted that the consideration of the spin-orbit
coupling (SOC) and noncollinear spin-polarization are essential to obtain the
DM interaction and magnetic anisotropic energy (MAE).

###  II.2 Spin-model Hamiltonian

A spin model Hamiltonian for a 2D magnetic hexagonal lattice in each layer is
utilized to provide spin-spin interactions in a magnetic 2D bilayer CrI3
system:

$\displaystyle{\cal H}=$ $\displaystyle\sum_{i,j}[\frac{1}{2}J_{ij}{\bf
S}_{i}\cdot{\bf S}_{j}+{\gamma}_{i}|S_{iz}|^{2}+$
$\displaystyle\frac{1}{2}{\bf D}_{ij}\cdot({\bf S}_{i}\times{\bf
S}_{j})+\sum_{i}\mu_{B}{{\bf B}_{ext}}\cdot{\bf{S}}_{i}],$ (1)

where $J_{ij}$ is the isotropic exchange coupling parameter between $i$ and
$j$ atoms, $\gamma_{i}$ denotes a magnetic anisotropy coefficient, ${\bf
D}_{ij}$ is the DM interaction vector and the last term shows the effect of an
external magnetic field, ${\bf B}_{ext}$, and $\mu_{B}$ is the magnetic Bohr.
The magnetic moment of ith Cr atom is indicated by ${\bf S}_{i}$ in Eq. (1).
It is intriguing to compare the strength of the interlayer and intralayer
exchange coefficient parameters which are respectively referred by $v$ and $t$
indices in this paper.

###  II.3 Spin Dynamics

The Landau-Lifshitz-Gilbert (LLG) equation Ellis  _et al._ (2015) is used to
atomistic model the spin dynamics of the system after the coefficients of the
spin Hamiltonian are established. We follow the method introduced by Evans et
al. Evans  _et al._ (2014) implemented in VAMPIRE code. The time evolution of
the spin direction is given by LLG equation according to

$\frac{\partial{\bf S}_{i}}{\partial t}=-\frac{\gamma}{1+\lambda^{2}}[{\bf
S}_{i}\times{\bf H}_{eff}^{i}+\lambda{\bf S}_{i}\times({\bf S}_{i}\times{\bf
H}_{eff}^{i})],$ (2)

where $\gamma$ is the gyromagnetic ration and $\lambda$ is microscopic damping
equal to unity for finding the equilibration states quickly. The effective
field applied on each spin is defined by ${\bf H}_{eff}^{i}=-{\partial\cal
H}/{\partial{\bf S}_{i}}+{\bf H}_{th}^{i}$ and the effective thermal field,
${\bf H}_{th}^{i}$, is counted by Langevin dynamic approach Brown (1979) to
consider thermal spin fluctuations.

In order to obtain the spin dynamics simulations, a 100 nm $\times$ 100 nm
square bilayer CrI3 was kept in thermal equilibrium above the critical (Néel
or Curie) temperature and then linearly cooled to 0 K in 2 ns. The magnetic
moments of atoms are randomized at a temperature above the critical
temperature, so we chose random initial spin directions for the spin dynamics
simulation. The time step of 1 fs is considered in the integration scheme of
Heun’s method García-Palacios and Lázaro (1998) for solving LLG equations,
suitable for materials with low critical temperature and large dimensions to
obtain stable results Evans  _et al._ (2014).

The long-range demagnetization field is also considered in the spin dynamics
calculations. The magnetic dipole-dipole interaction induces local spin-flip
in the magnetic materials and leading to the creation of complicated spin
texture even in the materials with a FM ground state Kawaguchi  _et al._
(2007); Ezawa (2010). We discretize the layer to macrocells with $1$ nm3
volume to avoid expenses of the demagnetization field computing at the
atomistic level. The magnetic moments of macrocells, ${\bf m}_{mc}$, are
defined by the sum of atomic magnetic moments within the macrocell and the
demagnetization field of ith macrocell is given by

${\bf H}_{d}^{mci}={\frac{{\mu}_{0}}{4\pi}}(\sum_{i\neq j}{\bf
M}_{ij}\cdot{\bf m}_{mc}^{j})-{\frac{\mu_{0}}{3}}{\frac{{\bf
m}_{mc}^{i}}{V_{mc}^{i}}},$ (3)

where $V_{mc}^{i}$ is the volume of ith macrocell, $\mu_{0}=4\pi\times
10^{-7}$ is the vacuum permeability and ${\bf M}_{ij}$ is the dipole-dipole
interaction matrix between $i$ and $j$ macrocells defined as Evans  _et al._
(2014),

${\bf
M}_{ij}=\begin{bmatrix}\frac{3r_{x}r_{x}-1}{r_{ij}^{3}}-\frac{1}{3}&3r_{x}r_{y}&3r_{x}r_{z}\\\
3r_{x}r_{y}&\frac{3r_{y}r_{y}-1}{r_{ij}^{3}}-\frac{1}{3}&3r_{y}r_{z}\\\
3r_{x}r_{z}&3r_{y}r_{z}&\frac{3r_{z}r_{z}-1}{r_{ij}^{3}}-\frac{1}{3}\end{bmatrix},$
(4)

where the positions of macrocells are obtained by the magnetic center of the
mass relation and the vector between the position of $i$ and $j$ macrocells is
shown by ${\bf r}_{ij}=r_{ij}(r_{x}\hat{i}+r_{y}\hat{j}+r_{z}\hat{k})$. We
consider that the sheet is linearly cooled in $2$ ns from an equilibrium
temperature, upper the Curie temperature, to zero Kelvin. To obtain
temperature dependent magnetization, we use the Monte Carlo Metropolis
algorithm with $10^{6}$ time steps for a $50$$\times 50$ nm2 sheet.

##  III Numerical Results and Discussions

The magnetic ground state, symmetric and asymmetric exchange coefficients, and
temperature-dependent magnetization of bilayer CrI3 in both the HT and LT
phases in the ground state are presented in this section. The effects of an
external electric field on the coefficients are then studied. Finally, we
investigate the spin dynamics of both phases in the presence of a
perpendicular electric field.

The HT phase of bilayer CrI3 is in the form of the monoclinic structure, while
the LT bilayer possesses a rhombohedral crystal structure as displayed in
Figs. 1(a)-(d). The different crystal structures leads to the unique magnetic
behaviors. Both the HT and LT phases are relaxed by spin-dependent DFT
calculations and their optimized lattice constant, bonding lengths and angles
are reported in Table 1. The results reveal that the bonding angles and
lengths of both layers are the same, and that the geometry of the bilayer is
unaffected by the electric field, which is in good agreement with those
described in Ref. Morell  _et al._ (2019). To obtain a formation of the
energy, we use the straight formula of $E_{\rm
F}=E_{t}(bilayer)-n_{Cr}E(Cr)-n_{I}E(I)$ where $E_{t}(bilayer)$ is the total
energy of the bilayer CrI3, $n_{Cr}(n_{I})$ is the number of the chromium
(iodine) atoms in the unit cell, and $E(Cr)$ and $E(I)$ are the energy of free
Cr and I atoms, respectively. The formation energies of the HT and LT bilayers
are given in Table 1 indicating that the both phases are energetically favored
to occur.

Table 1: Lattice constant, $a$, Cr-Cr (I-I) distance between the layers, $h_{Cr-Cr}(h_{I-I})$, Cr-I (Cr-Cr) bonding length, $d_{Cr-I}(d_{Cr-Cr})$, two different I-Cr-I bonding angles, $\beta_{I-Cr-I}$ and $\alpha_{I-Cr-I}$, and Cr-I-Cr bonding angle, $\theta_{Cr-I-Cr}$, and formation energy, $E_{f}$ of the HT and LT bilayer CrI3 structures. The bonding angles are shown in Fig. 1 (b).  Phase |  $a(\AA)$ | $h_{Cr-Cr}(\AA)$ | $h_{I-I}(\AA)$ | $d_{Cr-I}(\AA)$ | $d_{Cr-Cr}(\AA)$ | $\beta_{I-Cr-I}(^{o})$ | $\alpha_{I-Cr-I}(^{o})$ | $\theta_{Cr-I-Cr}(^{o})$ |  $E_{F}$ (eV/per Cr)
---|---|---|---|---|---|---|---|---|---
HT | 6.977 | 7.24 | 3.99 | 2.81 | 4.03 | 177.4 | 89.9 | 91.6 | -8.81
LT | 6.979 | 6.59 | 3.35 | 2.81 | 4.03 | 177.6 | 89.9 | 91.7 | -8.86

###  III.1 Exchange coefficients in the HT phase

In the HT phase of bilayer CrI3, there are four Cr atoms in the unit cell; two
in each layer. At zero external magnetic field, if the nearest neighbor
interlayer and intralayer interactions are considered, the extended first two
terms of the model Hamiltonian (Eq. (1)) in the unit cell can be written as

$\displaystyle{\cal H}=$ $\displaystyle\frac{1}{2}(\eta_{t}J_{t}[({\bf
S}_{1}\cdot{\bf S}_{2})+({\bf S}_{2}\cdot{\bf S}_{1})+({\bf S}_{3}\cdot{\bf
S}_{4})+({\bf S}_{4}\cdot{\bf S}_{3})]+$ $\displaystyle\eta_{1v}J_{1v}[({\bf
S}_{1}\cdot{\bf S}_{3})+({\bf S}_{3}\cdot{\bf S}_{1})+({\bf S}_{2}\cdot{\bf
S}_{4})+({\bf S}_{4}\cdot{\bf S}_{2})]+$ $\displaystyle\eta_{2v}J_{2v}[({\bf
S}_{1}\cdot{\bf S}_{4})+({\bf S}_{4}\cdot{\bf S}_{1})+({\bf S}_{2}\cdot{\bf
S}_{3})+({\bf S}_{3}\cdot{\bf S}_{2})])+$ $\displaystyle\eta\gamma|S_{z}|^{2}$
(5)

where $\eta_{t}=3$ is the number of the nearest neighbors in each layer and
$\eta_{1v}=\eta_{2v}=1$ are the number of type-$1$ (between atoms numbered by
$1$ ($2$) and $3$ ($4$) in Fig. 1 (b)) and type-$2$ (between atoms numbered by
$1$ ($2$) and $4$ ($3$) in Fig. 1 (b)) neighbors in the interlayer. $\eta=4$
is the number of Cr atoms in the unit cell. The exchange coupling coefficients
are obtained by mapping the model Hamiltonian onto the DFT-based calculations
which are computed for various magnetic configurations. For calculation of the
isotropic exchange coefficient, $4$-distinct spin configurations are
considered as; FM (${\bf S}_{1}=S\hat{k},{\bf S}_{2}=S\hat{k},{\bf
S}_{3}=S\hat{k},{\bf S}_{4}=S\hat{k}$), AFM1 (${\bf S}_{1}=S\hat{k},{\bf
S}_{2}=S\hat{k},{\bf S}_{3}=-S\hat{k},{\bf S}_{4}=-S\hat{k}$), AFM2 (${\bf
S}_{1}=S\hat{k},{\bf S}_{2}=-S\hat{k},{\bf S}_{3}=-S\hat{k},{\bf
S}_{4}=S\hat{k}$) and AFM3 (${\bf S}_{1}=S\hat{k},{\bf S}_{2}=-S\hat{k},{\bf
S}_{3}=S\hat{k},{\bf S}_{4}=-S\hat{k}$) in which all the magnetic moments of
Cr atoms are aligned in the $z$-direction and $S=3/2$. In Fig. 2, the
considered magnetic orders are represented in bilayer CrI3. In these
configurations, the DM terms of Eq. (1) are vanished owing to the parallel
magnetic moments of Cr atoms, and their total energies are obtained by

$\displaystyle
E_{FM}=2S^{2}(\eta_{t}J_{t}+\eta_{1v}J_{1v}+\eta_{2v}J_{2v})+\eta\gamma
S^{2},$ (6) $\displaystyle
E_{AFM1}=2S^{2}(\eta_{t}J_{t}-\eta_{1v}J_{1v}-\eta_{2v}J_{2v})+\eta\gamma
S^{2},$ $\displaystyle
E_{AFM2}=2S^{2}(-\eta_{t}J_{t}-\eta_{1v}J_{1v}+\eta_{2v}J_{2v})+\eta\gamma
S^{2},$ $\displaystyle
E_{FM3}=2S^{2}(-\eta_{t}J_{t}+\eta_{1v}J_{1v}-\eta_{2v}J_{2v})+\eta\gamma
S^{2}.$

We carry out spin-dependent DFT calculations to find the total energies of the
considered magnetic configurations and as a consequence the exchange coupling
parameters of the HT bilayer are given by

$\displaystyle
J_{t}=\frac{(E_{FM}+E_{AFM1})-(E_{AFM2}+E_{AFM3})}{8S^{2}\eta_{t}},$ (7)
$\displaystyle J_{v}=\frac{E_{FM}-E_{AFM1}}{4S^{2}\eta_{v}},$ (8)

where $J_{v}=J_{1v}+J_{2v}$ is the interlayer exchange coupling and
$\eta_{v}=\eta_{1v}=\eta_{2v}$. The intralayer exchange coupling is
$J_{t}=-7.00$ meV which presents the strong ferromagnetic coupling of Cr atoms
in the layers. $J_{v}=+0.08$ meV indicates the layers are coupled in a form of
antiferromagnetic. Furthermore, the weak strength of the isotropic interlayer
exchange coupling denotes the magnetic order of the layers can be tuned by the
external electrical field Morell  _et al._ (2019).

Figure 2: (Color online) The magnetic ordering of Cr atoms in (a) FM, (b)
AFM1, (c) AFM2 and (d) AFM3 spin configurations of HT bilayer CrI3. The dark
(light) balls represent the Cr atoms in the top (bottom) layer and arrows show
the magnetic moment vectors of Cr atoms in the z-direction. The dashed lines
connect the atoms with intralayer and interlayer exchange coupling which are
indicated, respectively, by t and v indices.

We observe the variation of the isotropic interlayer and intralayer exchange
coupling by DFT-based calculations, shown in Fig. 3. It is observed that the
sign of the interlayer exchange coupling is changed between $2$-$2.5$ V/nm,
hence a magnetic transition from the AFM to the FM occurs between layers while
the strength of the intralayer exchange coefficient is no longer altered by
the electric field. This means that the ferromagnetic coupling of the Cr atoms
in the layer is extremely strong.

Figure 3: (Color online) (a) The intralayer (squares) and interlayer
(lozenges) isotropic exchange coefficients as a function of the external
electric field, $E_{l}$. The change of the sign of the interlayer exchange
shows a magnetic transition from the AFM to a FM between the layers. The
ferromagnetic coupling of Cr atoms in the layer is extremely strong. (b) The
non-zero component of interlayer and intralayer DM vectors of the bilayer CrI3
structures as a function of $E_{l}$. The intralayer exchange coupling is
changed significantly by the electric field. Although the $z$ and $x$
components of the interlayer DM interactions are independent of the electric
field, the sign and the value of $D_{v}^{y}$ vary.

The extended Hamiltonian for the DM term is given by

$\displaystyle{\cal H}_{DM}=$ $\displaystyle\frac{1}{2}[{\bf
D}_{t12}\cdot({\bf S}_{1}\times{\bf S}_{2})+{\bf D}_{t12^{\prime}}\cdot({\bf
S}_{1}\times{\bf S}_{2^{\prime}})+{\bf D}_{t12^{\prime\prime}}\cdot$
$\displaystyle({\bf S}_{1}\times{\bf S}_{2^{\prime\prime}})+{\bf
D}_{t21}\cdot({\bf S}_{2}\times{\bf S}_{1})+{\bf D}_{t21^{\prime}}\cdot({\bf
S}_{2}\times{\bf S}_{1^{\prime}})$ $\displaystyle+{\bf
D}_{t21^{\prime\prime}}\cdot({\bf S}_{2}\times{\bf S}_{1^{\prime\prime}})+{\bf
D}_{t34}\cdot({\bf S}_{3}\times{\bf S}_{4})+{\bf D}_{t34^{\prime}}\cdot$
$\displaystyle({\bf S}_{3}\times{\bf S}_{4^{\prime}})+{\bf
D}_{t34^{\prime\prime}}\cdot({\bf S}_{3}\times{\bf S}_{4^{\prime\prime}})+{\bf
D}_{t43}\cdot({\bf S}_{4}\times{\bf S}_{3})$ $\displaystyle+{\bf
D}_{t43^{\prime}}\cdot({\bf S}_{4}\times{\bf S}_{3^{\prime}})+{\bf
D}_{t43^{\prime\prime}}\cdot({\bf S}_{4}\times{\bf S}_{3^{\prime\prime}})+$
$\displaystyle{\bf D}_{v13}\cdot({\bf S}_{1}\times{\bf S}_{3})+{\bf
D}_{v31}\cdot({\bf S}_{3}\times{\bf S}_{1})+$ $\displaystyle{\bf
D}_{v14}\cdot({\bf S}_{1}\times{\bf S}_{4})+{\bf D}_{v41}\cdot({\bf
S}_{4}\times{\bf S}_{1})+$ $\displaystyle{\bf D}_{v23}\cdot({\bf
S}_{2}\times{\bf S}_{3})+{\bf D}_{v32}\cdot({\bf S}_{3}\times{\bf S}_{2})+$
$\displaystyle{\bf D}_{v24}\cdot({\bf S}_{2}\times{\bf S}_{4})+{\bf
D}_{v42}\cdot({\bf S}_{4}\times{\bf S}_{2})],$ (9)

where ${\bf S}_{i^{\prime}}$ and ${\bf S}_{i^{\prime\prime}}$ are equal to the
magnetic moment of ith Cr atom, and ${\bf S}_{i^{\prime}}={\bf
S}_{i^{\prime\prime}}={\bf S}_{i}$. On the other hand, the DM vector between
two atoms can be obtained by ${\bf D}_{ij}=D_{ij}{\hat{\bf u}}_{ij}$, where
$D_{ij}=D_{ji}$ is a scalar and ${\hat{\bf u}}_{ij}$ is a unit-vector in the
direction of the line connecting $i$ and $j$ atoms Yang  _et al._ (2015).
Therefore, ${\bf D}_{ij}=-{\bf D}_{ji}$, ${\bf D}_{t12^{\prime}}={\bf
R}_{z}(120){\bf D}_{t12}$, ${\bf D}_{t12^{\prime\prime}}={\bf R}_{z}(240){\bf
D}_{t12}$ and ${\bf R}_{z}(\theta)$ is a $\theta$ degree rotation around the
$z$-axis. Three first terms of Eq. (9) can be written as

$\displaystyle{\cal H}$ ${}_{DM12}=$ (10) $\displaystyle({\bf D}_{t12}+{\bf
R}_{z}(120){\bf D}_{t12}+{\bf R}_{z}(240){\bf D}_{t12})\cdot({\bf
S}_{1}\times{\bf S}_{2}).$

If ${\bf D}_{t12}=(D_{t12}^{x}\hat{i}+D_{t12}^{y}\hat{j}+D_{t12}^{z}\hat{k})$,
then $({\bf D}_{t12}+{\bf R}_{z}(120){\bf D}_{t12}+{\bf R}_{z}(240){\bf
D}_{t12})=3D_{t12}^{z}\hat{k}$, hence only the $z$-component of the DM vector
between two Cr atoms in the layers of the bilayer CrI3 can be entered in the
model Hamiltonian of Eq. (1). In the similar way, we can write $({\bf
D}_{t34}+{\bf R}_{z}(120){\bf D}_{t34}+{\bf R}_{z}(240){\bf
D}_{t34})=3D_{t34}^{z}\hat{k}$ for the DM interaction between Cr atoms in the
top layer. Due to the symmetries between the layers, it is possible to
consider $D_{t12}^{z}=D_{t34}^{z}=D_{t}^{z}$. On the other hand, ${\bf
D}_{t12}=-{\bf D}_{t21}$ and ${\bf D}_{t13}=-{\bf D}_{t31}$, hence for the
first five lines of Eq. (9), we have the intralayer terms of the DM
interaction as,

${\cal H}_{DMt}=\frac{1}{2}[6D_{t}^{z}\hat{k}\cdot({\bf S}_{1}\times{\bf
S}_{2})+6D_{t}^{z}\hat{k}\cdot({\bf S}_{3}\times{\bf S}_{4})].$ (11)

Therefore, the symmetries of the HT phase dictate that just the $z$-component
of the intralayer DM interaction vector, $D_{t}^{z}$, has a non-zero value.

Figure 4: (Color online) (a) The normalized magnetization of the HT bilayers
as a function of temperature under various external electric field. (b) The
magnetic order temperature, $T^{*}$ ($T^{*}=T_{N}$ for $E_{l}\leqslant 2$ and
$T^{*}=T_{C}$ for $E_{l}\geqslant 2.5$), (lozenges) and magnetic anisotropy
energy (squares) of the HT bilayer CrI3 versus the external electric field.
The magnetic order temperature of the HT bilayer is increased by applying the
electric field that arises from the increasing of the intralayer isotropic
exchange coefficient. Furthermore, the magnetic moments of the Cr atoms are in
out-of-plane direction and it is unchangeable by the external electric field.

Moreover, the geometry of the bilayer crystal indicates that $\hat{\bf
u}_{13}=\hat{\bf u}_{24}$ and $\hat{\bf u}_{32}={\bf R}_{y}(180)\hat{\bf
u}_{14}$, then we can assume ${\bf D}_{v13}={\bf D}_{v24}$ and ${\bf
D}_{v32}={\bf R}_{y}(180){\bf D}_{v14}$. To obtain the interlayer contribution
of the DM vector, we have to focus on the last four lines of Eq. (9) which can
be summarized in the following form,

$\displaystyle{\cal H}_{DMv}$ $\displaystyle={\bf D}_{v13}\cdot({\bf
S}_{1}\times{\bf S}_{3})+{\bf D}_{v13}\cdot({\bf S}_{2}\times{\bf S}_{4})+{\bf
D}_{v14}\cdot$ $\displaystyle({\bf S}_{1}\times{\bf S}_{4})+{\bf
R}_{y}(180){\bf D}_{v14}\cdot({\bf S}_{3}\times{\bf S}_{2})$ (12)

where ${\bf D}_{v13}=(D_{v1}^{x},D_{v1}^{y},D_{v1}^{z})$, ${\bf
D}_{v14}=(D_{v2}^{x},D_{v2}^{y},D_{v2}^{z})$ and ${\bf R}_{y}(180){\bf
D}_{v14}=(-D_{v2}^{x},D_{v2}^{y},-D_{v2}^{z})$. Eight different spin
configurations are needed to obtain enough information about the interlayer
and intralayer DM vectors given in Table 2.

Table 2: Eight spin configurations for calculation of the interlayer and intralayer DM vectors of the HT bilayer.  Config. | $1^{st}$ | $2^{nd}$ | $3^{rd}$ | $4^{th}$ | $5^{th}$ | $6^{th}$ | $7^{th}$ | $8^{th}$
---|---|---|---|---|---|---|---|---
${\bf S}_{1}$ | $S\hat{i}$ | $S\hat{i}$ | $S\hat{i}$ | $S\hat{i}$ | $S\hat{k}$ | $S\hat{k}$ | $S\hat{k}$ | $S\hat{k}$
${\bf S}_{2}$ | $S(\cos(\frac{\pi}{6})\hat{i}+$ | $S(\cos(\frac{\pi}{6})\hat{i}-$ | $S\hat{i}$ | $S\hat{i}$ | $S\hat{k}$ | $S\hat{k}$ | $S\hat{k}$ | $S\hat{k}$
$\cos(\frac{\pi}{6})\hat{j})$ | $\sin(\frac{\pi}{6})\hat{j})$ |  |  |  |  |  |
${\bf S}_{3}$ | $S\hat{i}$ | $S\hat{i}$ | $S(\cos(\frac{\pi}{6})\hat{i}+$ | $S(\cos(\frac{\pi}{6})\hat{i}-$ | $S(\cos(\frac{\pi}{6})\hat{k}+$ | $S(\cos(\frac{\pi}{6})\hat{k}-$ | $S(\cos(\frac{\pi}{6})\hat{k}+$ | $S(\cos(\frac{\pi}{6})\hat{k}-$
|  | $\sin(\frac{\pi}{6})\hat{j})$ | $\sin(\frac{\pi}{6})\hat{j})$ | $\sin(\frac{\pi}{6})\hat{i})$ | $\sin(\frac{\pi}{6})\hat{i})$ | $\sin(\frac{\pi}{6})\hat{j})$ | $\sin(\frac{\pi}{6})\hat{j})$
${\bf S}_{4}$ | $S(\cos(\frac{\pi}{6})\hat{i}+$ | $S(\cos(\frac{\pi}{6})\hat{i}-$ | $S(\cos(\frac{\pi}{6})\hat{i}+$ | $S(\cos(\frac{\pi}{6})\hat{i}-$ | $S(\cos(\frac{\pi}{6})\hat{k}+$ | $S(\cos(\frac{\pi}{6})\hat{k}-$ | $S(\cos(\frac{\pi}{6})\hat{k}+$ | $S(\cos(\frac{\pi}{6})\hat{k}-$
$\sin(\frac{\pi}{6})\hat{j})$ | $\sin(\frac{\pi}{6})\hat{j})$ | $\sin(\frac{\pi}{6})\hat{j})$ | $\sin(\frac{\pi}{6})\hat{j})$ | $\sin(\frac{\pi}{6})\hat{i})$ | $\sin(\frac{\pi}{6})\hat{i})$ | $\sin(\frac{\pi}{6})\hat{j})$ | $\sin(\frac{\pi}{6})\hat{j})$

The DM interactions are calculated by mapping the considered spin
configurations on the model Hamiltonian and comparing the spin-dependent DFT-
obtained total energies. At this place, we want to focus on the DM terms
described in Eqs. (11) and (12).

$\displaystyle{\cal
H}_{DMt(1^{st})}=6S[(D_{t}^{z}\hat{k})\cdot(\sin(\pi/6)\hat{k})],$
$\displaystyle{\cal
H}_{DMt(2^{nd})}=-6S[(D_{t}^{z}\hat{k})\cdot(\sin(\pi/6)\hat{k})],$
$\displaystyle{\cal H}_{DMv(1^{st})}={\cal H}_{DMv(2^{nd})}=0.$ (13)

By imposing the above relations into Eq. (1), the total energy equations are
defined and the $z$-component of the intralayer DM vector is given by

$D_{t}^{z}=\frac{E_{1st}-E_{2nd}}{12S^{2}\sin(\pi/6)}.$ (14)

In a similar way, we can find the DM terms of model Hamiltonian for the other
configurations and the interlayer DM vector of the HT bilayer is achieved by
DFT-obtained total energies according to

$\displaystyle D_{v}^{z}$
$\displaystyle=D_{v1z}+D_{v2z}=\frac{E_{3rd}-E_{4th}}{4S^{2}\sin(\pi/6)},$
$\displaystyle D_{v}^{y}$
$\displaystyle=D_{v1y}=\frac{E_{5th}-E_{6th}}{4S^{2}\sin(\pi/6)},$
$\displaystyle D_{v}^{x}$
$\displaystyle=D_{v1x}+D_{v2x}=\frac{E_{8th}-E_{7th}}{4S^{2}\sin(\pi/6)}.$
(15)

We obtain $D_{t}^{z}=1.1$ $\mu$eV, $D_{v}^{z}=3.5$ $\mu$eV , $D_{v}^{y}=-38.8$
$\mu$eV and $D_{v}^{x}=3.1$ $\mu$eV for the HT bilayer. Accordingly, the
intralayer DM interaction is ignorable in comparison with the interlayer one.
The results show the strongest component of the interlayer DM vector is in the
$y$-direction; this means the magnetic moments of Cr atoms in the layers tend
to rotate along the $x$-$z$ plane, as long as they align in the $z$-direction
because of the high MAE of the bilayer.

The interlayer and intralayer DM interactions can be tuned by an external
electric field, very similar to the isotropic exchange coupling. In Fig. 3
(b), it is observed that the intralayer exchange coupling is changed
significantly by the electric field. Although the $z$ and $x$ components of
the interlayer DM interactions are independent of the electric field, the sign
and the value of $D_{v}^{y}$ vary. The sign of the $y$-component of the DM
changes from negative to positive values between $2$-$2.5$ V/nm where a
transition from the AFM to FM is achieved for the HT bilayer. In fact, a phase
transition needs a spin-canting between layers in the $x$-$z$ plane satisfied
by the non-zero $D_{v}^{y}$ which is in the same order of the $J_{v}$ in the
HT bilayer. The application of the electric field shifts the atomic positions
slightly, resulting in a change in the orbital couplings between atomic
neighbors. Indeed, if the $e_{g}-t_{2g}$ coupling, the resulting FM exchange
will larger than that causing the AFM exchange ($e_{g}-e_{g}$ orbital
coupling), the bilayer is in the FM ground states and vice versa Sivadas  _et
al._ (2018); Akram  _et al._ (2021). Upon increasing the electric field, an
AFM to FM transition is observed in the HT bilayer CrI3, showing the dominant
effect of $e_{g}-t_{2g}$ coupling. We expect these changes to be accompanied
by spin tilt, and consequently, the DM interactions representing the spin
tilting ability are tuned by the electric field.

To calculate the magnetic order temperature, Curie $T_{C}$, or Néel $T_{N}$
temperature of the FM or AFM bilayer, a sample of 50$\times$50 nm2 is
considered and the metropolis Monte Carlo algorithms are used for the
calculation of temperature-dependent magnetization (see Fig. 4 (a)). According
to Fig. 4 (b), the magnetic order temperature of the HT bilayer is increased
by applying the electric field that arises from the increasing of the
intralayer isotropic exchange coefficient.

The variation of $MAE=(E_{\perp}-E_{\parallel})/\eta S^{2}$ as a function of
the electric field is shown in Fig. 4 (b). It shows that the magnetic moments
of the Cr atoms are in out-of-plane direction, and it is unchangeable by the
external electric field.

###  III.2 Exchange coefficients in the LT phase

Now we turn our attention to another phase of the bilayer system. Four Cr
atoms in a unit cell of the LT-phase bilayer CrI3 prefer to determine a
parallel spin orientation in which the minimized total energy belongs to the
FM magnetic configuration Morell  _et al._ (2019). Here, similar to the HT
phase, we are interested in obtaining the symmetric and antisymmetric exchange
coefficients of the LT bilayer. We should note that the intralayer nearest
neighbors of the LT is equivalent to that in the HT phase, while there is only
an interlayer nearest neighbor bond between the Cr atoms numbered by $2$ and
$3$ in Fig. 1 (d). In the case of the LT bilayer CrI3, the first two terms of
Eq. (1) are extended to be as

$\displaystyle{\cal H}=$ $\displaystyle\frac{1}{2}(\eta_{t}J_{t}[({\bf
S}_{1}\cdot{\bf S}_{2})+({\bf S}_{2}\cdot{\bf S}_{1})+({\bf S}_{3}\cdot{\bf
S}_{4})+({\bf S}_{4}\cdot{\bf S}_{3})]+$ $\displaystyle\eta_{v}J_{v}[({\bf
S}_{2}\cdot{\bf S}_{3})+({\bf S}_{3}\cdot{\bf
S}_{2})])+\eta\gamma|S_{z}|^{2},$ (16)

where $\eta_{t}=3$ and $\eta_{v}=1$ are the number of the intralayer and
interlayer nearest neighbors. We consider four different magnetic
configurations of Cr atoms, as a FM (${\bf S}_{1}=S\hat{k}$, ${\bf
S}_{2}=S\hat{k}$, ${\bf S}_{3}=S\hat{k}$, ${\bf S}_{4}=S\hat{k}$), AFM1 (${\bf
S}_{1}=S\hat{k}$, ${\bf S}_{2}=S\hat{k}$, ${\bf S}_{3}=-S\hat{k}$, ${\bf
S}_{4}=-S\hat{k}$), AFM2 (${\bf S}_{1}=S\hat{k}$, ${\bf S}_{2}=-S\hat{k}$,
${\bf S}_{3}=-S\hat{k}$, ${\bf S}_{4}=S\hat{k}$) and AFM3 (${\bf
S}_{1}=S\hat{k}$, ${\bf S}_{2}=-S\hat{k}$, ${\bf S}_{3}=S\hat{k}$, ${\bf
S}_{4}=-S\hat{k}$), and thus able to find interlayer and intralayer symmetric
exchange coefficients. By applying the spin vectors on Eq. (16), the relations
for the total energies of the magnetic configurations are given by,

$\displaystyle E_{FM}=E_{0}+2S^{2}\eta_{t}J_{t}+S^{2}\eta_{v}J_{v}+\eta\gamma
S^{2},$ $\displaystyle
E_{AFM1}=E_{0}+2S^{2}\eta_{t}J_{t}-S^{2}\eta_{v}J_{v}+\eta\gamma S^{2},$
$\displaystyle
E_{AFM2}=E_{0}-2S^{2}\eta_{t}J_{t}+S^{2}\eta_{v}J_{v}+\eta\gamma S^{2},$
$\displaystyle
E_{AFM3}=E_{0}-2S^{2}\eta_{t}J_{t}-S^{2}\eta_{v}J_{v}+\eta\gamma S^{2}.$ (17)

The DFT-obtained total energies give us the symmetric exchange coefficients as

$\displaystyle J_{t}$
$\displaystyle=\frac{(E_{FM}+E_{AFM1})-(E_{AFM2}+E_{AFM3})}{8S^{2}\eta_{t}},$
$\displaystyle J_{v}$
$\displaystyle=\frac{(E_{FM}-E_{AFM1})}{2S^{2}\eta_{v}}.$ (18)

The results displayed in Fig. 5 (a) show the interlayer and intralayer
symmetric exchange couplings possess negative and consistent values with the
FM ground state of the LT bilayers. The sign of $J_{t}$ and $J_{v}$ remain
negative by applying the electric field, while their amounts are slightly
increased. In fact, the electric field develops the ferromagnetic coupling of
layers increase in the LT bilayer.

Figure 5: (Color online) (a) The intralayer (squares) and interlayer
(lozenges) isotropic exchange coefficients and (b) the $z$-component of
interlayer and intralayer DM vectors of the LT bilayer CrI3 as a function of
the external electric field, $E_{l}$. The sign of the exchange coefficients
remain negative by applying the electric field, while their amounts are
slightly increased showing that the electric field makes the ferromagnetic
coupling of layers stronger. Furthermore, the $z-$ component of the DM
interaction of the LT bilayer increases by the electric field meaning that the
spin-canting are in the plane of the Cr atoms.

The third term of the Hamiltonian in Eq. (1) can be written for the LT bilayer
as,

$\displaystyle{\cal H}_{DM}$ $\displaystyle=\frac{1}{2}[{\bf
D}_{t12}\cdot({\bf S}_{1}\times{\bf S}_{2})+{\bf D}_{t12^{\prime}}\cdot({\bf
S}_{1}\times{\bf S}_{2}^{\prime})+{\bf D}_{t12^{\prime\prime}}$
$\displaystyle\cdot({\bf S}_{1}\times{\bf S}_{2}^{\prime\prime})+{\bf
D}_{t21}\cdot({\bf S}_{2}\times{\bf S}_{1})+{\bf D}_{t21^{\prime}}\cdot({\bf
S}_{2}\times{\bf S}_{1}^{\prime})$ $\displaystyle+{\bf
D}_{t21^{\prime\prime}}\cdot({\bf S}_{2}\times{\bf S}_{1}^{\prime\prime})+{\bf
D}_{t34}\cdot({\bf S}_{3}\times{\bf S}_{4})+{\bf D}_{t34^{\prime}}\cdot({\bf
S}_{3}$ $\displaystyle\times{\bf S}_{4}^{\prime})+{\bf
D}_{t34^{\prime\prime}}\cdot({\bf S}_{3}\times{\bf S}_{4}^{\prime\prime})+{\bf
D}_{t43}\cdot({\bf S}_{4}\times{\bf S}_{3})+{\bf D}_{t43^{\prime}}$
$\displaystyle\cdot({\bf S}_{4}\times{\bf S}_{3}^{\prime})+{\bf
D}_{t43^{\prime\prime}}\cdot({\bf S}_{4}\times{\bf S}_{3}^{\prime\prime})+{\bf
D}_{v23}\cdot({\bf S}_{2}\times{\bf S}_{3})$ $\displaystyle+{\bf
D}_{v32}\cdot({\bf S}_{3}\times{\bf S}_{2})].$ (19)

Similar to the HT phase, the symmetric of atomic positions dictates to have
only the $z$-component of the intralayer DM interaction in the Hamiltonian
which is given by Eq. (11). For the interlayer DM interaction of the LT phase

${\cal H}_{DMv}={\bf D}_{v23}\cdot({\bf S}_{2}\times{\bf S}_{3}).$ (20)

It should be noticed that there is a $C_{3}$ rotation axis on the bonding
length between Cr atoms numbered by $2$ and $3$ in the unit cell of the LT
bilayers (see Fig. 1 (d)). According to the Moriya symmetry rules Moriya
(1960), the DM interaction between atoms $2$ and $3$ is parallel to the
$C_{3}$ axis which is along the $z$-direction. Therefore, we need to find the
$z$-component of the interlayer DM interaction which can be obtained by two
spin configurations (see Table 3). Finally, the $z$-component of intralayer
and interlayer DM interactions are respectively given by

$\displaystyle D_{t}^{z}=$
$\displaystyle\frac{(E_{1st}-E_{2nd})+(E_{3rd}-E_{4th})}{12S^{2}\sin(\pi/6)},$
$\displaystyle D_{v}^{z}=$
$\displaystyle\frac{E_{3rd}-E_{4th}}{2S^{2}\sin(\pi/6)}.$ (21)

Table 3: Four spin configurations for calculation of the interlayer and intralayer DM vectors of the LT bilayer CrI3. Config. | $1^{st}$ | $2^{nd}$ | $3^{rd}$ | $4^{th}$
---|---|---|---|---
${\bf S}_{1}$ | $S\hat{i}$ | $S\hat{i}$ | $S\hat{i}$ | $S\hat{i}$
${\bf S}_{2}$ | $S(\cos(\frac{\pi}{6})\hat{i}+$ | $S(\cos(\frac{\pi}{6})\hat{i}-$ | $S\hat{i}$ | $S\hat{i}$
$\sin(\frac{\pi}{6})\hat{j})$ | $\sin(\frac{\pi}{6})\hat{j})$ |  |
${\bf S}_{3}$ | $S\hat{i}$ | $S\hat{i}$ | $S(\cos(\frac{\pi}{6})\hat{i}+$ | $S(\cos(\frac{\pi}{6})\hat{i}-$
|  | $\sin(\frac{\pi}{6})\hat{j})$ | $\sin(\frac{\pi}{6})\hat{j})$
${\bf S}_{4}$ | $S(\cos(\frac{\pi}{6})\hat{i}+$ | $S(\cos(\frac{\pi}{6})\hat{i}-$ | $S(\cos(\frac{\pi}{6})\hat{i}+$ | $S(\cos(\frac{\pi}{6})\hat{i}-$
$\sin(\frac{\pi}{6})\hat{j})$ | $\sin(\frac{\pi}{6})\hat{j})$ | $\sin(\frac{\pi}{6})\hat{j})$ | $\sin(\frac{\pi}{6})\hat{j})$

The negligible values of $D_{t}^{z}=-0.4$ $\mu$eV and $D_{v}^{z}=-0.3$ $\mu$eV
are obtained using our spin-dependent DFT-based calculations. It should be
noted that there is an inversion symmetry in the center of the $2$ and $3$
atomic bond in the LT bilayer, so $D_{v}^{z}$ should be zero, and therefore,
there is no any DM interaction between the layers. We investigate the effect
of the external electric field on the interlayer and intralayer DM
interactions of the LT phase, and the results are reported in Fig. 5 (b). It
is observed that the interlayer and intralayer DM interactions of the LT phase
are increased by developing the electric field. It is interesting that the
order of DM interaction variation is similar for both LT and HT phases, but
their behavior is different. In the presence of an electric field, just the
$z$-component of the DM interaction of the LT bilayer increases; this means
that the spin-canting are in the plane of the Cr atoms. On the other hand, the
interlayer symmetric exchange coupling, $J_{v}$, is much stronger than
$D_{v}^{z}$, subsequently, the FM coupling survives in a competition between
FM coupling of layers and in-plane canting of spins.

In the HT bilayer, the amount of the interlayer DM vector and interlayer
symmetric exchange coupling are in the same order. In fact, they are in a
close competition with each other and it remains in the external electric
field. It is intriguing that the both of them leads to a transition from the
AFM to FM in the HT phase; the interlayer DM interaction is dominant in the
$y$-direction which leads to the spin canting in the $x$-$z$-plane and, as a
consequence, helps to the rotation of the spin-direction from the AFM to FM
and also the sign of $J_{v}$ is changed by applying the electric field.

It should be noticed that the interlayer and intralayer isotropic exchange
coupling in the Heisenberg Hamiltonian of the bilayer CrI3 have been
calculated in Ref. Lei  _et al._ (2021) which agrees well with our work. They
obtained $J_{t}=-6.49$ meV (by DFT) and $J_{v}=+0.10$ meV (by experimental
results from Ref. Huang  _et al._ (2017)) which are close to our findings
($J_{t}=-7.00$ meV and $J_{v}=+0.8$ meV) in order of magnitude and sign for HT
bilayer CrI3. They have also shown that the isotropic exchange couplings are
tuned by the electric field. There are DFT-based calculations Jiang  _et al._
(2019); Morell  _et al._ (2019); Lei  _et al._ (2021) and experimental results
Huang  _et al._ (2018); Jiang  _et al._ (2018b, a) showing that the interlayer
exchange coupling sign of HT phase of bilayer CrI3 is changed by the electric
field. One can find the DFT-obtained exchange couplings of different stacks of
bilayer CrI3 in the absence of an electric field in Refs. Sivadas  _et al._
(2018); Jiang  _et al._ (2019); Leon  _et al._ (2020); Akram  _et al._ (2021).
In contrast, the interlayer and intralayer DM interactions of the bilayer CrI3
have not been calculated in previous works, and we report here their electric
field tunable behavior.

The temperature dependence of the magnetic moments of the Cr atoms, $M(T)$ in
the LT bilayer is illustrated in Fig. 6 (a). It is shown that $M(T)$ is nearly
independent of the electric field. The variation of the Currie temperature,
$T_{C}$, by the electric field arises directly from the variation of $J_{t}$
and it is gradually increased in a more intense electric field (see Fig. 6
(b)). The MAE of the LT phase as a function of the electric field is shown in
Fig. 6 (b). The MAE remains negative by varying the electric field so the easy
axis of the LT bilayer lies in the $z$-direction and its value does not change
significantly.

We obtained a Néel temperature of 44 K for the AFM-HT bilayer and a Curie
temperature of 59 K for the FM-LT bilayer CrI3, which are tuned by the
electric field due to the variation in the parameters of the spin model
Hamiltonian. Increasing the number of layers from monolayer to bulk increases
the critical temperature between 45 K and 61 K Huang  _et al._ (2017); Kashin
_et al._ (2020). A Néel temperature of 45 K was reported experimentally for
the HT bilayer Huang  _et al._ (2018). The larger Curie temperature is due to
the stronger FM interlayer exchange coupling in LT bilayer CrI3 Wang and
Sanyal (2021).

Figure 6: (Color online) (a) The magnetization of the LT bilayers as a
function of temperature under different external electric field. (b) The Curie
temperature (lozenges) and magnetic anisotropy energy (squares) of LT bilayer
CrI3 versus external electric field. Notice that the MAE remains negative by
changing the electric field, therefore, the easy axis of the LT bilayer is in
the $z$-direction and its value does not change significantly.

Although the magnetic coupling coefficients of both the LT and HT bilayers are
tuned, the magnetic ground state of the HT bilayer can be manipulated by the
electric field. To perceive the possibility of manipulating magnetic textures
in the bilayer CrI3, we simulate their spin dynamics in the presence of the
external electric and magnetic fields and the emergence of the magnetic
domains and skyrmion patterns are explored.

To obtain the spin dynamics of the bilayer, the atomic position of the system
is constant while different atomic magnetic moments are appeared by time
evolutions Evans  _et al._ (2014). Here, we obtained that the formation energy
is 1000 times more than the exchange energy in the bilayer CrI3 which
guarantees the stability of the atomic structure under spin evolutions. In
fact, we want to mention that the spin evolution happens in a constant atomic
structure in any electric field. Under the electric field, the atomic
positions are relaxed by DFT-based calculations and leading variation of the
interband and intraband isotropic exchange and DM coefficients, so the spin
dynamics are calculated in the new structure.

###  III.3 Spin dynamics of the HT phase

In the HT bilayer, the AFM coupling between layers is dictated in the spin
dynamics of the system. The spin dynamics of the top and bottom layers are
plotted in Fig. 7. In $T=0.3$ K, it is observed that the bilayer is divided
into two areas in which layers are coupled in form of an anti-ferromagnetic.
While the bottom layer has a magnetic moment in the $+z$-direction, the top
layer has a magnetic moment in the $-z$-direction and vice versa; this means
that the AFM order is constant between the layers. This regularity is also
established on the domain walls, where the perpendicular components of the
magnetic moments, $M_{\perp}$, of two layers are in opposite directions. The
results of the simulations show that the created magnetic texture is stable in
size, but $M_{\perp}$ is slowly changed to be perpendicular to the domain wall
by the time evolution as shown in Figs. 7 (a)-(c). Thus, an ultrathin Néel-
type domain wall appears on both the top and bottom layers. It is clear that
the zero magnetic moments are obtained for the HT bilayer, but there is an
interesting microscopic magnetic pattern on each layer.

Figure 7: (Color online) The magnetic texture of the top and bottom layers of
the HT bilayer in $E_{l}=0$ V/nm, $B=0$, (a) $T=0.3$ K, ${\text{time}}=2$ ns,
(b) $T=0$, ${\text{time}}=2.2$ ns and (c) $T=0$, ${\text{time}}=3$ ns. The
blue (red) areas show where the magnetic moments of Cr atoms are in the
-z(+z)-direction. ${M}_{\perp}$ is represented by black arrows on the domain
walls. The created magnetic texture is stable in size, but ${M}_{\perp}$ is
slowly changed to be perpendicular to the domain wall by the time evolution.

We obtain the symmetric and asymmetric exchange coupling coefficients of the
spin Hamiltonian depending on the external electric field, therefore, we
explore the variation of the spin dynamics of the bilayer by the electric
field shown in Figs. 8 (a)-(c). In $E_{l}=1$ V/nm, the magnetic order of the
HT bilayer is the AFM and $J_{t}$ and $J_{v}$ are near to their values in the
absence of the electric field case, while the value of $D_{t}^{z}$ is
increased to $25.5$ $\mu$eV, therefore the new spin texture is created (see
Fig. 8 (b)). If we focus on one of the layers, we see that the wider magnetic
domains in $E_{l}=0$ case (Fig. 8 (a)) changes to two smaller quasi-circle
magnetic domains with finite diameters (nearly $25$ nm and $10$ nm). Skyrmion
patterns occur around the magnetic moments, while their chirality is different
for the top and bottom layers. The places and sizes of magnetic domains are
the same in the top and bottom layers, however, the direction of the magnetic
moments are exactly opposite for inside and outside of the domain walls. We
obtain these stable quasi-circle magnetic domains under the time evolution.

Figure 8: (Color online) The magnetic texture of the top and bottom layers of
the HT bilayer in (a) $E_{l}=0$, (b) $E_{l}=1$ V/nm and (c) $E_{l}=2$ V/nm,
$T=0.3$ K, $B=0$ and ${\text{time}}=2$ ns. The skyrmion patterns occur around
the magnetic moments, while their chirality is different for the top and
bottom layers. The places and sizes of magnetic domains are the same in the
top and bottom layers, however, the direction of the magnetic moments are
exactly opposite for inside and outside of the domain walls.

For the HT bilayer, the amount of the intralayer DM interaction reaches to
$41$ $\mu$eV by increasing the electric field to $E_{l}=2$ V/nm, and the
interlayer isotropic exchange coupling is equal to $-20$ $\mu$eV. Although the
value of $J_{v}$ is reduced, the negative sign shows the bilayer is still in
the AFM phase. These changes lead to a new spin texture in the bilayer. In
Fig. 8 (c), 30 nm-diameter quasi-circle magnetic domains are created in the
both top and bottom layers with opposite spin orientation, so the AFM
configuration of the bilayer is conserved. Most importantly, the diameter,
position, number of the magnetic domain and the domain wall chirality are
related to the external electric field. In fact, the magnetic texture of the
HT bilayer can be manipulated by an external electric field owing to the
variation of the DM and isotropic exchange coupling interactions. Therefore,
stable magnetic domains are constituted in the HT bilayer that possess
fascinating microscopic information, while the total magnetism is zero in the
macroscopic level.

Figure 9: (Color online) The magnetic texture of the top and bottom layers of
HT bilayer in $E_{l}=3$ V/nm, $T=0.3$ K, $B=0$ and ${\text{time}}=2$ ns.
Notice that the the magnetic domain in the FM bilayer is stable under the time
evolution.

In $E_{l}=3$ V/nm, the layers are coupled in a form of a ferromagnetic. This
phase transition is obviously observed in the spin dynamics simulation which
both layers have the same spin orientation (see Fig. 9) and originates from
the positive sign of $J_{v}$. On the other hand, the perpendicular-component
of the magnetic moments of the top and bottom layers are similar on the domain
walls. The spin dynamics simulations show that the obtained magnetic domain in
the FM bilayer is still stable under the time evolution.

Figure 10: (Color online) The spin dynamic results of the HT bilayer in
$E_{l}=0$ V/nm, (a) $B=0$, $T=0.3$ K and $time=2$ ns, (b) $B=0.5$ T, $T=0.3$ K
and ${\text{time}}=2$ ns and (c) $B=0.5$ T, $T=0$ K and ${\text{time}}=3$ ns.
the external magnetic field is in competition with the interlayer exchange
coupling and desire to bring the spin of the Cr atoms parallel to the external
magnetic field.

According to Figs. 10 (a)-(c), our simulation results show that the external
magnetic field changes the spin dynamic of the bilayer. In the absence of an
external electric field, Fig.10 (b) shows that there is an area with FM
configurations in the AFM background of the bilayer. Indeed, the external
magnetic field is in competition with the interlayer exchange coupling and
desire to accompany the spin of the Cr atoms parallel to the external magnetic
field. It is not completely winner in the competition when $B_{ext}$ is less
than $0.5$ T. By increasing the external magnetic field to $B_{ext}=2$ T, the
layers gain parallel magnetic moments in the direction of the applied magnetic
field.

Figure 11: (Color online) The spin dynamic results of the HT bilayer in
$E_{l}=2$ V/nm, $B=0.01$ T (a) $T=0.3$ K, $time=2$ ns, (b) $T=0$ K,
${\text{time}}=2.2$ ns and (c) $T=0$ K, ${\text{time}}=3$ ns. Notice, here,
the applied magnetic field can not create a phase transition in the magnetic
order of the HT bilayer.

Another example is the system when exposed to an external electric field,
$E_{l}=2$ V/nm. In $B_{ext}=0.5$ T, the AFM pattern is completely disappeared
and the magnetic moments of the Cr atoms are parallel. In the competition
between exchange coupling coefficients and the Zeeman energy, the Zeeman
energy of the strong magnetic field has conquered. In the presence of $0.01$ T
external magnetic field, the simulation results (Figs. 11(a)-(c)) show that
the quasi-circular magnetic domain form converts to a quasi-elliptical shape,
but the AFM configuration between the top and bottom layers is conserved
because the applied magnetic field is not large enough to create a phase
transition in the magnetic ordering of the bilayer. Interestingly, the
obtained magnetic domain and almost domain wall do not change by the time
evolution after cooling time (2 ns). Accordingly, a stable magnetic texture
can be obtained and tuned by a weak external magnetic field, while the strong
magnetic field leads to the magnetic phase transition in the bilayer.

###  III.4 Spin dynamic of the LT phase

In the LT bilayer, a strong negative interlayer exchange coupling leads to the
creation of the FM magnetic configuration. The spin dynamic simulation, in
Fig. 12 (a) shows two wide areas with opposite spin directions in each layer
while the magnetic moments of the layers are actually analogous. In the
boundary of the two areas, the perpendicular magnetic moments of the Cr atoms
are dominated and their chirality is the same in both layers. In fact, a 1D
spin-wave is created at the boundary and the results show that the chirality
of the spin-wave is changed smoothly in time (compare Fig. 12 (a) and (b)). On
the other hand, the created areas sizes are stable by the time evolution, so
we can consider that the 1D in-plane spin-wave is metastable (due to the
change of chirality). In the absence of an electric field, the ignorable DM
interaction cannot create a magnetic domain, and the observed spin-wave is due
to the demagnetization field. For clarity, we report the results of the spin
dynamic simulation in the absence of a demagnetization field in Fig. 12 (c).
It is obviously observed that the magnetic domain is disappeared and all the
Cr atoms of both layers have parallel magnetic moments due to strong
interlayer and intralayer isotropic exchange couplings of the LT bilayer CrI3.

Figure 12: (Color online) The spin dynamic results of the LT bilayer in
$E_{l}=0$ V/nm, $B_{ext}=0$ T (a) $T=0.3$ K and $time=2$ ns, (b) $T=0$ K and
${\text{time}}=3$ ns and (c) $T=0.3$ K and ${\text{time}}=2$ ns in the absence
of a demagnetization field. A 1D metastable spin-wave is created at the
boundary and the chirality of the spin wave is changed smoothly in time. In
this case, the Cr atoms have parallel magnetic moments due to strong
interlayer and intralayer isotropic exchange couplings.

The interlayer and intralayer DM interactions are increased by applying the
external electric field, while the sign of the exchange coupling and their
coefficients values remain practically constant in the LT bilayer CrI3. The
spin dynamic simulations show the existence of a similar quasi-circular
magnetic domain shape in the top and bottom layers in the presence of the
electric field. Hence, the meta-stable magnetic domain is disturbed by
increasing the DM interaction in $E_{l}=1$ V/nm and $E_{l}=2$ V/nm, see Figs.
13 (a)-(c). The results show that this quasi-circular magnetic domains are
disappeared by cooling the system and the magnetic moments of the Cr atoms
gained parallelized, see Fig. 14 for the LT bilayer under $E_{l}=1$ V/nm. We
can assume that the DM interactions leads to the creation of a quasi-circular
domain shape in the layers, but its value is extremely smaller than the
exchange coupling to stabilize it in the system, therefore, the obtained
magnetic domain is unstable and disappears rapidly in the LT bilayer.

Figure 13: (Color online) The magnetic texture of the top and bottom layers of
the LT bilayer in (a) $E_{l}=0$, (b) $E_{l}=1$ V/nm and (c) $E_{l}=2$ V/nm,
$T=0.3$ K, $B=0$ and ${\text{time}}=2$ ns. The existence of a similar quasi-
circle magnetic domain in the top and bottom layers in the presence of the
electric field are seen and the meta-stable magnetic domain is disturbed by
increasing the DM interaction. Figure 14: (Color online) The time evolution of
the top layer of the LT phase under $E_{l}=1$ V/nm, $T=0$ K and $B=0$. The
bottom layer has a similar spin dynamic results due to the FM coupling of the
layers.

The spin dynamics in the twisted bilayer CrI3 was recently calculated using an
arbitrary DM interaction and a local moiré field instead of an isotropic
exchange coupling between the layers Akram  _et al._ (2021) while ignoring the
demagnetization field. They showed that the twisted texture depends on the
twist angle. It is worth noting that applying an electric field is easier than
twisting the double layer CrI3 in experiments. It is experimentally shown that
the moiré magnetic domain in AFM twisted bilayer magnetic crystals can be
controlled by gate voltages Xu  _et al._ (2021). Despite a recent interest in
studying spin dynamics in 2D magnets Akram  _et al._ (2021); Xu  _et al._
(2021); Wahab  _et al._ (2021); Tong  _et al._ (2018); Khoshlahni  _et al._
(2019); Behera  _et al._ (2019); Abdul-Wahab  _et al._ (2021) and their
promising application in new generation magnetic domain-based devices Allwood
_et al._ (2005); Vélez  _et al._ (2019); Luo  _et al._ (2020), there is no
report of spin texture formation in the bilayer CrI3 by considering all DFT-
obtained isotropic and anisotropic exchange couplings in the presence of an
electric field. Here we calculate the spin dynamics of a large bilayer CrI3 as
an experimental domain by solving the LLG equations where the demagnetizing
field is important to find the best magnetic ground state of the system. The
results show that the magnetic domain in the HT bilayer CrI3 can be tuned by
the electric field, which is consistent with the experimental results Xu  _et
al._ (2021).

##  IV Conclusion

The interlayer and intralayer isotropic exchange coupling coefficients and DM
interactions are obtained by DFT-based calculations for bilayer CrI3. We have
shown that in the presence of an external electric field, interlayer DM
interactions have a critical impact on the magnetic phase transition of the HT
bilayer. The effects of the external electric field on the parameters of the
spin-model Hamiltonian and temperature-dependent magnetic moments of bilayers
are explored. Most importantly, a stable magnetic domain in the HT bilayer
CrI3 which can be manipulated by the electric and magnetic fields can be
constituted owing to the tunable interlayer and intralayer exchange coupling
and DM interactions. The spin dynamics simulations show that the HT bilayer
CrI3 can represent a promising candidate for spintronic and logic memory
applications. The interlayer exchange coupling, on the other hand, is only
weakly affected by the external electric field, hence there is no stable
magnetic domain in the LT bilayer CrI3. Current experiments can investigate
these findings. Furthermore, when disorder and magnetic impurity are taken
into account, our results can be extended in a system.

###### Acknowledgements.

We thank Z. Torbatian and S. Heydari for fruitful discussions. R. A thanks for
support from the Australian Research Council Centre of Excellence in Future
Low-Energy Electronics Technologies (project number CE170100039).

## References

  * Huang  _et al._ (2017) B. Huang, G. Clark, E. Navarro-Moratalla, D. R. Klein, R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A. McGuire, D. H. Cobden, _et al._ , “Layer-dependent ferromagnetism in a van der Waals crystal down to the monolayer limit,” Nature 546, 270 (2017).
  * Deng  _et al._ (2018) Y. Deng, Y. Yu, Y. Song, J. Zhang, N. Z. Wang, Z. Sun, Y. Yi, Y. Z. Wu, S. Wu, J. Zhu, _et al._ , “Gate-tunable room-temperature ferromagnetism in two-dimensional fe 3 gete 2,” Nature 563, 94 (2018).
  * Gong  _et al._ (2017) C. Gong, L. Li, Z. Li, H. Ji, A. Stern, Y. Xia, T. Cao, W. Bao, C. Wang, Y. Wang, _et al._ , “Discovery of intrinsic ferromagnetism in two-dimensional van der waals crystals,” Nature 546, 265 (2017).
  * Gibertini  _et al._ (2019) M. Gibertini, M. Koperski, A. F. Morpurgo, and K. S. Novoselov, “Magnetic 2d materials and heterostructures,” Nat. nanotechnol. 14, 408 (2019).
  * Gong and Zhang (2019) C. Gong and X. Zhang, “Two-dimensional magnetic crystals and emergent heterostructure devices,” Science 363, 4450 (2019).
  * Meseguer-Sánchez  _et al._ (2021) J. Meseguer-Sánchez, C. Popescu, J. L. García-Muñoz, H. Luetkens, G. Taniashvili, E. Navarro-Moratalla, Z. Guguchia, and E. J. Santos, “Coexistence of structural and magnetic phases in van der waals magnet cri3,” Nat. commun. 12, 1 (2021).
  * Lado and Fernández-Rossier (2017) J. L. Lado and J. Fernández-Rossier, “On the origin of magnetic anisotropy in two dimensional ${\mathrm{CrI}}_{3}$,” 2D Mater. 4, 035002 (2017).
  * Morell  _et al._ (2019) E. S. Morell, A. León, R. H. Miwa, and P. Vargas, “Control of magnetism in bilayer ${\mathrm{CrI}}_{3}$ by an external electric field,” 2D Mater. 6, 025020 (2019).
  * Torelli and Olsen (2020) D. Torelli and T. Olsen, “First principles heisenberg models of 2d magnetic materials: the importance of quantum corrections to the exchange coupling,” J. Condens. Matter Phys. 32, 335802 (2020).
  * Frey  _et al._ (2019) N. C. Frey, A. Bandyopadhyay, H. Kumar, B. Anasori, Y. Gogotsi, and V. B. Shenoy, “Surface-engineered mxenes: electric field control of magnetism and enhanced magnetic anisotropy,” ACS Nano 13, 2831 (2019).
  * Menichetti  _et al._ (2019) G. Menichetti, M. Calandra, and M. Polini, “Electronic structure and magnetic properties of few-layer cr2ge2te6: the key role of nonlocal electron–electron interaction effects,” 2d Mater. 6, 045042 (2019).
  * Hubert and Schäfer (2008) A. Hubert and R. Schäfer,  _Magnetic domains: the analysis of magnetic microstructures_ (Springer Science & Business Media, 2008).
  * Wahab  _et al._ (2021) D. A. Wahab, M. Augustin, S. M. Valero, W. Kuang, S. Jenkins, E. Coronado, I. V. Grigorieva, I. J. Vera-Marun, E. Navarro-Moratalla, R. F. Evans, _et al._ , “Quantum Rescaling, Domain Metastability, and Hybrid Domain-Walls in 2D CrI3 Magnets,” Adv. Mater. 33, 2004138 (2021).
  * Chen  _et al._ (2021) L. Chen, J.-H. Chung, M. B. Stone, A. I. Kolesnikov, B. Winn, V. O. Garlea, D. L. Abernathy, B. Gao, M. Augustin, E. J. G. Santos, and P. Dai, “Magnetic field effect on topological spin excitations in ${\mathrm{cri}}_{3}$,” Phys. Rev. X 11, 031047 (2021).
  * Akram  _et al._ (2021) M. Akram, H. LaBollita, D. Dey, J. Kapeghian, O. Erten, and A. S. Botana, “Moiré skyrmions and chiral magnetic phases in twisted crx3 (x= i, br, and cl) bilayers,” Nano Lett. 21, 6633 (2021).
  * Xu  _et al._ (2021) Y. Xu, A. Ray, Y.-T. Shao, S. Jiang, D. Weber, J. E. Goldberger, K. Watanabe, T. Taniguchi, D. A. Muller, K. F. Mak, _et al._ , “Emergence of a noncollinear magnetic state in twisted bilayer cri3,” arXiv preprint arXiv:2103.09850 (2021), https://doi.org/10.48550/arXiv.2103.09850.
  * Jiang  _et al._ (2019) P. Jiang, C. Wang, D. Chen, Z. Zhong, Z. Yuan, Z.-Y. Lu, and W. Ji, “Stacking tunable interlayer magnetism in bilayer CrI3,” Phys. Rev. B 99, 144401 (2019).
  * Kim  _et al._ (2019) H. H. Kim, B. Yang, S. Li, S. Jiang, C. Jin, Z. Tao, G. Nichols, F. Sfigakis, S. Zhong, C. Li, _et al._ , “Evolution of interlayer and intralayer magnetism in three atomically thin chromium trihalides,” PNAS 116, 11131 (2019).
  * Lei  _et al._ (2021) C. Lei, B. L. Chittari, K. Nomura, N. Banerjee, J. Jung, and A. H. MacDonald, “Magnetoelectric Response of Antiferromagnetic CrI3 Bilayers,” Nano Lett. 21, 1948 (2021).
  * Song  _et al._ (2019) T. Song, Z. Fei, M. Yankowitz, Z. Lin, Q. Jiang, K. Hwangbo, Q. Zhang, B. Sun, T. Taniguchi, K. Watanabe, _et al._ , “Switching 2D magnetic states via pressure tuning of layer stacking,” Nat. Mater. 18, 1298 (2019).
  * Li  _et al._ (2019) T. Li, S. Jiang, N. Sivadas, Z. Wang, Y. Xu, D. Weber, J. E. Goldberger, K. Watanabe, T. Taniguchi, C. J. Fennie, _et al._ , “Pressure-controlled interlayer magnetism in atomically thin CrI3,” Nat. Mater. 18, 1303 (2019).
  * Sivadas  _et al._ (2018) N. Sivadas, S. Okamoto, X. Xu, C. J. Fennie, and D. Xiao, “Stacking-dependent magnetism in bilayer ${\mathrm{CrI}}_{3}$,” Nano Lett. 18, 7658 (2018).
  * Jiang  _et al._ (2018a) S. Jiang, J. Shan, and K. F. Mak, “Electric-field switching of two-dimensional van der Waals magnets,” Nat. Mater. 17, 406 (2018a).
  * Jiang  _et al._ (2018b) S. Jiang, L. Li, Z. Wang, K. F. Mak, and J. Shan, “Controlling magnetism in 2d cri 3 by electrostatic doping,” Nat. nanotechnol. 13, 549 (2018b).
  * Kashin  _et al._ (2020) I. Kashin, V. Mazurenko, M. Katsnelson, and A. Rudenko, “Orbitally-resolved ferromagnetism of monolayer cri3,” 2D Mater. 7, 025036 (2020).
  * Huang  _et al._ (2018) B. Huang, G. Clark, D. R. Klein, D. MacNeill, E. Navarro-Moratalla, K. L. Seyler, N. Wilson, M. A. McGuire, D. H. Cobden, D. Xiao, _et al._ , “Electrical control of 2D magnetism in bilayer ${\mathrm{CrI}}_{3}$,” Nat. nanotechnol. 13, 544 (2018).
  * Wang and Sanyal (2021) D. Wang and B. Sanyal, “Systematic study of monolayer to trilayer cri3: Stacking sequence dependence of electronic structure and magnetism,” J. Phys. Chem. C . 125, 18467 (2021).
  * Chen  _et al._ (2018) L. Chen, J.-H. Chung, B. Gao, T. Chen, M. B. Stone, A. I. Kolesnikov, Q. Huang, and P. Dai, “Topological Spin Excitations in Honeycomb Ferromagnet ${\mathrm{CrI}}_{3}$,” Phys. Rev. X 8, 041028 (2018).
  * Lee  _et al._ (2020) I. Lee, F. G. Utermohlen, D. Weber, K. Hwang, C. Zhang, J. van Tol, J. E. Goldberger, N. Trivedi, and P. C. Hammel, “Fundamental spin interactions underlying the magnetic anisotropy in the Kitaev ferromagnet CrI3,” Phys. Rev. Lett. 124, 017201 (2020).
  * Vishkayi  _et al._ (2020) S. I. Vishkayi, Z. Torbatian, A. Qaiumzadeh, and R. Asgari, “Strain and electric-field control of spin-spin interactions in monolayer CrI3,” Phys. Rev. Mater. 4, 094004 (2020).
  * Jaeschke-Ubiergo  _et al._ (2021) R. Jaeschke-Ubiergo, E. S. Morell, and A. S. Nunez, “Theory of magnetism in the van der Waals magnet CrI3,” Phys. Rev. B 103, 174410 (2021).
  * Chen  _et al._ (2020) L. Chen, J.-H. Chung, T. Chen, C. Duan, A. Schneidewind, I. Radelytskyi, D. J. Voneshen, R. A. Ewings, M. B. Stone, A. I. Kolesnikov, _et al._ , “Magnetic anisotropy in ferromagnetic ${\mathrm{CrI}}_{3}$,” Phys. Rev. B 101, 134418 (2020).
  * Ellis  _et al._ (2015) M. O. Ellis, R. F. Evans, T. A. Ostler, J. Barker, U. Atxitia, O. Chubykalo-Fesenko, and R. W. Chantrell, “The Landau–Lifshitz equation in atomistic models,” Low Temp. Phys. 41, 705 (2015).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, _et al._ , “QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials,” J. Phys. Condens. Matter. 21, 395502 (2009).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation made simple,” Phys. Rev. Lett. 77, 3865 (1996).
  * Grimme (2006) S. Grimme, “Semiempirical GGA-type density functional constructed with a long-range dispersion correction,” J. Comput. Chem 27, 1787 (2006).
  * Cococcioni and De Gironcoli (2005) M. Cococcioni and S. De Gironcoli, “Linear response approach to the calculation of the effective interaction parameters in the LDA+ U method,” Phys. Rev. B 71, 035105 (2005).
  * Leon  _et al._ (2020) A. Leon, J. González, F. D. C. de Lima, J. Mejia, and E. S. Morell, “Strain-induced phase transition in ${\mathrm{CrI}}_{3}$ bilayers,” 2D Mater. 7, 035008 (2020).
  * Evans  _et al._ (2014) R. F. Evans, W. J. Fan, P. Chureemart, T. A. Ostler, M. O. Ellis, and R. W. Chantrell, “Atomistic spin model simulations of magnetic nanomaterials,” J. Phys. Condens. Matter 26, 103202 (2014).
  * Brown (1979) W. Brown, “Thermal fluctuation of fine ferromagnetic particles,” IEEE Trans. Magn. 15, 1196 (1979).
  * García-Palacios and Lázaro (1998) J. L. García-Palacios and F. J. Lázaro, “Langevin-dynamics study of the dynamical properties of small magnetic particles,” Phys. Rev. B 58, 14937 (1998).
  * Kawaguchi  _et al._ (2007) Y. Kawaguchi, H. Saito, and M. Ueda, “Can spinor dipolar effects be observed in Bose-Einstein condensates?” Phys. Rev. Lett. 98, 110406 (2007).
  * Ezawa (2010) M. Ezawa, “Giant skyrmions stabilized by dipole-dipole interactions in thin ferromagnetic films,” Phys. Rev. Lett. 105, 197202 (2010).
  * Yang  _et al._ (2015) H. Yang, A. Thiaville, S. Rohart, A. Fert, and M. Chshiev, “Anatomy of dzyaloshinskii-moriya interaction at Co/Pt interfaces,” Phys. Rev. Lett. 115, 267210 (2015).
  * Moriya (1960) T. Moriya, “Anisotropic superexchange interaction and weak ferromagnetism,” Phys. Rev. 120, 91 (1960).
  * Tong  _et al._ (2018) Q. Tong, F. Liu, J. Xiao, and W. Yao, “Skyrmions in the moiré of van der waals 2d magnets,” Nano Lett. 18, 7194 (2018).
  * Khoshlahni  _et al._ (2019) R. Khoshlahni, A. Qaiumzadeh, A. Bergman, and A. Brataas, “Ultrafast generation and dynamics of isolated skyrmions in antiferromagnetic insulators,” Phys. Rev. B 99, 054423 (2019).
  * Behera  _et al._ (2019) A. K. Behera, S. Chowdhury, and S. R. Das, “Magnetic skyrmions in atomic thin ${\mathrm{CrI}}_{3}$ monolayer,” Appl. Phys. Lett. 114, 232402 (2019).
  * Abdul-Wahab  _et al._ (2021) D. Abdul-Wahab, E. Iacocca, R. F. Evans, A. Bedoya-Pinto, S. Parkin, K. S. Novoselov, and E. J. Santos, “Domain wall dynamics in two-dimensional van der waals ferromagnets,” Appl. Phys. Rev. 8, 041411 (2021).
  * Allwood  _et al._ (2005) D. A. Allwood, G. Xiong, C. Faulkner, D. Atkinson, D. Petit, and R. Cowburn, “Magnetic domain-wall logic,” Science 309, 1688 (2005).
  * Vélez  _et al._ (2019) S. Vélez, J. Schaab, M. S. Wörnle, M. Müller, E. Gradauskaite, P. Welter, C. Gutgsell, C. Nistor, C. L. Degen, M. Trassin, _et al._ , “High-speed domain wall racetracks in a magnetic insulator,” Nat. Commun. 10, 1 (2019).
  * Luo  _et al._ (2020) Z. Luo, A. Hrabec, T. P. Dao, G. Sala, S. Finizio, J. Feng, S. Mayr, J. Raabe, P. Gambardella, and L. J. Heyderman, “Current-driven magnetic domain-wall logic,” Nature 579, 214 (2020).