The celebrated Mermin-Wagner theorem [1] points out that long-range magnetic order does not exist in low-dimensional (d ≤ 2) Heisenberg magnets at T > 0. Surprisingly, two-dimensional (2D) magnetic materials exhibit stable long range magnetic orders which have been discovered in experiments in the past few years. One of the most important reason is that some of the 2D magnets exhibit strong anisotropy [2], which open the Goldstone mode and thus stabilizes the long-range magnetic order. The 2D magnetic materials display various fascinating properties and important applications, such as topological superconductivity in monolayer transition metal dichalcogenides (TMDs) [3], quantum anomalous Hall (QAH) effect in layered materials of MnBi₂Te₄ [4], quantum spin liquid (QSL) state in metal–organic frameworks (MOFs) [5] and strained CrSiTe₃ [6], room temperature multiferroicity [7], and skyrmions in Fe₃GeTe₂ [8]. The emergence of 2D magnetic materials provides an ideal platform and extend playground for the study of low dimensional exotic magnetic phenomena.

In recent years, binary transition metal trihalides TMX₃ such as CrI₃ are extensively studied [9-12]. Similar to the CrI₃, the new vdW layered structure Ⅵ₃ has been paid some attentions recently. Kong [13], Tian [14] and Son et al. [15] synthesized bulk Ⅵ₃ experimentally, and observed that the Ⅵ₃ is a strongly anisotropic hard ferromagnetic semiconductor (band gap is about 0.7 eV). As the temperature increases, the Ⅵ₃ will undergo a magnetic phase transition and a structural phase transition. The Curie temperature T_C of the Ⅵ₃ is about 50 K, and when the temperature is higher than T_S = 80 K, a structural phase transition from R3¯ to C2m [14] will occur. Unlike the CrI₃, the magnetic ground state of the Ⅵ₃ is still controversial. One possible reason is that the t_2g manifold in the Ⅵ₃ is less than half-filled, which presumes that the t_2g manifold has multiple possible occupied states. So far, there have been abundant theoretical research progress on the Ⅵ₃, including monolayer [16-20], bilayer [21-23] and bulk phase [14, 24]. However, the theoretical calculations of the Ⅵ₃ not only contradict experimental results, but different theoretical calculation results also contradict each other. This is because the charge, spin, orbital, and lattice degrees of freedom of the transition metal compounds, especially in strongly correlated systems, are often highly entangled, resulting in computational complexity. The first urgent issue is why density functional theory (DFT) calculations often yield half-metallic metastable results in the Ⅵ₃ [16, 17, 21]. Interestingly, the DFT calculations also confirmed that the Ⅵ₃ is a ferromagnetic Mott insulator [18-20]. Another controversial issue is the magnetocrystalline anisotropy of the Ⅵ₃. Wang et al. showed that the Ⅵ₃ has a negligible magnetic anisotropy energy (MAE) as small as 0.15 meV/f.u. [19]. In contrast, Yang et al. [18] identified that the Ⅵ₃ is a 2D Ising-like ferromagnetic semiconductor with very strong uniaxial anisotropy and MAE is as high as 17 meV/f.u. In general, one possible reason for this large discrepancy is that the LDA + U method allows the presence of many metastable occupation states for the orbitals corrected with the Hubbard U parameter, even in hybrid functionals [25, 26]. Such numerous metastable states make it difficult to find the global minimum of the system through self-consistent calculations. Therefore, it is very important to study the electronic structure of the Ⅵ₃, determine the possible magnetic ground state, and reveal the origin of the strong anisotropy by controlling the orbital occupation.

The crystal structure and three sets of coordinate systems, namely the global coordinate system {xyz}, the right hand local coordinate system {x′y′z′} and the lattice vectors {abc} are shown in Fig. 1. The local coordinate chosen for the octahedron in the Ⅵ₃ monolayer, where the O-x′, O-y′ and O-z′ are along the three Ⅴ-Ⅰ bonds respectively. For the lattice vectors {abc}, the c axis is perpendicular to the ab plane, and the angle between a and b is 120°. In the global coordinate system {xyz}, the xy plane is parallel to the ab plane, the x and z-axes are parallel to the a and c axes respectively, and the c or z axis is exactly the [111] direction of the local coordinate system {x′y′z′}. The space group of the monolayer Ⅵ₃ is P3¯, and the corresponding point group is S₆. The center of the two Ⅴ atoms in each unit cell is the space inversion center, the Ⅴ atoms form a nearly planer honeycomb structure and the bond length of the Ⅴ-Ⅴ bond is 3.95 Å. There are six Ⅰ atoms around the Ⅴ atom, forming a Ⅴ-Ⅰ deformed octahedron. The crystal structure and three coordinate systems of the Ⅵ₃ are shown in Fig. 1. From group theory analysis (see details shown in Supporting information), five symmetry-adapted orbitals are identified, i.e., a_g, egπ+, egπ−, egπ′+ and egπ′−.

