Introduction

The intimate interplay between correlated electrons, lattice, and magnetism can result in a rich variety of interesting and important physical phenomena1. In two-dimensional (2D) materials with open d- or f-shells, additional quantum confinement caused by the reduced dimensionality suppresses the screening2,3 and thus may further enhance the electron correlation. The recent development of magnetic 2D van der Waals (vdW) materials adds a new magnetic functionality to the already vast appeal of the 2D materials family4,5,6,7. Their magnetism is very sensitive to, and can be controlled by, pressure8, stacking arrangement9, and external magnetic4 and electric10,11,12 fields. Such tunability offers opportunities to design and construct new energy-efficient spin-based devices. Even in their bulk form, the reduced coordination number in quasi-2D lattices constrains the electron hopping, thereby increasing the role of the Coulomb interaction. The resulting enhancement of electron correlations directly affects exchange interactions. The critical question is how the magnetic excitations in these confined systems can be understood, controlled, and exploited.

Despite considerable attention, accurate ab initio descriptions of spin excitations and a comprehensive understanding of magnetic interactions in magnetic 2D vdW materials (m2Dv) are still lacking. Some complexity is imparted by the role of spin–orbit coupling (SOC). For example, the presence and role of the Dzyaloshinskii–Moriya interaction (DMI), which is essential for materials to host topological magnons and other topological objects such as skyrmions13,14, remains puzzling. Specifically, 2D honeycomb ferromagnets can be viewed as the magnetic analog of graphene15,16. Earlier theoretical work17,18 demonstrated that introducing a next-nearest-neighboring DMI interaction, which breaks the inversion symmetry of the 2D honeycomb lattice, can induce a spinwave (SW) gap at the Dirac points and allow topological magnons. Recent inelastic neutron-scattering (INS) experiments have shown a gap opening along the high-symmetry lines in pristine CrI316, implying a sizable DMI or Kitaev interaction19 may indeed exist in CrI3. In contrast, ab initio investigations show that both DMI20,21,22,23 and Kitaev interactions23,24 in pristine CrI3 are insufficient to open up the sizable gap observed in experiments. The role of the relativistic magnetic interactions is still not fully understood. However, irrespective of this, the exchange interactions should be reinvestigated beyond density functional theory (DFT) to address the electron-correlation effects25,26,27.

In this work, using the most studied m2Dv—CrI3—as a prototype, we identify the role of electron correlations on the magnetic interactions and excitations in m2Dv. The central quantity that characterizes the spin excitations—the dynamic transverse spin susceptibility (DTSS)—is calculated and directly compared with the SW spectra measured by INS16,28. The electron-correlation effects beyond DFT are included via the quasiparticle self-consistent GW (QSGW) method29,30, as incorporated through its effective one-particle potential, which consists of both on-site and off-site nonlocal parts. We demonstrate that the explicit treatment of electron correlations is required to accurately describe the magnetic interactions, especially the interlayer interactions, in m2Dv. Furthermore, we made the remarkable discovery that a sizable magnon gap opens along the high-symmetry line in CrI3 even without relativistic exchanges. Instead, this gap is caused by a correlation-enhanced interlayer Cr-I-I-Cr super-superexchange coupling.

Results and discussion

Linear response theory

Starting from a self-consistent ab initio band structure, we first calculate the bare transverse spin susceptibility \({\chi }_{0}({\bf{r}},{\bf{r}}^{\prime} ,{\bf{q}},\omega )\) using a linear response method31,32,33. Then the full transverse susceptibility χ is calculated, within the random phase approximation (RPA), as χ = χ0 + χ0Iχ, where I is the exchange-correlation kernel. Two-particle quantities χ0, χ, and I are functions of coordinates r and \({\bf{r}}^{\prime}\) within the unit cell. As magnetic moments and excitations are nearly completely confined within the Cr sites, we calculate χ0 on a product basis30,34 and then project it onto the local spin densities of Cr pairs. This projection discretizes \({\chi }_{0}({\bf{r}},{\bf{r}}^{\prime} {\bf{q}},\omega )\) into a matrix χ0(i, j, q, ω), where i and j index the Cr sites in the unit cell. Such discretization allows us to (1) determine the kernel I using a sum rule31; (2) map χ−1 into the Heisenberg model to extract pair exchange parameters; and (3) greatly reduce the computational effort. As for the electronic structures, we employ the QSGW method29,30, wherein nonlocal potentials, both on-site and long-range off-site, are explicitly calculated. The widely used DFT + U methods35, which provide a simplistic correction of on-site Coulomb interactions, are also employed for comparison. Besides the nonlocal on-site potential, the off-site potential in QSGW could also be critical, as it directly affects the relative positions of cation-3d and anion-p bands26, and thus the indirect exchange interactions of Cr ions via one (superexchange) or multiple (super-superexchange) Iodine sites. Further details of the QSGW method and implementation29,30,36,37, and applications on CrI326 can be found in the Supplementary Methods.

