Microscopic origin of ferromagnetism in trihalides CrCl$_3$ and CrI$_3$

TL;DR

This study investigates the microscopic origins of ferromagnetism in CrCl3 and CrI3 using density functional theory and modeling, highlighting the crucial role of ligand states and comparing multiple computational approaches.

Contribution

It provides a comparative analysis of exchange interactions in Cr trihalides using various first-principles methods, emphasizing the importance of ligand states in ferromagnetism.

Findings

Ligand states significantly influence ferromagnetic coupling.

Different computational methods yield varying estimates of exchange interactions.

The constraint method best captures the role of ligand relaxation in magnetic interactions.

Abstract

Microscopic origin of the ferromagnetic (FM) exchange coupling in CrCl$_3$ and CrI$_3$, their common aspects and differences, are investigated on the basis of density functional theory combined with realistic modeling approach for the analysis of interatomic exchange interactions. We perform a comparative study based on the pseudopotential and linear muffin-tin orbital methods by treating the effects of electron exchange and correlation in GGA and LSDA, respectively. The results of ordinary band structure calculations are used in order to construct the minimal tight-binding type models describing the behavior of the magnetic Cr $3d$ and ligand $p$ bands in the basis of localized Wannier functions, and evaluate the effective exchange coupling ($J_{\rm eff}$) between two Cr sublattices employing four different technique: (i) Second-order Green's function perturbation theory for infinitesimal spin rotations of the LSDA (GGA) potential at the Cr sites; (ii) Enforcement of the magnetic force theorem in order to treat both Cr and ligand spins on a localized footing; (iii) Constrained total-energy calculations with an external field, treated in the framework of self-consistent linear response theory. We argue that the ligand states play crucial role in the ferromagnetism of Cr trihalides, though their contribution to $J_{\rm eff}$ strongly depends on additional assumptions, which are traced back to fundamentals of adiabatic spin dynamics. Particularly, by neglecting ligand spins in the Green's function method, $J_{\rm eff}$ can easily become antiferromagnetic, while by treating them as localized, one can severely overestimate the FM coupling. The best considered approach is based on the constraint method, where the ligand states are allowed to relax in response to each instantaneous reorientation of the Cr spins, controlled by the external field.

Tables

Table 1. Table 1: Nearest-neighbor exchange interaction ( J 1 subscript 𝐽 1 J_{1} , in meV) in the d ​ p 𝑑 𝑝 dp model for CrCl 3 and CrI 3 as obtained in the Green’s function method for infinitesimal rotations of the Cr 3 ​ d 3 𝑑 3d spins (results of the Cr 3 ​ d 3 𝑑 3d model for CrCl 3 are shown for comparison in parentheses). The model was constructed in the LMTO method and using maximally localized Wannier functions technique (MLWF) in the QE method.

	CrCl3	CrI3
MLWF	$-1.63$ ($-2.75$)	$2.37$
LMTO	$4.07$ ($2.47$)	$2.98$

Table 2. Table 2: Parameters of the exchange interaction (in meV) as obtained using the QE (MLWF) method, supplemented with the Green’s function method for the infinitesimal rotations of the Cr 3 ​ d 3 𝑑 3d spins. Notations of exchange interactions are shown in Fig.  6 .

	$J_1$	$J_2$	$J_3$	$J_4$	$J_5$	$J_6$
CrCl3	$-1.63$	$-0.10$	$0.21$	$0.09$	$0.07$	$-0.28$
CrI3	$2.37$	$0.25$	$0.84$	$0.46$	$0.68$	$-0.26$

Table 3. Table 3: Parameters of the effective exchange coupling (in meV) as obtained using the QE (MLWF) method and Eq. ( 6 ) for the energies of infinitesimal rotations of spins, which enforces the magnetic force theorem. J eff A superscript subscript 𝐽 eff A J_{\rm eff}^{\rm A} and J eff F superscript subscript 𝐽 eff F J_{\rm eff}^{\rm F} are the values obtained for the antiferromagnetic and ferromagnetic rotations of spins in two Cr sublattices, as shown in Fig.  8 a and b, respectively, and J eff = J eff A − J eff F subscript 𝐽 eff superscript subscript 𝐽 eff A superscript subscript 𝐽 eff F J_{\rm eff}=J_{\rm eff}^{\rm A}-J_{\rm eff}^{\rm F} is the total coupling. The results of LMTO calculations are shown in parentheses for comparison.

	$J_{\mathrm{eff}}^A$	$J_{\mathrm{eff}}^F$	$J_{\mathrm{eff}}$
CrCl3	$207.74$ ($286.10$)	$95.98$ ($145.11$)	$111.76$ ($140.99$)
CrI3	$193.23$ ($203.14$)	$124.63$ ($144.78$)	$68.60$ ($58.36$)

Table 4. Table 4: Local magnetic moments (in μ B subscript 𝜇 B \mu_{\rm B} ) at the Cr and ligand sites as obtained for the d ​ p 𝑑 𝑝 dp model in the QE (MLWF) method. The results of LMTO calculations are shown in parenthesis for comparison.

	Cr	ligand
CrCl3	$2.951$ ($3.145$)	$0.016$ ($-0.048$)
CrI3	$3.105$ ($3.370$)	$-0.035$ ($-0.123$)

Table 5. Table 5: Parameters of effective exchange coupling (in meV) as obtained using QE (MLWF) method and Eq. ( 11 ) for the constrained energies obtained as self-consistent response to the external field applied at the Cr sites.

	$J_{\mathrm{eff}}$
CrCl3	$27.71$
CrI3	$41.55$

Table 6. Table 6: Parameters of exchange interaction (in meV) obtained from solution of model ( 12 ) in the mean-field Hartree-Fock approximation. Notations of exchange interactions are shown in Fig.  6 .

	$J_1$	$J_2$	$J_3$	$J_4$	$J_5$	$J_6$
CrCl3	$-0.97$	$-0.14$	$-0.26$	$-0.02$	$-0.03$	$-0.16$
CrI3	$1.12$	$-0.27$	$-0.36$	$0.02$	$0.07$	$-0.45$

Equations

$$E_{\mathrm{xc}}^{\mathrm{LSDA}}[n,\mathbf{m}]=\int d\mathbf{r}n(\mathbf{r}){\epsilon }_{\mathrm{xc}}[n(\mathbf{r}),|\mathbf{m}(\mathbf{r})|]$$

$${\mathcal{H}}_S=-\frac{1}{2}\sum _{ij}{J}_j{\mathbf{e}}_i\cdot {\mathbf{e}}_{i+j},$$

$$J_i=\frac{1}{2\pi }\mathrm{Im}{\int }_{-\infty }^{{\epsilon }_{\mathrm{F}}}d\epsilon {\mathrm{Tr}}_L\left\{{\hat{b}}_0^z{\hat{G}}_{0i}^{\uparrow }(\epsilon ){\hat{b}}_i^z{\hat{G}}_{i0}^{\downarrow }(\epsilon )\right\},$$

$$J_1\approx \left(\frac{1}{{\Delta }_{te}^{\uparrow \uparrow }}-\frac{1}{{\Delta }_{te}^{\uparrow \downarrow }}\right)t_{te}^2-\frac{1}{{\Delta }_{tt}^{\uparrow \downarrow }}{t}_{tt}^2,$$

