Magnetic ground state of monolayer CeI₂: Occupation matrix control and DFT+U calculations

The two-dimension(2D) magnets exhibit magnificent potential in spintronic devices and magnetic storage devices due to their unique properties and the nanoscale size. Although Mermin-Wagner theorem Mermin denies the possibility of long-range magnetic order in 2D systems at non-zero temperature, a series of experimental studies Jiang2021 ; Huang2017 ; Tian2019 ; Cai2019 ; Zhang2019 ; Banerjee2017 ; Zhou2019 ; Bonilla2018 ; Li2018 ; Hara2018 ; Gong2017 ; Wang2016 ; Kim2019 ; Flem1982 ; Fei2018 ; May2019 have confirmed the presence of magnetic orders in 2D materials. These findings encourage people to further explore 2D magnetic materials with superior properties, such as robust dynamic stability, high magnetic transition temperature, and unique magnetic anisotropy.

The magnetic ground state (GS) of a 2D material is significant. It not only dominates the behaviour about magnetism at low temperature, but decides the low-energy excitation of the spin waves at finite temperature. For magnetic crystals, a comprehensive description of the magnetic GS should contain the information of two aspects: the specific magnetic structure and the electronic states of the magnetic sites. For the former part, spin-polarized scanning tunneling microscope (SP-STM), Ramam spectroscopy, and the observation of anomalous Hall effect have been adopted to examine the magnetic order of 2D magnetic insulators and 2D magnetic metals Huang2017 ; Gong2017 ; May2019 ; Lee2016 ; Tian2016 ; Deng2018 . For the latter part, x-ray magnetic circular dichroism (XMCD) and x-ray magnetic linear dichroism(XMLD) are two techniques which help to obtain the local electronic configuration Liu2016 ; Laan1999 . Nonetheless, exploring magnetic ground states of 2D materials experimentally is still challenging. The limitations include the acquirement of high quality samples, strict experimental conditions, and even expensive costs on economy and time. From this point of view, exploring 2D magnetic materials by using density functional theory (DFT) calculations is another feasible scheme. The DFT tool gives the prediction of the magnetic ground state without synthesizing actual samples. It has been applied to study a series of 2D magnetic materials such as layered $\alpha$-RuCl₃, VSe(Te)₂, CrGe(Si)Te₃, etc. Sarikurt2018 ; Kim20161 ; Sandilands2015 ; Li2014 ; Ataca2012 ; Chen2020 ; Fuh2016 ; Pan2014 ; Kan2014 ; Lin2017 ; Lin2016 For 2D CrI₃ and CrBr₃, DFT calculations give quite consistent prediction of magnetic features with the experiments, including the magnetic structure and Curie temperature ($T_{\text{C}}$) Wang20162 ; Zhang2015 ; Liu20162 ; Jiang20212 . With the assistance of DFT calculations, the studying progress of 2D magnetic materials has been accelerated significantly.

In 2D magnetic insulators and semiconductors, the spin magnetic moments are commonly contributed by polarized electrons from partially filled $d$ or f atomic shells. These electrons are localized around the lattice sites and need extra corrections for strong electronic correlation. Therefore, the density functional theory plus Hubbard $U$ (DFT+U) scheme Liechtenstein ; Dudarev is widely applied to give reasonable electronic structures of the correlated systems. In fact, the conventional DFT+U computational scheme faces challenges in locating the ground state of the localized electrons due to the presence of metastable states, especially when the d and f shells are less filled. In other words, although a specific magnetic ordered state can be determined by DFT+U calculations, however, the obtained electronic state of the localized electrons may be not identical for different researchers. This phenomenon has occurred in the DFT+U studies of Ti-based systems Watson , Ce-based systems Watson ; Shick , U-based systems Miskowiec ; Dorado ; Dorado1 ; Dorado2 ; Matthew ; Claisse , cubic fluorite PuO₂ Jomard ; Hou2022 , and monolayer VX₃ (X=Cl, Br, I) hou2024 etc., and the metastable states can originate from 3$d$, 4$f$, or 5$f$ localized electrons (not found in 4$d$ or 5$d$ electronic systems yet, to our knowledge). Moreover, these works have also shown that the occurrence of trapping into metastable states is not relevant to the form of adopted exchange-correlation functional. Local density approximation (LDA) Shick , generalized gradient approximation (GGA) Watson ; Miskowiec ; Dorado ; Dorado1 , strongly constrained appropriate normed (SCAN) type meta-GGA Hou2022 , and HSE06 type hybrid Jollet ; Ratcliff functionals have been adopted in the studies and this accident still happens constantly. To solve the problem, a few techniques have been applied including U-ramping Claisse ; Meredig , quasi-annealing Geng , occupation matrix control (OMC) Watson ; Dorado ; Claisse ; Amadon ; Zhou2011 ; Zhou2022 , DFT+DMFT Amadon2012 and even the combination of them Rabone . These techniques help to reduce the probability of trapping into metastable states to some extent. However, it remains an open question as to how the ground state of the localized electrons can be strictly obtained.

