Environmental Screening and Ligand-Field Effects to Magnetism in CrI$_3$ Monolayer

We present a detailed study on the microscopic origin of magnetism in suspended and dielectrically embedded CrI$_3$ monolayer. To this end, we down-fold two distinct minimal generalized Hubbard models with different orbital basis sets from \emph{ab initio} calculations using the constrained random phase approximation. Within mean-field approximation, we show that these models are capable of describing the formation of localized magnetic moments in CrI$_3$ and of reproducing electronic properties of full \emph{ab initio} calculations. We utilize the magnetic force theorem to study microscopic magnetic exchange channels between the different orbital manifolds. We find a multi-orbital super-exchange mechanism as the origin of magnetism in CrI$_3$ resulting from a detailed interplay between effective ferro- and anti-ferromagnetic Cr-Cr $d$ coupling channels, which is decisively affected by the ligand (I) $p$ orbitals. We show how environmental screening such as resulting from encapsulation with hexagonal boron nitride (hBN) of the CrI$_3$ monolayer affects the Coulomb interaction in the film and how this successively controls its magnetic properties. Driven by a non-monotonic interplay between nearest and next-nearest neighbour exchange interactions we find the magnon dispersion and the Curie temperature to be non-trivially affected by the environmental dielectric screening.


I. INTRODUCTION
Ferromagnetic (FM) layered materials hold high promises for becoming one of the key ingredients in future spintronic nanodevices 1-4 based on van-der-Waals heterostructures. Recent observations of magnons in thin layers of chromium trihalides via inelastic electron tunneling spectroscopy (IETS) 1,5 and magneto-Raman spectroscopy 6 have also opened new ways to explore magnon-based low-consumption spintronics at the twodimensional (2D) limit 7 . Whether or how the magnetic properties are modified or stabilized in these heterostructures is, however, still under debate. Recently reported theoretical predictions suggest, for instance, a strong dependence of the magneto-optical response on the film thickness and the spin-orbit coupling of CrI 3 8 , while recent experimental studies point towards the possibility to tune the magnetic properties of monolayer and bilayer CrI 3 using electrostatic gating [9][10][11][12] . These observations highlight the importance of addressing how the environment of layered magnetic systems modify their electronic and magnetic properties via proximity, gating, or screening effects.
For the case of chromium trihalides, CrX 3 (X = I, Br Cl), the nearest neighbour magnetic exchange couplings have been theoretically reported to be between 1 and 3.2 meV 13-18 depending on the material and calculation scheme, which renders rigorous theoretical descriptions absolutely necessary. This holds especially with respect to the Coulomb interactions in layered material, which are enhanced due reduced polarization in the surrounding. At the same time changes to the environmental polarization and thus to the environmental screening can significantly modify the Coulomb interaction within the layered material. Thus, properties of layered materials affected by the Coulomb interaction, such as the formation of excitons and plasmons or the sta-bilization of charge and superconducting order, are in principle strongly dependent on the material's environment allowing for Coulomb engineering of these manybody properties [19][20][21][22][23][24] . So far it is not known how efficient Coulomb engineering can be for tailoring the magnetic ground state or magnonic excitations of layered magnetic systems. For the case of CrI 3 both can be described by effective spin Hamiltonians describing exchange interactions between the Cr atoms as the spin carriers. Thus any changes from environmental screening to the exchange couplings of this effective Hamiltonian can give rise to changes in the magnetic transition temperature or the magnon spectrum opening new routes towards the control and potential tailoring of magnetic properties of chromium trihalides.
In the following, we study in detail the impact of Coulomb interactions on the magnetic properties of suspended and dielectrically embedded monolayer CrI 3 . We start from spin-unpolarized ab initio band structure calculations for monolayer CrI 3 using density functional theory (DFT), which we use to down-fold two minimal models defined by localized Wannier orbitals describing the Cr d states only and a combined (d+p)-basis also taking the I p states into account. After constructing the corresponding Coulomb tensors from constrained Random Phase Approximation (cRPA) calculations we solve the resulting multi-orbital Hubbard models within meanfield Hartree-Fock (HF) theory. The resulting interacting and spin-resolved quasi-particle band structures are used afterwards to evaluate the orbitally resolved magnetic exchange interactions using the magnetic force theorem (MFT). Using our Wannier Function Continuum Electrostatics (WFCE) approach we additionally include the environmental screening effects to the Coulomb interactions, which finally allows us to study how the microscopic exchange coupling channels, the magnetic transition temperature, and the magnon dispersions are af-fected by different parts of the Coulomb tensor and/or by the dielectric encapsulation of the CrI 3 monolayer. In the crystal structure of monolayer CrI 3 the Cr atoms are arranged on a honeycomb lattice surrounded by the I ligands forming an edge-sharing octahedral environment around each metal ion, as depicted in Fig. 1 (a). The ligand field splits the Cr d-orbitals into two sets, namely, the t 2g (d xy , d xz , d yz ) and the e g (d x 2 −y 2 , d z 2 ) 13 . Spin unpolarized first-principles calculations in the Generalized Gradient Approximation (GGA) level predict a metallic ground state with half-filled bands of predominantly Cr t 2g character separated by sizable gaps from fully occupied bands of mostly I p character and unoccupied bands of Cr e g type, as shown in the orbital resolved density of states in Fig. 1 (b). Spin-polarized calculations in the GGAs+U approximation 25 shift down the majority spin t 2g bands which hybridize with the ligand p-bands. The unoccupied majority spin e g bands appear at approx. 1.5 eV above the valence states followed by the mi-nority spin t 2g and e g bands, respectively, resulting in a FM insulating state, as depicted in Fig. 1 (c). Although LSDA+U approaches seem to reasonably describe the magnetic properties of CrI 3 , its validity still needs to be benchmarked with higher-level theories such as dynamical mean field theory, as the local Coulomb repulsion U is rather large compared to the t 2g band width as discussed in detail later. Thus, to fully microscopically understand each part of the problem, ranging from hybrdizations effects, to the Coulomb tensor, or the choice of basis set we proceed in the following with a down-folding procedure to generate reliable material-realistic minimal models.