By far, except for a few studies23,38,39 using the magnetic force theorem (MFT)40, most of the theoretical investigations of the magnetic interactions in m2Dv are based on mapping the total energies of various collinear spin configurations into the Heisenberg model (referred to hereafter as the energy-mapping method), often employed with DFT + U. However, the applicability and accuracy of such an approach in CrI3 are not clear. Specifically, using such an energy-mapping method, DFT overestimates the exchange couplings by 50%16,41; the additional U correction on Cr sites in DFT + U further increase coupling, only worsening the agreement with experiments. In principle, the MFT approach is more suitable to describe the small spin deviation from the ground state, such as the SW excitations; moreover, it also allows one to resolve coupling into orbital contributions and elucidate the underlying exchange mechanism39,42. However, the resulting values still vary and are inconclusive; even the opposite trend of exchange-coupling dependence on on-site Cr U-values have been obtained23. The discrepancy likely lies in the details of the constructions of the TB Hamiltonian and Green’s function43, which are often assisted by the Wannier function approach44. On the other hand, DTSS is challenging to compute in practice. Studies to date have been mostly limited to simple systems, likely due to its computationally demanding nature and other complications45. For example, the evaluation of kernel I is not explicit in the case beyond the mean-field scheme and the Goldstone theorem, which ensures SW at q = 0 and ω = 0 in the absence of magnetocrystalline anisotropy, is often not guaranteed. In this work, we calculate DTSS without including SOC; the Goldstone magnon mode at q = 0 and ω = 0 is ensured as the kernel I is determined to satisfy the sum rule (Eq. (2) in Supplementary Methods). Finally, to investigate the detailed SW dispersion near the Dirac point, high-resolution SW spectra are needed. This is achieved by calculating the real-space bare susceptibility χ0(R, ω) on a R-mesh first and then Fourier transforming χ0(R, ω) to obtain χ0(q, ω) with a dense set of q points along high-symmetry paths.

CrI3 crystallizes in either the low-temperature rhombohedral structure (R-CrI3) or high-temperature monoclinic structure (M-CrI3)46. A honeycomb Cr monolayer is sandwiched between two I layers; then, the blocks of I-Cr-I triple layers are stacked along the z direction, held together by a weak vdW force. M-CrI3 has a slightly distorted honeycomb Cr lattice and, more importantly, a different stacking arrangement, resulting in the A-type antiferromagnetic (AFM) ordering9,47,48,49 instead of the ferromagnetic (FM) ordering as in R-CrI3. The sensitivity of interlayer Cr coupling on stacking arrangement reflects the changes of super-superexchange pathways across the vdW gap. A thorough understanding requires an explicit treatment of interlayer exchanges, instead of using a single effective interlayer exchange parameter as is often employed to describe the system. Thus, the exchange couplings and SWs in R-CrI3 need to be considered in the context of rhombohedral symmetry. In this work, we first focus on the magnetism in FM R-CrI3 and then discuss stacking effects on magnetism using M-CrI3.

Figure 1 shows the crystal structure of R-CrI3 and the corresponding Brillouin zone (BZ). The 2D honeycomb Cr sublattice and corresponding BZ are also shown for better comparison with previous studies, in which the SWs are often discussed in the hexagonal notation. The rhombohedral primitive unit cell contains two formula units, whereas the conventional hexagonal cell includes six. Correspondingly, the former has a larger BZ than the latter. As shown in Fig. 1d, along the kx direction, its BZ boundary M3 is three times further from Γ than that of the hexagonal structure, Γ–M1. Along the ky direction, M2 (equivalent to M3) is at the BZ boundary for both rhombohedral and hexagonal cells.