Recently, the monolayer CeI₂ was predicted to be an in-plane ferromagnetic (FM) semiconductor with a high $T_{\text{C}}$ of $374$ K by DFT+U calculations Sheng . The dynamic stability of the monolayer CeI₂ was supported by both phonon calculations and a molecular dynamics simulation at room temperature. In the monolayer CeI₂, two 6s electrons of the Ce atom move away to form the Ce-I bond, leaving a $4f^{1}5d^{1}$ electronic configuration for the Ce²⁺ ion. According to above analysis, a number of metastable states originated from the highly localized $4f$ electrons should be expected. Owing to the importance for both basic physics and potential applications, it is necessary to provide a systematic analysis of the magnetic GS of monolayer CeI₂.

In this study, we perform DFT+U calculations combined with OMC to investigate the magnetic GS of monolayer CeI₂. The magnetic structure is once again confirmed to be FM. The FM GS as well as twenty-three FM metastable states of monolayer CeI2 with different local electronic states are identified. Our results show that the computed isotropy exchange parameter and magnetic anisotropy energy are distinguishable for the GS and metastable states. Based on the exchange parameters of the magnetic GS, the $T_{\text{C}}$ is estimated to be about $128$ K by using Monte Carlo simulations. To reasonably analyze the contributions to the total energies, the Coulomb potential energies of the 4f electron with different occupied electronic orbitals in the crystal field are discussed. Furthermore, we take spin-orbit coupling (SOC) into consideration. The magnetic structure of monolayer CeI₂ then becomes FM with both in-plane and off-plane components of magnetic moments and the easy axis of magnetization is coupled with the crystal structure. The single Ce-$4f$ electron occupies another distinguishable electronic orbital compared to the situation without SOC, which illustrates that SOC is essential to describe the magnetic GS of monolayer CeI₂. A fully off-plane FM state is confirmed to be a metastable state of monolayer CeI₂ whose total energy is just slightly higher than the magnetic GS. We anticipate that achieving saturation magnetization in the off-plane direction for monolayer CeI₂ is feasible under appropriate external conditions.

The DFT+U calculations are carried out by employing the Vienna ab-initio simulation package (VASP) code Kresse . OMC is implemented by applying the patch of Allen et al. tailored for VASP Watson . The adopted atomic pseudo potentials are constructed by the project augmented wave (PAW) method Blochl ; Kresse2 ; Blochl2 . The cut-off energy of the plane wave bases is set to $600$ eV. The first Brillouin Zone is sampled by a $\Gamma$-center $11\times 11\times 1$ $k$-point mesh. The convergence of both cut-off energy and $k$-point mesh has been carefully examined. The crystal structures are relaxed until the Feynman-Hellman force on each atom is smaller than $0.005$ eV/Å. The optimizing of the charge density is ended when the energy difference is smaller than $10^{-6}$ eV between the current two steps of the iteration. We make use of the DFT+U method of Dudarev et al. Dudarev with a form of

where $n_{m_{1},m_{2}}^{\sigma}$ is the matrix element of the occupation matrix with spin $\mathit{\sigma}$. The standard DFT part here is handled by the Perdew-Burke-Ernzerhof (PBE) type GGA exchange-correlation functional. Perdew The Hubbard $U$ and Hund $J$ for the Ce-$4f$ electron are set to be $7.47$ eV and $0.99$ eV, respectively, which have been applied in previous study Sheng ; Larson . To be cautious, a discussion is also made to demonstrate the rationality of the $U$ parameter (see Section A of the Supporting Information SI ). Based on the computed Bloch states with plane waves bases, Wannier90 code Mostofi is used to construct the tight-binding Hamiltonian. TB2J code He is used to calculate the Green’s function and then the isotropic exchange parameter for magnetic properties is obtained. The $T_{\text{C}}$ is gained by Monte Carlo simulations based on the heat bath algorithms Miyatake . The visualizations of the crystal structure and spin densities in this work is achieved by VESTA Momma . The VASPKIT code is used to process the data of the DFT calculations Wang2021 .