B. Model Hamiltonian
In the following we will construct model Hamiltonians of the form while the Coulomb term takes local density-density U mm and Hund's exchange J H mm interactions on the Cr atoms into account and is defined by wheren σ = d † mσ d mσ is the number operator and σ is the spin index. These Hamiltonians are solved within a mean-field HF ansatz and analyzed in terms of its microscopic magnetic properties as described in the Methods section.

C. Basis Sets and Non-Interacting Hamiltonians
Motivated by the orbital-resolved ab initio nonmagnetic density of states depicted in Fig. 1 (b), we utilize two minimal basis sets: one including only effective Cr d Wannier orbitals and one also including I p states, which are referred to as d-only and (d+p) model, respectively, in the following. To this end, we start with conventional DFT calculations utilizing the Perdew-Burke-Ernzerhof GGA exchange correlation functional 26 within a PAW basis 27,28 as implemented in the Vienna Ab initio Simulation Package (vasp) 29,30 for CrI 3 monolayers with a lattice constant of a 0 ≈ 6.97Å embedded in a supercell with a height of 35Å. We use (16 × 16 × 1) k-meshes together with an energy cutoff of 230 eV and apply a Methfessel-Paxton smearing of σ = 0.02 eV. For the d-only model we project the Kohn-Sham states onto Cr-localized d orbitals (with a rotated axis parallel to the tetragonal main axis) and maximally localize them afterwards using the wannier90 package 31 . For the (d+p)-model, we start with the same Cr-centered d orbitals and add also I-centered p orbitals to the list of initial projects. In this case we do not perform a maximal localization since it would increase the localization of the I p orbitals at the cost of a delocalization of Cr d orbitals. The resulting localized Wannier orbitals are used in a subsequent step to calculate the needed hopping matrix elements for the definition of H W . The resulting Wannier models perfectly interpolate the mostly Cr d Kohn-Sham states in the d-only model as well as also the I p states in the (d+p)-model.