$$J_1^{\mathrm{DE}}=\frac{1}{2\pi }\mathrm{Im}{\int }_{-\infty }^{{\epsilon }_{\mathrm{F}}}d\epsilon {\mathrm{Tr}}_L\left\{{\hat{G}}_{01}^{\uparrow }(\epsilon ){\hat{H}}_{10}^{\downarrow }\right\}.$$

$$E(\vec{\mathbf{m}})=E_{\mathrm{sp}}(\vec{\mathbf{b}}+\vec{\mathbf{h}})-\frac{1}{2}{\mathrm{Tr}}_L\{\vec{\mathbf{h}}\cdot \vec{\mathbf{m}}\},$$

$$J_{\mathrm{eff}}^{\mathrm{A}}=\frac{1}{3}\frac{{\partial }^{2}E(\theta ,-\theta )}{\partial (2\theta {)}^2},$$

$$E_{\mathrm{xc}}=-\frac{1}{4}\sum _i{\mathrm{Tr}}_L\left\{{\hat{\mathbf{m}}}_i{\hat{I}}_i{\hat{\mathbf{m}}}_i\right\},$$

$${\hat{\mathbf{b}}}_i=-\frac{1}{2}\left({\hat{\mathbf{m}}}_i{\hat{I}}_i+{\hat{I}}_i{\hat{\mathbf{m}}}_i\right).$$

$$\delta {\vec{b}}^{\alpha }=\left({\left[1-\mathbb{I}{\mathbb{R}}^{\alpha }\right]}^{-1}-1\right){\vec{h}}^{\alpha },$$

$$\delta E=-\frac{1}{4}{\mathrm{Tr}}_L\{\vec{\mathbf{h}}\cdot \delta \vec{\mathbf{m}}\}.$$

$$\hat{\mathcal{H}}=\sum _{ij}\sum _{\sigma }\sum _{ab}{H}_{ij}^{ab}{\hat{c}}_{ia\sigma }^{†}{\hat{c}}_{jb\sigma }^{}+\frac{1}{2}\sum _i\sum _{\sigma {\sigma }^{\prime }}\sum _{abcd}{U}_{abcd}{\hat{c}}_{ia\sigma }^{†}{\hat{c}}_{ic{\sigma }^{\prime }}^{†}{\hat{c}}_{ib\sigma }^{}{\hat{c}}_{id{\sigma }^{\prime }}^{},$$

$$\hat{b}=-\frac{1}{2}\left(\hat{m}\hat{I}+\hat{I}\hat{m}\right),$$

$$\hat{\tilde{b}}=-\frac{1}{2}\left(\hat{\tilde{m}}\hat{\tilde{I}}+\hat{\tilde{I}}\hat{\tilde{m}}\right),$$

$${\tilde{I}}^{ab}=-2{\tilde{b}}^{ab}/({\tilde{m}}^a+{\tilde{m}}^b).$$

Full text

I Introduction

Cr trihalides, Cr$X_{3}$ (where $X=$ Cl, Br, or I), and other van der Waals compounds with the rhombohedral structure have attracted recently much attention due to discovery of two-dimensional (2D) ferromagnetism CrGeTe3_Nature ; CrI3_Nature , which has tremendous importance for the development of ultra-compact spintronics. Furthermore, there is a general fundamental interest in microscopic origin of this phenomenon. Basically, there are two main question: (i) the existence of ferromagnetic (or any other long-range magnetic) order at finite temperatures in the 2D case, which formally contradicts to Mermin-Wagner theorem MerminWagner ; (ii) the origin of ferromagnetic (FM) interactions themselves. The first restriction is known to be removed by anisotropic interactions BrunoPRL2001 , as was recently confirmed for van der Waals magnets CrGeTe3_Nature ; CrI3_Lui ; CrI3_Lado . Regarding the second issue, the ferromagnetism of Cr trihalides is basically a bulk property. For instance, bulk samples of CrBr3 and CrI3 are known to be ferromagnetic with the Curie temperature of $33$ K and $68$ K, respectively Tsubikawa1960 ; Bene1969 , while CrCl3 is a layered antiferromagnet (i.e., also consisting of FM layers) with the Néel temperature of $17$ K Bene1969 .

Then, why are Cr trihalides ferromagnetic? An easy answer would be the following: Cr is trivalent, located in the octahedral environment, and the Cr-$X$-Cr angle is close to $90^{\circ}$. Therefore, the Cr-Cr coupling is expected to be ferromagnetic according to Goodenough-Kanamori-Anderson (GKA) rules Kanamori_GKA , mainly due to intraatomic exchange interaction responsible for Hund’s first rule at the ligand ($X$) sites. However, the GKA rule for the $90^{\circ}$-degree exchange is not very conclusive, as it relies on specific paths and processes for the exchange coupling, while in reality there can be several competing mechanisms supporting either ferromagnetism or antiferromagnetism, so that the final answer strongly depends on the ratio of relevant physical parameters Chaloupka . In fact, Kanamori himself admitted that there are exceptions from this rule in the $90^{\circ}$-degree case Kanamori_GKA , and in our work we will show that the exchange coupling in Cr trihalides can easily become antiferromagnetic, depending on the considered model.

Another important issue is how these GKA rules work in practical first-principles calculations based on density functional theory (DFT), and how to map properly this or any other theory of interacting electrons onto effective spin model and not to lose the contributions that can be responsible for the GKA rules JHeisenberg ; Stepanov ; review2008 ; CrO2PRB2015 . The problem is indeed very complex and there are many tricky issues in it. Particularly, the mapping onto the spin model frequently assumes the local character of exchange-correlation (xc) interactions, following a similar concept of local spin density approximation (LSDA), which often supplements DFT calculations and is also “local”. However, the main FM mechanism of the $90^{\circ}$ exchange is basically non-local: the magnetic $3d$ states spread from two Cr sites to the neighboring ligand site and there interact via Hund’s rule coupling Kanamori_GKA , which has the same origin as the canonical direct exchange interaction proposed by Heisenberg almost a centaury ago Heisenberg . In this respect, it is interesting to note that CrBr3 was indeed regarded as a rare example of materials (together with CrO2) where ferromagnetism arises from the direct exchange interaction Yosida .

There is quite a common tendency in the electronic structure calculations: if (typically due to additional approximations) such calculations underestimate the FM coupling, the missing effect is ascribed to the direct exchange CrO2PRB2015 ; Ku ; Mazurenko2007 . Nevertheless, this step looks somewhat irrational, because the effect of direct exchange interactions is already included in LSDA, which is used as a starting point for many electronic structure calculations, even despite it is formally “local”. Furthermore, LSDA also includes the important effects of screening for the direct exchange interactions due to electron correlations. Indeed, the xc energy in LSDA is given by

[TABLE]

($n$ and $\boldsymbol{m}$ being total electron and spin magnetization density, respectively). Thus, at each point, the xc energy is given by $n$ and $\boldsymbol{m}$ at the same point, which is the definition of locality of the LSDA functional. However, $n$ and $\boldsymbol{m}$ can be equivalently expanded in terms of a complete set of localized Wannier functions for the occupied states WannierRevModPhys . In this representation, $\boldsymbol{m}(\boldsymbol{r})$ is the sum of Wannier function contributions coming from different magnetic sites. Therefore, although $E_{\rm xc}^{\rm LSDA}$ is local with respect to the total magnetization density, it includes non-local effects caused by overlap of individual contributions to $\boldsymbol{m}$ centered at different magnetic sites in the Wannier representation.

