Electron correlation effects on exchange interactions and spin excitations in 2D van der Waals materials

Despite serious effort, the nature of the magnetic interactions and the role of electron-correlation effects in magnetic two-dimensional (2D) van der Waals materials remains elusive. Using CrI3 as a model system, we show that the calculated electronic structure including nonlocal electron correlations yields spin excitations consistent with inelastic neutron-scattering measurements. Remarkably, this approach identifies an unreported correlation-enhanced interlayer super-superexchange, which rotates the magnon Dirac lines off, and introduces a gap along the high-symmetry Γ-K-M path. This discovery provides a different perspective on the gap-opening mechanism observed in CrI3, which was previously associated with spin–orbit coupling through the Dzyaloshinskii–Moriya interaction or Kitaev interaction. Our observation elucidates the critical role of electron correlations on the spin ordering and spin dynamics in magnetic van der Waals materials and demonstrates the necessity of explicit treatment of electron correlations in the broad family of 2D magnetic materials.


INTRODUCTION
The intimate interplay between correlated electrons, lattice, and magnetism can result in a rich variety of interesting and important physical phenomena 1 . In two-dimensional (2D) materials with open d-or f-shells, additional quantum confinement caused by the reduced dimensionality suppresses the screening 2,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 family [4][5][6][7] . Their magnetism is very sensitive to, and can be controlled by, pressure 8 , stacking arrangement 9 , and external magnetic 4 and electric [10][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 skyrmions 13,14 , remains puzzling. Specifically, 2D honeycomb ferromagnets can be viewed as the magnetic analog of graphene 15,16 . Earlier theoretical work 17,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 CrI 3 16 , implying a sizable DMI or Kitaev interaction 19 may indeed exist in CrI 3 . In contrast, ab initio investigations show that both DMI [20][21][22][23] and Kitaev interactions 23,24 in pristine CrI 3 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 effects [25][26][27] . In this work, using the most studied m2Dv-CrI 3 -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 INS 16,28 . The electroncorrelation effects beyond DFT are included via the quasiparticle self-consistent GW (QSGW) method 29,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 highsymmetry line in CrI 3 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 χ 0 ðr; r 0 ; q; ωÞ using a linear response method [31][32][33] . Then the full transverse susceptibility χ is calculated, within the random phase approximation (RPA), as χ = χ 0 + χ 0 Iχ, where I is the exchange-correlation kernel. Two-particle quantities χ 0 , χ, and I are functions of coordinates r and r 0 within the unit cell. As magnetic moments and excitations 1 are nearly completely confined within the Cr sites, we calculate χ 0 on a product basis 30,34 and then project it onto the local spin densities of Cr pairs. This projection discretizes χ 0 ðr; r 0 q; ωÞ 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 rule 31 ; (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 method 29,30 , wherein nonlocal potentials, both on-site and longrange off-site, are explicitly calculated. The widely used DFT + U methods 35 , 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 bands 26 , 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 implementation 29,30,36,37 , and applications on CrI 3 26 can be found in the Supplementary Methods. By far, except for a few studies 23,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 CrI 3 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 mechanism 39,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 obtained 23 . The discrepancy likely lies in the details of the constructions of the TB Hamiltonian and Green's function 43 , which are often assisted by the Wannier function approach 44 . 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 complications 45 . 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. CrI 3 crystallizes in either the low-temperature rhombohedral structure (R-CrI 3 ) or high-temperature monoclinic structure (M-CrI 3 ) 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-CrI 3 has a slightly distorted honeycomb Cr lattice and, more importantly, a different stacking arrangement, resulting in the A-type antiferromagnetic (AFM) ordering 9,47-49 instead of the ferromagnetic (FM) ordering as in R-CrI 3 . 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-CrI 3 need to be considered in the context of rhombohedral symmetry. In this work, we first focus on the magnetism in FM R-CrI 3 and then discuss stacking effects on magnetism using M-CrI 3 . Figure 1 shows the crystal structure of R-CrI 3 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 k x direction, its BZ boundary M 3 is three times further from Γ than that of the hexagonal structure, Γ-M 1 . Along the k y direction, M 2 (equivalent to M 3 ) is at the BZ boundary for both rhombohedral and hexagonal cells.
Dynamic transverse spin susceptibility First, we calculate within DFT the bare and full transverse susceptibilities, χ 0 q; ω ð Þ and χ q; ω ð Þ, which characterize the single-particle Stoner excitations and collective SW excitations, respectively. The intensities of Im½χ 0 q; ω ð Þ and Im½χ q; ω ð Þ 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 minorityspin 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 Im½χ q; ω ð Þ, 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.
We now compare the DFT SW spectra with INS experiments. Along the Γ-M 2 path, Fig. 2c shows that two SW branches cross near the point K, before reaching the BZ boundary M 2 . For other directions, a gap exists between two magnon branches. This is consistent with previous studies 16 without considering DMI. Along the Γ-M 3 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 CrI 3 . The SW energy at the Zpoint 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 INS 16 . 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 systems 50 .
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 CrI 3 . More rigorous and comprehensive frameworks, such as the dynamical mean-field theory (DFMT) + QSGW 52 , 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 J ij for a Heisenberg model H ¼ À P i≠j J ijêi Áê j , whereê i is the unit vector of the atomic spin moment at site i. The J ij parameters are obtained from the inverse of susceptibility matrix, [χ(q, ω = 0)] −1 , with a subsequent Fourier transform 31,53-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 J ij parameters. The obtained spectra agree well with those determined by the peaks of Im½χ q; ω ð Þ (see Supplementary Fig. 2). 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 values 41 derived from total energy mapping and~50% larger than those extracted from INS 16 . In comparison to DFT, DFT + U decreases J 1 by~19% and the overall interlayer (FM) couplings PJ by~27%, whereas QSGW decreases  The Heisenberg model is defined as H ¼ À P i≠j J ijêi Áê j , whereê i is the unit vector of the atomic spin moment at site i. The degeneracy (No.) and distance (R ij ) of J ij are also provided. Symbol * denotes intra-sublattice couplings, which do not contribute the magnon gap between acoustic and optical modes. Positive (negative) J ij values correspond to FM (AFM) couplings. U = 3 eV and U = 1.36 eV are used in DFT + U with QSGW + U calculations, respectively.J 20 vanishes in all calculations. J 1 by~14% and PJ 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 multisite exchange interactions [56][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 systems 57 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 calculated 40,56 . On the other hand, mapping the total energy of spin-spiral configurations may provide a better description of J ij than the collinear configurations.
Interlayer coupling induced gap opening along Γ-K-M Remarkably, besides the improvement of SW energies, correlations beyond DFT also open up a gap along the Γ-K-M 2 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 Γ-K-M 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-superexchangeJ 2 . If one only considers the Cr sublattice itself, which has a higher symmetry (R3m),J 20 andJ 2 should be equivalent. The presence of the Iodine sublattice breaks the mirror and twofold rotational symmetry, lifting the degeneracy ofJ 20 andJ 2 . Thus, althoughJ 20 andJ 2 connect similar Cr pairs with the same distance of 9.517 Å, their exchange paths are not the same, which allows for J 2 and J 20 to adopt different values. However, up until now, all previous work has considered that J 2 ¼ J 20 ¼ 0, which indeed is supported, to a great extent, by DFT. Surprisingly, the inclusion of correlations beyond DFT results in a sizable AFMJ 2 . As shown in Table 1,J 20 vanishes in all calculations, whereas AFMJ 2 is negligible in DFT but becomes stronger in QSGW and DFT + U, reachingJ 2 ¼ 10%J 1 in QSGW + U. As shown in Fig. 3b,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 20 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 outof-phase precessions of two Cr-spin sublattices, and depends on the inter-sublattice couplings 2|B(q)| (see Eq. (5) and Eq. (6) of the "Methods" section). Within the considered exchange range, we have BðqÞ ¼ J 1 ðqÞ þ J 3 ðqÞ þJ 0 ðqÞ þJ 1þ ðqÞ þJ 20 ðqÞ þJ 2 ðqÞ, where J i (q) is the corresponding Fourier component of J i (R). Couplings J 1 (q), J 3 (q), andJ 1þ ðqÞ are real functions along the Γ-K-M path and vanish at the K-point, resulting a bandcrossing at K if other terms are ignored. The interlayer couplingJ 0 is along the z direction;J 0 ðqÞ ¼J 0 is a real constant when q is in the basal plane, shifting the bandcrossing along the Γ-K-M path. In real space, the connecting vectors (in-plane components) ofJ 20 andJ 2 are rotated by π/6 with respect to those of J 1 , J 3 , andJ 1þ . Correspondingly, in reciprocal space,J 20 ðqÞ andJ 2 ðqÞ are complex functions along the Γ-K-M path (see Supplementary Fig. 4). With J 20 ¼ 0,J 2 ðqÞ itself results in a non-vanishing |B(q)| and thus a gap along the Γ-K-M path. However, ifJ 20 ¼J 2 , then ðJ 20 ðqÞ þJ 2 ðqÞÞ is a real function when q is in the basal plane, shifting the magnon crossing asJ 0 does along the Γ-K-M path. Thus, the combination of vanishingJ 20 and correlation-enhancedJ 2 will induce the magnon gap along the Γ-K-M path. To illustrate, we calculate the SW spectra in a simple J 1 -J 20 -J 2 model withJ 20 ¼J 2 ¼ J 1 =12 and a J 1 -J 2 model withJ 2 ¼ J 1 =6, respectively. As shown in Fig. 3c, the latter case opens a gap along the Γ-K-M path.
However, unlike the global DMI-induced gap, the gap induced by the nonequivalence of the exchange interactionsJ 20 and 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 Γ-K-M lines (at the k z = 0 plane). Using the J ij parameters obtained within QSGW + U, we calculate, within the linear SW theory, the in-plane SW spectra at various k z 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 k z . 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-CrI 3 has different stacking than R-CrI 3 , which dramatically modifies the interlayer supersuperexchange paths and results in A-type AFM ordering in M-CrI 3 . 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-CrI 3 , whereas DFT + U calculations 47,48 have shown that AFM interlayer configurations can be stabilized in M-CrI 3 , depending on the U-value. Within QSGW, we calculate χ q; ω ð Þ in M-CrI 3 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-CrI 3 . 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.
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 CrI 3 . An unexpected correlation-enhanced interlayer super-superexchange induces a gap along the Γ-K-M path in the absence of previously proposed relativistic exchanges or magnon-phonon interactions 23 . 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 CrI 3 . Identifying the existence of the magnon gap along the Γ-K-M path in monolayer CrI 3 , 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, used 61 . Our work suggests the necessity of an explicit treatment of electron correlations to accurately describe the magnetism in the broad family of layered magnetic materials 62,63 , including magnetic topological vdW materials 64,65 , where the interaction between magnetization and the topological surface state is essential.

Pair exchange parameters for Heisenberg model
In adiabatic approximation, we map spin susceptibility into a classical Heisenberg model whereê i is the unit vector of magnetic moment on site i. Effective pair exchange parameters are calculated as where m i (r) is the density of magnetic moment on Cr site i and J(r 1 , r 2 ) is related to χ +− (r 1 , r 2 ) and satisfies Z Ω dr 2 Jðr 1 ; r 2 Þχ þÀ ðr 2 ; r 3 ; ω ¼ 0Þ ¼ δðr 1 À r 3 Þ: Thus, the effective pair exchange parameters can be obtained from the matrix elements of the inverse of spin susceptibilities matrix 31,53,54 with a subsequent Fourier transform, Here, the χ þÀ ij matrices are obtained by projecting onto the functions {m i (r)/||m i ||} representing normalized local spin densities on each magnetic Cr site, which gives a matrix χ ij (q, ω) in magnetic basis site indices 31 . This projection corresponds to the rigid spin approximation, which is a modest approximation for CrI 3 . A 8 × 8 × 8 q-mesh was employed to calculate J ij (R).
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: where A(q) and B(q) are sum of J(q) over the intra-sublattice and intersublattice pairs, respectively: AðqÞ ¼ J 2 ðqÞ þJ 1À ðqÞ BðqÞ ¼ J 1 ðqÞ þ J 3 ðqÞ þJ 0 ðqÞ þJ 1þ ðqÞ þJ 20 ðqÞ þJ 2 ðqÞ: Here, is the Fourier transform of real-space exchange parameters J i with corresponding connecting vectors δ. As shown in Eq. (5), the gap between two modes is ∝ |2B(q)|. Along the Γ-K-M path, J 1 (q), J 3 (q), andJ 1þ ðqÞ are real cosine functions of q and vanish at point K, whereasJ 0 ðqÞ is a real constant andJ 2 ðqÞ a complex number.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.