Interplay between interlayer exchange and stacking in CrI 3 bilayers
IInterplay between interlayer exchange and stacking in CrI bilayers D. Soriano , , C. Cardoso , , J. Fern´andez-Rossier , QuantaLab, International Iberian Nanotechnology Laboratory (INL),Av. Mestre Jos´e Veiga, 4715-330 Braga, Portugal Radboud University, Institute for Molecules and Materials, NL-6525 AJ Nijmegen, the Netherlands CNR-Nanoscience Institute, S3 Center, 41125 Modena, Italy and Departamento de F´ısica Aplicada, Universidad de Alicante, 03690, Spain
We address the interplay between stacking and interlayer exchange for ferromagnetically orderedCrI , both for bilayers and bulk. Whereas bulk CrI is ferromagnetic, both magneto-optical andtransport experiments show that interlayer exchange for CrI bilayers is antiferromagnetic. BulkCrI is known to assume two crystal structures, rhombohedral and monoclinic, that differ mostlyin the stacking between monolayers. Below 210-220 Kelvin, bulk CrI orders in a rhombohedralphase. Our density functional theory calculations show a very strong dependence of interlayerexchange and stacking. Specifically, the ground states of both bulk and free-standing CrI bilayersare ferromagnetic for the rhombohedral phase. In contrast, the energy difference between bothconfigurations is more than one order of magnitude smaller for the monoclinic phase, and eventuallybecomes antiferromagnetic when either positive strain or on-site Hubbard interactions ( U ≥
3) areconsidered. We also explore the interplay between interlayer hybrydization and stacking, usinga Wannier basis, and between interlayer hybrydization and relative magnetic alignment for CrI bilayers, that helps to account for the very large tunnel magnetoresistance obvserved in recentexperiments. PACS numbers:
I. INTRODUCTION
The recent discovery of several 2D ferromagneticmaterials is expanding the research horizons in 2DMaterials. These new materials, and more specificallyCrI , open new venues for the fabrication of low dimen-sional spintronic and optoelectronic devices basedon multi-layer structures, and are being the object ofstrong interest . Importantly, some of these appli-cations rely on the antiparallel interlayer alignment atzero field, which can be reverted by the application of amagnetic field or, intriguingly, electric fields.The family of CrX (X = Cl, Br, I) magnetic insu-lators is representative of this type of 2D ferromagneticmaterials. In the single layer limit, magnetic order isvery sensitive to magnetic anisotropy, that is governed byanisotropic superexchange in the case of CrI . Ex-periments in bulk show that CrCl is the only one show-ing an anti-ferromagnetic (AF) order , while CrBr andCrI are bulk ferromagnets with Curie temperatures T c = 37 and 61 K, respectively. In contrast with bulk,interlayer coupling for few layer CrI is found to be an-tiferromagnetic, based on optical , transport and, more recently, microscopic probes. This providesa first motivation for this work.The second motivation arises from the following ob-servation. Bulk CrI undergoes a structural transition at 210-220 Kelvin, between a rhombohedral phase at lowtemperature and monoclinic phase at higher tempera-ture. The layer stacking in these structures is shownin Figure 1(a,b). Interestingly, the differential mag-netic susceptibility, dχdT , presents a kink at the structuraltransition , which is consistent with a variation of the interlayer exchange interaction.In monolayer CrI , the intralayer FM coupling, thatultimately drives the long-range ordering between Cratoms, can be anticipated by the Goodenough-Kanamorirules of single-ligand superexchange (M-L-M), on ac-count of the almost perpendicular alignment between theCr-I-Cr bonds. The extension of these rules for morethan one ligand, in order to predict the interlayer ex-change coupling in van der Waals structures (M-L—L-M), is not straightforward. In the following, we assumea different approach aiming to address the interplay be-tween stacking and interlayer exchange coupling for CrI bilayer, combining density functional theory and an ef-fective interlayer coupling model. II. METHODOLOGY
Our calculations are based on density functional the-ory. For each CrI stacking shown in Fig.1(a), we firstperform a geometry relaxation starting from previouslyreported experimental crystal structure ( a = b = 6 . .
866 ˚A for themonoclinic one). The relaxation is carried out usingthe plane-wave based code PWscf as implemented in theQuantum-Espresso ab-initio package . For the self con-sistent calculations, we use a 8 × × k -point grid for thebilayer calculations and a 12 × × forthe exchange-correlation functional are used for Cr and Iatoms. Van der Waals interactions are included throughthe Grimme-D2 model. Spin-orbit coupling is not con-sidered in these calculations. a r X i v : . [ c ond - m a t . m e s - h a ll ] M a r FIG. 1: (Color online) Atomic structure of bilayer CrI with(a) rhombohedral and (b) monoclinic crystal structure. Thesestructures are found in bulk CrI at low and high temperature.Orange and gray atoms correspond to Cr and I, respectively.(c) and (d) show the AB and AA / stacking of the hexag-onal Cr lattices between the layers. The green and red linescorrespond to the top and bottom hexagonal Cr lattices re-spectively. The green and red dots denote the Cr atoms in thetop and bottom layers in the elementary unit cell. (e) Detailof the atomic structure of the I atoms at the bilayer interfacefor the rhombohedral stacking. (f) Idem for the monoclinicone. III. INTERLAYER EXCHANGE ANDSTACKING
Figures 1(c) and 1(d) show the two Cr hexagonallattices (in red and green) and the stacking details ofthe bilayer. In the rhombohedral case, the two Cr lat-tices follows an AB or Bernal stacking, similar to bilayergraphene. The monoclinic case can be obtained by start-ing with an AA stacking and displacing the top layer a /3 along one of the in-plane lattice vectors (cid:126)a or (cid:126)b (wehave labeled this stacking as AA / ). The two atomicstructures induce a different arrangement of the I atomsat the interface which is shown in Figure 1(e,f). Theposition of the I atoms has an important effect on theinterlayer coupling since it affects the Cr-Cr interlayerdistance through steric effects (see Table I).We have performed first-principles calculations of bi-layer (bulk) CrI using both structures, namely AB(rhombohedral) and AA / (monoclinic), with differentinterlayer magnetic order (FM vs AF). Our results aresummarized in table (I). We find that, for both bilayerand bulk CrI the energy difference between the FM andAF configurations is dramatically reduced in the case ofmonoclinic stacking. As we discuss in section V, two different perturbations lead to an antiferromagnetic in-terlayer interaction for the AA / stacking of the bilayer:addition of a Hubbard U correction, keeping the samegeometry obtained without U, and modification of theinterlayer distance.We now discuss how to relate the DFT results to theaverage interlayer exchange coupling ( J ). For thatmatter, we assume that the interlayer exchange can bedescribed by a classical Heisenberg model: H inter = − (cid:88) i ∈ ,j ∈ J inter ij (cid:126)S i · (cid:126)S j (1)where 1 and 2 label the two CrI layers and J ij is arethe interlayer exchange interactions. The sign conven-tion we take is such that J ij > J ij <
0) stands forferromagnetic (antiferromagnetic) interaction. Assum-ing that all spins are parallel, and have a length S , theenergy for configurations where all spins in a given layerare parallel, and collinear with those of the other layer,is U = ∓ S J where J = (cid:88) i ∈ ,j ∈ J inter ij (2)is the average interlayer coupling and the sign − (+)corresponds to the ferromagnetic (antiferromagnetic) in-terlayer alignment. It is self-evident that J is an in-creasing function of the number of atoms in the unit cell.For CrI , there are two atoms per unit cell and plane.We now break down this average exchange, and thecorresponding total interlayer exchange, as a sum overthe contribution coming from each unit cell , J = N j ,where N is the number of unit cells and j is split asthe sum of intracell and intercell contributions. j = (cid:88) i ∈ ( I, ,j ∈ ( I, J intra ij + (cid:88) ∈ ( I, ,j ∈ ( I (cid:48) , J inter ij (3)Since J inter ij decays very rapidly with distance, j con-verges. Therefore, in the case of the bilayer, the totalenergy per unit cell, that can be compared with DFTcalculations, reads as an Ising model for a dimer: U in ( s , s ) = − j S s s (4)where s , s = ± E ≡ U in (+ , − ) − U in (++) = 2 j S (5)We now carry out the same analysis for the case ofbulk, the unit cell has 3 planes. Therefore, the effectivemodel has to keep track of the magnetization of 6 layers: U in ( s , s , s , s , s , s ) = − j S (cid:88) i =1 , s i s i +1 (6) TABLE I: Summary of the calculations for bilayer (AA / vs AB) and bulk (Monoclinic vs Rhombohedral) CrI . The equilibriuminterlayer distance ( d PW in ) is obtained by relaxing the geometry starting from the experimental crystal structure . We keepthe same geometry for the two functionals, namely PBE-D2 and PBE+U-D2. D2 stands for the Grimme-D2 van der Waalscorrection . The energy difference is defined as ∆ E = E AF − E FM . The values of the exchange interlayer coupling areobtained from Equations 5 and 7 for bilayer and bulk respectively. The last column shows the values of j /N at , where N at isthe number of atoms in the unit cell. d PW in (˚A) ∆ E PW (meV) j in ( µ eV) j in ( µ eV/u.c.)Bilayer AA / (PBE-D2) 6.621 0.21 46.7 11.7Bilayer AB (PBE-D2) 6.602 9.43 2095.6 523.9Bilayer AA / (PBE+U-D2) 6.621 -0.36 -80.7 -20.2Bilayer AB (PBE+U-D2) 6.602 17.82 3959.7 989.9Bulk Monoclinic (PBE-D2) 6.621 0.44 16.3 1.4Bulk Rhombohedral (PBE-D2) 6.602 4.50 166.7 13.9 where we assume s = s , to account for the periodicboundary conditions in the off-plane direction. Thus, forthe bulk calculations we have:∆ E ≡ U in (+ − + − + − ) −U in (++++++) = 12 j S (7)Equations (5) and (7) permit to relate our DFT cal-culations with the average interlayer exchange. By sodoing, we find that the interlayer exchange shows alwaysa stronger ferromagnetic character than the AB (rhombo-hedral) phase (see Table I) than the AA / (monoclinic).This clearly indicates that there is a correlation betweenstacking geometry and the interlayer exchange. This isexpected, since different stacking imply both different in-terlayer Cr-Cr distances and Cr-I bond angles (see figure1), which are the structural variables that control ex-change. IV. RELATION WITH EXPERIMENTS
Our results for bulk are consistent with the ferromag-netic interlayer interaction observed experimentally. Inaddition, our results might help to understand the kinkin the magnetic susceptibility observed by McGuire etal. at the structural phase transition observed in bulkCrI at 220 Kelvin, between a low temperature rhom-bohedral and a high temperature monoclinic structures.The connection can be established as follows. At hightemperature, CrI is paramagnetic. The susceptibility ofa ferromagnet in the paramagnetic regime is describedby the Curie law, χ CW = S ( S + 1)( gµ B ) k B T − T CW (8)where S is the total spin of the Cr atoms, g is the g-factor, µ B is the Bohr magneton, k B is the Boltzmann constantand T CW is the Curie temperature, that logically dependson both the interlayer and intralayer couplings throughthe relation 3 k B T CW = S ( S + 1) J (9) where J = (cid:80) j ( J inter ij + J intra ij ) is the sum of all the ex-change interactions for a given spin i . We are assuminghere that the average magnetization of all spins is thesame. The calculated variation of the interlayer couplingfor the two different stackings, shown in the table, willlead to a modification of T CW , ∆ T CW at the tempera-ture of the structural transition, that leads to an addi-tional contribution to dχdT ∝ ∂χ∂T CW ∆ T CW .We now discuss the relation of our results withexperimental results for CrI bilayers and thinfilms , showing antiferromagnetic interlayerinteraction. The application of off-plane magnetic fieldsof B c (cid:39) . T and 0 . T revert the interlayer magnetiza-tion in CrI bilayers as recently proven by transport andoptical measurements respectively. We can estimate theinterlayer exchange by equating the Zeeman energy perunit cell, E Z = 4 × gµ B SB c to ∆ E in eq. (4). We thusobtain j = − gS µ B B c (cid:39) − µeV and − µeV .Our DFT results show that interlayer exchange is muchsmaller for the AA / stacking, although still weakly fer-romagnetic. Other density functional calculations, ap-peared after a first version of our work was posted in thearXiv, show that interlayer exchange can indeed becomeantiferromagnetic for the AA / stacking, using function-als different from GGA . The common point in allDFT calculations is that interlayer exchange has a weakerferromagnetic character for the AA / stacking than forthe AB.Our DFT calculations still predict that the AB stack-ing is the ground state structure for the freestanding bi-layer. However, the energy difference between these twostacking configurations is ∆ E struct ≈ .
25 meV/Cr atom,much smaller than its bulk counterpart, ∆ E struct ≈ bilayers deposited on top of substrates, such asgraphene and silicon oxide, it can be that these favourthe AA / stacking, leading to an antiferromagnetic in-teraction. It is also possible that stacking energetics isdifferent in bulk and in very thin films, on account ofthe different contributions coming from long-range dis-persive forces in both cases. Recent experimental work provides evidence that this might be indeed the case. V. EFFECT OF INTERLAYER DISTANCE ANDON-SITE HUBBARD INTERACTION ( U ) ONINTERLAYER EXCHANGE We now consider two types of perturbations that couldfurther reduce the ferromagnetic interlayer exchange andeventually yield an antiferromagnetic coupling, namelya modification of the interlayer distance and the inclu-sion of an on-site Hubbard interaction ( U ) using the socalled PBE+U functional, in the spirit of the LDA+Uapproximation . The first one could be driven by thecoupling to the substrate, whose effect is missing in ourDFT calculations. For instance, charge transfer is pre-dicted to occur in the graphene/CrI interface , thatcould modify interlayer separation.Figure 2(a) shows the energy difference ∆ E = E AF − E F M for different interlayer distances and for both struc-tures. For ∆
E < / stacking for d − d > . remains FM. The maximal value for the AF exchange, j AF = − µeV is obtained for d − d = 1 ˚A.We now discuss the scaling of ∆ E as we change theon-site Hubbard inderaction U , keeping the same ge-ometry obtained for U = 0. Our results are shown infigure 2(c,d), for the AA / and AB stacking respec-tively. We observe that, in contrast to the AB stacking,the interlayer exchange in the AA / case scales non-monotonically with U and, for U ≥ U correction, follows the same pattern. First,none of these perturbations drive the system to the AFinterlayer interaction in the case of the AB stacking. Sec-ond, both perturbations drive the interlayer interactionAF in the AA / stacking. Third, the dependence of ∆ E on both interlayer distance and U is non-monotonic onlyin the case of the AA / stacking.Given that all known mechanisms for exchange lead tomonotonic dependence with distance, the non-monotonicbehaviour of ∆ E , for the AA / stacking, that includesa change of sign, clearly shows that interlayer exchangeinteraction is the result of at least two contributions withopposite signs: J in = J FMin − J AFin . (10)The first contribution, which favours a ferromagneticcoupling, arises both from interlayer superexchange path-ways and direct exchange. The second contributionfavours antiferromagnetic exchange. FIG. 2: (Color online) (a) Energy difference ∆ E = E AF − E FM of bilayer CrI for AB (orange) and AA / (blue) stack-ing. ∆ E < d − d > . / bilayer becomesAF. (b) Scaling of ∆ E for different values of the Coulombrepulsion ( U ) in AA / stacking. For U ≥ / (mon-oclinic) stacking becomes AF. (c) Idem for the AB (rhombo-hedral) stacking. In contrast to the AA / case, the AB casenever becomes AF when increasin the Coulomb interaction. VI. INTERPLAY BETWEEN INTERLAYERHYBRIDIZATION AND STACKING
We now explore if interlayer antiferromagnetic ex-change could be accounted for by the theory of kineticexchange of Anderson (see also Hay et al. ), for elec-trons occupying otherwise degenerate orbitals that be-come weakly hybridized by an interlayer hopping γ . Theinterlayer hopping leads to the formation of bonding-antibonding states that delocalize the states among thetwo layers. In the limit of on-site Hubbard repulsion U much larger than γ , the low energy states of this Hub-bard dimer are described by a spin Heisenberg modelwith antiferromagnetic exchange J AFin = γ ˜ U . Here, weuse a Hubbard ˜ U to differentiate it from the Hubbard FIG. 3: (Color online) Spin unpolarized band structures ofbilayer CrI for (a) AA / and (b) AB stacking configurations.The half-filled bands correspond to the t g bands of Cr. Thered and blue lines on top of the t g bands corresponds to theWannier bands. U used in LDA+U calculations. The former it is alwayspresent even for U = 0 calculations, and stands for theenergy that should be paid to doubly occupy an atomicorbital. The second one is the extra energy that shouldbe paid when the orbitals are strongly localized in orderto account for correlation effects.In order to explore whether the interlayer hybdridiza-tion γ is significantly different for the two stacking ge-ometries for CrI bilayer, we calculate the hybridizationbetween crystal-field split t g . For that matter we obtaina tight-binding model from our DFT calculations, using arepresentation of the DFT hamiltonian in a basis of max-imally localized Wannier orbitals. To do so, we use DFTas implemented in the plane-wave based PWscf code (see
Methodology section), with spin-unpolarized solutions , toensure that the band splitting comes only from the inter-layer hybridization. In the non-magnetic solutions, the t g bands are half filled, in contrast to the spin-polarizedcase, where the 3 electrons with same spin fill the 3 t g bands in one spin channel. The spin-unpolarized t g bands of bilayer CrI for AA / and AB cases are shownin Figure 3(a,b).In order to obtain a representation of the Hamiltonianin a basis of atomic-like orbitals, we transform our plane-wave basis into a localized Wannier one using wannier90 code . The representation of the Hamiltonian in thatbasis allows us to extract the interlayer hopping ampli-tudes Γ ij from the Wannier Hamiltonian. We choose aprojection over the subspace spanned by the t g manifold,namely { d xy , d xz , d yz } centered in the Cr atoms. Red andblue lines on top of the t g bands in Figure 3(a,b) cor-respond to the Wannier bands obtained for each stack-ing configuration. The Wannier Hamiltonian for intracellatoms takes the form H W = (cid:32) E t g Γ Γ E t g (cid:33) (11) where E and E are 6 × t g orbitals in layer 1 and layer 2. The Γmatrices contain the hopping terms connecting both lay-ers. Equations 12 and 13 show the detailed structure ofthe coupling matrices. The ball and stick models close tothe matrices indicate the unit cell atoms for both stack-ing configurations. Atoms with different colors belong todifferent layers.Γ AB = d xy d xz d yz d xy d xz d yz d xy d xz −
22 0 0 0 0 d yz − d xy −
21 0 0 0 d xz
22 0 22 0 0 0 d yz −
21 22 0 0 0 −
22 (12)Γ AA = d xy d xz d yz d xy d xz d yz d xy − d xz
23 0 0 0 0 0 d yz −
25 0 0 d xy d xz d yz d -type orbitals together with the contribution of the io-dine atoms at the interface makes difficult to comparethis system with a typical single-orbital based model ofbilayer honeycomb lattice. Comparing both hopping ma-trices, we observe that the higher contribution to the an-tiferromagnetic kinetic exchange in the AB case comesonly from the interaction between d yz orbitals in atoms1 and 4 (Γ AB = −
29 meV). In contrast, the AA / inter-layer hopping matrix shows two important contributions(Γ AA = Γ AA = 29 meV) between orbitals d xy and d xz .However, from this analysis, we can not conclude thatthe average interlayer hybridization is very different forthe two stackings. Therefore, the mechanism that ac-counts for the different interlayer exchange interactionmust arise from other exchange mechanism. VII. SPIN POLARIZED ENERGY BANDS ANDIMPLICATIONS FOR VERTICAL TRANSPORT
We now discuss the spin-polarized band structures ob-tained from first-principles calculations. In Figure 4, weshow the spin polarized bands of the AB (top panels) andAA / (bottom panels) bilayer CrI . For the spin polar-ized calculations, the spin majority t g bands are fullyoccupied and the first set of empty bands is made of spinmajority e g states. The FM cases (a,c) show a clear band FIG. 4: (Color online) Spin-polarized band structures of bi-layer CrI for (a) AB(FM), (b) AB(AF), (c) AA / (FM) and(d) AA / (AF). splitting of the majority e g (red bands above the Fermienergy) and minority t g bands (blue bands above theFermi energy) coming from the interlayer coupling. Thisis similar to the unpolarized calculations in Figure 3.In contrast, in the antiparallel interlayer (AF) cases(b,d), the interlayer splitting is absent. The reason isthat e g and t g for a given spin channel in one layer aredegenerate with the same bands with opposite spin in theother layers. Since interlayer coupling is spin conserv-ing, the resulting hybridization is dramatically reduced.This difference of interlayer coupling in the FM and AFconfigurations definitely contributes to explain the verylarge magnetoresistance observed in vertical transportwith CrI bilayers in the barriers. Given that the e g states are the lowest energy channels available for tunnel-ing electrons in the barrier, they can only be transferredelastically between adjacent CrI layers when their rela-tive alignment is not antiparallel. VIII. DISCUSSION AND CONCLUSIONS.-
Our results provide a plausible explanation for two dif-ferent experimental observations. First, the kink in the differential spin susceptibility dχdT observed in bulk CrI at the structural phase transition . Second, the anti-ferromagnetic interlayer coupling observed for few-layerCrI bilayers, at low temperatures . Given that inthese experiments the CrI layers are either deposited ona substrate or embedded in a circuit, we conjecture thatthe stacking might be different than in bulk. However,the confirmation of this hypothesis will require furtherexperimental and computational work.To summarize, we have computed the interlayer ex-change for CrI bilayers in two different stacking, thatcorrespond to the rhombohedral and monoclinic struc-tures observed for bulk CrI . We find that the interlayercoupling shows a much weaker ferromagnetic characterfor the monoclinic than the rhombohedral phase, andeventually undergoes an antiferromagnetic transition un-der shear strain. We claim that this provides a possibleexplanation for two different experimental observations.First, the kink observed in the differential susceptibilityat the structural transition in bulk. Second, the factthat CrI bilayers deposited on graphene are known tohave antiferromagnetic interlayer coupling. Note:
During the final completion of this manuscript,we became aware of the work of four other groups addressing the relation between interlayer exchange andstacking in CrI using different approximations and ob-taining results compatible with ours. Acknowledgements.-
We acknowledge Efr´en Navarro-Moratalla for pointing out the different interlayer cou-pling for bilayer and bulk CrI . We thank Jos´e LuisLado, Francisco Rivadulla , A. H. MacDonald and JeilJung for fruitful discussions. DS thanks NanoTRAIN-forGrowth Cofund program at INL and the financialsupport from EU through the MSCA Individual Fel-lowship program at Radboud University. J. F.-R. ac-knowledge financial support from FCT for the P2020-PTDC/FIS-NAN/4662/2014, the P2020-PTDC/FIS-NAN/3668/2014 and the UTAPEXPL/NTec/0046/2017projects, as well as Generalitat Valenciana fundingPrometeo2017/139 and MINECO Spain (Grant No.MAT2016-78625-C2). CC and JFR acknowledge FEDERproject NORTE-01-0145-FEDER-000019. The authorsthankfully acknowledge the computer resources at Cae-saraugusta and the technical support provided by the In-stitute for Biocomputation and Physics of Complex Sys-tems (BIFI) (RES-QCM-2018-2-0032). Part of this workwas carried out on the Dutch national e-infrastructurewith the support of SURF Cooperative. C. Gong, L. Li, Z. Li, H. Ji, A. Stern, Y. Xia, T. Cao,W. Bao, C. Wang, Y. Wang, et al., Nature , 265 (2017). B. Huang, G. Clark, E. Navarro-Moratalla, D. R. Klein,
R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A.McGuire, D. H. Cobden, et al., Nature , 270 (2017). M. A. McGuire, Crystals , 121 (2017). Z. Fei, B. Huang, P. Malinowski, W. Wang, T. Song,J. Sanchez, W. Yao, D. Xiao, X. Zhu, A. May, et al., Na-ture materials , 778 (2018). D. J. O’Hara, T. Zhu, A. H. Trout, A. S. Ahmed, Y. K.Luo, C. H. Lee, M. R. Brenner, S. Rajan, J. A. Gupta,D. W. McComb, et al., Nano Letters , 3125 (2018). T. Song, X. Cai, M. W.-Y. Tu, X. Zhang, B. Huang, N. P.Wilson, K. L. Seyler, L. Zhu, T. Taniguchi, K. Watanabe,et al., Science , 1214 (2018). D. R. Klein, D. MacNeill, J. L. Lado, D. Sori-ano, E. Navarro-Moratalla, K. Watanabe, T. Taniguchi,S. Manni, P. Canfield, J. Fern´andez-Rossier, et al., Science , 1218 (2018). Z. Wang, I. Guti´errez-Lezama, N. Ubrig, M. Kroner,M. Gibertini, T. Taniguchi, K. Watanabe, A. Imamo˘glu,E. Giannini, and A. F. Morpurgo, Nature Communications , 2516 (2018). K. L. Seyler, D. Zhong, D. R. Klein, S. Gao, X. Zhang,B. Huang, E. Navarro-Moratalla, L. Yang, D. H. Cobden,M. A. McGuire, et al., Nature Physics , 277 (2018). D. Zhong, K. L. Seyler, X. Linpeng, R. Cheng, N. Sivadas,B. Huang, E. Schmidgall, T. Taniguchi, K. Watanabe,M. A. McGuire, et al., Science Advances , e1603113(2017). B. Huang, G. Clark, D. R. Klein, D. MacNeill, E. Navarro-Moratalla, K. L. Seyler, N. Wilson, M. A. McGuire, D. H.Cobden, D. Xiao, et al., Nature Nanotechnology , 544(2018). S. Jiang, J. Shan, and K. F. Mak, Nature materials ,406 (2018). S. Jiang, L. Li, Z. Wang, K. F. Mak, and J. Shan, NatureNanotechnology , 549 (2018). J. L. Lado and J. Fern´andez-Rossier, 2D Materials ,035002 (2017). J. Liu, M. Shi, J. Lu, and M. P. Anantram, Phys. Rev. B , 054416 (2018). Y. Liu and C. Petrovic, Phys. Rev. B , 014420 (2018). N. Richter, D. Weber, F. Martin, N. Singh, U. Schwingen-schl¨ogl, B. V. Lotsch, and M. Kl¨aui, Phys. Rev. Materials , 024004 (2018). J. Zhang, B. Zhao, T. Zhou, Y. Xue, C. Ma, and Z. Yang,Phys. Rev. B , 085401 (2018). K. Zollner, M. Gmitra, and J. Fabian, New J. Phys. ,073007 (2018). P. Jiang, L. Li, Z. Liao, Y. Zhao, and Z. Zhong, Nanoletters (2018). J. Liu, M. Shi, P. Mo, and J. Lu, AIP Advances , 055316(2018). C. Cardoso, D. Soriano, N. A. Garc´ıa-Mart´ınez, andJ. Fern´andez-Rossier, Phys. Rev. Lett. , 067701 (2018). Q. Tong, F. Liu, J. Xiao, and W. Yao, Nano Letters ,7194 (2018). I. Lee, F. G. Utermohlen, K. Hwang, D. Weber, C. Zhang,J. van Tol, J. E. Goldberger, N. Trivedi, and P. C. Hammel,arXiv e-prints arXiv:1902.00077 (2019), 1902.00077. O. Besbes, S. Nikolaev, and I. Solovyev, arXiv e-printsarXiv:1901.09525 (2019), 1901.09525. C. Xu, J. Feng, H. Xiang, and L. Bellaiche, npj Computa-tional Mathematics , 57 (2018), 1811.05413. B. Morosin and A. Narath, The Journal of ChemicalPhysics , 1958 (1964). L. L. Handy and N. W. Gregory, Journal of the AmericanChemical Society , 891 (1952). M. A. McGuire, H. Dixit, V. R. Cooper, and B. C. Sales,Chemistry of Materials , 612 (2015). D. R. Klein, D. MacNeill, Q. Song, D. T. Larson, S. Fang,M. Xu, R. A. Ribeiro, P. C. Canfield, E. Kaxiras,R. Comin, et al., arXiv e-prints arXiv:1903.00002 (2019),1903.00002. L. Thiel, Z. Wang, M. A. Tschudin, D. Rohner,I. Guti´errez-Lezama, N. Ubrig, M. Gibertini, E. Gian-nini, A. F. Morpurgo, and P. Maletinsky, arXiv e-printsarXiv:1902.01406 (2019), 1902.01406. J. B. Goodenough, Phys. Rev. , 564 (1955). J. B. Goodenough, Journal of Physics and Chemistry ofSolids , 287 (1958). J. Kanamori, Journal of Physics and Chemistry of Solids , 87 (1959). S. Feldkemper and W. Weber, Phys. Rev. B , 7755(1998). P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car,C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni,I. Dabo, et al., Journal of Physics: Condensed Matter ,395502 (2009). J. P. Perdew, K. Burke, and M. Ernzerhof, Physical reviewletters , 3865 (1996). S. Grimme, Journal of computational chemistry , 1787(2006). P. Jiang, C. Wang, D. Chen, Z. Zhong, Z. Yuan, Z.-Y. Lu,and W. Ji, ArXiv e-prints (2018), 1806.09274. N. Sivadas, S. Okamoto, X. Xu, C. J. Fennie, and D. Xiao,Nano Letters , 7658 (2018). S. W. Jang, M. Y. Jeong, H. Yoon, S. Ryee, and M. J.Han, ArXiv e-prints (2018), 1809.01388. C. Lei, B. Lingam Chittari, K. Nomura, N. Baner-jee, J. Jung, and A. H. MacDonald, arXiv e-printsarXiv:1902.06418 (2019), 1902.06418. V. I. Anisimov, J. Zaanen, and O. K. Andersen, PhysicalReview B , 943 (1991). P. W. Anderson (Academic Press, 1963), vol. 14 of
SolidState Physics , pp. 99 – 214. P. J. Hay, J. C. Thibeault, and R. Hoffmann, Journal ofthe American Chemical Society , 4884 (1975). A. A. Mostofi, J. R. Yates, G. Pizzi, Y. S. Lee, I. Souza,D. Vanderbilt, and N. Marzari, Comput. Phys. Commun.185