In this work we will address this problem in details by considering two characteristic examples of Cr trihalides: CrCl3 and CrI3. Particularly, how should one calculate the interatomic exchange interactions within DFT in order to take into account all relevant processes, responsible for the FM coupling in the $90^{\circ}$ case? To this end, we will argue that it may not be enough to consider the reorientation of only transition-metal spins (for instance, within the commonly used Green’s function perturbative approach for infinitesimal spin rotations JHeisenberg ), as there is a large contribution to the exchange energy associated with the ligand sites Oguchi , which should be also taken into account. Nevertheless, it does not automatically guarantee the “right” solution of the problem because the answer strongly depends on how the ligand spins are treated, tracing back to fundamentals of adiabatic spin dynamics. For instance, by enforcing the magnetic force theorem in order to treat the Cr and ligand spins on the same localized footing, one can severely overestimate the FM coupling. It turns out that the best considered approach is to fully relax the ligand states for each instantaneous configuration of the Cr spins: it correctly describes the FM coupling, both in CrCl3 and CrI3, and yields the effective exchange coupling constant comparable with results of brute force total energy calculations. We will show how this problem can be efficiently solved in the framework of self-consistent linear response theory SCLR .

The rest of the paper is organized as follows. In Sec. II we will describe our method based on construction and analysis of the minimal tight-binding type models for the Cr $3d$ and ligand $p$ bands of CrCl3 and CrI3 in the basis of localized Wannier functions. All our calculations are based on either LSDA or generalized gradient approximation (GGA) as implemented in, respectively, the linear muffin-tin orbital (LMTO) and pseudopotential Quantum-ESPRESSO (QE) methods. The on-site Coulomb repulsion $U$ should not play a decisive role in the origin of ferromagnetism in Cr trihalides, as argued in Appendix A. Then, in Sec. III we will consider different approaches for calculations of effective intersublattice exchange coupling: the Green’s function technique (Sec. III.1), magnetic force theorem (Sec. III.2), and constrained calculations with an external magnetic field (Sec. III.3). Finally, in Sec. IV, we will give a brief summary on our work.

II Method

We use the $R\overline{3}$ structure for both CrCl3 and CrI3 (shown in Fig. 1) with the experimental parameters reported in Refs. CrCl3str and CrI3str , respectively. As for the electronic structure calculations, we perform a comparative study based on the LMTO method in the atomic-spheres approximation (ASA) LMTO and the pseudopotential QE method realized in the plane-wave basis qe . Furthermore, all LMTO calculations have been performed in LSDA, while in QE we employed the generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof exchange-correlation (xc) functional PBE . We use the mesh of $8\times 8\times 8$ ($10\times 10\times 10$) ${\bf k}$-points in the Brillouin zone in the case of QE (LMTO) and the energy cutoff of $110$ Ry in the QE calculations.

Corresponding electronic structure, obtained for the FM state, is shown in Fig. 2.

Owing to the octahedral environment of the Cr sites, the magnetic $3d$ states are split into the low-lying $t_{2g}$ ($t$) and high-lying $e_{g}$ ($e$) groups. In CrCl3, the Cr $3d$ bands are well separated from the Cl $3p$ one, so that even for the majority ($\uparrow$) spin channel there is a finite gap between these two groups of bands. Therefore, in this case one can consider two types of models: the $dp$ one, which explicitly includes the $3d$ states of Cr as well as the $p$ states of the ligand atoms, and the pure $d$ one, consisting only of the Cr $3d$ states. Nevertheless, in CrI3, due to the additional upward shift of the I $5p$ states, the $\uparrow$-spin Cr $t_{2g}$ band merges into the I $5p$ one. Therefore, in this case we are able to consider only the more general $dp$ model. The electronic structure of CrCl3 and CrI3 obtained in the QE method is in excellent agreement with the results of previous GGA calculations Wang2011 . The differences between LMTO and QE results will be discussed below.

Next, we construct the Wannier functions spanning the subspace of the Cr $3d$ and ligand $p$ bands and define the model $dp$ Hamiltonian, $\hat{H}^{\uparrow,\downarrow}=[H_{ij}^{ab}]^{\uparrow,\downarrow}$, separately for the majority ($\uparrow$) and minority ($\downarrow$) spin states, as the matrix elements of the original LSDA (GGA) Hamiltonian in the basis of Wannier functions, which are numbered by indices $a$ and $b$ at sites $i$ and $j$. In order to construct the Wannier functions, we use the projector operator method and the maximally localized Wannier functions (MLWF) technique (implemented in the Wannier90 package) in the case of, respectively, LMTO and QE review2008 ; WannierRevModPhys ; wannier90 . For the trial orbitals in LMTO we use basis functions for the Cr $3d$ and ligand $p$ states. The completeness of the Wannier basis guarantees that the obtained band structure in the region of the Cr $3d$ and ligand $p$ bands fully coincides with the original LSDA (GGA) one. After that, the $d$ model for CrCl3 was constructed by starting from the $dp$ model and eliminating the Cr $3p$ band by means of the projector operator technique, both in LMTO and QE methods.

The atomic level splitting, which is obtained by diagonalizing the site-diagonal part of $\hat{H}^{\uparrow,\downarrow}$, is shown in Fig. 3.

The $\uparrow$-$\downarrow$ spin splitting of the Cr $3d$ states, driven by the xc interactions in LSDA (GGA), is larger than the $t_{2g}$-$e_{g}$ crystal-field splitting within each projection of spin. In the ASA-LMTO method, the $\uparrow$-$\downarrow$ level splitting is nearly rigid and well described by the phenomenological expression $I_{\rm Cr}M_{\rm Cr}$, where $I_{\rm Cr}\approx 0.85$ eV is the Stoner parameter and $M_{\rm Cr}=3\mu_{\rm B}$ is the local spin moment of Cr3+, as expected in LSDA Gunnarsson . In the QE (MLWF) method, however, this simple picture is no longer valid and there is an additional splitting between the $t_{2g}$ and $e_{g}$ orbitals, which also depends on spin. Particularly, the $t_{2g}$-$e_{g}$ splitting for spin $\downarrow$ is substantially smaller than that for spin $\uparrow$, which reflects similar behavior of the density of states in Fig. 2. Nevertheless, this is to be expected because the QE method takes into account the asphericity of the Kohn-Sham potential lda , which is totally ignored in ASA. Furthermore, this asphericity is additionally amplified in GGA (in comparison with LSDA). The spin splitting of the ligand $p$ levels is small in comparison with the crystal field splitting within each projection of spin. The latter effect is again larger in the QE calculations.

Fig. 4 shows the behavior of averaged transfer integrals, $\bar{t}_{ij}=\sqrt{\sum_{ab}H_{ij}^{ab}H_{ji}^{ba}}$, for CrCl3.

Despite substantial difference in details of QE and LMTO calculations, the obtained parameters are very similar in two considered methods. The transfer integrals around Cr sites are basically limited to the neighboring Cr-Cl and Cr-Cr bonds. The transfer integrals around ligand sites involve more ligand-ligand bonds. Nevertheless, all transfer integrals are short-ranged and practically vanish at the distance $d\sim 4.5$ Å. Very similar tendency was found for CrI3, which is not shown here.

Parameters of the $d$ model for CrCl3 are summarized in Fig. 5.

The elimination of the Cl $3p$ bands considerably increases the crystal-field splitting between the $t_{2g}$ and $e_{g}$ orbitals, being consistent with an old conjecture by Kamanori that the $10Dq$ splitting results mainly from the hybridization of the transition-metal $3d$ states with the ligand $p$ states Kanamori . Nevertheless, the structure of the atomic Cr $3d$ levels remains quite different in LMTO and QE (MLWF). As we shall see below, this difference strongly affects the behavior of nearest-neighbor (nn) exchange interactions. As expected, the transfer integrals become more long-ranged in comparison with the ones in the $dp$ model PRB06 . However, in all other respects their behavior is very similar in LMTO and QE.