D. Constrained Random Phase Approximation
The Coulomb interaction matrix elements of H U are evaluated using the Cr d Wannier orbitals from these two models within the cRPA scheme 32 according to Here U represents the partially screened Coulomb interaction defined by with V being the bare Coulomb interaction and π rest the partial or rest-polarization from all states expect those of the correlated sub-space defined by the Cr d states. π rest thus describes screening from all other Cr as well as all I states including a significant amount of empty states from the full Kohn-Sham basis. In detail, we use in total 128 bands and define π rest by explicitly excluding all Kohn-Sham states between −0.5 and 3 eV from the full polarization. To this end, we use a recent cRPA implementation by Kaltak   We start with analyzing the bare and cRPA screened density-density Coulomb matrix elements as obtained from the d-only basis. The full matrices are shown in Tab. I together with the corresponding Hund's exchange elements. In all cases the resulting density-density matrices are approximately of the form where v t (j H t ) and v e (j H e ) are intra-orbital densitydensity (inter-orbital Hund's exchange) matrix elements within the t 2g and e g manifolds, respectively, v et = 1/2(v e + v t ), and j H 1...4 are inter-orbital Hund's exchange elements between the two manifolds. This form of the density-density Coulomb matrix is similar to the one obtained for a fully rotational-invariant d shell which is here, however, perturbed due to the ligand-induced crystal-filed splitting. Therefore, instead of five (U 0 , J H 1...4 ) or even just three (U 0 , F 2 , F 4 34 ) parameters to represent the full density-density matrix, we need here eight. As a result, the t 2g and e g channels themselves are easily parameterized using two Hubbard-Kanamori parameters (v and j H ). The inter-channel elements show, however, a significant orbital dependence which can be represented by the four j H 1...4 Hund's exchange elements. Since the Wannier orbital spread Ω α = w α |r 2 |w α − w α |r|w α 2 31 of the t 2g wave functions (Ω t2g ≈ 3Å 2 ) is smaller then the corresponding e g one (Ω eg ≈ 5.3Å 2 ) the bare intra-orbital density-density elements v t ≈ 15.5 eV are larger than the v e ≈ 12.3 eV elements and also represent the largest elements in v mm . This is also reflected in the intra-channel j H t ≈ 0.54 eV and j H e ≈ 0.51 eV elements which are, however, still rather similar. The interchannel j H 1...4 vary between about 0.25 eV and 0.47 eV. We note that these bare Hund's exchange interactions are significantly smaller than the approximated values in bulk Cr of J H ≈ 0.75 eV 35 , which we attribute here to the enhanced Wannier function spread of the effective Cr d-only basis.
The cRPA density-density matrix elements U mm are significantly reduced by factors between 1/4 to nearly 1/6 by the rest-space screening from the other Cr and I states. Although this effective screening is strongly orbital dependent, it does not change the overall orbital structure of the Coulomb matrix depicted in Eq. (6). We find U t ≈ 3.6 eV and U e ≈ 3.1 eV, which are still rather large compared to the band-width of the half-filled Cr t 2g band of about 1 eV. Thus even taking cRPA screening into account correlation effects can be expected to play an important role. The screened Hund's exchange interactions, J H t ≈ 0.49 eV and J H e ≈ 0.44 eV, are reduced by no more than 10% in comparison to the bare values, which also holds for the inter-channel ones.
We proceed with analyzing how these local cRPAscreened Coulomb interactions affect the band structure of the Cr d-only model within mean-field theory. As depicted in Fig. 2 the Coulomb interactions drive the systems into an insulator with a sizeable band gap and the same band ordering as known from LSDA(+U ) calculations (i.e. fully occupied t 2g followed by completely empty e g,↑ , t 2g,↓ , and e g,↓ manifolds). On a qualitative level this minimal basis thus seems to be capable of reproducing the Cr d band structure of full ab initio calculations. To analyze to what extend this minimal model is also capable of reproducing the super-exchange mechanism responsible for the FM ordering in CrI 3 (see Methods) we have applied the MFT to these HF solutions in a subsequent step (for details see Methods). The resulting orbitally-resolved exchange couplings are given in Tab. II and show a strong AFM total coupling of J = −1.964 meV for nearest neighbours and J = −0.042 meV for next nearest neighbours, driven by a strong AFM J t2g-t2g interaction. This is obviously in contradictions to LSDA(+U ) calculations and to all available experimental data. This wrong prediction can be understood from the Kugel-Khomskii (KK) formalism 15,36 , which describes the total magnetic exchange as the sum of the t 2g -t 2g and t 2g -e g contributions as Here,t t2g-t2g (r) andt t2g-eg (r) are effective hopping matrix elements between the different orbital channels, ∆ represents the crystal field splitting, andŨ andJ H are the averaged Coulomb and Hund's exchange parameters. We immediately understand that J t2g-t2g is by definition of AFM nature and controlled only by the effective t 2g -t 2g hopping and the density-density channel of the Coulomb interaction (Ũ ). In contrast, the FM channel, J eg-eg , is controlled by the effective t 2g -e g hopping and bothŨ and J H . By using the hopping values of our Wannier construction, given in Tab. VI of the Supplement, and the averaged electron-electron interactions from the cRPA Coulomb tensor (Ũ ≈ 3.4 andJ H ≈ 0.4 eV), we obtain an AFM intralayer exchange coupling in line with the MFT results. We carefully checked that this wrong KK prediction also holds in the case of the extended (d+p) basis by using effective hoppings from this models with and without p-mediated hopping, as summarized in Tab. VI of the Supplement. In both cases, using only the direct hoppings as well as those after integrating our the I p contributions, we find that the effective n.n.t t2g-t2g is significantly larger than thet t2g-eg one so that the AFM coupling channels always dominates. Thus the correct microscopic origin of the FM coupling in CrI 3 cannot be described in a Cr d-only basis and can neither be modelled within the KK approach.

