Exploration of the two-dimensional Ising magnetic materials in the triangular prismatic crystal field

Since the discovery of two-dimensional (2D) ferromagnetism (FM) in atomically thin layers CrI₃¹ and Cr₂Ge₂Te₆², two-dimensional magnetic materials have been the subject of extensive research due to their unique electronic structure and rich tunable properties^{3, 4, 5, 6, 7}. The Mermin-Wagner theorem establishes that 2D isotropic Heisenberg spin systems lack long-range magnetic order at any non-zero temperature⁸. As a result, magnetic anisotropy (MA) becomes critical for stabilizing 2D magnetic systems. Moreover, a huge MA offers greater resistance to thermal perturbations, which is beneficial for ensuring stable data storage and enhancing the reliability of spintronic devices^{1, 2}. Both the CrI₃ monolayer¹ and Cr₂Ge₂Te₆² have a closed $t_{2g}^{3}$ shell for the octahedral Cr³⁺ $S$ = 3/2 ion. As a result, no single-ion anisotropy (SIA) is produced from the orbital singlet. Instead, the finite perpendicular MA in these materials arises from exchange anisotropy, which is induced by spin-orbit coupling (SOC) involving the ligand heavy $p$ orbitals and their hybridization with Cr 3$d$ orbitals^{9, 10}. The search for materials with huge MA, especially 2D Ising magnetic materials, is pivotal for the advancement of spintronics. Their resilience to thermal perturbations enables their use in room-temperature magnetic applications at the 2D limit. This breakthrough is key in developing next-generation spintronic devices, promising to revolutionize information storage and processing with improved efficiency, speed, and miniaturization ^{9, 10, 11, 12, 13, 14, 15, 16}.

Recently, MSi₂N₄ monolayer, where M represents either Mo or W, has been successfully synthesized using chemical vapor deposition techniques¹⁷. These materials show semiconductor characteristics and are notable for their excellent electronic properties as well as their robust environmental stability. Building on this achievement, there has been heightened research interest in the MX₂Y₄ family of layered materials ^{18, 19}. Noteworthy research includes the discovery of unique electronic effects in CrSi₂Y₄ (Y = N, P)²⁰, the electronic and excitonic properties in MSi₂Y₄ (M = Mo, W; Y = N, P, As, Sb) films²¹, the exploration of new stable Janus MX₂Y₄ (M = Sc$\sim$Zn, Y$\sim$Ag, Hf$\sim$Au; X = Si, Ge; Y = N, P) monolayers for potential applications²², the investigation of the electronic structure of MX₂Y₄ (M = Ti, Cr, Mo; X = Si; Y = N, P) bilayers under vertical strain²³, an approach to construct MX₂Y₄ family monolayers with a septuple-atomic-layer structure²⁴, and the prediction of an indirect band-gap semiconductor Janus MSiGeN₄ (M = Mo, W) monolayers²⁵.

In the case of MSi₂N₄, where M adopts a +4 charge state within a unique triangular prismatic crystal environment¹⁷. The five $d$ orbitals split into three distinct sets: a lowest-energy singlet labeled $3z^{2}-r^{2}$, a middle doublet denoted as $xy$/$x^{2}-y^{2}$, and the highest-energy doublet, $xz$/$yz$, as shown in Fig. 1(c). When the SOC is considered, these initially degenerate doublets ($xy$/$x^{2}-y^{2}$ and $xz$/$yz$) further separate into complex $L_{\pm 2}$ and $L_{\pm 1}$ orbitals, respectively. This further orbital splitting, induced by SOC, leads to a huge perpendicular SIA^{9, 10}. Therefore, the unique triangular prismatic crystal field in these MSi₂N₄ monolayers is conducive to achieving huge MA.