III Magnetic interactions

In this Section we explore abilities of different models and techniques for the evaluation of interatomic exchange interactions in CrCl3 and CrI3. Our ultimate goal is to elucidate the microscopic origin of the FM coupling in these compounds and we will show that the problem is indeed very nontrivial. The basic idea is to map the total energy difference associated with the reorientation of the Cr spins onto the model

[TABLE]

where $\boldsymbol{e}_{i}$ is the direction of spin at site $i$, and $J_{j}$ is the corresponding exchange coupling between (central) site $i$ and the one located at the lattice point $i+j$. Unless specified otherwise, we will concentrate on the behavior of the effective coupling between two Cr sublattices $I$ and $II$ (see Fig. 6): $J_{\rm eff}=\sum_{j\in II}J_{j}$ (provided that the central site belongs to the sublattice $I$). If interactions with other Cr sites, except the nn ones, are negligibly small, we will have a trivial relation $J_{\rm eff}\approx J_{1}$, which holds approximately for CrCl3 and CrI3.

The most straightforward estimate for $J_{\rm eff}$ is to use the total energy difference between the antiferromagnetic (AFM) and FM states in the QE method, $\Delta E=E_{\rm AFM}-E_{\rm FM}$, which yields $21.35$ and $34.93$ meV per one formula unit of CrCl3 and CrI3, respectively, corresponding to $J_{\rm eff}=\Delta E/3=7.12$ and $11.64$ meV. This result clearly indicates that the coupling is indeed ferromagnetic and stronger in the case of CrI3. Thus, this result can be used as a reference for other models and approximations.

III.1 Green’s function method

The first model we consider is based on the Green’s function method JHeisenberg . The main idea here is to define the xc field $\hat{\boldsymbol{b}}_{i}=(0,0,\hat{b}_{i}^{z})$, associated with an arbitrary site $i$, where $\hat{b}_{i}^{z}=\hat{H}_{ii}^{\uparrow}-\hat{H}_{ii}^{\downarrow}$, and consider the energy change caused by infinitesimal rotations of this field, \hat{\boldsymbol{b}}_{i}\rightarrow\hat{\boldsymbol{b}}_{i}=(\theta_{i},0,1$$-$$\theta_{i}^{2}/2)\hat{b}_{i}^{z} ($\theta_{i}$ being the azimuthal angle), as a perturbation. Then, to 2nd order in $\theta$, the corresponding energy change will be given by Eq. (2) with the parameters

[TABLE]

where $\hat{G}^{\uparrow,\downarrow}_{0i}(\varepsilon)=[\varepsilon-\hat{H}^{\uparrow,\downarrow}]^{-1}_{0i}$ is the one-electron Green function between sites [math] and $i$, $\varepsilon_{\rm F}$ is the Fermi energy, and ${\rm Tr}_{L}$ denotes the trace over orbital indices JHeisenberg . Importantly, Eq. (3) describes the change of single-particle (Kohn-Sham) energies for the occupied states, while other contributions to the total energy are supposed to cancel out in the 2nd order of $\theta$ as a result of the magnetic force theorem JHeisenberg . This technique is widely used in the LMTO method where, owing to the pseudoatomic basis, one can easily define local magnetic moments at each site. The advantage of the Wannier basis is that it allows us to transfer this technique to other methods, including those originally formulated in the basis of plane waves, as in the QE method.

In principle, Eq. (3) can also be used for the calculation of the exchange coupling between Cr and ligand sites. Nevertheless, we begin here with the standard procedure and employ this equation only for the Cr-Cr interactions, while assuming that the main role of the intermediate ligand states is to assist these interactions by connecting two Cr sites. The results of these calculations for $J_{1}$ are summarized in Table 1.

The LMTO and QE methods provide quite consistent description for $J_{1}$ in the case of CrI3, where the FM character of this interaction is consistent with experimental data. Nevertheless, the situation is totally different in CrCl3, where even the sign of $J_{1}$ is different in LMTO and QE, thus putting the latter data in clear disagreement with the experiment. This is the first puzzle we need to solve and below we will argue that such behavior is related to the electronic structure of CrCl3 and CrI3 in LMTO and QE.

Since the $dp$ and $d$ models provide qualitatively similar results and encounter the same problem for CrCl3 ($J_{1}>0$ in LMTO, while $J_{1}<0$ in QE), it is convenient to start our analysis with the simplest $d$ model. An analytical expression for $J_{1}$ can be found by employing superexchange theory PWA ; KugelKhomskii , which relies on additional approximations but provides some insight into behavior of $J_{1}$. It yields

[TABLE]

where $t^{2}_{te}=\sum_{a\in t,b\in e}\left(H_{01}^{ab}\right)^{2}$, $t^{2}_{tt}=\sum_{a,b\in t}\left(H_{01}^{ab}\right)^{2}$, and $\Delta$ is the intraatomic energy splitting between the $t_{2g}$ and either $e_{g}$ ($te$) or $t_{2g}$ ($tt$) levels with the same ($\uparrow\uparrow$) or opposite ($\uparrow\downarrow$) directions of spin. Then, it is clear that smaller $\Delta^{\uparrow\uparrow}_{te}$ and larger $\Delta^{\uparrow\downarrow}_{te}$ in LMTO (see Fig. 5a) favor the FM coupling. On the other hand, smaller $\Delta^{\uparrow\downarrow}_{tt}$, also in LMTO, will tend to stabilize the AFM coupling. However, since $t^{2}_{tt}<t^{2}_{te}$ ($t^{2}_{te}=0.044$ eV2, while $t^{2}_{tt}=0.028$ eV2, according to MLWF calculations), the last term plays a less important role. These tendencies clearly explain why $J_{1}>0$ in LMTO, while $J_{1}<0$ in QE. One of the lessons learned from this analysis is that the FM coupling in the $d$ model can be stabilized only if $\Delta^{\uparrow\uparrow}_{te}$ is relatively small. This is the main reason why the ferromagnetism totally disappears, even in LMTO, if one tries to correct this superexchange picture by adding the on-site Coulomb repulsion, which additionally increases $\Delta^{\uparrow\uparrow}_{te}$, as discussed in Appendix A. Similar conclusion can be arrived by considering the so-called “double exchange” (DE) contribution to $J_{1}$, obtained in the 2nd order expansion for $\hat{G}_{i0}^{\downarrow}(\varepsilon)$ with respect to $(\hat{b}_{0}^{z})^{-1}=(\hat{b}_{i}^{z})^{-1}$ in Eq. (3) PRL99 ; CrO2PRB2015 , which yields

[TABLE]

In comparison with Eq. (4), this expression takes into account the contributions proportional to $1/\Delta^{\uparrow\uparrow}$ and neglects the contributions proportional to $1/\Delta^{\uparrow\downarrow}$, assuming that the latter are small. As expected PRL99 , $J_{1}^{\rm DE}$ is ferromagnetic and can be estimated in LMTO and QE as $18.85$ and $16.23$ meV, respectively. Thus, we clearly observe again the excess of ferromagnetism in the case of LMTO, which is closely related to the structure of intraatomic level splitting.