Fig. 1: Crystal structure of R-CrI3 and corresponding BZ in comparison to hexagonal lattice.
figure 1

a The rhombohedral primitive unit cell and hexagonal conventional unit cell of R-CrI3 (\(R\overline{3}\), BiI3-type, space group no. 148). b The primitive unit cell of 2D honeycomb Cr sublattice. c The first BZ and the reciprocal lattice vectors of the rhombohedral primitive cell. Along the \({\hat{k}}_{z}\) direction, the rhombohedral BZ boundary Z is denoted as (0.5, 0.5, 0.5) in rhombohedral notation or (0, 0, 1.5) in hexagonal notation. Dashed lines denote the BZ of the 2D hexagonal lattice. d The top view of c. Vectors b1 and b2 are the reciprocal lattice vectors of the 2D hexagonal lattice shown in b.

Dynamic transverse spin susceptibility

First, we calculate within DFT the bare and full transverse susceptibilities, \({\chi }_{0}\left({\bf{q}},\omega \right)\) and \(\chi \left({\bf{q}},\omega \right)\), which characterize the single-particle Stoner excitations and collective SW excitations, respectively. The intensities of \({\rm{Im}}[{\chi }_{0}\left({\bf{q}},\omega \right)]\) and \({\rm{Im}}[\chi \left({\bf{q}},\omega \right)]\) along high-symmetry paths are shown in Fig. 2a, c, respectively. The corresponding special q points are denoted in Figs 1c and 2d. The energy scale of the Stoner excitations, starting at ~1.7 eV and peaking at ~3 eV, are two orders of magnitude larger than SW excitations, leaving the latter with no damping. This means that we can likely use the physical picture of well-defined local magnetic moments on Cr atoms and map the low-lying spin dynamics onto a purely localized-spin Hamiltonian as will be described in detail below. The threshold energy of ~1.7 eV corresponds to the gap size of the spin-flip transition from the top of majority-spin Cr valence states to the bottom of minority-spin Cr conduction states, as shown in Fig. 2b, whereas the peak energy corresponds to the spin splitting of the Cr-d states. The SW energies, defined by the peaks of \({\rm{Im}}[\chi \left({\bf{q}},\omega \right)]\), are solely determined by the RPA poles of 1 − Iχ0 = 0. With two Cr atoms in the primitive cell, we find two poles for each q, resulting in two magnon branches, as shown in Fig. 2c–f. Experimental SW energies extracted from previous INS work by Chen et al.16,28 are plotted to compare with the calculated SW spectra.

Fig. 2: Spin excitations calculated from \({\rm{Im}}[{{\boldsymbol{\chi }}}_{0}\left({\bf{q}},\omega \right)]\) and \({\rm{Im}}[{\boldsymbol{\chi }}\left({\bf{q}},\omega \right)]\) in R-CrI3.
figure 2

The special q points along the high-symmetry path ZΓM2M3ΓL are denoted in Fig. 1c, d. a \({\rm{Im}}[{\chi }_{0}\left({\bf{q}},\omega \right)]\) calculated in DFT. b The density of states calculated in DFT. The horizontal dashed-and-dotted line indicates the top of valence band. c–f \({\rm{Im}}[\chi \left({\bf{q}},\omega \right)]\) calculated in DFT (c), DFT + U (d), QSGW (e), and QSGW + U (f). The intensity of \({\rm{Im}}[\chi \left({\bf{q}},\omega \right)]\) is shown in log scale. Experimental SW energies, adopted from INS work by Chen et al.16,28, are denoted by open circles.