In this study, we employ crystal field theory, spin-orbital state analyses, and density functional calculations to examine the MSi₂N₄ (M = V, Cr) monolayers. Our findings indicate that the transition metals (M) exist in a +4 valence state within the triangular prismatic crystal field. However, the strong hybridization between neighboring M⁴⁺ ions interferes with the $d$ orbital splitting, leading to distinct magnetic properties. Specifically, pristine VSi₂N₄ exhibits a magnetic state characterized by the V⁴⁺ ion with a $S$ = 1/2 charge-spin state. Our findings reveal that VSi₂N₄ features an orbital singlet, which results in a small in-plane MA of approximately 2 $\mu$eV per V atom. In contrast, the pristine CrSi₂N₄ monolayer is nonmagnetic, characterized by the Cr⁴⁺ 3$d^{2}$ $S$ = 0 state. However, when the nonmagnetic Cr⁴⁺ is substituted by Si⁴⁺, the Cr${}_{\frac{1}{3}}$Si${}_{\frac{8}{3}}$N₄ monolayer transforms into an antiferromagnetic insulator with Cr⁴⁺ 3$d^{2}$ $S$ = 1 state. Moreover, the Cr⁴⁺ ions display a large orbital moment of –1.06 $\mu_{B}$, oriented along the $z$-axis, as well as a huge perpendicular MA energy of 18.63 meV per Cr atom in the Cr${}_{\frac{1}{3}}$Si${}_{\frac{8}{3}}$N₄ monolayer. Therefore, we demonstrate that the triangular prismatic crystal field serves as a promising experimental platform for discovering novel 2D Ising materials.

The density functional theory (DFT) calculations are carried out using the Vienna ab initio Simulation Package (VASP) ²⁶. The exchange-correlation effect is described by the generalized gradient approximation (GGA) using the functional proposed by Perdew, Burk, and Ernzerhof (PBE) ²⁷. The kinetic energy cutoff for plane-wave expansion is set to 450 eV. A $3\times 3\times 1$ supercell is chosen for MSi₂N₄. The crystal structure of the VSi₂N₄ (CrSi₂N₄) monolayer is illustrated in Fig. 1(a). It can be understood as the incorporation of a 2H phase VN₂ (CrN₂) layer into an $\alpha$-InSe-type Si₂N₂ framework. The resulting composite forms a hexagonal lattice with the $P\overline{6}m2$ space group.

The local V ion has a triangular prismatic crystal field, splitting the degenerate five $d$ orbitals into the lowest $3z^{2}-r^{2}$ singlet, middle $xy$/$x^{2}-y^{2}$ doublet and the highest $xz$/$yz$ doublet. Taking into account the SOC, the degenerate $x^{2}-y^{2}$/$xy$ ($xz$/$yz$) doublet can further split into $L_{\pm 2}$ ($L_{\pm 1}$) orbitals (see Fig. 1 (c)). The complex orbital wave functions can be written as

Thus, the utilization of the triangular prismatic crystal field emerges as a promising avenue for the achieving of 2D Ising magnetization with huge MA.

To better describe the on-site Coulomb interactions of V (Cr) 3$d$ electrons, the typical value of the Hubbard $U_{\rm eff}$ = 4.0 eV is used in the GGA+$U$ calculations ²⁸. The spin-orbit coupling (SOC) is also included in our GGA+SOC+$U$ calculations to study MA. The method for controlling the density matrix is implemented using open-source software developed by Watson ²⁹, and this method has been widely used in previous studies ^{30, 12, 31, 32, 33}. In addition, using the site-projected wave function character of $d$ orbitals in DFT calculations, we can project the density of states (DOS) based on different eigen wavefunctions, which are linear combinations of the five common $d$ orbitals ($xz$, $yz$, $xy$, $x^{2}-y^{2}$, $3z^{2}-r^{2}$). As seen in equation (1), this method can be used in the complex crystal field and the complex orbitals due to SOC, which is essential in our calculations.

We first perform the spin-polarized GGA calculations to study the electronic and magnetic structure of the VSi₂N₄ monolayer. In the FM state, the calculated V local spin moment of 0.86 $\mu_{\rm B}$ refers to the nominal V⁴⁺ 3$d^{1}$ $S$ = 1/2 state. The reduction of V⁴⁺ spin moment is due to the strong covalency with the nitrogen ligands. The total spin moment of 0.92 $\mu_{\rm B}$ per V atom agrees well with the V⁴⁺, Si⁴⁺ and N^3- charge states in the VSi₂N₄ monolayer. We plot in Fig. 2(a) the orbitally resolved DOS for the FM state. We find that the $xz$/$yz$ doublet has the highest energy about 3 eV above the Fermi level. In contrast, the lowest $3z^{2}-r^{2}$ singlet and middle $x^{2}-y^{2}$/$xy$ doublet have a strong hybridization, both of which are on average 1/3 filled to fulfill the formal V⁴⁺ 3$d^{1}$ state. Due to the strong V-N covalency in the VN₂ layer, the adjoining nitrogen becomes negatively spin-polarized and has a local spin moment of –0.04 $\mu_{\rm B}$. The Si₂N₂ layer is not affected by spin-polarized V⁴⁺ ions and remains nonmagnetic.