Then, why does not this problem occur in CrI3? Regarding the intraatomic level splitting, the situation is pretty much similar in CrCl3 and CrI3 (see Fig. 3): in both cases, the LMTO and QE methods provide two different schemes of atomic level splitting (but with clear similarity between CrCl3 and CrI3, if one compares separately the results of either LMTO or QE calculations). Nevertheless, the LMTO and QE methods provide rather consistent description for $J_{1}$ in CrI3 (see Table 1). The reason is again related to the electronic structure of CrI3, which can be classified as the charge-transfer insulator ZSA : the $\uparrow$-spin I $5d$ and Cr $3d$ bands strongly overlap so that these two groups of states become entangled with each other (see Fig. 2). In such situation, the distribution of the Cr $3d$ states over the energy, which mainly controls the value of $J_{1}$ in Eq. (3), is determined by the $dp$ hybridization, while the atomic level splitting is less important. This naturally explains similarity of LMTO and QE data for CrI3.

In order to further elucidate the impact of the electronic structure on interatomic exchange interactions in CrCl3 and CrI3, it is instructive to consider the band-filling dependence of $J_{1}$, which is again calculated using Eq. (3), but varying the position of the Fermi level (7).

In CrCl3, $J_{1}(\varepsilon)$ practically vanishes after integration over the Cl $3p$ band, spreading from $-5$ to $-1$ eV (see Fig. 2), even though there is a substantial variation of $J_{1}(\varepsilon)$ within the band. Thus, the contribution of the Cl $3p$ band to $J_{1}$ is small, that justifies the use of the $d$ model for these purposes, as was done above. Then, in the QE method, the shape of $J_{1}(\varepsilon)$ in the Cr $t_{2g}$ band (spreading from $-1$ to [math] eV) is nearly symmetric relative to the band center so that the total contribution of this band to $J_{1}$ vanishes. This means that the FM and AFM interactions originating from virtual excitations from the Cr $t_{2g}$ band to unoccupied states cancel each other. In the LMTO method, the shape of $J_{1}(\varepsilon)$ in the Cr $t_{2g}$ region is slightly deformed, leading to some asymmetry, which is in turn responsible for the appearance of uncompensated FM interactions, as was discussed above. In CrI3, however, the I $5p$ and Cr $t_{2g}$ states form a common band. The magnetic interactions in this band obey some general principles: namely, the ferromagnetism is expected at the beginning and the end of the band filling, while the antiferromagnetism is expected in the middle of the band filling Heine1 ; Heine2 . Obviously, if the band were fully isolated, we would have $J_{1}(\varepsilon_{\rm F})=0$. Nevertheless, the interaction with unoccupied states makes in ferromagnetic, both in LMTO and QE calculations.

Parameters of the main exchange interactions, spreading up to six nearest neighbors (see Fig. 6), are reported in Table 2. One can see that the nn interaction clearly dominates, so that $J_{\rm eff}\approx J_{1}$ holds for CrCl3 and to lesser extent for CrI3. Furthermore, using the obtained parameters $J_{i}$ and employing Tyablikov’s random-phase approximation tyab , the Curie temperature for CrI3 can be estimates as $83$ K, which is in fair agreement with the experimental value of $68$ K Bene1969 .

Summarizing this part, the Green’s function based expression, Eq. (3), correctly describes the FM coupling in CrI3, but not in CrCl3. The FM character of exchange interactions obtained for CrCl3 in the ASA-LMTO method is probably fortuitous as it does not properly treat the asphericity of the Kohn-Sham potential, as was demonstrated by more rigorous full-potential QE calculations. Thus, we are still missing some important mechanism that stabilizes the FM ground state in CrCl3 (and probably plays some role also in CrI3). In the remaining part of this section we will go beyond Eq. (3) for the analysis of interatomic exchange interactions and argue that missing FM contributions are related to the ligand $p$ states, though the answer strongly depends on how these contributions are treated.

III.2 Enforcement of the magnetic force theorem

First, we note that the expression (3) suffers from a systematic error: although it follows from the magnetic force theorem, which allows us to replace the total energy change by the change of single-particle energies, it formally violates this theorem. The magnetic force theorem can be proven if a certain rotation of the xc field $\hat{\boldsymbol{b}}_{i}\rightarrow\tensor{R}_{\theta_{i}}\hat{\boldsymbol{b}}_{i}$ rotates the spin magnetization by the same angle: $\hat{\boldsymbol{m}}_{i}\rightarrow\tensor{R}_{\theta_{i}}\hat{\boldsymbol{m}}_{i}$. Note that in these notations, $\hat{\boldsymbol{b}}_{i}$ is the vector $\hat{\boldsymbol{b}}_{i}=(\hat{b}_{i}^{x},\hat{b}_{i}^{y},\hat{b}_{i}^{z})$, where each of its components is the matrix in the Wannier basis, while $\tensor{R}_{\theta_{i}}$ is the $3\times 3$ tensor acting on $\hat{b}_{i}^{x}$, $\hat{b}_{i}^{y}$ , and $\hat{b}_{i}^{z}$. Then, since the xc energy is invariant under spin rotations, the total energy change is described by the single-particle part given by Eq. (3) PRB98 . This property holds in LSDA, unrestricted Hartree-Fock approximation, and for any xc functional satisfying the condition of gauge invariance PRB98 . However, the rotation $\hat{\boldsymbol{b}}_{i}\rightarrow\tensor{R}_{\theta_{i}}\hat{\boldsymbol{b}}_{i}$ (which is typically used in practical calculations) does not necessarily lead to the rotation $\hat{\boldsymbol{m}}_{i}\rightarrow\tensor{R}_{\theta_{i}}\hat{\boldsymbol{m}}_{i}$. More specifically, considering an iterative process inherent to solving the Kohn-Sham equations, the rotation $\hat{\boldsymbol{b}}_{i}\rightarrow\tensor{R}_{\theta_{i}}\hat{\boldsymbol{b}}_{i}$ will indeed lead to the rotation $\hat{\boldsymbol{m}}_{i}\rightarrow\tensor{R}_{\theta_{i}}\hat{\boldsymbol{m}}_{i}$ at the input of iteration. However, at the output of the same iteration, the magnetization will tend to relax to the ground state configuration and additionally rotate away from $\tensor{R}_{\theta_{i}}\hat{\boldsymbol{m}}_{i}$. This problem was noticed, for instance, by Stocks et al. Stocks and considered in details by Bruno BrunoPRL2003 , who also proposed to adjust the magnetization by applying a constraining field fixing the direction of the spin magnetization along $\tensor{R}_{\theta_{i}}\hat{\boldsymbol{m}}_{i}$ for a given rotation of the xc field, $\hat{\boldsymbol{b}}_{i}\rightarrow\tensor{R}_{\theta_{i}}\hat{\boldsymbol{b}}_{i}$. Below, we closely follow this strategy, also using for these purposes the linear response theory SCLR .