F. Extended (d+p)-basis
Motivated by the wrong predictions of the minimal Cr d-only model, we expanded the basis to also include the so that the resulting bare matrix elements are significantly enhanced as compared to the d-only case. Also the density-density matrix elements in the e g channel are now larger than in the t 2g channel. The cRPA screening due to all other Cr and all I states, again significantly reduce all matrix elements in an orbital-dependent manner. The final intra-orbital density-density interactions are now on the order of 4 eV with inter-orbital elements of the order of 3 eV. Notably, the cRPA screened Hund's exchange interactions are also enhanced by the increased localization of the Cr d states. They are, however, still on the order of 0.5 to 0.7 eV and thus still significantly smaller than the common approximation of 0.9 eV. As before in the d-only basis the interaction term H U of our Hamiltonian from Eq. (1) acts on the Cr d states only. The kinetic part, H W , however now also includes I p contributions. Thus, to counter-act double counting effects we subtract here the double counting potential defined in Eq. (8) from the Cr d states using a nominal Cr 3+ d occupation of N σ imp = 3. This corrects for the relative positioning of the interacting Cr d bands with respect to the uncorrelated I p states.
The resulting HF density of states is shown in Fig. 3. In contrast to the d-only model we now find a band structure which is vastly reminiscent of the full ab initio GGAs+U results shown in Fig. 1. Next to the spin-splitting and -ordering also the full and sub band gaps are in good agreement with GGAs+U predictions.
The orbital-resolved intralayer exchange couplings calculated using the MFT are given in Tab explore how the different Coulomb interaction channels (U , U , J H ) affect the electronic band structure and the magnetic properties, we present in Fig. 4 the band widths and positions together with the orbitally-resolved exchange interactions upon individual reductions of each Coulomb channel by 20%. Reducing the intra-orbital Coulomb interactions (∆U ) results in a lower spin splitting of the t 2g and e g bands, similar to the single-orbital Hubbard model. The other two cases, namely, ∆U and ∆J H , are more difficult to understand, since these terms mix different types of orbitals. When reducing U , we observe a strong reduction of the splitting between the t 2g,↑ /e g,↑ manifolds, accompanied by a strong gap reduction of about 1 eV. In contrast, the splitting between the t 2g,↓ /e g,↓ manifolds is weakly affected by U compared to the initial full cRPA case. Finally, when we reduce Hund's exchange J H , the t 2g,σ /t 2g,σ and the e g,σ /e g,σ splittings are both reduced, while the splitting between the majority spin occupied t 2g,↑ bands and empty e g,↑ bands increases with respect to the unperturbed case, leading to a small increase in the bandgap.
These band-structure renormalizations have nontrivial effects to the microscopic exchange interactions, which we depict in Fig. 4 (bottom panel) and which can be just partially understood within the KK formalism. In line with KK a reduction of U or U yields an enhancement of the AFM J t2g-t2g channel, as shown in the ∆U and ∆U columns of Fig. 4 (bottom panel), while modifications to J H do not affect J t2g-t2g at all. Also in qualitative agreement with KK we find that a reduction of J H reduces the FM J t2g-eg . In contrast to the simple model predictions from KK we find from MFT that J t2g-t2g is slightly reduced upon reducing U while modifications to U yields no modifications. This underlines the need for the full microscopic MFT in combination with an extended basis to quantitatively understand magnetism in CrI 3 . We expect that just an extended model combining the super-exchange mechanism from the Goodenough-Kanamori description and the multi-orbital KK mechanism might yield a qualitatively understanding in line with our MFT findings, which will be studied separately. Here, we conclude that the extended (d+p)-basis together with the corresponding cRPA Coulomb matrix elements and using a HF solver result in a reliable band-structure and microscopically correct magnetic properties.