To account for electronic correlations of the V ions and to incorporate SOC effect, we carry out the GGA+SOC+$U$ calculations with the out-of-plane magnetization axis. As shown in Table 1, the enhanced electron localization gives rise to an increasing local spin moment of 1.12 $\mu_{\rm B}$ for V ion, and a small orbital moment of –0.01 $\mu_{\rm B}$ induced by SOC. The total spin moment of 1.00 $\mu_{\rm B}$ per V atom suggests the formal V⁴⁺ 3$d^{1}$ $S$ = 1/2 state again. A clear 0.2 eV semiconductor energy gap is observed, as seen in Fig. 2(b). The electron of V⁴⁺ ion mainly occupies the lowest $3z^{2}-r^{2}$ singlet, while the higher $x^{2}-y^{2}$/$xy$ doublet remains partially occupied. This indicates that the $3z^{2}-r^{2}$ and $x^{2}-y^{2}$/$xy$ orbitals have a strong hybridization even in localized V ions, and the phenomenon is similarly observed in hexagonal MoS₂^{34, 35, 36, 37, 38, 39}. Then, we assume the magnetization axis along the in-plane and run self-consistent calculations till the energy difference converges within 1 $\mu$eV per V atom. Comparing the total energies of the different magnetizations (out-of-plane and in-plane), we find the VSi₂N₄ has the easy in-plane MA, and the out-of-plane has a higher energy by only 2 $\mu$eV per V atom, which is insufficient to stabilize a 2D Ising magnetization effectively.

To realize a large orbital moment and giant MA, the single electron in the V⁴⁺ 3$d^{1}$ configuration should occupy the $x^{2}-y^{2}$/$xy$ doublet, which makes the SOC active and may eventually determine the 2D Ising magnetization. To address this, we perform the GGA+SOC+$U$ calculations and initialize the system in two distinct electronic configurations: one with the single electron occupying the $3z^{2}-r^{2}$ state and the other in the $L_{-2}$ state. Our results show that the $L_{-2}$ state is unfavorable in the VSi₂N₄ monolayer and ultimately transitions to the $3z^{2}-r^{2}$ and $x^{2}-y^{2}$/$xy$ hybridized orbital state, as visualized in Fig. 2. Therefore, the VSi₂N₄ monolayer has the V⁴⁺ 3$d^{1}$ $S$ = 1/2 configuration favoring a hybridized orbital state over the energetically unfavorable $L_{-2}$ state, resulting in a small in-plane MA of 2 $\mu$eV per V atom.

As we mentioned above, the potential for a large orbital moment and huge MA arises from its unique triangular prismatic crystal field. However, this potential is significantly limited due to the strong hybridization between the lower-energy $3z^{2}-r^{2}$ orbital and the higher $xy$/$x^{2}-y^{2}$ orbitals. In the VSi₂N₄, each vanadium (V) atom possesses five $d$ orbitals that are initially orthogonal within the triangular prismatic crystal field. Yet, the hybridization occurs due to the influence of adjacent V ions, which themselves have a high coordination number of 6 in the triangular lattice. This results in VSi₂N₄ exhibiting only a minimal MA.

In order to explore the feasibility of achieving the $L_{-2}$ state characterized by a large orbital moment and huge MA, we employ calculations within an isolated triangular prismatic crystal field. For this purpose, a $3\times 3\times 1$ VSi₂N₄ supercell comprising 9 V ions is utilized. To prevent magnetic interactions, 8 of the V ions are substituted with Si ions named V${}_{\frac{1}{9}}$Si${}_{\frac{26}{9}}$N₄, which share the same +4 valence state but are nonmagnetic, as seen in Fig. 3(a). Additionally, Si’s conductive band is its 3$p$ orbitals, which cannot be hybridized with the adjacent V 3$d$ orbitals. Using the corresponding occupation number matrix of the $3z^{2}-r^{2}$ and $L_{-2}$ states, our LSDA+SOC+$U$ calculations show that both insulating solutions can be stabilized, as seen in Fig. 4. The $3z^{2}-r^{2}$ state is more stable than the $L_{-2}$ state by 756 meV per V atom, agreeing well with expectations based on the isolated triangular prismatic crystal field. This energy difference is the result of the crystal field energy difference between the $3z^{2}-r^{2}$ and the degenerate $xy$ and $x^{2}-y^{2}$ orbitals, minus the energy contribution from spin-orbit coupling in the degenerate orbitals. Our calculations demonstrate that the $L_{-2}$ state represents a local minimum in the energy landscape, and we have successfully stabilized this state. In contrast to the prominent hybridization between the $3z^{2}-r^{2}$ and $x^{2}-y^{2}$/$xy$ orbitals observed in VSi₂N₄ (refer to Fig. 2), this interaction is significantly reduced in the isolated V${}_{\frac{1}{9}}$Si${}_{\frac{26}{9}}$N₄ composition.