Namely, let $\vec{\boldsymbol{b}}$ be the column vector composed of $\hat{\boldsymbol{b}}_{i}$ at different sites of the unit cell (including ligands), $(\vec{\boldsymbol{b}})^{T}=(\ \dots\ ,\hat{\boldsymbol{b}}_{i},\ \dots\ )$, and $\vec{\boldsymbol{m}}$ be a similar vector for the magnetization. Then, knowing the $x$ component $(\vec{m}^{x})^{T}=(\ \dots\ ,\theta_{i}\hat{m}_{i}^{z},\ \dots\ )$ of the rotated magnetization, one can find the constraining field $\vec{\boldsymbol{h}}=(\vec{h}^{x},0,0)$, which reproduces $\vec{m}^{x}$, using the linear response theory as $\vec{h}^{x}=\left[\left(\mathcal{R}^{\uparrow\downarrow}\right)^{-1}+\left(\mathcal{R}^{\downarrow\uparrow}\right)^{-1}\right]\vec{m}^{x}-\vec{b}^{x}$, where $\mathcal{R}^{\uparrow\downarrow}$ and $\mathcal{R}^{\downarrow\uparrow}$ are the matrix elements of the response tensor, whose explicit expression can be found in Ref. SCLR . Note, that throughout this work we use the notations where the spin magnetization is related to the density matrix $\hat{n}$ as $\hat{\boldsymbol{m}}={\rm Tr}_{S}\{\hat{\boldsymbol{\sigma}}\hat{n}\}$ ($\hat{\boldsymbol{\sigma}}$ being the vector of Pauli matrices and ${\rm Tr}_{S}$ being the trace over spin variables) and the magnetic field (both constraining and xc one) enters the Kohn-Sham equations as $\frac{1}{2}\hat{\boldsymbol{h}}\cdot\hat{\boldsymbol{\sigma}}$, which holds at each site of the system.

Using $\vec{\boldsymbol{h}}$ and $\vec{\boldsymbol{b}}$, one can calculate the total energy corresponding to the “rotated” configuration of spins as is typically done in the constraint method:

[TABLE]

where $E_{\rm sp}(\vec{\boldsymbol{b}}+\vec{\boldsymbol{h}})$ is the sum of the occupied Kohn-Sham energies calculated for the field $\vec{\boldsymbol{b}}+\vec{\boldsymbol{h}}$, and $\vec{\boldsymbol{h}}\cdot\vec{\boldsymbol{m}}$ stands for the dot product of two vectors with the summation over atomic sites. For $\vec{\boldsymbol{h}}=0$, Eq. (6) is reduced to Eq. (3).

In order to obtain $J_{\rm eff}$, we first consider the magnetic configuration shown in Fig. 8a, where the spins in two Cr sublattices are rotated by $\theta$ and -$$\theta relative to the FM axis $z$. The corresponding energy (per one unit cell including two Cr sites) is denoted as $E(\theta,-\theta)$.

We call these rotations “antiferromagnetic”, referring to the alignment of transversal components of the rotated spins. Corresponding exchange coupling can be evaluated as

[TABLE]

where the prefactor $1/3$ takes into account three nn bonds in the $xy$ plane. The obtained energies are shown in Fig. 9 and the corresponding parameters are summarized in Table 3.

The constrained energy reveals perfect quadratic dependence on the angle $\theta$. Moreover, after enforcing the magnetic force theorem, Eq. (6) correctly reproduces the FM coupling for both CrI3 and CrCl3. However, the value of this coupling is too large compared to typical energies of magnetic excitations in transition-metal oxides and related compounds.

In the following, we will show that the reason is related to a very fundamental problem of how to treat the contribution of ligand sites to the exchange coupling between Cr sites. Indeed, the $dp$ hybridization induces feasible magnetization at the ligand sites (Table 4).

The difference between the QE (MLWF) and LMTO data are quite expected: the values of local magnetic moments depend on the crystal-field splitting, which is very different in QE and LMTO, as was discussed above (see Fig. 3). Nevertheless, the QE and LMTO data show the same tendency: the Cr moment increases from CrCl3 to CrI3 (due to the increase of $dp$ hybridization and overlap of the I $5p$ and Cr $3d$ bands), and this change is compensated by the ligand sites as so to keep the total moment equal to $3$ $\mu_{\rm B}$ in both FM semiconductors CrCl3 and CrI3. $J_{\rm eff}$ may of course depend on the magnetization at the ligand sites. Nevertheless, one should clearly understand how this magnetization is treated and what is the relevant physical picture behind Eq. (6). In this respect, we would like to note that, in order to use the magnetic force theorem, the constraining field should be applied at all sites of the system, including the ligand ones. Therefore, Eq. (6) implies that the directions of the ligand spins are also rigidly fixed and there is no conceptual difference between the Cr and ligand spins: both groups of spins are treated as localized in the sense that the rotation of the Cr spins has absolutely no effect on the ligand spins and vice versa. Moreover, Fig. 8a implies that the directions of the ligand spins are fixed to be the same as in the FM state. Certainly, this is an approximation and apparently very crude one.

The simplest correction to $J_{\rm eff}^{\rm A}$ can be obtained by noting that, besides the interaction between the Cr spins, it also includes interaction with the FM background of the ligand spins. The latter can be corrected by considering the “ferromagnetic” rotations, where the spins in both Cr sublattices are rotated by the same angle $\theta$ (Fig. 8b). The corresponding parameter $J_{\rm eff}^{\rm F}=\frac{1}{3}\frac{\partial^{2}E(\theta,\theta)}{\partial(2\theta)^{2}}$ describes the interaction with the FM background of the ligand spins, and the proper coupling between Cr spins can be found as $J_{\rm eff}=J_{\rm eff}^{\rm A}-J_{\rm eff}^{\rm F}$. All these parameters are listed in Table 3. One can see that the interaction $J_{\rm eff}^{\rm F}$ is indeed very strong and substantially reduces the total value of $J_{\rm eff}$. Nevertheless, the absolute values of $J_{\rm eff}\sim 100$ meV are still much too high from the viewpoint of magnetic excitation energies, meaning that we can be facing another serious problem and the magnetic excitations described by Eq. (6) may be not the lowest energy ones.

III.3 Constrained calculations with the external field

Then, what is wrong and what is missing in the corrected Green’s function theory, described by Eq. (6)? In this respect, it is important to recall that the theory of spin dynamics, underlying Eqs. (3) and (6), is based on the so-called adiabaticity concept spindynamics1 ; spindynamics2 , which states that one can distinguish two types of variables, slow magnetic and fast electronic, so that for each instantaneous magnetic configuration the electronic variables have sufficient time to relax and reach equilibrium. In practical terms, this means that for each constraining field, which controls the direction of magnetization, one should solve self-consistently the set of Kohn-Sham equations and calculate the total energies, which then can be used for the description of low-energy magnetic excitations. This property is implied in derivation of Eqs. (3) and (6), where the total energies are additionally replaced by the sum of Kohn-Sham energies, owing to magnetic force theorem. Then, it is clear that the magnetization at the Cr sites, to a good approximation, can be treated as “slow” variable. However, what is the nature of the ligand states? Should they be treated as “slow” magnetic or “fast” electronic variables? The question is rather nontrivial. On the one hand, the ligand magnetization is induced by the Cr one by means of the $dp$ hybridization. From this point of view, it can be viewed as an “electronic” variable. On the other hand, since it is a part of the total magnetization, it is always tempting to treat it on the same footing as its Cr counterpart, i.e. as the “magnetic” variable. This is precisely what was done in Eq. (6), where we rigidly constrained the direction of magnetization at all sites of the system, including the ligand ones. Furthermore, $\vec{\boldsymbol{h}}$ in Eq. (6) has a matrix form in the subspace of orbital variables, which is required in order to constrain similar matrix of the magnetization density $\vec{\boldsymbol{m}}$. This is certainly a very strong constraint imposed on the magnetic system. In reality, however, we do not need to constrain all matrix elements of $\vec{\boldsymbol{m}}$. What we need is to fix only the directions of magnetic moment, $\boldsymbol{M}_{i}={\rm Tr}_{L}\left\{\hat{\boldsymbol{m}}_{i}\right\}$, while other constraints can be released. This is also related to how we view the adiabaticity concept. Eq. (6) implies that all variables of $\vec{\boldsymbol{m}}$ are slow and magnetic, while the reality can be different: only $\boldsymbol{e}_{i}=\boldsymbol{M}_{i}/|\boldsymbol{M}_{i}|$ are slow variables, while other elements of $\vec{\boldsymbol{m}}$ are fast and have sufficient time to adjust an instantaneous change of $\boldsymbol{e}_{i}$.