A quantum state of a single electron can be express as a linear combination of a group of complete and orthogonal bases in the Hilbert space. In the limitation of localized electrons of an isolated atom, the dimension of the space is equal to $2l+1$, where $l$ is the angular quantum number of that electron. In monolayer CeI₂, $l$ is no longer a good quantum number for the $5d$ electron of Ce²⁺ since the $5d$ electron is not well localized. Additionally, the delocalized $5d$ electron has a more plane-wave-like charge density rather than the atom-orbital-like charge density of the $4f$ electron. Hence only the $4f$ electron is accompanied by high risks of falling into metastable states in DFT+U calculations. Its electronic orbital requires careful control and observation.

The occupation matrix is the density matrix in particle-number representation and a $14\times 14$ occupation matrix is sufficient to represent the electronic state of the single $4f$ electron exactly. Here the general set of 4f atomic orbitals is used as the bases, which are $f_{y(3x^{2}-y^{2})^{\uparrow}}$, $f_{xyz^{\uparrow}}$, $f_{yz^{2\uparrow}}$, $f_{z^{3\uparrow}}$, $f_{xz^{2\uparrow}}$, $f_{z(x^{2}-y^{2})^{\uparrow}}$, $f_{x(x^{2}-3y^{2})^{\uparrow}}$, $f_{y(3x^{2}-y^{2})^{\downarrow}}$, $f_{xyz^{\downarrow}}$, $f_{yz^{2\downarrow}}$, $f_{z^{3\downarrow}}$, $f_{xz^{2\downarrow}}$, $f_{z(x^{2}-y^{2})^{\downarrow}}$, and $f_{x(x^{2}-3y^{2})^{\downarrow}}$, respectively. The spin-quantization axis is chosen to be along $z$ direction of the Cartesian coordinate system. We firstly do not take SOC into consideration. Since there is only one localized electron in the $4f$ shell of Ce²⁺, any quantum state of the single electron with both spin up and spin down components is unphysical. We thus neglect the spin index and a quantum state of the single $4f$ electron can be written as

where $c_{i}$ ($i=1,2,\cdots,7$) is the complex expansion coefficient of the basis. There should be

For simplicity, we mark the quantum state as ($c_{1}$, $c_{2}$, $c_{3}$, $c_{4}$, $c_{5}$, $c_{6}$, $c_{7}$) for the following part of this paper. The eigenvectors and eigenvalues of the occupation matrix are actually quantum states and the corresponding occupation numbers, respectively.

An OMC procedure is to set an initial occupation matrix artificially for the localized electrons in the DFT+U calculations. The initial occupation matrix has also been called the starting point Zhou2022 . Different starting points may access to different final states of the localized electrons. In our calculations, we make our starting points for the Ce-$4f$ electron all diagonal with a average occupation on several orbital bases. The down spin block and the two off-diagonal blocks of the $14\times 14$ occupation matrix are null matrices. The up spin block of the occupation matrix has a form of

with

where $n$ is the number of orbital bases whose occupation number is non-zero. The trace of the up spin block is the number of occupying electrons, which is equal to $1$ for the Ce-$4f$ shell here. Hence there is

Based on Eqs. (4)-(6), we have $\sum\nolimits_{i=1}^{7}C_{7}^{i}=127$ different starting points for the DFT+U calculations. After the self-consistent iteration processes based on different starting points, we obtain twenty-four final states of monolayer CeI₂, whose $4f$ electronic states are different with each other. A detailed discussion of the OMC procedure for the SOC situation is provided in Section B of the Supporting Information. SI

