Using fully first-principles non-collinear self-consistent field density functional theory (DFT) calculations with relativistic spin-orbital coupling effects, we show that, by applying an out-of-plane electrical field on a free-standing two-dimensional chromium tri-iodide (CrI_{3}) ferromagnetic monolayer, the Néel-type magnetic Skyrmion spin configurations become more energetically-favorable than the ferromagnetic spin configurations. It is revealed that the topologically-protected Skyrmion ground state is caused by the breaking of inversion symmetry, which induces the non-trivial Dzyaloshinskii-Moriya interaction (DMI) and the energetically-favorable spin-canting configuration. Combining the ferromagnetic and the magnetic Skyrmion ground states, it is shown that 4-level data can be stored in a single monolayer-based spintronic device, which is of practical interests to realize the next-generation energy-efficient quaternary logic devices and multilevel memory devices.

## I. INTRODUCTION

According to the Mermin-Wagner theorem, the long-range magnetic order and intrinsic ferromagnetism are forbidden in two-dimensional materials. Notably, recent experimental and theoretical studies have revealed that the magnetocrystalline anisotropy, which stems from the relativistic spin-orbital coupling (SOC) effects, could remove this restriction and permit intrinsic ferromagnetism in two-dimensional van der Waals monolayers.^{1–3} These important discoveries have offered tremendous new opportunities to create novel monolayer-based spintronic memory and logic devices.^{4}

The traditional spintronic devices are mainly based on multilayer heterostructures. For instance, the operations of the magnetoresistive devices (e.g., the spin transfer torque memory devices) rely on the quantum mechanical tunneling through the sandwich-like magnetic tunnel junctions;^{5,6} the operations of magnetic Skyrmion electronic devices rely on the breaking of inversion symmetry at the interface between the magnetic films and the heavy metal materials.^{7,8} These traditional multilayer-based devices are primarily binary (spin up/down), i.e., only one bit (“0”/“1”) can be stored per device. Furthermore, it is extremely challenging to downscale the feature size of these traditional multilayer-based spintronic devices into the sub-10 nm range.^{5–8}

By contrast, the emerging monolayer-based spintronic devices rely on the strong SOC effects in the ligand atoms inside the monolayer.^{3} In a ferromagnetic CrI_{3} monolayer, for instance, the strong SOC effects lead to magnetocrystalline anisotropy, as shown in the approximate XXZ Hamiltonian^{3,9}

where *S*_{i} is a three-dimensional vector to represent the spin on the *i*^{th} Cr site; |*S*_{i}|=3/2; *J* is the Heisenberg isotropic symmetric exchange coefficient; *γ* is the anisotropic exchange coefficient; ** D** is the Dzyaloshinskii–Moriya interaction (DMI) vector which represents the anti-symmetric exchange effects; and the summation is performed for the first nearest neighbor Cr-Cr pairs. As shown in Fig. 1(a), the magnetocrystalline anisotropy energy (MAE) is much lower along the

*z*-axis than in the

*xy*-plane, which is defined by the CrI

_{3}monolayer plane. The MAE counteracts the thermal fluctuation to align the spins along the

*z*-axis, establishing long-range magnetic order and intrinsic ferromagnetism.

^{1–3}

This new physical mechanism, as pointed out by Lado et al,^{3} may lead to electric-field-induced skyrmionic ground state whose magnonic Hamiltonian is topologically non-trivial, which is potentially useful to realize novel devices. Inspired by these findings,^{3} we focus on analyzing the topology of spin configurations of the CrI_{3} monolayers in this paper. In Section II, the simulation methodologies and results are shown. In Section III, the technological importance of the results is discussed. In Section IV, a conclusion is made.

## II. SIMULATIONS

When computing the angle-dependence of MAE, all spins on the Cr sites are restricted to a certain direction (*ϕ*, *θ*) in the non-collinear self-consistent field DFT calculations with the relativistic SOC effects (NC-SOC-DFT). The relative energy difference between the in-plane spin configuration (*ϕ* = 90°) and the out-of-plane spin configuration (*ϕ* = 0°) is the easy-axis MAE value, which is computed to be about 0.82 meV per Cr atom. The calculations are performed using the OpenMX package,^{10} based on the linear combination of pseudo-atomic orbital (PAO) with s2p2d1 basis set and 6.0/7.0 a.u. cutoff radius for Cr/I atoms. Both the relativistic SOC and the non-collinear spin effects are included. The Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional is used.^{11} The plane wave cutoff is set as 200 Ry. The supercell size in the *z*-direction is set as 4 nm, leaving enough vacuum to separate the CrI_{3} monolayers due to the artificial periodic boundary condition in the *z*-direction. As shown in Fig. 1(a), 2 Cr atoms and 6 I atoms are included. The 8×8×1 ** k**-points are used to sample the Brillouin Zone, which is adequate for Brillouin Zone integration.