Therefore, as the next step, we explore the opposite strategy by assuming that all magnetic variables are restricted by the Cr sites, while the ligand magnetization simply follows the Cr one as any other electronic variable. Namely, like in the constraint formalism, we rotate the Cr spins by applying the external field $\vec{\boldsymbol{h}}=(\vec{h}^{x},0,0)$. The rotated configuration of the Cr spins is the same as in Fig. 8a. However, now this field acts only at the Cr sites. Furthermore, each $\hat{\boldsymbol{h}}_{i}$ is constant in the space of orbital variables, $\boldsymbol{h}_{i}^{ab}=\boldsymbol{h}_{i}\delta^{ab}$ ($\delta^{ad}$ being Kronecker delta), and acts only at $\boldsymbol{M}_{i}$, while other degrees of freedom are found from the self-consistent solution of Kohn-Sham equations. In principle, this is well known numerical technique, which is widely used in electronic structure calculations Stocks . In this work, we propose an analytical solution of this problem based on the self-consistent linear response theory SCLR , which allows us to find the self-consistent xc field $\vec{\boldsymbol{b}}$ to first order in $\vec{\boldsymbol{h}}$, the corresponding magnetization $\vec{\boldsymbol{m}}$, and the total energy to second order in $\vec{\boldsymbol{h}}$, which are sufficient to evaluate $J_{\rm eff}$.

First, we approximate the xc energy in LSDA (GGA) as Gunnarsson :

[TABLE]

where $\hat{I}_{i}$ is taken in the matrix form in order to describe asphericity of the xc potential. The corresponding xc field, $\hat{\boldsymbol{b}}_{i}=2\delta E_{\rm xc}/\delta\hat{\boldsymbol{m}}_{i}$, at site $i$ is given by

[TABLE]

By applying this equation for $\hat{b}_{i}^{z}$ and $\hat{m}_{i}^{z}$, obtained in collinear calculations for the FM state, one can find the Stoner matrix $\hat{I}_{i}$ (see Appendix B). We need this step in order to include self-consistency effects for noncollinear magnetic structures, like the one shown in Fig. 8a. Namely, by introducing the rank-4 tensor ${\mathbb{I}}=\left[{\mathcal{I}}^{abcd}\right]$ with ${\mathcal{I}}^{abcd}=-\frac{1}{2}\left(I^{bc}\delta^{ad}+I^{ad}\delta^{bc}\right)$, Eq. (9) can be rewritten in the compact form $\vec{\boldsymbol{b}}={\mathbb{I}}\vec{\boldsymbol{m}}$, which implies summation over last two indices of ${\mathbb{I}}$ and both indices of $\vec{\boldsymbol{m}}$. Then, the change of the xc field caused by $\vec{\boldsymbol{h}}$ can be found as SCLR :

[TABLE]

where ${\mathbb{R}}^{\alpha}$ is related to spin-dependent elements of the response tensor as ${\mathbb{R}}^{x}={\mathbb{R}}^{y}=\frac{1}{2}\left({\mathcal{R}}^{\uparrow\downarrow}+{\mathcal{R}}^{\downarrow\uparrow}\right)$ and ${\mathbb{R}}^{z}=\frac{1}{2}\left({\mathcal{R}}^{\uparrow\uparrow}-{\mathcal{R}}^{\downarrow\downarrow}\right)$ SCLR . Then, using $\vec{\boldsymbol{h}}$ and $\delta\vec{\boldsymbol{b}}$, one can find the corresponding magnetization change as $\delta\vec{m}^{\alpha}={\mathbb{R}}^{\alpha}(\vec{h}^{\alpha}+\delta\vec{b}^{\alpha})$ and evaluate the constrained energy (i.e., including the penalty term) to 2nd order in $\vec{\boldsymbol{h}}$ as SCLR :

[TABLE]

Finally, using $\delta\vec{\boldsymbol{m}}$ one can find the canting angle $\theta$ and evaluate the effective coupling constant $J_{\rm eff}$ from Eq. (7). The obtained energies are also shown in Fig. 9, which reveals an excellent quadratic dependence $\delta E=6J_{\rm eff}\theta^{2}$. The corresponding parameters $J_{\rm eff}$ are summarized in Table 5.

One can clearly see that the new values of $J_{\rm eff}$ are substantially reduced and become comparable with the ones derived from the total energy difference between the FM and AFM states. Furthermore, $J_{\rm eff}$ is larger in CrI3, being again in total agreement with the trend obtained in the total energy calculations. We would also like to emphasize that this $J_{\rm eff}$ characterizes the local stability of the FM state. It should be comparable with the total energy difference between AFM and FM states but does not necessarily need to reproduce this difference exactly. The latter may include some other effects, also related to the magnetic polarization of the ligand states, which can be very different in the AFM and FM structures.

IV Conclusions

Microscopic origin of FM coupling in quasi-2D van der Waals compounds Cr$X_{3}$ (where $X=$ Cl and I) was investigated on the basis of ab initio electronic structure calculations within density functional theory. Although FM coupling in CrCl3 and CrI3 is formally expected from GKA rules for the nearly $90^{\circ}$ Cr-$X$-Cr exchange path, the realization of this rule in DFT calculations is somewhat nontrivial, which requires special attention in the selection of practical methods and approximation for calculations of interatomic exchange coupling. In the considered Cr trihalides, the ligand states play a crucial role in the origin of FM coupling. However, the value of the exchange coupling strongly depends on “philosophy” of how these ligand states should be treated, which is traced back to fundamentals of adiabatic spin dynamics. To certain extend, the exchange coupling depends on the approximation, which is used for the exchange-correlation potential and whether it is treated on the level of LSDA or GGA. We have revealed the microscopic origin of this difference, which is mainly related to the intraatomic Cr $3d$ level splitting. More importantly, this coupling depends on the approximations used for the ligand states. Depending on the approximation, each of which has certain logic behind, the interatomic exchange coupling can vary from weakly antiferromagnetic to strongly ferromagnetic. The best considered approximation, which is consistent with brute force total energy calculations, is to explicitly include the effect of ligand states, but to treat them as electronic degrees of freedom and allow to relax for each instantaneous change of the directions of Cr spins. This strategy can be realized in constrained total energy calculations with the external magnetic field, which can be efficiently treated in the framework of self-consistent linear response theory.

Appendix A

$d$ model decorated with on-site Coulomb interactions

In this appendix we explore the effect of the on-site Coulomb and exchange interactions on the interatomic exchange coupling, following the general strategy of construction and solution of the effective low-energy model for the Cr $3d$ bands,

[TABLE]

constructed in the basis of Wannier orbitals review2008 , where $\hat{c}^{\phantom{\dagger}}_{ia\sigma}$ ($\hat{c}^{\dagger}_{ia\sigma}$) stands for the creation (annihilation) of an electron in the Wannier orbital $a$ at site $i$ with spin $\sigma$. We start with a non-magnetic band structure, obtained in the LMTO method in the local density approximation (LDA). Since Cr $3d$ bands are well isolated from the ligand bands, both in CrCl3 and CrI3 (see Fig. 10), the construction of such model is rather straightforward.