The optimized crystal structure of monolayer CeI₂ is shown in Fig. 1(a). Each Ce²⁺ ion is surrounded by six nearest I^- ions which form a regular triangle prism (RTP). The distance between two nearest Ce²⁺ ions is $4.28$ Å which makes the direct exchange from different Ce sites negligible. The main reason that causes the magnetic order should be the superexchange interaction mediated by I-$p$ electrons. Since both Ce-$4f$ and Ce-$5d$ electrons contribute to the superexchange interaction, the total inter-site exchange strength can be decomposed to $J^{d-d}$, $J^{d-f}$, and $J^{f-f}$, as shown in Fig. 1(b). These are the exchange parameter between $5d$ shell and $5d$ shell, $5d$ shell and $4f$ shell, and $4f$ shell and $4f$ shell from different sites, respectively. Meanwhile, with the help of coordinated I^- ions, the $d-f$ hybridization on each Ce site should be considered. This on-site $d-f$ exchange parameter is denoted by $J^{\text{site}}$, as shown in Fig. 1(b). For the inter-site exchange parameters, we only consider the first nearest exchange parameter $J_{1}$ and second nearest exchange parameter $J_{2}$ as shown in Fig. 1(a). Hence the total Hamiltonian of CeI₂ can be express as

with

where $\mathit{\delta}$ is the spin-independent part of the total energy. $i$ and $j$ denote different sites in the unit cell.

All the following discussions for monolayer CeI₂ are based on the fact that the magnetic structure of the GS is FM no matter which orbital the localized electron is occupying. This feature of CeI₂ can be verified by our computed total energies of the magnetic structures, as well as the computed values of the exchange parameters, which will be discussed later. With OMC applied, the relative energies of the GS and the metastable states of FM CeI₂ are listed in Table I. It is noteworthy that the FM state studied by the previous study Sheng is confirmed to be a metastable state (S₂₀) in our work.

By using a model Hamiltonian of spin interactions like Eq. (7) in combination with the DFT total energies of different magnetic structures, one can fit the values of the exchange parameters. However, this scheme is invalid for monolayer CeI₂ because the exchange parameter is not identical for different magnetic structures. A detailed discussion about it is given in Section C of the Supporting Information. SI Here, we make use of the obtained Bloch states by a self-consistent calculation to construct tight-binding Hamiltonian and calculate the exchange parameter by using Green’s function method He . The values of $J_{1}$ and $J_{2}$ for monolayer CeI₂ with different local electronic states are shown in Fig. 1(c). One can find that $J_{1}$ is always positive and $J_{2}$ is always negative for all the states. The absolute value of $J_{1}$ is far greater than $J_{2}$, which favors FM order of monolayer CeI₂. Note that $J_{1}$ is quite sensitive to the $4f$ electronic states and S₀ has the largest $J_{1}$.

The exchange parameters of the ground state S₀ are used to simulate the $T_{\text{C}}$ of monolayer CeI₂. There are $J_{1}=10.7$ meV and $J_{2}=-0.34$ meV. Monte-Carlo method is used and $T_{\text{C}}$ is estimated to be about $128$ K as shown in Fig. 1(d). It is important to note that the presence of metastable states of monolayer CeI₂ can lead to different values of the inter-site exchange parameters, potentially affecting the reliability of the predicted $T_{\text{C}}$ for the material. Therefore, it is also recommended to exclude the metastable states before conducting further calculations for other local-electron systems, in addition to monolayer CeI₂.

In the case of weak SOC (common in $3d$ electronic systems), the orbital state of localized electrons is primarily decided by on-site Coulomb interactions and the crystal field. Thus, SOC can be treat as an additional perturbation. In other words, the SOC energy can be safely gained based on the simulated charge density which does not include SOC effect. The SOC Hamiltonian is written as

where $\mathit{\lambda}$ parameter denotes the strength of SOC. Applying this approximation to FM monolayer CeI₂, the relative energies of different spin directions are computed. As shown in Table I, all the magnetic states have the highest energy when the spin is along $z$ direction. The easy axis of spin polarization varies depending on the specific orbital state of the 4f electron. For S₀ (the GS), both $x$ and $y$ directions are the easy axes which are $390$ $\mu$eV lower in energy than the $z$ direction. Hence, monolayer CeI₂ is recognized as an isotropic in-plane FM material within the theoretical framework of weak SOC.