G. Substrate Tunability
In the following we proceed with the analysis of dielectric screening effects to the magnetic properties of CrI 3 monolayer. We investigate how a dielectric encapsulation, as depicted in the inset of Fig. 5 (a), modifies the Coulomb interactions in CrI 3 and how this affects its bands structure and eventually its microscopic magnetic properties. Fig. 5 (a) summarizes how the orbital averaged local intra-and inter-orbital density-density as well as Hund's exchange matrix elements, scale upon increasing the environmental screening ( ). Since affects the macroscopic monopole-like interactions only, U and U are equally reduced by , while J H is not modified at all. The resulting screening-induced effects will thus be a combination of the ∆U and ∆U columns from Fig. 4. In Fig. 5 (b) we show the HF density of states (DOS) for three different values of . As discussed in the previous section, a decreasing U results in a reduction of the dbands splitting ∆ t2g,σ/t2g,σ and ∆ eg,σ/eg,σ and a decrease of the inter-orbital Coulomb terms U reduces the gap between t 2g,↑ and e g,↑ . The overall effect of increasing is thus to decrease all gaps between all sub-bands.
In Fig. 5 (c) we summarize the resulting effects to the magnetic properties. Upon increasing we find the total nearest-neighbour magnetic exchange decreasing from about J ≈ 1.1 meV at = 1 to J ≈ 0.2 meV at = 20, while the total next-nearest-neighbour exchange interaction slightly increases from J ≈ 0.3 meV to J ≈ 0.7 meV. From Tab. IV we understand that the decreasing trend in J is mostly driven by the enhancement of the AFM coupling within the J t2g-t2g channel, while the FM J t2g-eg channel is barely affected. The nextnearest-neighbour J is affected simultaneously by both, increasing FM and AFM microscopic exchange interactions. The FM channel grows, however, slightly faster so that the overall FM J is enhanced. From this we expect non-trivial effects of the environmental dielectric screening to the magnetic properties of CrI 3 monolayer relevant for most experimental setups dealing with supported and/or encapsulated films.
We start with the analysis of the magnon dispersion which reacts to the environmental screening differently in different parts of the Brillouin zone, as shown in Fig. 5 (d). At the Γ point the optical (high energy) branch is continuously lowered in energy upon increasing the screening, while the Dirac point at K is nonmonotonously affected. Upon increasing from 1 to 8 the magnon energy at K first increases before it starts to decrease for larger . This behaviour becomes clear from the spin-wave Hamiltonian Eq. (12) evaluated at K yielding the degenerated magnon-energies E(K) = 3S(J + 3J + λ). As one can see from Fig. 5 (c), J + 3J is indeed a non-monotonic function of with a maximum around = 6. For the van-Hove singularities of both magnon branches at the M point we find similar nontrivial and non-monotonic behaviours with yielding a partially extended flat dispersion of the optical branch for intermediate . These non-trivial modifications to the magnon dispersion due to changes in are also reflected in the total magnon spectrum. Thereby we can relate each (partial) maximum in the magnon spectrum with either the Γ or the M point. As a result, these most prominent spectral features either monotonously decrease in energy or follow the non-monotonous trends from the M point.
In Figure Fig. 5 (e), we additionally show the Curie temperature (T C ) for the same . Again, we find a nonmonotonic behavior with an initially increasing T C upon increasing the screening, a maximal T C around = 6, and a subsequently strongly decreasing trend. This trend approximately follows the spectral peak arising from the optical magnon branch at the M point. Due to its similarity with the K-point behaviour we conclude that the initial increasing trend in T C is driven by the increasing trend of J while the final decreasing trend is driven by J. The non-monotonic behavior of T C is thus due to the non-monotonic interplay between nearest and nextnearest neighbor exchange interactions as a function of the environmental screening.