A key factor affecting the order of the singlet a_g and the doublet egπ± energy levels is trigonal distortions – elongation or contraction of TMX₆ octahedra along the c axes [27, 28]. In our calculations, the octahedra is the trigonal contraction. Fig. 2 shows the energy level order and orbital occupation of the triangular crystal field (CF) under the elongation and compression, the t_2g levels are here split into a singlet a_g and a doublet egπ± by trigonal distortions. In the ion limit, Ⅴ³⁺ has two d electrons occupying the low-lying t_2g manifold. When the octahedron undergoes trigonal contraction, one of d² electrons preferentially occupy the a_g orbital, and another electron will simultaneously occupy the doublet egπ± orbital, resulting in the so-called half metallic metastable state. On the contrary, the trigonal elongation of the octahedron lead to two d electrons to preferentially occupy the doublet egπ±, and the system will open a band gap. We noticed that Yang et al. reported that the Ⅵ₃ is a robust out-of-plane anisotropic ferromagnetic semiconductor with an orbital magnetic moment (MM) of -1.05 µ_B [18]. Notice that the orbital MM of the doublet egπ± is ∓ 1 μ_B, the possible electronic structure of the Ⅵ₃ is that the d² electrons fully occupy the a_g and the egπ+ orbitals. Therefore, how the two d electrons occupy the t_2g manifold not only determines whether the Ⅵ₃ is a metal or a semiconductor, but also determines the anisotropy of the Ⅵ₃. It is extremely important to study the electronic structure of the Ⅵ₃ by controlling the orbital occupation. Note that the orbits in the local coordinate system {x′y′z′} and the global coordinate system {xyz} connect through a unitary transformation (see details shown in Supporting information). For clarity, orbits in the global coordinate system {xyz} are marked with symbol ˜, and the expression of the density matrix is in the global coordinate system {xyz}.

First, we do not consider spin-orbit coupling (SOC) effect. We tested a series of Hubbard U values. When U = 3.68 eV, the calculated band gap value is consistent with the experimental measurement value [14]. By controlling the orbital occupation, the initial density matrix was fully relaxed by the VASP software package and the electronic structures and total energies of different occupied states were calculated [29-33]. The spin-polarized band structures of the four different occupied states are shown in Fig. 3. It can be seen from the band structures that the three constrained occupation states are all direct band gap semiconductors, and the band gaps are about 0.8 eV. Since the d² electrons fully occupy the two orbitals in the t_2g manifold, the remaining one empty orbital is pushed up to the conduction band, and thus opening up a band gap. But for the unconstrained state, the DFT + U calculation results show a Dirac-like half-metallic band structure [16, 17]. As shown in Table 1, comparing the energy of the four states, we find that the lowest energy occupied state is the a˜ge˜2 state (notice that e˜1,2 are real orbitals obtained from linear combination of e˜gπ± which are complex orbitals for non-SOC calculations). In addition, the total energy of the metallic state is about 0.7 eV higher than that of the a˜ge˜2 state, indicating that the metallic state is a metastable state. We notice that FeO is similar to the Ⅵ₃, that is, it is difficult to open up a band gap by merely increasing the value of Hubbard U [34]. To determine the true ground state of the FeO, one needs to specify the occupation of the electrons in the d orbitals.

The spin-resolved projected density of states (PDOS) of the 4 different occupied states are plotted in Fig. 4. We note that the fully occupied spin-up lower Hubbard band is mainly located around -2 eV below the Fermi level, and the spin-down empty upper Hubbard band is located around 2 eV above the Fermi level in Fig. 4. It is estimated that the Hubbard U is about 4 eV, very close to the 3.68 eV we set in the DFT + U approach. The contribution of the VBM of the a˜ge˜1 and a˜ge˜2 states is mainly from the dz2 orbital of V; while the VBM of the e˜1e˜2 state is the contribution of the dx2−y2 orbital. For the e˜1e˜2 occupancy, d_xy and dx2−y2 are doubly degenerate; d and d are doubly degenerate. Note that the wave functions are linear combination of the d orbitals, and the diagonal elements of the density matrix are the sum of the squares of the expansion coefficients of a certain d orbital in all occupied state wave functions. For example, for the a˜ge˜1 occupied state, ideal contributions of dz2, dx2−y2 and d orbitals are 1, 23 and 13 respectively; the corresponding self-consistent results are 0.87, 0.78 and 0.47, respectively. Since the Ⅴ-Ⅰ bond has both ionic and covalent contributions, a small amount of d and d orbitals also participate in the a˜ge˜1 state [35], and the corresponding matrix elements in the a˜ge˜1 state are 0.18 and 0.27. For the other two occupied states (a˜ge˜2 and e˜1e˜2), the calculated results are also in line with expectations. That is to say, the initial guessed occupied states are very close to the self-consistent results, indicating that the calculation results are reliable.