In order to better understand the origin of the FM order in monolayer CeI₂, the components of $J_{1}$ are computed and listed in Table I. It is observed that the values of $J^{d-d}$, $J^{d-f}$, and $J^{f-f}$ are always in the order of $J^{d-d}>J^{d-f}>J^{f-f}$ for different states. $J_{1}$ is primarily contributed by $J^{d-d}$, which is similar to the situation in GdI₂ (Gd²⁺: $4f^{7}5d^{1}$) Wang2020 . According to Eq. (7), $J^{\text{site}}$ can be obtained from the energy difference of the two spin configurations $4f^{\uparrow}5d^{\uparrow}$ and $4f^{\uparrow}5d^{\downarrow}$ for each state, which is denoted as $J^{\text{site}}=1/2\left(E_{\text{tot}}^{\uparrow\uparrow}-E_{\text{tot}}^{\uparrow\downarrow}\right)$. The computed values of $J^{\text{site}}$ for different states are listed in Table I. The positive values of $J^{\text{site}}$ helps to stabilize the saturated spin moment of Ce²⁺, i.e. $4f^{\uparrow}5d^{\uparrow}$ spin configuration is favorable in energy for monolayer CeI₂ rather than $4f^{\uparrow}5d^{\downarrow}$ spin configuration, although the spin singlet state (¹G${}_{0}^{4}$) has been recommended by laser spectroscopy for Ce atoms Worden . The magnetic moment is 2 $\mu_{B}$ per Ce coming from the two parallel spins. Furthermore, S₀ has the largest value of $J^{\text{site}}$, which is one of the reason why it leads to the lowest total energy referring to Eq. (7).

The spin densities and electronic structures of S₂₀ and S₀ are illustrated in Fig. 2. The significant difference between the two states can be noticed by inspecting the spin densities. Meanwhile, S₀ exhibits more low-energy contributions to the bands from $4f$ orbitals than that of S₂₀. The $E-k$ dispersion of the $4f$ electron is much more weak for S₀, almost forming a flat band at $-4.1$ eV. Since there is $d-f$ hybridization on each Ce site mediated by I^- ions, the $5d$ electronic states are slightly affected by the $4f$ electronic states as well. The energy gaps of S₀ and S₂₀, which are caused by the spin splitting of the $5d$ electrons, are $386$ meV and $336$ meV, respectively. In fact, there exists another degenerate ground state S${}_{0}^{+}$, whose $4f$ electronic state is (0.66, 0, 0, 0.74, 0, 0, 0). This state is recognized to show exactly the same magnetic properties as S₀. The detailed discussion about S${}_{0}^{+}$ is given in Section F of the Supporting Information. SI The spin densities and electronic structures of all the other metastable states are also shown in Section F of the Supporting Information. SI

In order to have a better understanding of the ground state of FM monolayer CeI₂, we calculate the electronic orbital energies of the single $4f$ electron when occupying different orbitals in the RTP crystal field. This so-called orbital energy is originated from pure Coulomb interactions. As shown in Fig. 3(a), the single $4f$ electron is located at $O$, with six surrounding anions located at $R_{i}$ ($x_{i}$, $y_{i}$, $z_{i}$) ($i=1,2,\cdots,6$) in the Cartesian coordinate system. The height of the RTP is denoted as $l_{0}$ and the edge length of the bottom surface is denoted as $a_{0}$. $l_{0}$ and $a_{0}$ are taken to be $4.24$ Å and $4.28$ Å, respectively, which are based on the relaxed monolayer CeI₂ by DFT calculations. The orbital energy can be expressed as

where $\psi_{4f}$ is the wave function of the single $4f$ electron and $V_{\text{CF}}$ is the crystal field potential energy with a from of

in which the constant term has been set to be unity for simplicity in the calculations. The charge densities of the six I^- ions are treated as point charges. Since the single $4f$ electron of CeI₂ is localized, the numerical integration for Eq. (10) is calculated within a integral area $\Omega$:

where $\Omega$ denotes the spherical area $\sqrt{x^{2}+y^{2}+z^{2}}\leq R$. The integral radius $R$ is set to be $1.2$ times of the Wigner-Seitz radius of Ce to reasonably describe both the local character and the $f\mathtt{-}p$ hybridization of the $4f$ electron.