III. CONCLUSIONS & OUTLOOK
By combining state-of-the-art cRPA-based ab initio down-folding with our WFCE approach and the MFT method, we were able to study on a microscopic level how magnetism in CrI 3 monolayer builds up and how it is controlled by environmental screening properties. We showed that a mean-field description within the HF approximation to treat local Cr Coulomb interactions is sufficient to reproduce all characteristics of full ab initio GGAs+U calculations. From a thorough investigation of different minimal models we understood that only an I p-based super-exchange mechanism together with the full multi-orbital t 2g -e g structure of the Cr atoms allows for a realistic description of CrI 3 magnetism. A minimal model thus needs to involve all Cr d and I p states with sizable Coulomb interactions acting on the Cr d states. We also showed how environmental screening significantly reduces the local Coulomb interactions, which decisively modifies the electronic band structure and which finally yields non-monotonic changes to the microscopic magnetic exchange interactions. In detail we found that dielectric encapsulation of the CrI 3 monolayer strongly reduces the nearest-neighbour exchange interaction, while the next-nearest-neighbour interaction is just slightly enhanced, which leaves remarkable footprints in the magnon spectral function and the Curie temperature.
These findings point to a variety of new questions and problems to tackle in the future: On the modeling side we found sizeable local Coulomb interactions as compared to the non-interacting band width, possibly rendering dynamical mean field theory or similar approaches necessary to capture all correlation effects [37][38][39] . Our extended (d+p)-model is an optimal starting ground for studies like these. Together with sizeable magnonphonon couplings 40 , light-matter interactions, and longrange Coulomb interaction effects 8,41 , we can expect a plethora of correlation effects in this material to be found in the future.
In light of our findings on the environmental-screeninginduced modifications to the magnetic properties of monolayer CrI 3 we expect that magnetism in multilayer CrI 3 is more involving than expected. Next to electronic inter-layer hybridization effects layer-dependent changes to Coulomb interactions seem to be important as well to gain a full understanding. Also the role of anisotropic magnetic interactions, including both symmetric and asymmetric forms, needs to be studied in further detail. Substrate effects might be especially relevant in the context of Dzyaloshinskii-Moriya interactions (DMI), as they could be considerably enhanced by breaking/lowering inversion symmetries. At the same time, DMI appear to be promising in for the stabilization of topological magnons 42 and skyrmionic phases 43,44 in chromium trihalides. Coulomb engineering of DMI could be considered using a computational scheme similar to the one proposed in Ref. 45.
Finally, our findings render CrI 3 monolayer based heterostructures with spatially structured environments a possibly fascinating fundamentally new playground to non-invasively structure the magnetic properties of layered magnetic materials similar to what has been discussed for correlation effects in layered semiconductors 20,21,23,46,47 . Together with the recent discovery of other 2D ferromagnets and antiferromagnets we thus expect that magnetic van der Waals heterostructures are most promising platforms to engineer and design next-generation magnetic and opto-magnetic devices. The encapsulation-mediated tunability of the magnetic exchange has important implications in the future application of 2D ferromagnets as spintronic devices. The possibility to combine strong and weak FM regions on 2D ferromagnets using different substrates may find application as memory storage devices. Also, real-space manipulation of the magnon dispersion can open new possibilities for low-consumption magnon-based devices.