We now compare the DFT SW spectra with INS experiments. Along the Γ–M2 path, Fig. 2c shows that two SW branches cross near the point K, before reaching the BZ boundary M2. For other directions, a gap exists between two magnon branches. This is consistent with previous studies16 without considering DMI. Along the Γ–M3 path, the SW minimum occurs at [300] (hexagonal notation), instead of [100] (Γ1-point), reflecting the symmetry of the rhombohedral structure. The maximum energy of the optical mode measured in INS is about 20 meV, whereas DFT gives 30 meV and overestimates it by 50%. As shown in Fig. 2c, DFT also overestimates the acoustic magnon energies, most severely (by a factor of ~4) along the confined z direction. The overestimation of interlayer coupling also affects the in-plane SW energies, because nearly all interlayer couplings have in-plane components in their connecting vectors. Hence, in bulk materials, an accurate description of interlayer exchanges is also essential to describe the in-plane SW accurately.

Next, we investigate the effects of electronic correlations on magnetic interactions within DFT + U and QSGW. We found that increasing the U-value in DFT + U increases the on-site Cr moment (Supplementary Table I) and lowers the energies of both magnon branches (Supplementary Fig. 3). Figure 2d shows the SW spectra for U = 3 eV, a typical value used for CrI3. The SW energy at the Z-point is decreased to 5.8 meV, about 2/3 of the DFT value; however, it is still about a factor of three larger than that found in INS experiments. Moreover, the optical mode becomes flatter, in contrast to INS16. Surprisingly, electron-correlation effects also introduce a SW gap at the Dirac point, which will be discussed in detail later. Figure 2e shows the SW spectra calculated within QSGW. The optical magnon is centered at around 24 meV, similar as in DFT + U, but recovers some dispersion and agrees better with INS experiments. Interestingly, the acoustic SW energies are reduced in comparison to DFT + U. The SW energy at the Z-point drops to ~3.4 meV, showing a much weaker interlayer FM coupling and better agreement with experiments. Thus, the elaborate treatment of electron interactions in QSGW improves the description of spin excitations in these systems. Considering GW methods may underestimate the on-site electron interactions, they have been applied on top of DFT + U for various systems50. Here we also apply QSGW on top of DFT + U (referred to hereafter as QSGW + U)51, to roughly mimic the additional on-site interactions. Figure 2f shows the SW spectra calculated with U = 1.36 eV. Although the in-plane acoustic SW is still somewhat overestimated, the overall spectra compare well with experiments, suggesting that additional on-site interactions beyond QSGW may be needed to best describe electronic structures and SW in CrI3. More rigorous and comprehensive frameworks, such as the dynamical mean-field theory (DFMT) + QSGW52, could prove valuable for future research.

Exchange parameters

To develop a quantitative understanding of how electron correlations affect magnetic interactions and excitations, we calculate the effective pair exchange parameters Jij for a Heisenberg model \(H=-{\sum }_{i\ne j}{J}_{ij}\ {\hat{{\bf{e}}}}_{i}\cdot {\hat{{\bf{e}}}}_{j}\), where \({\hat{{\bf{e}}}}_{i}\) is the unit vector of the atomic spin moment at site i. The Jij parameters are obtained from the inverse of susceptibility matrix, [χ(q, ω = 0)]−1, with a subsequent Fourier transform31,53,54,55. As shown in ref. 54, the exchange parameters determined in this way coincide with those calculated via the MFT. The exchange parameters between the first few neighbors, as depicted in Fig. 3a, are listed in Table 1; couplings beyond 12 Å are negligible. Using the linear SW theory, we recalculate the SW spectra with the extracted Jij parameters. The obtained spectra agree well with those determined by the peaks of \({\rm{Im}}[\chi \left({\bf{q}},\omega \right)]\) (see Supplementary Fig. 2).

Fig. 3: Exchange parameters, SW spectra, and nodal lines calculated using linear spinwave theory.
figure 3

a Pair exchange parameters for the first few neighbors in R-CrI3. b Atomic configurations around the exchange paths of \({\tilde{J}}_{{2}_{0}}\) and \({\tilde{J}}_{2}\). AFM \({\tilde{J}}_{2}\) corresponds to a Cr-I-I-Cr super-superexchange coupling. c SW spectra along the Γ–KM path calculated within a J1-\({\tilde{J}}_{2}\)-\({\tilde{J}}_{{2}_{0}}\) (\({\tilde{J}}_{2}={\tilde{J}}_{{2}_{0}}={J}_{1}/12\)) model and a J1-\({\tilde{J}}_{2}\) (\({\tilde{J}}_{2}={J}_{1}/6\)) model. d Dirac nodal lines, where the magnon bands cross, wind around the K-point and along the \({\hat{k}}_{z}\) direction; the exchange parameters calculated in QSGW + U are used.