The orbital energies of all the identified $4f^{1}$ electronic states listed in Table I are computed. As shown in Fig. 3 (b), the $4f$ electronic state of S₀ has the lowest orbital energy in the simulated RTP crystal field environment of the relaxed monolayer CeI₂. Now besides the largest value of $J^{site}$, the RTP crystal field in monolayer CeI₂ also helps to stabilize the ground state S₀ by attaching the lowest $4f$ electronic orbital energy. In the RTP crystal field with D${}_{{}^{3h}}$ point group symmetry, the original seven-fold $4f$ orbitals split into three non-degenerate levels $A_{1}^{\prime}$, $A_{2}^{\prime}$, $A_{2}^{\prime\prime}$ and a pair of two-fold levels $E^{\prime}$ and $E^{\prime\prime}$. The degeneracy of the symmetry-protected orbitals is well-reflected in Fig. 3(b). The $4f$ orbital states are $f_{xyz}$, $f_{z(x^{2}-y^{2})}$, and $-0.92f_{xyz}+0.4f_{z(x^{2}-y^{2})}$ for S₁₅, S₁₆ and S₁₇, respectively. As $f_{xyz}$ and $f_{z(x^{2}-y^{2})}$ are belong to degenerate $E^{\prime\prime}$, the numerical results of crystal orbital energy are the same for these two orbitals along with their linear combination. The similar degeneracy occurs in S₄ and S₅, or S₁₂ and S₁₃, which originates from the degenerate $E^{\prime}$: $f_{yz^{2}}$ and $f_{xz^{2}}$. What’s more, the degenerate $4f$ crystal field orbitals also access to exactly equal on-site $d-f$ exchange energy and total energy for monolayer CeI₂ referring to Table I. This feature reflects the strong constrain from orbital symmetry for the system.

We explore how the orbital energies are influenced by biaxial strain of monolayer CeI₂. We examine the seven general $4f$ atomic orbitals and the $4f$ electronic orbital of S₀, as shown in Fig. 4 (c). The degeneracy of two groups of orbitals $E^{\prime}$: $f_{yz^{2}}$&$f_{xz^{2}}$ and $E^{\prime\prime}$: $f_{xyz}$&$f_{z(x^{2}-y^{2})}$ are maintained under biaxial strain, since biaxial strain does not change the D_3h point group symmetry of monolayer CeI₂. The $4f$ electronic orbital of S₀: $0.66f_{y(3x^{2}-y^{2})}-0.74f_{z^{3}}$ keeps its lowest orbital energy until the compressive strain reaches 4% of the lattice parameter. Fig.3 (c) reveals the change of relative energies among different $4f$ atomic orbitals caused by strain of the crystal structure. Although the $4f$ orbital energy itself does not dominate the local-electron ground state of monolayer CeI₂, we expect potential strain-caused phase transition of local electronic state in magnetic $3d$ transition systems, where the crystal field plays a more important role than $4f$ electronic systems.

Since OMC is applied, the DFT total energy of the system with the $4f$ electron occupying arbitrary orbital is able to be gained as long as the corresponding $4f$ electronic state is metastable. Due to the high symmetry of the crystal field in monolayer CeI₂, the seven general $4f$ atomic orbitals are almost metastable. The corresponding DFT total energies are excerpted from Table I into Fig. 3 (d). Fig.3 (d) shows the feasibility of simulating crystal field excitation by using DFT plus OMC for crystals. This approach includes electronic interactions from the whole crystal field environment, such as non-nearest ligand field and orbital characters of the ligands, which is not that of point charges like in a simple crystal field model. Thus, the simulated crystal field excitation energy is expected to agree with experiments better than model calculations.

So far three kinds of electronic interactions have been verified to have influence on the FM ground state of monolayer CeI₂ in this paper. They are inter-site exchange interaction, on-site $d-f$ exchange interaction and $4f$ crystal field Coulomb interaction. In fact, it is natural to take the crystal field energy of the $5d$ electron into account. However, with both local and itinerant characters, the simulations for the $5d$ electron are complicated and are beyond our study. A qualitative discussion about the $5d$ electronic states in monolayer CeI₂ is given in Section E of the supportting information.

In the more realistic situation, the electronic state of the localized 4f electron in monolayer CeI₂ is not only decided by electronic interactions, but also decided by strong SOC. Thus, the acquirement of the magnetic GS requires fully unconstrained noncollinear self-consistent calculations Worden with SOC included in the DFT+U energy functional. Here OMC is once again applied to exclude metastable states. With SOC being considered, the magnetic ground state of monolayer CeI₂ still favors FM order, but the magnetic moments are with components of both $y$ and $z$ directions, as shown in Fig. 4(a). The easy axis of magnetization is coupled with the crystal field, which makes the magnetic moments point toward certain directions of the crystal structure. The ratio of the off-plane ($z$) and in-plane ($y$) components of the total magnetic moment is $0.65:1$. The net magnetic moment for each Ce is about $1.92$ $\mu_{B}$ which is contributed by both spin magnetic moment and orbital moment.