As in Se. II, the one-electron parameters $H_{ij}^{ab}$ are identified with the matrix elements of LDA Hamiltonian in the Wannier basis. The only difference is that, since now we are dealing with the non-magnetic band structure, these parameters do not depend on spin indices. The parameters of screened Coulomb interactions, $U_{abcd}$, have been evaluated in the framework of constrained random-phase approximation (cRPA), as described in Ref. review2008 . The matrix $\hat{U}=\left[U_{abcd}\right]$ can be fitted in terms of the on-site Coulomb repulsion $U=F^{0}$, the intraatomic exchange interaction $J=(F^{2}+F^{4})/14$, and “nonsphericity” $B=(9F^{2}-5F^{4})/441$ ($F^{0}$, $F^{2}$, and $F^{4}$ being radial Slater’s integrals), responsible for the charge stability, and first and second Hund’s rules, respectively. This fitting yields (in eV): $U=1.79$ ($1.15$) , $J=0.85$ ($0.78$), and $B=0.09$ ($0.07$) for CrCl3 (CrI3). We note that the screened $U$ is not particularly large. Furthermore, $U$ is smaller in CrI3 due to proximity of the I $5d$ band and, therefore, more efficient screening of Coulomb interactions in the Cr $3d$ band by the ligand $p$ band review2008 .

Then, the model (12) can be solved in the mean-field Hartree-Fock approximation and the exchange parameters evaluated using Eq. (3) review2008 . Alternatively, one can use superexchange theory PWA ; KugelKhomskii , which yields similar parameters $J_{i}$. The results are summarized in Table 6.

In CrCl3 all interactions are antiferromagnetic. Therefore, the FM order is clearly unstable, both in and between the $xy$ planes. At the first sight, the situation in CrI3 looks somewhat better: at least, the nn coupling $J_{1}$ remains ferromagnetic. However, closer analysis shows that $J_{1}$ is counterbalanced by other AFM interactions, particularly $J_{3}$ and $J_{6}$ in the $xy$ plane as well as $J_{2}$ between the plane. For instance, the parameter $\sum_{i}J_{i}$, which is the measure of Curie temperature in the mean-field approximation and also the Curie-Weiss temperature, appears to be negative. Therefore, the FM state is unstable also in CrI3, contrary to the experimental data. This is consistent with qualitative analysis based on Eq. (4): the Coulomb $U$ increases the $t_{2g}$-$e_{g}$ level splitting $\Delta^{\uparrow\uparrow}_{te}$ and thus suppresses the FM coupling.

Thus, to conclude this Appendix we would like to stress the following: (i) The simple $d$ model (12), incorporating on-site Coulomb and exchange interactions, fails to describe properly the FM ordering in CrI3 (and FM ordering in the $xy$ plane of CrCl3). The proper model should include the effects of polarization of the ligand states; (ii) The screened Coulomb interaction $U$ is relatively small in CrCl3 and particularly in CrI3. Thus, plain LSDA (GGA) should provide a good starting point for the analysis of magnetic properties of these compounds, at least on a semi-quantitative level. To certain extent, large $t_{2g}$-$e_{g}$ splitting in QE calculations, based on GGA, mimics the effect of small on-site Coulomb repulsion $U$. A similar conclusion has been reached recently for CrO2, which is another important FM compound, in the joint experimental-theoretical study combining soft-x-ray angle-resolved photoemission spectroscopy and first-principles electronic structure calculations CrO2ARPES .

Appendix B Calculation of Stoner Matrix

In this Appendix, we briefly discuss the solution of matrix equation,

[TABLE]

for the matrix $\hat{I}$. Here, $\hat{b}$ and $\hat{m}$ stand for the $z$-component of, respectively, the exchange field and magnetization density, derived from LSDA (GGA) calculations. For simplicity, we drop the site index $i$.

The basic idea is to use the representation that diagonalizes $\hat{m}$: $\hat{\tilde{m}}=\hat{O}^{T}\hat{m}\hat{O}$, where $\hat{\tilde{m}}=[\tilde{m}^{a}\delta^{ab}]$. Then, Eq. (13) can be rewritten as

[TABLE]

where $\hat{\tilde{b}}=\hat{O}^{T}\hat{b}\hat{O}$ and $\hat{\tilde{I}}=\hat{O}^{T}\hat{I}\hat{O}$, for which the matrix elements of $\hat{\tilde{I}}$ can be found as

[TABLE]

Finally, $\hat{I}$ is be obtained from $\hat{\tilde{I}}$ as $\hat{I}=\hat{O}\hat{\tilde{I}}\hat{O}^{T}$.

For further analysis, it is convenient to use the representation, which diagonalizes $\hat{I}$. The results are summarized in Fig. 11.

At the Cr-sites, the eigenvalues of $\hat{I}$ split in two groups of levels: threefold nearly degenerate $t_{2g}$ and twofold degenerate $e_{g}$. The $t_{2g}$ levels have practically the same energies for CrCl3 and CrI3 (about $3.3$ eV), while the position of $e_{g}$ levels is different. However, this is to be expected: according to Eq. (15), the matrix elements $\tilde{I}^{aa}$ are inversely proportional to $\tilde{m}^{a}$. The value of $\tilde{m}^{a}$ for the $e_{g}$ states depend on the $dp$ hybridization, which controls the degree of admixture of formally unoccupied Cr $e_{g}$ states into occupied $p$ band of ligands (see Fig. 2). This admixture is substantially larger in CrI3, which yields larger $\tilde{m}^{a}$ and, therefore, smaller $\tilde{I}^{aa}$ for the $e_{g}$ levels.

The averaged value of $\hat{I}$, which enters the regular Stoner model, can be estimated from $b_{\rm av}=I_{\rm av}M$, where $b_{\rm av}=\frac{1}{N}{\rm Tr}_{L}\{\hat{b}\}$, $M={\rm Tr}_{L}\{\hat{m}\}$, and $N$ is the number of orbitals ($N=5$ for Cr $3d$ and $N=3$ for ligand $p$ states) as $I_{\rm av}=\frac{1}{N^{2}}{\rm Tr}_{L}\{\hat{I}\}$. For the Cr $3d$ states, it yields $I_{\rm av}=1.55$ ($1.08$) eV in the case of CrCl3 (CrI3), which is pretty close to typical atomic values of the Stoner parameter for the $3d$ elements Gunnarsson . At the ligand sites, some of the eigenvalues of $\hat{I}$ are negative, so as $I_{\rm av}$. The situation is certainly different from the atomic limit. However, it should be understood that $\hat{b}$ and $\hat{m}$ at the ligand sites have a different microscopic origin: $\hat{b}$ is the difference of site-diagonal elements of $\hat{H}^{\uparrow,\downarrow}$ between spins $\uparrow$ and $\downarrow$. This difference is small (see Fig. 3) and does not play a decisive role in the formation of local moments at the ligand sites. These moments are induced by the hybridization with the Cr $3d$ sites and these processes are controlled by off-diagonal elements of $\hat{H}^{\uparrow,\downarrow}$ with respect to the site indices. Thus, the direction of $\hat{\boldsymbol{m}}$ at the ligand sites does not necessary coincide with the one of $\hat{\boldsymbol{b}}$, as it follows from the above analysis.