Table 1 Pairwise intralayer (Ji) and interlayer (\({\tilde{J}}_{i}\)) exchange parameters in R-CrI3calculated within various methods.

We find that the on-site Hubbard U corrections on Cr-d states included in DFT + U have a stronger effect on decreasing the nearest-neighboring coupling, while the explicit nonlocal potentials in QSGW have a more substantial effect on the longer-range interlayer couplings. Within DFT, the calculated in-plane exchanges are similar to values41 derived from total energy mapping and ~50% larger than those extracted from INS16. In comparison to DFT, DFT + U decreases J1 by ~19% and the overall interlayer (FM) couplings \(\sum \tilde{J}\) by ~27%, whereas QSGW decreases J1 by ~14% and \(\sum \tilde{J}\) by ~60%. QSGW not only reduces the optical SW energies but also significantly lowers the acoustic ones, especially along the interlayer direction.

Interestingly, regarding the dependence of exchange and SW energy on the U parameter, the energy-mapping method using multiple collinear states gives the opposite trend from the linear response method. Using the energy-mapping method, exchange and SW energy increase with U, worsening the agreement between theory and INS measurements. This suggests that non-Heisenberg interactions, such as biquadratic exchange or multi-site exchange interactions56,57,58,59,60, might be important in this system, especially when additional electron correlations are taken into account. By analogy with transition-metal oxide systems57 one can assume that the induced magnetic polarization on iodine ligands can play a decisive role in this non-Heisenberg behavior. Contrary to the total energy differences, the linear response method describes accurately the small spin deviations from the given (ground) state and thus the SW spectra. The corresponding exchange integrals thus depend on the initial equilibrium spin configuration in which they are calculated40,56. On the other hand, mapping the total energy of spin-spiral configurations may provide a better description of Jij than the collinear configurations.

Interlayer coupling induced gap opening along Γ–KM

Remarkably, besides the improvement of SW energies, correlations beyond DFT also open up a gap along the Γ–KM2 path in both DFT + U and QSGW, as shown in Fig. 2d–f. Although the calculated gap size (~1.8 meV in QSGW + U) is smaller than the experimental value of ~4 meV observed in INS, the existence of such a gap along the Γ–KM path is unexpected in the absence of DMI or Kitaev interaction. As we show later, this gap can be caused by the correlation-enhanced interlayer super-superexchange \({\tilde{J}}_{2}\). If one only considers the Cr sublattice itself, which has a higher symmetry (\(R\overline{3}m\)), \({\tilde{J}}_{{2}_{0}}\) and \({\tilde{J}}_{2}\) should be equivalent. The presence of the Iodine sublattice breaks the mirror and twofold rotational symmetry, lifting the degeneracy of \({\tilde{J}}_{{2}_{0}}\) and \({\tilde{J}}_{2}\). Thus, although \({\tilde{J}}_{{2}_{0}}\) and \({\tilde{J}}_{2}\) connect similar Cr pairs with the same distance of 9.517 Å, their exchange paths are not the same, which allows for J2 and \({J}_{{2}_{0}}\) to adopt different values. However, up until now, all previous work has considered that \({J}_{2}={J}_{{2}_{0}}=0\), which indeed is supported, to a great extent, by DFT. Surprisingly, the inclusion of correlations beyond DFT results in a sizable AFM \({\tilde{J}}_{2}\). As shown in Table 1, \({\tilde{J}}_{{2}_{0}}\) vanishes in all calculations, whereas AFM \({\tilde{J}}_{2}\) is negligible in DFT but becomes stronger in QSGW and DFT + U, reaching \({\tilde{J}}_{2}=10 \% {J}_{1}\) in QSGW + U. As shown in Fig. 3b, \({\tilde{J}}_{2}\) corresponds to a Cr-I-I-Cr super-superexchange with a Cr-I-I angle of 159°, giving the major AFM contribution to interlayer couplings, whereas vanishing \({J}_{{2}_{0}}\) has no obvious exchange path with intervening I anions.