A. Hartree-Fock Solver and Double Counting Corrections
We approximately solve the Hamiltonian from Eq. (1) utilizing a variational single-particle wave function which allows for the breaking of the spin-symmetry. The variational energy is computed by decoupling the interaction terms in the conventional way. The resulting intraatomic HF Hamiltonian takes the form of an effective single-particle one with local occupations n i which need to be self-consistently evaluated. It can be divided into spin-conserving (h ↑↑ , h ↓↓ ) and spin-mixing terms (h ↑↓ , h ↓↑ ). In a general form, the full effective single-particle Hamiltonian reads as: where the i, j are orbital indices. The explicit expressions for each of the terms in the Hamiltonian are given below: where ρ is the self-consistent density matrix containing the local occupations: ρ σσ ii =n σ i and ρ σσ ij = d † i,σ d j,σ . In order to minimize double-counting errors in the (d+p)-model due to interaction effects in the hopping matrix elements from the ab initio calculations, we subtract a double-counting potential based on the fully localized limit (FLL) 25 approximation: are the mean Coulomb and Hund exchange interactions obtained from the cRPA tensors, N imp is the total occupancy of the Cr d-orbitals, and N σ imp is the occupancy per spin (N σ imp = N imp /2 in the paramagnetic ground state).

B. Magnetic Force Theorem
Magnetism in CrI 3 and related compounds results from a detailed interplay between local and non-local kinetic and local Coulomb interactions terms yielding an effective magnetic exchange between neighbouring Cr atoms. Generally speaking it can be understood as a super-exchange mechanism mediated by the ligand atoms which follows approximately the Goodenough-Kanamori mechanism 13 . With an approximate 90 • angle between neighbouring Cr-I-Cr atoms we thus expect a FM coupling. On a fully microscopic level, Kashin et al. 15 have recently shown that the total FM exchange interaction results from an interplay between an AFM coupling channel between Cr t 2g orbitals and a FM channel between t 2g and e g orbitals. In both, LSDA and LSDA+U calculations this interplay is dominated by the FM t 2g -e g channel so that the total exchange interaction is also FM. To understand how this microscopic picture is affected by different choices of the target space, the different orbital channels of the Coulomb tensor, and by environmental screening effects, we follow Kashin et al. 15 and analyze the results of our HF calculation by means of the MFT 48 , which allows us to calculate the orbitally resolved exchange interaction matrix elements via the second variations of the total energy with respect to infinitesimal rotations of the magnetic moments, leading to the expression 49,50 (ε)e ik·Rij is the real-space Green's function, whose k-space matrix representation reads where I is the unity matrix, η → 0 + is a numerical smearing parameter, R ij is the translation vector connecting atoms i and j, and H σ k is the reciprocal-space Hamiltonian for spin σ =↑, ↓, whose matrix elements are obtained in the basis of Wannier functions from our HF calculations. Positive and negative J ij correspond here to FM and AFM couplings, respectively.

C. Magnetic Properties
The spin Hamiltonian describing the exchange interactions between Cr atoms in CrI 3 can be written as: where the first term A describes single-ion anisotropy, J ij is the isotropic Heisenberg exchange, and λ ij is the anisotropic symmetric exchange.
To calculate the spin-wave spectrum, we transform Eq. (11) into a bosonic Hamiltonian using the linearized Holstein-Primakoff transformation with A = 0 and λ = 0.09 meV (see Ref. 13): where J and J correspond to nearest and next-nearest neighbor isotropic exchange couplings, and λ is the nearest neighbor anisotropic symmetric exchange. In k-space, the Hamiltonian Eq. (12) for a honeycomb lattice takes the form where 0 = 3J + 6J + 3λ, and f 1 (q) = R e −iq·R and f 2 (q) = R e −iq·R are form factors with R and R running over cells with nearest and next-nearest neighbor atoms, respectively. The diagonalization of this Hamiltonian yields the magnon spectrum, which reads (14) To calculate the magnetic T C , we use the Tyablikov's decoupling approximation (also known as RPA) 51 through the expression