^{3,9}The atomic lattice relaxation continues until all residue atomic forces are smaller than 10

^{−4}Hartree/Bohr (about 5 meV/Å).

To validate the NC-SOC-DFT methodology mentioned above, the calculation results are compared against the experimental measurements and simulation results in literature. The equilibrium lattice constant is computed to be *a*_{0}=6.880Å, which agrees with the value 6.867Å measured in Ref. 12. The Heisenberg isotropic symmetric exchange coefficient is computed to be *J*=−2.03 meV, which is close to the calculated value −2.2 meV in Ref. 3. The MAE is computed to be around 0.82 meV per Cr atom, which is close to the values 0.6−0.7 meV calculated in Refs. 3 and 13. The anisotropic exchange coefficient |*γ*| is computed to be about 0.2 meV per pair of Cr atoms, which is much smaller than |*J*| and in line with the calculation conclusions shown in Ref. 3. By using the Metropolis Monte Carlo method described in our prior work^{9} and the computed coefficients mentioned above, the Curie temperature is computed to be about 40K, which is close to the value 45K measured in Ref. 1; and the coercive magnetic field is computed to be 0.14-0.16 T, which is within the range 0.05-0.2 T measured in Ref. 1.

By using this validated DFT methodology, the lattice and electronic responses of the off-plane electrical field, ** F**, are studied below by adding a saw-like potential along the

*z*-direction. Under the electric field, the redistribution of electron density and the related dipole effects are automatically accounted for in the NC-SOC-DFT calculations. As shown in Fig. 1(b), within the range of our interests,

**has very small impact on the MAE profile. As shown in Fig. 2,**

*F***also has very small impact on the nearest Cr-Cr distance,**

*F**d*

_{Cr-Cr}, and the lattice constant

*a*(

*E*). For example, when |

**|=0.2 V/nm, |(**

*F**a*(

*E*)–

*a*

_{0})/

*a*

_{0}| is only about 0.06%. However,

**has significant impacts on the Cr-I bond length. When**

*F***is applied along the +**

*F**z*direction, the negatively charged I atoms move toward the −

*z*direction and the positively charged Cr atoms move toward the +

*z*direction. This significantly alters the distances between the Cr-plane and the two I-planes, i.e.,

*d*

_{i}(

*E*) where

*i=*1,2. The relative changes |(

*d*

_{i}(

*E*)–

*d*

_{i}(0))/

*d*

_{i}(0)| are very large (approximately 3.6% when |

**|=0.2 V/nm), indicating the**

*F***-induced breaking of inversion symmetry.**

*F*When |** F**|=0,

*d*

_{1}=

*d*

_{2}so the relativistic SOC effects from the two I-planes cancel each other, leading to a vanishing DMI in Eq. (1). When |

**|≠0,**

*F**d*

_{1}≠

*d*

_{2}so the DMI in Eq. (1) becomes non-trivial. As quantitatively analyzed in our prior work,

^{9}the DMI magnitude in Eq. (1) is |

**| ≈**

*D**k*

_{DMI}× cos

^{−1}

*α*× |

**|, where**

*F**k*

_{DMI}≈ 0.4 meV/(V/nm) as shown in Fig. 4 of Ref. 9 and

*α*≈ 58° is the angle between

**vector and**

*D**z*-axis as shown in Fig. 1 of Ref. 9. This corresponds to 0.5 ×

*c*

_{1}×

*c*

_{2}× |