The gap between acoustic and optical modes, at an arbitrary q point, gives the energy difference between the in-phase and out-of-phase precessions of two Cr-spin sublattices, and depends on the inter-sublattice couplings 2B(q) (see Eq. (5) and Eq. (6) of the “Methods” section). Within the considered exchange range, we have \(B({\bf{q}})={J}_{1}({\bf{q}})+{J}_{3}({\bf{q}})+{\tilde{J}}_{0}({\bf{q}})+{\tilde{J}}_{1+}({\bf{q}})+{\tilde{J}}_{{2}_{0}}({\bf{q}})+{\tilde{J}}_{2}({\bf{q}})\), where Ji(q) is the corresponding Fourier component of Ji(R). Couplings J1(q), J3(q), and \({\tilde{J}}_{1+}({\bf{q}})\) are real functions along the Γ–KM path and vanish at the K-point, resulting a bandcrossing at K if other terms are ignored. The interlayer coupling \({\tilde{J}}_{0}\) is along the z direction; \({\tilde{J}}_{0}({\bf{q}})={\tilde{J}}_{0}\) is a real constant when q is in the basal plane, shifting the bandcrossing along the Γ–KM path. In real space, the connecting vectors (in-plane components) of \({\tilde{J}}_{{2}_{0}}\) and \({\tilde{J}}_{2}\) are rotated by π/6 with respect to those of J1, J3, and \({\tilde{J}}_{1+}\). Correspondingly, in reciprocal space, \({\tilde{J}}_{{2}_{0}}({\bf{q}})\) and \({\tilde{J}}_{2}({\bf{q}})\) are complex functions along the Γ–KM path (see Supplementary Fig. 4). With \({\tilde{J}}_{{2}_{0}}=0\), \({\tilde{J}}_{2}({\bf{q}})\) itself results in a non-vanishing B(q) and thus a gap along the Γ–KM path. However, if \({\tilde{J}}_{{2}_{0}}={\tilde{J}}_{2}\), then \(({\tilde{J}}_{{2}_{0}}({\bf{q}})+{\tilde{J}}_{2}({\bf{q}}))\) is a real function when q is in the basal plane, shifting the magnon crossing as \({\tilde{J}}_{0}\) does along the Γ–KM path. Thus, the combination of vanishing \({\tilde{J}}_{{2}_{0}}\) and correlation-enhanced \({\tilde{J}}_{2}\) will induce the magnon gap along the Γ–KM path. To illustrate, we calculate the SW spectra in a simple J1\({\tilde{J}}_{{2}_{0}}\)\({\tilde{J}}_{2}\) model with \({\tilde{J}}_{{2}_{0}}={\tilde{J}}_{2}={J}_{1}/12\) and a J1\({\tilde{J}}_{2}\) model with \({\tilde{J}}_{2}={J}_{1}/6\), respectively. As shown in Fig. 3c, the latter case opens a gap along the Γ–KM path.

However, unlike the global DMI-induced gap, the gap induced by the nonequivalence of the exchange interactions \({\tilde{J}}_{{2}_{0}}\) and \({\tilde{J}}_{2}\) does not persist through the whole BZ because a solution of B(q) = 0 can be found near the K-point. The Dirac nodal lines in the SW spectra still form but do not cross the Γ–KM lines (at the kz = 0 plane). Using the Jij parameters obtained within QSGW + U, we calculate, within the linear SW theory, the in-plane SW spectra at various kz planes. Figure 3d shows the helical Dirac nodal lines form around the edges of the hexagonal BZ; each line crosses only twice the face of the first BZ. It would be interesting to see whether future INS experiments can confirm the small displacement of the Dirac point off the Γ–K line or at finite kz. However, such measurement may be challenging due to the complexity from the modulation of the dynamic structure factor close to the gap, the requirements of high instrumental resolution and good sample mosaic, and the finite SOC effects.

Stacking-dependent magnetic ordering