The DFT calculations also recommend a 12-fold degenerate magnetic GS for monolayer CeI₂ with SOC included, as displayed in Fig. 5(b). All thedegenerate states can be related to the FM GS shown in Fig. 4(a) by the symmetry operations of D_3h point group, i.e. {$E$, $2C_{3}$, $3C_{2}^{\prime}$, $\sigma_{h}$, $2S_{3}$ and $3\sigma_{v}$}. Due to SOC, the symmetry operations act on both electronic orbitals and spins. Thus these degenerate states can be clearly distinguished by the different orientations of magnetization. Except for the original FM GS, the $C_{3}$ and $C_{2}^{\prime}$ rotation operations give birth to another five degenerate states with identical chirality of the electronic orbital. This identical chirality is marked by “$+$” in Fig. 5(b). The $\sigma_{v}$ inversion operations then double the number of degenerate states, and produce six more states with opposite chirality of the electronic orbital. This opposite chirality is marked by “$-$” in Fig. 5(b). We expect this 12-fold degenerate magnetic GS of monolayer CeI₂ could be experimentally verified by magnetic measurements.

The SOC bands of monolayer CeI₂ is shown in Fig. 4(c). The two spins from $4f$ and $5d$ shells are still parallel to benefit the total energy. No significant and qualitative changes are detected for the electronic structure near the fermi level when compared with the situation without SOC. In spite of this, SOC is essential to reasonably describe the magnetic anisotropy in monolayer CeI₂. As shown in Fig. 4(d), with the strong enough external magnetic field along $x$, $y$, and $z$ direction, respectively, the total energies of the corresponding magnetized states are $5.53$ meV, $8.60$ meV, and $0.78$ meV higher per Ce than the magnetic GS. These energy differences are one order larger than the situation that SOC is treated as a perturbation. The variation of the electronic orbital contributes to these energy differences, in addition to the SOC energy introduced by Eq. (9). Hence the conclusion that monolayer CeI₂ is an isotropic in-plane FM material is no longer satisfied when getting away from the framework of weak SOC which is discussed in Section IV.

The fully off-plane ($z$ direction) FM state is concerned due to its potential application in spintronic devices. A comparison between this FM metastable state and the FM GS is shown in Table II. The included angle of the magnetic moments from the two FM states is 57^∘. The fully off-plane FM state is 0.78 meV higher per Ce in total energy that the FM GS. Apparently an appropriate external magnetic field can help to achieve saturated magnetization in the off-plane direction from the GS. We also expect other measures that can tune the magnetic characters of monolayer CeI₂ such as electric field or strain.

In this work, the ab-initio DFT+U calculations combined with OMC are performed to investigate the magnetic GS of monolayer CeI₂. The FM GS and FM metastable states are identified. It is shown that the metastable states have different magnetic properties with the magnetic GS. We recommend that something ought to be done to prevent the metastable states from producing errors when predicting the magnetic properties of such local-electron systems by DFT calculations. The $T_{\text{C}}$ of monolayer CeI₂ is simulated to be about $128$ K based on the magnetic ground state, not above room temperature as mentioned before. The calculations of on-site $d-f$ exchange parameter and $4f$ electronic orbital energy imply that the GS of monolayer CeI₂ with FM order requires the harmony of multiple types of electron-electron interactions to minimize the total energy. With SOC included, the easy axis of magnetization is found to be coupled with the crystal structure. The FM GS with both in-plane and off-plane components of magnetic moments, and a fully off-plane FM metastable state of monolayer CeI₂ are confirmed to be close in total energy. Our work shows it is necessary to absorb SOC into the energy functional to allow SOC to affect the electronic state of the local electrons during the self-consistent calculations. Only in this way could the reliable magnetic GS of monolayer CeI₂ be obtained and the magnetic anisotropy be reliably described.

For other systems with local-electron magnetism, the identification of the local electronic states are also significant. However, the phenomenon of trapping into metastable states during the calculations is quite common and it causes confusions if no extra information is given. Here it is recommended to give the specific representation of the computed local electronic state in one’s research. This action helps to enhance the repeatability and normalization of the computational studies for magnetic materials. (All the representations of the used starting points and corresponding final states in our work have been shown in Section G of the Supporting Information. SI )