**|/**

*D**A*≈ 0.36 mJ/m

^{2}when |

**| = 0.2 V/nm and 1.78 mJ/m**

*F*^{2}when |

**| = 1.0 V/nm, where each Cr atom is associated with**

*F**c*

_{1}=3 Cr-Cr bonds; each unit cell has

*c*

_{2}= 2 Cr atoms and an in-plane area

*A*= 0.5 ×

*a*

^{2}× sin

*θ*;

*a*≈ 6.86 Å is lattice constant; and

*θ*= 120° is the angle between the two in-plane lattice vectors. While the MAE tends to align the spins in a parallel manner along the

*z*-axis, which favors the ferromagnetic ground state, the DMI tends to tilt the spins, which favors the spin-canting ground state. As |

**| increases, the MAE is largely kept unchanged (Fig. 1), but the DMI is enhanced due to breaking of inversion symmetry (Fig. 2).**

*F*^{9}The consequent impacts on the topological spin configuration deserve investigation. It worth noting that the first nearest neighbor approximation and the XXZ model approximation are applied in Eq. (1), which may not be accurate enough to model the energetically-subtle spin-canting configurations.

^{13}So, here the fully first-principles NC-SOC-DFT calculations are used to replace the Monte Carlo method,

^{9}to compute the ground states. The NC-SOC-DFT calculations are performed on a big supercell that contains a 240-atom CrI

_{3}monolayer, where the Γ-point is used to sample the Brillouin Zone, while all other NC-SOC-DFT settings are kept the same as mentioned above.

It is not a trivial task to accomplish these NC-SOC-DFT calculations, mainly due to two reasons. Firstly, to account for the non-collinear spin-canting effects, the wave function has to be expressed using a two-component spinor, making the solution of Schrödinger equation computationally demanding. In addition, the non-collinear spin-canting effects make the SCF iterations difficult to converge, especially for the big systems. Secondly, the electronic structure problem to be solved here spans multiple energy scales: (i) the inter-atom bonding energy is roughly at the order of 10^{0} eV; (ii) the exchange effects as expressed in Eq. (1) are approximately at the order of 10^{-4} to 10^{−3} eV; and (iii) the energy cost of spin rotation is typically at the order of 10^{-6} to 10^{-4} eV as shown in Fig. 1. If not treated properly, the energetically-subtle non-collinear spin-canting effects in the scale (iii) could be easily destroyed by the SCF iteration error while converging the energies in the scales (i) and (ii).

In order to reliably capture the energetically-subtle spin-canting effects in the scale (iii), a two-step methodology is adopted to perform the NC-SOC-DFT calculations. To facilitate discussion below, the spin angle of the Cr atom located at the circle center of Fig. 3(a) is denoted as (*ϕ*_{c}, *θ*_{c}); the spin angles of the Cr atoms outside the circle of Fig. 3(a) are denoted as (*ϕ*_{o}, *θ*_{o}); the spin angles of the other Cr atoms are denoted as (*ϕ*_{a}, *θ*_{a}). In the first step, *ϕ*_{c} is constrained as 0°, *ϕ*_{o} is constrained as 180°, and (*ϕ*_{a}, *θ*_{a}) are allowed to change freely. The NC-SOC-DFT calculation is run until the SCF iteration error is reduced below 10^{-6} eV/atom to converge the energies in the scales (i) and (ii). In the second step, all constraints are removed and all spin vectors are allowed to change freely. Starting from the results of the first step, the NC-SOC-DFT iteration is performed until the SCF iteration error is reduced below 10^{-9} eV/atom, which is much smaller than the spin rotation energy cost in the scale (iii). The converged spin configuration is shown in Fig. 3(b), which is a Néel-type magnetic Skyrmion. Since *ϕ*_{c} ≈ 0°, this spin configuration is denoted as “upSkm”. To explore the symmetry properties, *ϕ*_{c} is constrained as 180° and *ϕ*_{o} is constrained as 0° in the first step to run this two-step NC-SOC-DFT calculation again, while keeping all other settings and procedures the same as those in the “upSkm” calculation. The converged spin configuration is shown in Fig. 3(c), which is denoted as “dwSkm” below.

To compute the ferromagnetic spin configurations, the initial condition is set such that all spins in the 240-atom supercell point in the +*z* direction. Then, without constraining spin directions, the NC-SOC-DFT calculation is performed until the SCF iteration error is reduced below 10^{-9} eV/atom. The converged spin configuration is denoted as “upFM”. By setting the initial conditions such that all spins point in the –*z* direction, the similar NC-SOC-DFT calculation is run again. The converged spin configuration is denoted as “dwFM”. As shown in Fig. 4(a), the “upFM” and “dwFM” states are degenerate. So, the “dwFM” total energy, *E*_{dwFM}, is used as the reference energy below. Compared to the “upFM” and “dwFM” states, the antiferromagnetic (“AFM”) state always has a higher energy, as shown in Fig. 4(a). Here, the “AFM” state is computed using the NC-SOC-DFT, by setting the initial spin directions as anti-parallel along the *z*-direction between any pair of the nearest Cr atoms.

As shown in Fig. 4(a) and Fig. 4(b), when |** F**|<

*F*

_{0}, the degenerate “upFM” and “dwFM” states are more energetically-favorable than the degenerate “upSkm” and “dwSkm” states, indicating the ferromagnetic ground state. When |

**|>**

*F**F*

_{0}, however, the “upSkm” and “dwSkm” states become more energetically-favorable than the “upFM” and “dwFM” states, indicating the magnetic Skyrmion ground state. The

*F*

_{0}, which is computed to be about 0.14 V/nm as shown in Fig. 4(a), is the threshold electrical field intensity to switch between the ferromagnetic ground state and the magnetic Skyrmion ground state. The energy splitting between the “upFM”/“dwFM” state and the “upSkm”/“dwSkm” state is large. For instance, when |

**|=0.2 V/nm, Δ**

*F**E*

_{upSkm}≈Δ

*E*

_{dwSkm}≈3 meV/Cr, which is much larger than the energy cost of spin rotation (smaller than 1 meV/Cr as shown in Fig. 1). This indicates that the Néel-type magnetic Skyrmions are topologically-protected spin configurations, which is desirable to store information in a stable manner.

## III. DISCUSSION

Due to the quantum confinement and tunneling effects in the ultra-scaled semiconductor devices, the further downscaling of traditional integrated circuit technologies in the sub-10 nm technology nodes is becoming more and more difficult and the “Moore’s Law” is approaching the insurmountable fundamental limit.^{7,14,15} In the “post-Moore era”, it is a still open question that which alternative technology would become the next torch bearer of the information revolution. The electrical-field-induced magnetic Skyrmion ground state of the emerging ferromagnetic monolayers unraveled in this early-stage research offers possible opportunities to develop next-generation monolayer-based spintronic devices, provided that several key properties could be improved. For example, in a CrI_{3} monolayer, the Curie temperature *T*_{c} ≈ 45 K and the threshold electrical field *F*_{0} ≈ 0.14 V/nm. To enable room temperature devices and to alleviate requirements on external circuits, it is necessary to research and develop new monolayer materials with lower *F*_{0} and higher *T*_{c}, e.g., by controlling crystalline structure^{12} and doping.^{16,17} This research direction deserves future attention.

It worth noting that the Skyrmion size hinges on the strength of Heisenberg isotropic symmetric exchange coefficient, MAE, and DMI, etc.,^{18} which are all electrical-field-dependent as quantitatively analyzed in our prior work.^{9} To rigorously quantify the Skyrmion size, enlarged supercells are required to minimize the impacts of the artificial periodic boundary conditions (PBC) in the *xy*-plane. However, as discussed in the two-step NC-SOC-DFT methodology, it is challenging to further enlarge the supercell due to the difficulties to achieve SCF convergence. Therefore, though the NC-SOC-DFT calculations here show that the Skyrmion size is around 1.5 nm, it should only be interpreted as an indication that the Skyrmions can exist in the CrI_{3} monolayers under out-of-plane electrical fields. The electrical-field dependence and the temperature-dependence of the Skyrmion size deserve future research attention, too.

## IV. CONCLUSION

Using the non-collinear self-consistent field density functional theory calculations with the relativistic spin-orbital coupling effects (NC-SOC-DFT), it is shown that the magnetic Skyrmion ground state can exist in the recently discovered ferromagnetic chromium tri-iodide (CrI_{3}) monolayers, under an external out-of-plane electrical field. It is shown that the magnetic Skyrmion is Néel-type. It is revealed that the magnetic Skyrmion is topologically-protected, which is beneficial to enable stable information storage. It is demonstrated that, by combining the degenerate ferromagnetic states and the degenerate magnetic Skyrmion states, 4-level per cell magnetization may be realized which could be of interests for enabling quaternary devices. Despite the above-mentioned interesting properties of the emerging ferromagnetic monolayers, significant future research attention is needed to realize practical monolayer-based spintronic device applications. Specifically, it is important to search new monolayer materials with higher Curie temperature and lower threshold electrical field to create Skyrmion ground states.

## ACKNOWLEDGMENTS

We gratefully acknowledge Dr. Xiaodong Xu and Dr. M. P. Anantram of the University of Washington for the helpful discussions.^{1,4} This research is supported by the “Fundamental Research Funds for Central Universities” from Hunan University of People’s Republic of China. The calculations are performed at the National Supercomputing Center in Changsha, People’s Republic of China. The authors gratefully thank Dr. Xianzhong Duan, Dr. Min He, Dr. Kenli Li, and Dr. Can Leng of Hunan University for their kind help and support.