Finally, we demonstrate that explicit treatments of electron correlations can correctly describe the dependence of interlayer interaction on stacking order. M-CrI3 has different stacking than R-CrI3, which dramatically modifies the interlayer super-superexchange paths and results in A-type AFM ordering in M-CrI3. This intimate interplay between stacking order and magnetic ordering plays a crucial role in manipulating the magnetism in these materials. DFT total energy calculations predict the wrong FM ground state for M-CrI3, whereas DFT + U calculations47,48 have shown that AFM interlayer configurations can be stabilized in M-CrI3, depending on the U-value. Within QSGW, we calculate \(\chi \left({\bf{q}},\omega \right)\) in M-CrI3 starting from both the FM and the A-type AFM ground-state configurations. The corresponding SW spectra along the high-symmetry paths are shown in Fig. 4a, b, respectively. The acoustic SW calculated with FM configuration, as shown in Fig. 4a, is negative along the Γ–Z path (normal to the basal plane), suggesting the instability of the FM interlayer configuration in M-CrI3. In contrast, the SW spectra calculated with the AFM configuration, as shown in Fig. 4b, are positive through the whole BZ. Thus, by taking into account the explicit electron correlations in QSGW, the magnetic ordering dependence on stacking can be correctly described in a parameter-free fashion.

Fig. 4: SW spectra in FM and A-type AFM M-CrI3 calculated within QSGW.
figure 4

Wavevectors q = hb1 + kb2 + lb3 is denoted as q = (h, k, l) in reciprocal lattice units (r.l.u.). SW are plotted along X–Γ–Y–Γ–Z for the FM structure and X–Γ–Y–Γ–Z–Γ for the AFM structure. High-symmetry q points X, Y, and Z are at BZ boundaries and denoted as (0.5, 0.5, 0), (−0.5, 0.5, 0), and (0, 0, 0.5), respectively. a \({\rm{Im}}[\chi \left({\bf{q}},\omega \right)]\) calculated in FM configuration. b \({\rm{Im}}[\chi \left({\bf{q}},\omega \right)]\) calculated in AFM configuration. c The primitive cell of AFM M-CrI3 structure, in which the lattice vector a3 is doubled, in comparison to the FM structure. d The first BZ of M-CrI3. Reciprocal lattice vectors b1 and b2 are the same for both FM and AFM structures, while \({{\bf{b}}}_{3}^{(\,\text{FM}\,)}=2{{\bf{b}}}_{3}^{(\,\text{AFM}\,)}\).

In conclusion, we have demonstrated that an ab initio description of electron correlations beyond DFT can accurately describe the SW spectra and the stacking-dependent magnetism in a parameter-free fashion, and reveals the gap-opening mechanism in CrI3. An unexpected correlation-enhanced interlayer super-superexchange induces a gap along the Γ–KM path in the absence of previously proposed relativistic exchanges or magnon–phonon interactions23. To experimentally verify the true physical mechanism of the gap opening and the nature of magnetic interactions, future INS experiments may be used to search for the bandcrossings off the high-symmetry line in bulk CrI3. Identifying the existence of the magnon gap along the Γ–KM path in monolayer CrI3, in which the interlayer exchange is absent, can also help illuminate the responsible interactions. Of course, INS studies of single-layer materials seem to be impossible due to a small number of atoms. However, other techniques such as electron energy loss spectroscopy can be, in principle, used61. Our work suggests the necessity of an explicit treatment of electron correlations to accurately describe the magnetism in the broad family of layered magnetic materials62,63, including magnetic topological vdW materials64,65, where the interaction between magnetization and the topological surface state is essential.

Methods

Pair exchange parameters for Heisenberg model

In adiabatic approximation, we map spin susceptibility into a classical Heisenberg model

$$H=-\sum _{i\ne j}{J}_{ij}\ {\hat{{\bf{e}}}}_{i}\cdot {\hat{{\bf{e}}}}_{j},$$
(1)

where \({\hat{{\bf{e}}}}_{i}\) is the unit vector of magnetic moment on site i. Effective pair exchange parameters are calculated as

