Magnetic domain wall dynamics under external electric field in bilayer CrI₃

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 CrI₃ 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 CrI₃ 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).

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 CrI₃ 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 i^th 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.

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 CrI₃ 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$ nm³ 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 i^th 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 i^th 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$ nm² sheet.

In the HT phase of bilayer CrI₃, 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 CrI₃. 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).

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.

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 i^th 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 CrI₃ 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.

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.

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 CrI₃, 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 nm² 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.

Now we turn our attention to another phase of the bilayer system. Four Cr atoms in a unit cell of the LT-phase bilayer CrI₃ 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 CrI₃, 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.

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)

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 CrI₃ 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 CrI₃. 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 CrI₃ is changed by the electric field. One can find the DFT-obtained exchange couplings of different stacks of bilayer CrI₃ 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 CrI₃ 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 CrI₃, 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 CrI₃ Wang and Sanyal (2021).

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 CrI₃, 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 CrI₃ 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.

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.

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.

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.

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.

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.

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.

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 CrI₃.

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 CrI₃. 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.

The spin dynamics in the twisted bilayer CrI₃ 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 CrI₃ 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 CrI₃ 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 CrI₃ 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 CrI₃ can be tuned by the electric field, which is consistent with the experimental results Xu et al. (2021).