Abstract
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.
Similar content being viewed by others
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.
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.
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).
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 Γ–K–M
Remarkably, besides the improvement of SW energies, correlations beyond DFT also open up a gap along the Γ–K–M2 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-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 2∣B(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 Γ–K–M 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 Γ–K–M 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 Γ–K–M 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 Γ–K–M 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 Γ–K–M path. Thus, the combination of vanishing \({\tilde{J}}_{{2}_{0}}\) and correlation-enhanced \({\tilde{J}}_{2}\) will induce the magnon gap along the Γ–K–M 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 Γ–K–M 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 Γ–K–M 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.
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 Γ–K–M 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 Γ–K–M 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
where \({\hat{{\bf{e}}}}_{i}\) is the unit vector of magnetic moment on site i. Effective pair exchange parameters are calculated as
where mi(r) is the density of magnetic moment on Cr site i and J(r1, r2) is related to χ+− (r1, r2) and satisfies
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,
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:
where A(q) and B(q) are sum of J(q) over the intra-sublattice and inter-sublattice pairs, respectively:
Here,
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 Γ–K–M 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.
Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.
References
Kotliar, G. & Vollhardt, D. Strongly correlated materials: insights from dynamical mean-field theory. Phys. Today 57, 53–59 (2004).
Wehling, T. O. et al. Strength of effective coulomb interactions in graphene and graphite. Phys. Rev. Lett. 106, 236805 (2011).
van Schilfgaarde, M. & Katsnelson, M. I. First-principles theory of nonlocal screening in graphene. Phys. Rev. B 83, 081409(R) (2011).
Gong, C. et al. Discovery of intrinsic ferromagnetism in two-dimensional van der Waals crystals. Nature 546, 265 (2017).
Huang, B. et al. Layer-dependent ferromagnetism in a van der Waals crystal down to the monolayer limit. Nature 546, 270–273 (2017).
Burch, K. S., Mandrus, D. & Park, J.-G. Magnetism in two-dimensional van der Waals materials. Nature 563, 47–52 (2018).
Gibertini, M., Koperski, M., Morpurgo, A. F. & Novoselov, K. S. Magnetic 2D materials and heterostructures. Nat. Nanotechnol. 14, 408–419 (2019).
Song, T. et al. Switching 2D magnetic states via pressure tuning of layer stacking. Nat. Mater. 18, 1298–1302 (2019).
Klein, D. R. et al. Probing magnetism in 2D van der Waals crystalline insulators via electron tunneling. Science 360, 1218–1222 (2018).
Jiang, S., Shan, J. & Mak, K. F. Electric-field switching of two-dimensional van der Waals magnets. Nat. Mater. 17, 406–410 (2018).
Jiang, S., Li, L., Wang, Z., Mak, K. F. & Shan, J. Controlling magnetism in 2D CrI3 by electrostatic doping. Nat. Nanotechnol. 13, 549–553 (2018).
Deng, Y. et al. Gate-tunable room-temperature ferromagnetism in two-dimensional Fe3GeTe2. Nature 563, 94–99 (2018).
Behera, A. K., Chowdhury, S. & Das, S. R. Magnetic skyrmions in atomic thin CrI3 monolayer. Appl. Phys. Lett. 114, 232402 (2019).
Park, T.-E. et al. Néel-type skyrmions and their current-induced motion in van der Waals ferromagnet-based heterostructures (2019).
Pershoguba, S. S. et al. Dirac magnons in honeycomb ferromagnets. Phys. Rev. X 8, 011010 (2018).
Chen, L. et al. Topological spin excitations in honeycomb ferromagnet CrI3. Phys. Rev. X 8, 041028 (2018).
Owerre, S. A. A first theoretical realization of honeycomb topological magnon insulator. J. Phys. Condens. Matter 28, 386001 (2016).
Kim, S. K., Ochoa, H., Zarzuela, R. & Tserkovnyak, Y. Realization of the haldane-kane-mele model in a system of localized spins. Phys. Rev. Lett. 117, 227201 (2016).
Lee, I. et al. Fundamental spin interactions underlying the magnetic anisotropy in the kitaev ferromagnet cri3. Phys. Rev. Lett. 124, 017201 (2020).
Liu, J., Shi, M., Lu, J. & Anantram, M. P. Analysis of electrical-field-dependent Dzyaloshinskii-Moriya interaction and magnetocrystalline anisotropy in a two-dimensional ferromagnetic monolayer. Phys. Rev. B 97, 054416 (2018).
Ghosh, S., Stojić, N. & Binggeli, N. Structural and magnetic response of CrI3 monolayer to electric field. Phys. B 570, 166–171 (2019).
Ghosh, S., Stojić, N. & Binggeli, N. Comment on “Magnetic skyrmions in atomic thin CrI3 monolayer” [Appl. Phys. Lett. 114, 232402 (2019)]. Appl. Phys. Lett. 116, 086101 (2020).
Kvashnin, Y. O., Bergman, A., Lichtenstein, A. I. & Katsnelson, M. I. Relativistic exchange interactions in CrX3 (X=Cl, Br, I) monolayers. Phys. Rev. B 102, 115162 (2020).
Xu, C., Feng, J., Xiang, H. & Bellaiche, L. Interplay between Kitaev interaction and single ion anisotropy in ferromagnetic CrI3 and CrGeTe3 monolayers. npj Comput. Mater. 4, 57 (2018).
Menichetti, G., Calandra, M. & Polini, M. Electronic structure and magnetic properties of few-layer Cr2Ge2Te6: the key role of nonlocal electron-electron interaction effects. 2D Mater. 6, 045042 (2019).
Lee, Y., Kotani, T. & Ke, L. Role of nonlocality in exchange correlation for magnetic two-dimensional van der Waals materials. Phys. Rev. B 101, 241409 (2020).
Lee, Y., Antonov, V. N., Harmon, B. N. & Ke, L. X-ray spectra in magnetic van der Waals materials Fe3GeTe2, CrI3, and CrGeTe3: a first-principles study. Preprint at https://arxiv.org/abs/2010.15880 (2020).
Chen, L. et al. Magnetic anisotropy in ferromagnetic CrI3. Phys. Rev. B 101, 134418 (2020).
van Schilfgaarde, M., Kotani, T. & Faleev, S. Quasiparticle self-consistent GW theory. Phys. Rev. Lett. 96, 226402 (2006).
Kotani, T., van Schilfgaarde, M. & Faleev, S. V. Quasiparticle self-consistent GW method: a basis for the independent-particle approximation. Phys. Rev. B 76, 165106 (2007).
Kotani, T. & van Schilfgaarde, M. Spin wave dispersion based on the quasiparticle self-consistent GW method: NiO, MnO and α-MnAs. J. Phys. Condens. Matter 20, 295214 (2008).
Karlsson, K. & Aryasetiawan, F. Spin-wave excitation spectra of nickel and iron. Phys. Rev. B 62, 3006–3009 (2000).
Ke, L., van Schilfgaarde, M., Pulikkotil, J., Kotani, T. & Antropov, V. Low-energy coherent stoner-like excitations in CaFe2As2. Phys. Rev. B Rapid Commun. 83, 060404 (2011).
Aryasetiawan, F. & Gunnarsson, O. Product-basis method for calculating dielectric matrices. Phys. Rev. B 49, 16214–16222 (1994).
Anisimov, V. I., Aryasetiawan, F. & Lichtenstein, A. I. First-principles calculations of the electronic structure and spectra of strongly correlated systems: the LDA+U method. J. Phys. Condens. Matter 9, 767–807 (1997).
Pashov, D. et al. Questaal: a package of electronic structure methods based on the linear muffin-tin orbital technique. Comput. Phys. Commun. 249, 107065 (2020).
Kotani, T. Quasiparticle self-consistent GW method based on the augmented plane-wave and muffin-tin orbital method. J. Phys. Soc. Jpn. 83, 094711 (2014).
Besbes, O., Nikolaev, S., Meskini, N. & Solovyev, I. Microscopic origin of ferromagnetism in the trihalides CrCl3 and CrI3. Phys. Rev. B 99, 104432 (2019).
Kashin, I. V., Mazurenko, V. V., Katsnelson, M. I. & Rudenko, A. N. Orbitally-resolved ferromagnetism of monolayer CrI3. 2D Mater. 7, 025036 (2020).
Liechtenstein, A. I., Katsnelson, M. I., Antropov, V. P. & Gubanov, V. A. Local spin density functional approach to the theory of exchange interactions in ferromagnetic metals and alloys. J. Magn. Magn. Mater. 67, 65–74 (1987).
Zhang, W.-B., Qu, Q., Zhu, P. & Lam, C.-H. Robust intrinsic ferromagnetism and half semiconductivity in stable two-dimensional single-layer chromium trihalides. J. Mater. Chem. C. 3, 12457–12468 (2015).
Yoon, H., Jang, S. W., Sim, J.-H., Kotani, T. & Han, M. J. Magnetic force theory combined with quasi-particle self-consistent GW method. J. Phys. Condens. Matter 31, 405503 (2019).
Ke, L. Intersublattice magnetocrystalline anisotropy using a realistic tight-binding method based on maximally localized wannier functions. Phys. Rev. B 99, 054418 (2019).
Mostofi, A. A. et al. An updated version of WANNIER90: a tool for obtaining maximally-localised Wannier functions. Comput. Phys. Commun. 185, 2309–2310 (2014).
Okumura, H., Sato, K. & Kotani, T. Spin-wave dispersion of 3d ferromagnets based on quasiparticle self-consistent GW calculations. Phys. Rev. B 100, 054419 (2019).
McGuire, M. A., Dixit, H., Cooper, V. R. & Sales, B. C. Coupling of crystal structure and magnetism in the layered, ferromagnetic insulator CrI3. Chem. Mater. 27, 612–620 (2015).
Sivadas, N., Okamoto, S., Xu, X., Fennie, C. J. & Xiao, D. Stacking-dependent magnetism in bilayer CrI3. Nano Lett. 18, 7658–7664 (2018).
Soriano, D., Cardoso, C. & Fernández-Rossier, J. Interplay between interlayer exchange and stacking in CrI3 bilayers. Solid State Commun. 299, 113662 (2019).
Soriano, D. & Katsnelson, M. I. Magnetic polaron and antiferromagnetic-ferromagnetic transition in doped bilayer CrI3. Phys. Rev. B 101, 041402 (2020).
Jiang, H., Rinke, P. & Scheffler, M. Electronic properties of lanthanide oxides from the GW perspective. Phys. Rev. B 86, 125115 (2012).
Sponza, L. et al. Self-energies in itinerant magnets: a focus on Fe and Ni. Phys. Rev. B 95, 041112 (2017).
Tomczak, J. M., van Schilfgaarde, M. & Kotliar, G. Many-body effects in iron pnictides and chalcogenides: Nonlocal versus dynamic origin of effective masses. Phys. Rev. Lett. 109, 237010 (2012).
Szczech, Y. H., Tusch, M. A. & Logan, D. E. Finite-temperature magnetism in the Hubbard model. Phys. Rev. Lett. 74, 2804–2807 (1995).
Katsnelson, M. I. & Lichtenstein, A. I. Magnetic susceptibility, exchange interactions and spin-wave spectra in the local spin density approximation. J. Phys. Condens. Matter 16, 7439 (2004).
Ke, L., Belashchenko, K. D., van Schilfgaarde, M., Kotani, T. & Antropov, V. P. Effects of alloying and strain on the magnetic properties of Fe16N2. Phys. Rev. B 88, 024404 (2013).
Turzhevskii, S. A., Lichtenstein, A. I. & Katsnelson, M. I. Degree of localization of magnetic moments and the non-heisenberg nature of exchange interactions in metals and alloys. Sov. Phys. Solid State 32, 1138–1142 (1990).
Logemann, R., Rudenko, A. N., Katsnelson, M. I. & Kirilyuk, A. Exchange interactions in transition metal oxides: the role of oxygen spin polarization. J. Phys. Condens. Matter 29, 335801 (2017).
Wysocki, A. L., Belashchenko, K. D. & Antropov, V. P. Consistent model of magnetism in ferropnictides. Nat. Phys. 7, 485–489 (2011).
Wysocki, A. L., Belashchenko, K. D., Ke, L., van Schilfgaarde, M. & Antropov, V. P. Biquadratic magnetic interaction in parent ferropnictides. J. Phys. Conf. Ser. 449, 012024 (2013).
Kartsev, A., Augustin, M., Evans, R. F. L., Novoselov, K. S. & Santos, E. J. G. Biquadratic exchange interactions in two-dimensional magnets. npj Comput. Mater. 6, 150 (2020).
Senga, R. et al. Position and momentum mapping of vibrations in graphene nanostructures. Nature 573, 247–250 (2019).
Cenker, J. et al. Direct observation of two-dimensional magnons in atomically thin CrI3. Nat. Phys. https://doi.org/10.1038/s41567-020-0999-1 (2020).
Zhang, X. X. et al. Gate-tunable spin waves in antiferromagnetic atomic bilayers. Nat. Mater. 19, 838–842. https://doi.org/10.1038/s41563-020-0713-9 (2020).
Li, B. et al. Competing magnetic interactions in the antiferromagnetic topological insulator MnBi2Te4. Phys. Rev. Lett. 124, 167204 (2020).
Li, B. et al. Two-dimensional ferromagnetism with long-range interactions in the layered magnetic topological insulator MnBi2Te4. https://arxiv.org/abs/2007.08468 (2020).
Acknowledgements
We are indebted to B. Harmon, R.J. McQueeney, A. Kressig, B. Li, J.R. Morris, and M. van Schilfgaarde for fruitful discussions. We are also grateful to L. Chen and P. Dai for providing their INS data. This work was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division, and Early Career Research Program. Ames Laboratory is operated for the U.S. Department of Energy by Iowa State University under Contract number DE-AC02-07CH11358. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract number DE-AC02-05CH11231. The work by MIK is supported by European Research Council via Synergy Grant 854843-FASTCORR.
Author information
Authors and Affiliations
Contributions
L.K. conceived and implemented the research. L.K. and M.K. analyzed the results and wrote the manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Ke, L., Katsnelson, M.I. Electron correlation effects on exchange interactions and spin excitations in 2D van der Waals materials. npj Comput Mater 7, 4 (2021). https://doi.org/10.1038/s41524-020-00469-2
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41524-020-00469-2
This article is cited by
-
Nonlinear Optical Properties of Plasmon–Exciton Nanostructures Based on PbS Quantum Dots and Gold Nanorods in the Field of Nanosecond Laser Pulses
Journal of Russian Laser Research (2024)
-
Control of magnetic states and spin interactions in bilayer CrCl3 with strain and electric fields: an ab initio study
Scientific Reports (2023)
-
A self-adaptive first-principles approach for magnetic excited states
Quantum Frontiers (2023)
-
High-throughput computation and structure prototype analysis for two-dimensional ferromagnetic materials
npj Computational Materials (2022)
-
Environmental screening and ligand-field effects to magnetism in CrI3 monolayer
npj Computational Materials (2021)