D. Substrate Screening
Next to the CrI 3 intrinsic properties we aim to also understand the role of external screening properties such as resulting for substrate materials or capping dielectrics. To this end, we utilize our WFCE approach 19 , which realistically modifies the CrI 3 Coulomb interaction tensor according to dielectric environmental screening. In this way we will be able to understand how the electronic band structure and the microscopic magnetic properties are affected, e.g., by encapsulating CrI 3 with hBN or under the influence of bulk substrates.
We start with the non-local bare Coulomb interaction of the CrI 3 monolayer as obtained from our cRPA calculations in momentum space. Within a matrix representation v αβ (q) using a product basis α, β = {n, m} we can easily diagonalize the Coulomb tensor with v ν (q) and |v ν (q) being the corresponding eigenvalues and eigenvectors of the Coulomb tensors. Assuming that the eigenbasis does not drastically change upon the effects of the cRPA screening, we can thus represent the full cRPA Coulomb tensor as where ε ν (q) are the corresponding pseduo-eigenvalues of the dielectric tensor describing the different screening channels. In Fig. 6 we show the first three ε ν (q) as a function of momentum q and find that only one shows a significant dispersion, which we refer to as the "leading" eigenvalue in the following. This behaviour becomes clear upon investigating the corresponding eigenvectors in the basis of the two Cr atoms in the long-wavelength limit, i.e. for q → 0: The leading eigenvalue v 1 (q) thus renders Coulomb penalties for mono-pole-like perturbations (all orbitally resolved electronic densities on both Cr atoms are in phase), while the second eigenvalue v 2 (q) corresponds to Cr-Cr dipole-like Coulomb penalties. ε 1 (q) and ε 2 (q) thus correspondingly describe mono-and Cr-dipole-like screening. While the dipole-like screening from the environment is negligible, the mono-pole-like screening as rendered by ε 1 (q) is strongly affected. This classical electrostatic screening can be modeled by solving the Poisson equation for a dielectric slab of height h embedded in some different dielectric environment 19,52-54 yielding withε For ε above sub = ε below sub = 1 this describes the leading dielectric function of a free-standing monolayer, which we can fit perfectly to the cRPA data, as shown in Fig. 6 and yielding h ≈ 5.2Å and ε (0) 1 ≈ 8.7 (for the case of the (d+p)-model). With these fitting parameters we can now modify the full cRPA Coulomb tensor to describe environmental screening rendered by finite ε above sub and ε below sub . Due to the monopole-like character of this environmental screening we only affect density-density Coulomb matrix elements in the very same way. Therefore Coulomb exchange elements, such as J H are not affected by the substrate screening. In the Supplemental Material we benchmark this approach to full cRPA calculations taking the screening from additional (strained) hBN layers into account. FIG. 6. Leading Screening Channels in the Wannier Function Continuum Electrostatics Approach. Red, green, and blue dots correspond to ε1(q), ε2(q), and ε3(q), respectively. ε1(q) is fitted using the expression from Eq. (19). The shown data is for the (d+p)-model.

ACKNOWLEDGMENTS
We thank M. Kaltak for sharing his cRPA implementation with us. MIK and ANR acknowledge support by European Research Council via Synergy Grant 854843 -FASTCORR. D.S. thanks financial support from EU through the MSCA project Nr. 796795 SOT-2DvdW. Part of this work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative.

Appendix A: WFCE benchmark
To benchmark the validity of the WFCE approach for CrI 3 we calculated the local cRPA Coulomb tensors for freestanding CrI 3 and hBN-embedded CrI 3 . For the latter we used a heterostructure consisting of a CrI 3 monolayer between two hBN monolayers. To keep the number of involved atoms minimal we strained the hBN lattice constant to 2.32Å (compared to 2.51Å) allowing us to use a primitive unitcell of CrI 3 embedded in 2 × 2 hBN monolayer supercells with 44 atoms in total. The distance between the outer I atoms and the hBN layers was set to 4Å. For the cRPA calculations we used 8 × 8 × 1 k and q meshes, 256 (384) Kohn-Sham states for the freestanding (embedded) monolayer, and utilized the donly model. The resulting leading screening channels are shown in Fig. 7 and the corresponding local Coulomb matrices are given in Tab. V. From Fig. 7 we see that indeed just the "leading" screening channel is affected by the dielectric environment formed by the hBN layers. From the Coulomb matrix elements in Tab. V we additionally see that indeed just the density-density elements are modified by the dielectric environment, while  Red, green, and blue dots /crosses correspond to ε1(q), ε2(q), and ε3(q), respectively. The shown data is for the d-only model using an interlayer-distance of hvac = 20Å and 8×8×1 k/q meshes.