When the SOC effect is included, reorientation of the spin and the orbital MM will cause a little change in the total energy. In the Ⅵ₃ monolayer, the orbital MM is constrained along the z or -z direction, and the orientation of the spin is not restricted. The interaction of the SOC and CF leads to single ion anisotropy [36]. We adopt the DFT + U + SOC combined with the constrained density matrix method to calculate the anisotropy properties. Consider that there are three possible occupations for the orbital part, namely electrons occupying a˜ge˜gπ+, a˜ge˜gπ− and e˜gπ+e˜gπ−. On the other hand, considering that the spin-up and down along the x, y and z directions respectively, there are 18 possible occupied states. For comparison, we also calculate the case without the density matrix constraints. The energy difference, the spin MM and the orbital MM of 18 possible occupied states are shown in Table 2, Table 3. The bold marked '0.000' in Table 2 represents the energy of the reference state, and the energies of other occupied states respect to the reference state are all positive values. The DFT + U + SOC calculations show that the orbital MM is about ∓0.75 µ_B/f.u., which is very close to the results of the ∓1 µ_B [18]. It can be seen from Table 2, Table 3 that when the spin MM and the orbital MM are antiparallel, Ⅵ₃ monolayer has the lowest energy. The easy axis is along the z direction, and the giant MAE is about 11.49 meV/f.u., which exhibits strong anisotropy. Interestingly, when the spin MM and the orbital MM are parallel, the Ⅵ₃ exhibits in-plane anisotropy and the total energy is about 26 meV higher when the spin is along the z direction. When the electrons occupy the e˜gπ+ and e˜gπ− states, the orbital angular momentum or orbital MM cancels out, there is nearly no net orbital MM and the system is almost isotropic. For reference, the calculation results without any occupancy constraints shown in Table 4. The calculated results show in-plane anisotropy with a higher total energy than that of the occupied constraints.

Additionally, we can use atomic picture to understand the origin of the out-of-plane MAE. The SOC effect of the Ⅴ³⁺ can be explained through the Hund's third rule that the level with the lowest value of the total angular momentum quantum number J lies lowest in energy. Each of the Ⅴ³⁺ ion has only two d electrons, that is, less than half-filled of the t_2g manifold. The orbital angular momentum L is antiparallel to that of the spin angular momentum S resulting the ground state. Taking the two occupied states of the agegπ+ and the agegπ− as examples, all the spin MMs are along the z direction. For each of the Ⅴ³⁺ ion, the DFT + U + SOC calculations show that E_SOC < 0 for the former and E_SOC > 0 for the latter, indicating that the spin MM and the orbital MM are arranged antiparallel. Owing to the orbital MMs of the egπ+ and the egπ− are opposite signs, the SOC effect splitting of the doublet egπ± with an energy of 2λ. This explains why the total energy of the system increases by ~49 meV when the spin MM and the orbital MM are parallel. Our study shows that the DFT calculations yield the metallic results with nearly quenched orbital angular momentum if the orbital occupation is not constrained. Therefore, the orbital occupation constrained calculations is a key skill to show that the Ⅵ₃ is a Mott insulator, and to reveal that the strong magnetic anisotropy.

In summary, by controlling orbital occupancy, we systematically studied the electronic structure and magnetic anisotropy of the Ⅵ₃. Our results suggest that the possible ground state of monolayer Ⅵ₃ is a ferromagnetic semiconductor with strong uniaxial magnetic anisotropy. In the Ⅴ³⁺ ion, the d² electrons occupy the a_g and the egπ+ orbitals resulting in an orbital MM of –0.75 µ_B, which is the origin of the strong magnetic anisotropy. Although we studied the Ⅵ₃ with fixed orbital occupation, this approach is universal and mainly depends on the symmetry of the system; it is easy to extend this approach to study the novel magnetic materials that contain 4d, 5d and f electrons.

Declaration of competing interest

Ke Xu: Writing – original draft, Software, Methodology, Investigation, Formal analysis, Conceptualization. Shulai Lei: Writing – review & editing, Supervision, Software, Resources, Project administration, Conceptualization. Panshuo Wang: Formal analysis. Weiyi Wang: Visualization. Yuan Feng: Writing – review & editing, Visualization. Junsheng Feng: Writing – review & editing.

CRediT authorship contribution statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work is partially supported by the Natural Science Foundation of Hubei Province (No. 2022CFC030), the Science and Technology Research Project of Hubei Provincial Department of Education (No. D20212603) and Hubei University of Arts and Science (No. 2020kypytd002). W. Y. Wang acknowledges the support from National Natural Science Foundation of China (No. 22303098). J. S. Feng acknowledges the support from Anhui Provincial Natural Science Foundation (No. 1908085MA10). We thank Prof. H. J. Xiang and Dr. J. Y. Ni for useful discussions.

Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.cclet.2024.110257.