In Fig. 4(a), we observe that the single electron of V⁴⁺ in the 3$d^{1}$ configuration mainly fills the lowest $3z^{2}-r^{2}$ orbital, leaving the higher $x^{2}-y^{2}$/$xy$ ($xz$/$yz$) orbitals unoccupied. The calculated magnetic anisotropy energy (MAE), at 59 $\mu$eV per V atom with easy in-plane magnetization, indicates that the $3z^{2}-r^{2}$ state is not ideal for sustaining stable 2D magnetism. Conversely, as shown in Fig. 4(b), the electron in the V⁴⁺ ion occupies the higher-energy $x^{2}-y^{2}$/$xy$ orbitals when in the $L_{-2}$ state, resulting in a large orbital moment of –1.32 $\mu_{\rm B}$. The calculated MAE value of 23 meV per V atom suggests that the $L_{-2}$ state could support strong Ising-type magnetism.

To enhance magnetic coupling, the coordination number of the V atom is increased in our calculations. We examine two configurations: V${}_{\frac{1}{3}}$Si${}_{\frac{8}{3}}$N₄ and V${}_{\frac{2}{3}}$Si${}_{\frac{7}{3}}$N₄, as displayed in Figs. 3(b) and 3(c). Our results show that the V${}_{\frac{1}{3}}$Si${}_{\frac{8}{3}}$N₄ structure stabilizes both $3z^{2}-r^{2}$ and $L_{-2}$ electronic states. Conversely, in the V${}_{\frac{2}{3}}$Si${}_{\frac{7}{3}}$N₄ structure, the $L_{-2}$ state is unstable and converges to the $3z^{2}-r^{2}$ state. Additionally, when examining the V${}_{\frac{1}{3}}$Si${}_{\frac{8}{3}}$N₄ structure, the MAE drops to 1 $\mu$eV per V atom for the $3z^{2}-r^{2}$ state and to 12 meV per V atom for the $L_{-2}$ state. This reduction is attributed to the interaction between neighboring V ions and is lower than the MAE values found in the V${}_{\frac{1}{9}}$Si${}_{\frac{26}{9}}$N₄ configuration.

After careful analysis of the spin-orbital states and Si-substituted levels in the triangular prismatic crystal field of VSi₂N₄, we find that maximum MAE, large orbital moment, and optimal exchange coupling are attainable in the V${}_{\frac{1}{3}}$Si${}_{\frac{8}{3}}$N₄ configuration. However, this configuration exists in a $3z^{2}-r^{2}$ ground state, which is contrary to our expectations. To address this, introducing Cr ions with an extra electron could stabilize the desired $L_{-2}$ ground state, which is critical for Ising magnetism.

We first calculated the pristine CrSi₂N₄ monolayer. Our results show that it is a non-magnetic material, which has been studied in valley physics. Comparing with the V⁴⁺ 3$d$¹ $S$ = 1/2 state, the extra electron of the Cr⁴⁺ occupies the hybridization orbitals of the lowest $3z^{2}-r^{2}$ singlet and middle $x^{2}-y^{2}$/$xy$ doublet, resulting in a nonmagnetic Cr⁴⁺ with $S$=0 state, as seen in Fig. 5(a).

Then, we calculated the Cr${}_{\frac{1}{3}}$Si${}_{\frac{8}{3}}$N₄, a similar structure of Fig. 3(b). The Cr⁴⁺ exhibits a local spin moment of 2.25 $\mu_{\rm B}$ and large orbital moment of –1.06 $\mu_{\rm B}$. The total spin moment of 6.00 $\mu_{\rm B}$/supercell agrees with the expected Cr⁴⁺ 3$d^{2}$ $S$ = 1 spin configuration. As seen in Fig. 5(b), the lowest $3z^{2}-r^{2}$ is fully occupied, while the higher $x^{2}-y^{2}$/$xy$ orbitals are half-filled. This results in an active SOC and leads to a nominal Cr⁴⁺ $(3z^{2}-r^{2})^{1}(L_{-2})^{1}$ spin configuration. The presence of sizable orbital moment and the DOS together suggest that Cr${}_{\frac{1}{3}}$Si${}_{\frac{8}{3}}$N₄ has the potential for exhibiting 2D Ising magnetism, making it an attractive candidate for applications in spintronics. To explore the stability of Cr${}_{\frac{1}{3}}$Si${}_{\frac{8}{3}}$N₄ under finite temperatures, ab initio molecular dynamics (AIMD) simulations are conducted. Comparisons are also made with CrSi₂N₄ and VSi₂N₄. The results indicate that Cr${}_{\frac{1}{3}}$Si${}_{\frac{8}{3}}$N₄ remains thermodynamically and dynamically stable at 300 K, as seen in Fig. 6.