$${J}_{ij}=\frac{1}{4}{\int}_{{\Omega }_{i}}\ {\rm{d}}{{\bf{r}}}_{1}{\int}_{{\Omega }_{j}}\ {\rm{d}}{{\bf{r}}}_{2}\ {m}_{i}({{\bf{r}}}_{1})J({{\bf{r}}}_{1},{{\bf{r}}}_{2}){m}_{j}({{\bf{r}}}_{2}),$$
(2)

where mi(r) is the density of magnetic moment on Cr site i and J(r1, r2) is related to χ+ (r1, r2) and satisfies

$${\int}_{\Omega }\ {\rm{d}}{{\bf{r}}}_{2}\ J({{\bf{r}}}_{1},{{\bf{r}}}_{2}){\chi }^{+-}({{\bf{r}}}_{2},{{\bf{r}}}_{3},\omega =0)=\delta ({{\bf{r}}}_{1}-{{\bf{r}}}_{3}).$$
(3)

Thus, the effective pair exchange parameters can be obtained from the matrix elements of the inverse of spin susceptibilities matrix31,53,54 with a subsequent Fourier transform,

$$\begin{array}{lll}&&{J}_{ij}({\bf{R}})=\frac{1}{{\Omega }_{{\rm{BZ}}}}\int \ {\rm{d}}{\bf{q}}\ {e}^{i2\pi {\bf{q}}\cdot {\bf{R}}}{J}_{ij}({\bf{q}})\\ &=&\frac{1}{{\Omega }_{{\rm{BZ}}}}\int \ {\rm{d}}{\bf{q}}\ {e}^{i2\pi {\bf{q}}\cdot {\bf{R}}}| | {m}_{i}| | {\left[{\left({\chi }^{+-}({\bf{q}},\omega = 0)\right)}^{-1}\right]}_{ij}| | {m}_{j}| | .\end{array}$$
(4)

Here, the \({\chi }_{ij}^{+-}\) matrices are obtained by projecting onto the functions {mi(r)/mi} representing normalized local spin densities on each magnetic Cr site, which gives a matrix χij(q, ω) in magnetic basis site indices31. This projection corresponds to the rigid spin approximation, which is a modest approximation for CrI3. A 8 × 8 × 8 q-mesh was employed to calculate Jij(R). Using a larger 12 × 12 × 12 q-mesh barely changes the SW spectra (see Supplementary Fig. 1).

The detailed calculation of transverse spin susceptibility using the linear response method and preparation of the single-particle Hamiltonian, including both DFT and QSGW, can be found in Supplementary Methods.

Linear SW theory

Starting from the Heisenberg Hamiltonian, Eq. (1), and using the Holstein–Primakoff transformation, the SW energies of the acoustic and optical modes are written as:

$$E({\bf{q}})=\frac{4}{m}\left(A(0)+B(0)-A({\bf{q}})\pm | B({\bf{q}})| \right),$$
(5)

where A(q) and B(q) are sum of J(q) over the intra-sublattice and inter-sublattice pairs, respectively:

$$\begin{array}{l}A({\bf{q}})={J}_{2}({\bf{q}})+{\tilde{J}}_{1-}({\bf{q}})\\ B({\bf{q}})={J}_{1}({\bf{q}})+{J}_{3}({\bf{q}})+{\tilde{J}}_{0}({\bf{q}})+{\tilde{J}}_{1+}({\bf{q}})+{\tilde{J}}_{{2}_{0}}({\bf{q}})+{\tilde{J}}_{2}({\bf{q}}).\end{array}$$
(6)

Here,

$${J}_{i}({\bf{q}})={J}_{i}\sum _{\delta }{e}^{-i2\pi {\bf{q}}\cdot \delta }$$
(7)

is the Fourier transform of real-space exchange parameters Ji with corresponding connecting vectors δ. As shown in Eq. (5), the gap between two modes is 2B(q). Along the Γ–KM path, J1(q), J3(q), and \({\tilde{J}}_{1+}({\bf{q}})\) are real cosine functions of q and vanish at point K, whereas \({\tilde{J}}_{0}({\bf{q}})\) is a real constant and \({\tilde{J}}_{2}({\bf{q}})\) a complex number.