The ground state of Cr${}_{\frac{1}{3}}$Si${}_{\frac{8}{3}}$N₄ features a local spin moment of 2.25 $\mu_{\rm B}$ for Cr⁴⁺, accompanied by an antiparallel orbital moment of –1.06 $\mu_{\rm B}$ along the $z$-axis. This indicates that the SOC, facilitated by the robust SIA, orients the magnetic moment along the $z$-axis, leading to perpendicular MA and giving rise to Ising magnetism. Here we assume the spin Hamiltonian

where the first term describes the Heisenberg isotropic exchange (AF when $J$>0), the second term is the SIA with the easy magnetization $z$-axis (when $D$>0), and the last term $J^{\prime}$ refers to the anisotropic exchange. To determine the magnetic parameters $J$, $D$, and $J^{\prime}$ in Cr${}_{\frac{1}{3}}$Si${}_{\frac{8}{3}}$N₄, we use a 3$\times$3$\times$1 supercell of MSi₂N₄ containing three Cr ions, as seen in Fig. 3(b). This allows us to flip the spin of the central Cr ion, creating an AF state. However, it is important to note that this AF state can be described as a ferrimagnetic state due to the non-zero net magnetic moment in which two Cr ions have spins aligned in one direction and one in the opposite, leading to a net magnetic moment of –2 $\mu_{\rm B}$ per supercell. We then compute the magnetic properties for four distinct states: FM and AF, each with perpendicular and in-plane magnetization (as detailed in Table 2). Counting $JS^{2}$ for each pair of Cr⁴⁺ $S$ = 1 ion (positive $J$ refers to AF exchange), the magnetic exchange energies of the four states per supercell are written as follows:

Based on our total energy calculations, we estimate the magnetic parameters as follows: $J$ = 1.676 meV, $D$ = 18.634 meV, $J^{\prime}$ = 0.705 meV. One of the key findings is that the Ising-type, represented by the $D$ term, plays a dominant role in establishing perpendicular MA in the Cr${}_{\frac{1}{3}}$Si${}_{\frac{8}{3}}$N₄ monolayer. Specifically, the $D$ is about twenty times stronger than $J^{\prime}$ in stabilizing the 2D Ising magnetism. The positive $J$ refers to the AF coupling of the adjacent Cr⁴⁺-Cr⁴⁺ ions. Therefore, the Cr${}_{\frac{1}{3}}$Si${}_{\frac{8}{3}}$N₄ monolayer can be a 2D AF Ising magnetic material. Furthermore, we study a biaxial strain effect on Cr${}_{\frac{1}{3}}$Si${}_{\frac{8}{3}}$N₄ monolayer. Our results show that the $(3z^{2}-r^{2})^{1}(L_{-2})^{1}$ ground state remains robust against the strains on the optimized lattice. The SIA strength rises in the feasible strain, as seen in Fig. 7.

In summary, we propose Cr${}_{\frac{1}{3}}$Si${}_{\frac{8}{3}}$N₄ monolayer can be a 2D AF Ising magnetic material, using crystal field theory, spin-orbital state analyses, and density functional calculations. Our results indicate that the strong $d$ orbital hybridization between adjacent M⁴⁺ ions in the MSi₂N₄ (M = V, Cr) monolayers disrupts the $d$ orbital splitting in this triangular prismatic crystal field. Through the Si⁴⁺-substituted, the Cr${}_{\frac{1}{3}}$Si${}_{\frac{8}{3}}$N₄ monolayer can achieve the huge perpendicular MA of 18.63 meV per Cr atom with a large orbital moment of –1.06 $\mu_{\rm B}$ along the $z$-axis. Our research emphasizes the importance of investigating the degrees of freedom in spin-orbital states as a fruitful avenue for discovering new 2D Ising materials.

This work was supported by National Natural Science Foundation of China (Grants No. 12104307). S. Chen and W. Xu contributed equally to this work.