Linear response theories for interatomic exchange interactions

The linear response is a perturbation theory establishing the relationship between given physical variable and the external field inducing this variable. A well-known example of the linear response theory in magnetism is the susceptibility relating the magnetization with the magnetic field. In 1987, Liechtenstein et al came up with the idea to formulate the problem of interatomic exchange interactions, which would describe the energy change caused by the infinitesimal rotations of spins, in terms of this susceptibility. The formulation appears to be very generic and, for isotropic systems, expresses the energy change in the form of the Heisenberg model, irrespectively on which microscopic mechanism stands behind the interaction parameters. Moreover, this approach establishes the relationship between the exchange interactions and the electronic structure obtained, for instance, in the first-principles calculations based on the density functional theory. The purpose of this review is to elaborate basic ideas of the linear response theories for the exchange interactions as well as more recent developments. The special attention is paid to the approximations underlying the original method of Liechtenstein et al in comparison with its more recent and more rigorous extensions, the roles of the on-site Coulomb interactions and the ligand states, and calculations of antisymmetric Dzyaloshinskii–Moriya interactions, which can be performed alongside with the isotropic exchange, within one computational scheme. The abilities of the linear response theories as well as many theoretical nuances, which may arise in the analysis of interatomic exchange interactions, are illustrated on magnetic van der Walls materials CrX 3 ( X= Cl, I), half-metallic ferromagnet CrO2, ferromagnetic Weyl semimetal Co3Sn2S2, and orthorhombic manganites AMnO3 ( A= La, Ho), known for the peculiar interplay of the lattice distortion, spin, and orbital ordering.


Introduction
On many occasions our image of magnetism rests on the picture interacting spins (e i and e j ) attached to the atomic sites (i and j). If the system is isotropic, such interactions have a form of the scalar product e i · e j [1][2][3][4][5]. The interactions is ferromagnetic (FM) if e i · e j = 1 and antiferromagnetic (AFM) if e i · e j = −1. If i and j are no longer connected by the spacial inversion, the spins tend to align neither ferromagnetically nor antiferromagnetically. 1 The corresponding interaction, which is called Dzyaloshinskii-Moriya (DM) interaction, is given by the cross product [e i × e j ] and driven by the relativistic spin-orbit (SO) coupling [6,7]. Thus, it is always nice to have a transparent toy model, which would explain that certain material has a particular magnetic structure because some interactions are strong or weak, ferromagnetic or antiferromagnetic, etc. The experimental inelastic neutron scattering data are typically fitted to extract parameters of such physically meaningful model. In theory, the spin model can be constructed by averaging the energies of interatomic interactions over non-magnetic degrees of freedom [1][2][3][4][5]7]. In this article we will explain how the interatomic exchange interactions can be generally derived starting from the electronic structure obtained in the first-principles calculations.
Let us consider the simplest possible model: where J ij is the isotropic exchange, d ij = (d x ij , d y ij , d z ij ) is DM vector, and the spin moments e i are normalized to the unity: |e i | = 1. Although we will be primarily interested in the behavior of isotropic interactions, it is appears to be possible to consider J ij in the combination with d ij within one computational scheme. The reason will become clear in a moment. Our goal is to find parameters of this model using the information about the electronic structure. Of course, the model (1) is an approximation as there is no reason why the energy of a general magnetic system should have such a simple form and be described exclusively by bilinear interactions. There are only few microscopic mechanisms, which are consistent with the form of Equation (1). These are the direct Heisenberg exchange [1], Anderson's superexchange [2], and long-range exchange interactions by Ruderman, Kittel, Kasuya, and Yosida (RKKY) [3-5, 8, 9]. In the first example, this is the property of exchange integral, which is finite only between states with the same spin. In the last two examples, this is the consequence of the 2nd order perturbation theory with respect to, respectively, transfer integrals and intraatomic exchange interactions, that couples localized spins to the conduction electrons.
Nevertheless, there is one more, very special case, where the magnetic energy can be also described by Equation (1). These are the infinitesimal rotations of spins near the equilibrium, as was realized by Liechtenstein et al [10]. Considering the perturbation e i = (θ cos qR i , θ sin qR i , 1 − θ 2 2 ) near the ground state e GS = (0, 0, 1), one can find the energy change (per one unit cell) caused by interactions between the transversal components of spins. For the model (1), this energy change is given by where X q = j X 0j exp(−iq · R j ) is the Fourier image of X ij . Then, the basic idea is to extract the same energy change from the electronic structure calculations, typically within spin-density functional theory (SDFT) [11][12][13], and map it on Equation (2). This should give us the parameters of exchange interactions J q and id z q (or J ij and d z ij after the Fourier transform to the real space). Since the perturbation of e i is chosen in the form of the conical spin spiral (which is compatible with the DM interactions), J ij can be considered in the combination with d z ij in the one computational scheme, where the energy change is uniquely specified by q [14]. Basically, this is a perturbation theory, which can be formulated in terms of the response function (or the susceptibility).
One of the most attractive points of the infinitesimal rotations of spins is that the bilinear form of Equation (1) remains valid irrespectively on which microscopic mechanism stands behind the interaction parameters. It can be superexchange [2], RKKY [3][4][5], double exchange [15] or any other mechanism, provided that the rotations are small. Even biquadratic exchange [16] for small θ can be reformulated in the bilinear form (1). Without SO coupling, the energy change caused by the infinitesimal rotations of spins can be always described by the Heisenberg model and in this sense the method is very universal.
Another important point of the work of Liechtenstein et al [10] is that they have proposed a practical scheme for calculating the exchange parameters and proved for these purposes the magnetic force theorem [17][18][19], which justifies the use of the singleparticle energies, obtained from the Kohn-Sham (KS) equations in SDFT [12,13], for evaluating the energy change caused by the infinitesimal rotations of spins. The theorem greatly simplifies the calculations and improves the numerical accuracy.
The basic variable of SDFT is the magnetization densitym. In the ground state,m is controlled by the exchange-correlation (xc) fieldb. Therefore, instead of rotatingm, Liechtenstein et al [10] have proposed to rotateb, assuming thatm will automatically rotate by the same angle. This leads to the commonly used expression for J ij : which is nothing but the 2nd order perturbation theory for the single-particle energy, formulated in terms of the single-particle Green's functionsĜ σ ij with the spins σ = ↑ or ↓ (ε F being the Fermi energy, Tr L stands for the trace over orbital indices, and Im denotes the imaginary part). Equation (3) is known for almost 40 years and was successfully applied for the analysis of interatomic magnetic interactions in various substances [20][21][22][23][24][25][26][27].
Nevertheless, there are also open questions. Particularly, several authors have raised doubts that the true energy change caused by the infinitesimal rotations of spins can be described by rotating onlyb and suggested that it should include the additional contribution steaming from the external magnetic field, which is needed to control the direction of the magnetization [28][29][30]. Thus, Equation (3) may be incomplete. Presumably, the most persuasive arguments were given in 2003 by Bruno [29], who proposed how Equation (3) should be corrected. Surprisingly, however, that even 20 years later after this publication there is no systematic analysis of the problem: Equation (3) is widely used, but little is known how good it is. The problem is complicated by rather common misunderstanding putting the equality between Equation (3) and more fundamental magnetic force theorem.
Another questions is what is the right object to rotate? The spin model (1) is typically formulated on the lattice. Therefore, the magnetization should be also associated with the atomic sites, in some basis of atomic-like orbitals. Then, one can define and rotate the local moments, which are scalars. Alternatively, if there are several atomic orbitals per magnetic site, one can define the magnetization matrix and rotate it. For instance, Equation (3) implies such matrix form. Nevertheless, which construction is more suitable, based on rotations of the scalar moments or the magnetization matrices, is absolutely unclear.
Then, what shall we do with the ligand sites, which can hardly be the source of the magnetism, but frequently carry an appreciable magnetization due to the hybridization with the magnetic transition-metal sites? The contributions of such ligand states are typically ignored, and the interactions (3) are computed only between the transition-metal sites without the justification. On the other hand, there are well-known Goodenough-Kanamori-Anderson (GKA) rules [31][32][33][34], which state, among others, that in certain circumstances the exchange interactions between the transition-metal sites can be controlled by the effective Stoner coupling on the intermediate ligand sites. Such coupling is typically added empirically to correct J ij given by Equation (3) between the transition-metal sites [35][36][37]. However, if the theory is general enough, it should include all such contributions automatically.
The aim of this article is to review some basic ideas as well as more recent developments related to the use of the linear response methods for the analysis of interatomic exchange interactions and give clear answers to all above questions. The general theory is discussed in Section 2. Then, Sections 3-5 deal with practical examples for several types of compounds, where our main idea is to explain the theoretical nuances, which may arise in various parts of calculations of the interatomic exchange interactions: the applicability of Equation (3) and its refinements, the role of on-site Coulomb correlations, the contributions of the ligand sites, the merging of correlated and uncorrelated bands, etc. Short Section 6 outline other developments beyond the main scopes of this review. The article is summarized in Section 7. Two Appendices deal with the construction of the tight-binding (TB) Hamiltonians using the band structure obtained in first-principles calculations and the evaluation of magnetic transition temperature for the spin model in the random phase approximation (RPA).

Basic idea, notations, and conventions
In the magnetic equilibrium, the 1st derivative of the total energy with respect to a small change of the magnetization m(r) vanishes and the energy changes is described by the 2nd derivative. In a general sense, the magnetic force theorem states that not only the total energy but also its 2nd derivative with respect to the infinitesimal rotations of the magnetization is the ground-state properties as it can be expressed via the eigenvalues and eigenfunctions of the ground state. This statement can be traced back to fundamentals of the quantum mechanics, where the system can be measured only via perturbations. Therefore, it is logical that the 2nd derivative can be connected to the properties of the ground state. The response function (or the susceptibility) is the useful tool, which establishes such connection by means of the perturbation theory.
In practical terms, we will deal mainly with SDFT [11][12][13], where the groundstate magnetization density and the total energy are described with the help of the single-particle spin-dependent KS Hamiltonian H σ (r) with some local self-consistent potential incorporating all effects of exchange and correlations [12,13]. This locality implies that the change of the potential in certain point r depends only on the change of magnetization in the same point r. As a consequence, the energy change caused by the infinitesimal rotations of the magnetization can be presented in the form of pairwise interactions. Without SO coupling it corresponds to the isotropic bilinear Heisenberg model. This is a general property of the 2nd order perturbation theory with the local potentials.
From the viewpoint of analysis and interpretation, it is more convenient to adopt the TB representation, which deals with the atomically resolved properties emerging from the solution of some lattice model. Another advantage of the TB representation is the on-site Coulomb interactions, which can be easily incorporated into the model. The purpose of these interactions is to correct limitations of the local density approximation (LDA) or the generalized gradient approximation (GGA), which are derived in the limit of homogeneous electron gas and typically used to describe the effects of exchange and correlations in SDFT. Mathematically, this can be done by constructing the orthonormal basis of localized Wannier functions centered at the atomic sites [38]. The KS Hamiltonian in this basis is the matrix specified by the lattice (i, j) and orbital (a, b) indices,Ĥ σ ≡ [H σ ia,jb ], so as all other matrices constructed fromĤ σ . For instance, the Green function isĜ σ (ε) = [ε −Ĥ σ ] −1 and the magnetization is given bym = Θ(ε F −Ĥ ↑ ) − Θ(ε F −Ĥ ↓ ), in terms of the Heaviside function Θ. We assume that the xc potential in this TB representation remains local in the sense that on each atomic site i it depends only on the magnetizationm i on the same site. It is a common practice to relate the magnetic part of such potential, the so-called xc fieldb, with the site-diagonal elements of the TB Hamiltonian:b i =Ĥ ↑ ii −Ĥ ↓ ii . However, suchb is illdefined because it ignores the non-local contributions steaming from the off-diagonal part ofĤ σ ij [25]. A more consistent definition of the localb in terms of the response function and the ground-state magnetization will be given in Section 2.6.
Other conventions can be formulated as follow: • For periodic systems, [Ĥ σ ia,jb ] can be Fourier transformed to [H σ µa,νb (k)], with µ and ν denoting the atomic positions within the primitive cell. Furthermore, the analysis throughout this paper assumes the use of the periodic gauge H σ µa,νb (k + G) = H σ µa,νb (k) for any reciprocal lattice translation G [39]; • The n×n matrixâ µ , specified by the n orbital indices, can be viewed as the column vector a µ of the length n 2 . The scalar product a † µ · b µ is the shorthand notation for Tr L {â µbµ } (where a † µ is the row vector corresponding to the column vector a µ ). In the case of Cartesian vectors (such as the magnetization matrix or the magnetic field interacting with the magnetization), the notation a † µ · b µ stands for the regular scalar product with the summation over the orbital indices. The notation a † · b implies the summation over the atomic indices as well; • The n×n×m×m tensor A = [A ab,cd ], with first two orbitals (ab) residing on the site µ and last two orbitals (cd) residing on the site ν, can be viewed as the n 2 × m 2 matrixÂ µν . The constructionÂ µν b ν implies the summation over the two orbital indices on the site ν; • The 2nd derivative is a local probe and the interaction parameter depends on the point in which it is calculated. For instance, considering the simplest interaction energy E = −J cos ϕ between two spins in the bond, the 2nd derivative near FM (ϕ = 0) and AFM (ϕ = π) configurations of spins will be, respectively, J and −J. Nevertheless, as it is typically done, we will additionally change the sign of interaction parameters for the antiferromagnetically coupled bonds, thus adopting the universal definition where J > 0 and < 0 stands for the FM and AFM interactions, respectively.

Spin spirals, spin-orbit coupling, and Dzyaloshinskii-Moriya interactions
The SO interaction is known to include the spin-diagonal, ξ 2L zσz , as well as off-diagonal, ξ 2 (L xσx +L yσy ), parts, where ξ is the SO coupling parameter,L = (L x ,L y ,L z ) is the vector of angular momenta, andσ = (σ x ,σ y ,σ z ) is that of the Pauli matrices. The antisymmetric DM interaction d = (d x , d y , d z ) emerges in the 1st order of ξ and generally one should be able to calculate all three vector projections onto x, y, and z (unless they are related by the symmetry properties). Nevertheless, by proper rotations of the coordinate frame, which transform xyz to zyx and zxy, d x and d y can be viewed as d z in the new coordinate frame. Therefore, we need the numerical procedure only for calculating d z . The important point in this respect was realized by Sandratskii [14], who suggested that in order to calculate d z , it is sufficient to consider only spin-diagonal part of the SO coupling. His idea was based on a simple observation that DM interactions give rise to spiral magnetic structures [40], which can be regarded as "eigenstates" of the spin model (1). Therefore, the energies of the model (1) can be uniquely specified by the vectors q, describing propagation of the spin spiral. Then, the same should hold for the electronic model, which is used for the mapping onto the spin one, and the spin spirals should be amount possible magnetic solutions of such model. In practical terms, this means that the electronic states should obey the generalized Bloch theorem, which combines translations with the SU(2) rotations of the spiral spin texture [41]. Nevertheless, this theorem can be applied only if, in the local coordinate frame where all spins are parallel to z, the spin is the good quantum number so that the Hamiltonian H σ remains diagonal with respect to the spin indices σ. Therefore, suchĤ σ can include the diagonal part of the SO coupling, but not the off-diagonal one.
The method is suitable for DM interactions, but not for the magnetic anisotropy, which emerges in the 2nd order of the SO coupling and typically include both diagonal and off-diagonal contributions. This is again in line with the idea of the spin-spiral approach: the DM interactions give rise to the spin spirals, while the magnetic anisotropy acts against them, by deforming the spin spirals and locking them to the crystallographic lattice [42,43]. The alternative to the spin-spiral technique is to work in the real space separately for each magnetic bond [44][45][46][47]. Such methods, which have certain limitations, will be briefly considered in Section 6.1.

General expression for the energy change
As was already pointed out, our basic idea is to "excite" the spin spiral, rotating the ground-state magnetizationm GS = (0, 0,m) as (see Figure 1) and evaluate the interactions between the transversal "fluctuations" of the magnetizationm ⊥ q,i = (cos qR i , sin qR i , 0) θm for small θ. In order to induce sucĥ m ⊥ q,i , we have to apply the external field h q,i = (cos{qR i + α q }, sin{qR i + α q }, 0)ĥ (5) in the direction perpendicular to the ground state magnetization. If the SO interaction is included,ĥ q,i is not necessarily parallel tom ⊥ q,i as the latter can experience the effect of the interaction d z , which tend to additionally rotate the magnetization in the xy plane. Therefore, the phases α q are needed to compensate the effect of DM interactions (see Figure 1). For small α q ,ĥ q,i can be written aŝ whereĥ 0 q,i = (cos qR i , sin qR i , 0)ĥ and n z = (0, 0, 1).  Figure 1. (a) Conical spin spiral, which is assumed in calculations of interatomic exchange interactions: q is the propagation vector, h 0 is the constraining field inducing the transversal magnetization m ⊥ i without spin-orbit coupling, and m i is the rotated magnetization on the site i. (b) Configuration of constraining field in the xy plane: h 0 is required to induce given transversal magnetization m ⊥ , while the additional perpendicular field, αn z × h 0 , is required to compensate the additional rotation of m ⊥ caused by DM interaction d z .
The corresponding energy can be evaluated in the framework of constrained SDFT [29,48] as: where T and E xc are, respectively, the kinetic and exchange-correlation (xc) energies, while the last term controls the size of the transversal magnetization. For simplicity, we drop here all irrelevant dependencies of E on the charge density. Then, the kinetic energy can be expressed as the sum of the occupied KS singleparticle energies, E sp , calculated for the external field h q and the xc field minus the interaction energy of m q with h q and b q [12,13], yielding The important property of the xc energy in this respect, being the consequence of fundamental gauge invariance of the density functional theory [49], is that rotations of the spin magnetizationm GS →m q,i do not change E xc [50][51][52]. Furthermore, since E xc [m q,i ] = E xc [m GS ], the rotation of the magnetization will rotate the xc field by the same angle:b Therefore, b † q · m q does not change either and Equation (8) will lead to the following energy change: Then, δE sp can be evaluated by treating h q + b q as a perturbation, to the 2nd order in h q + b ⊥ q and the 1st order in the longitudinal change of the xc field, − 1 2b θ 2 . The details are elaborated in ref. [53], leading to the simple but general expression: In fact, this result is well anticipated. On the one hand, δE q should be proportional to h q as there would be no energy change without the external field. On the other hand, there only possible interaction of h q with the constrained magnetization m ⊥ q is the scalar product given by Equation (10). It may look incomplete as δE q does not seem to know anything about α q and the DM interactions. Nevertheless, all necessary information is in Equation (10) and in Section 2.5 we will show how it should be used to derive practical expressions for the isotropic exchange and DM interactions.
In addition to the rotations given by (4), the magnetization can experience the longitudinal change, which is caused by these rotations. It will affectm, resulting in an additional change of each of the terms in Equation (8). Nevertheless, these contributions can be shown to cancel out in the lowest order of θ [10, 51].

Response tensor
The response theory is basically the perturbation theory relating the small change of the potential v with the induced density n: n =R v. Spin-dependent v can be generally specified by four elements: Then, each v σσ ′ induces the corresponding change n σσ ′ as where the rank-4 tensor R σσ ′ can be found in terms of the 1st-order perturbation theory for the wave functions [54]. In our case, the perturbation is h q + b ⊥ q and our goal is to evaluate m ⊥ q . Then, it is convenient to use the local coordinate frame wherê h i = (1, α q , 0)ĥ 0 , which is obtained by rotatingĥ q,i about z by the angles −qR i , and employ the generalized Bloch theorem, combining lattice translations with the SU(2) rotations of spins [41]. This will lead to the additional shift of the k-mesh for the states with σ =↑ relative to those with σ =↓. Moreover, sinceĥ i = (1, α q , 0)ĥ 0 corresponds to v h = 1 2 we have to consider onlyR ↑↓ andR ↓↑ . Then, the perturbation theory yields where ε σ mk and |C σ lk = [C aσ lk ] are, respectively, the eigenvalues and eigenvectors ofĤ σ , and f σ mk ≡ Θ(ε F − ε σ mk ) is the Fermi distribution function. The orbital indices in each of the pairs ab and cd belong to the same atomic sites in the unit cell. Then, R ↑↓ can be obtained from R ↓↑ using the proeprty IfĤ σ remains invariant under the time reversal (e.g., without SO interaction), Equation with the summation running over the 1st Brillouin zone (BZ). Then, using the definition (11), one can find: In the local coordinate frame, they should correspond to the magnetization m ⊥ = n ↑↓ + n ↓↑ along x: requires that the perpendicular to it magnetization along y, i( n ↑↓ − n ↓↑ ), should vanish (so as the y component of the xc field) according to our constraint conditions. These are the equations for h 0 and α q for given m ⊥ = θ m. Their meaning is very straightforward. For instance, in Equation (16), the isotropic paty of the magnetization α qR + q h 0 , which is induced by α q h 0 along y, is compensated by the one, which is induced due to the DM interaction by the field h 0 + b ⊥ acting in the perpendicular direction x. The same is with Equation (15), where m ⊥ also has two components: the isotropic one, induced by h 0 + b ⊥ along x, and the one caused by the DM interaction, transferring the effect of the magnetic field α q h 0 , applied along y, to the magnetization along x. This explains how one can naturally separate the contributions of the isotropic and DM interactions in Equation (10).

Exchange interactions
The next step is the mapping of the total energy change (10) onto the spin model: where we explicitly consider the possibility of having several magnetic sublattices. In the local coordinate frame, Equation (10), can be rearranged as where we have added and subtracted the xc field b ⊥ . Our strategy is to start with the expression for m ⊥ without the SO coupling, which is given by the 1st term in Equation (15), and then consider the corrections arising in the 1st order of the SO coupling, which are given by the 2nd term. 3 Then, noting that , one immediately finds the expression for the isotropic exchange interactions: Considering the 2nd term in Equation (15), the construction i( h 0 + b ⊥ ) † ·R − q α q h 0 describes the interaction between x and y components of the magnetic field caused by the DM interactions. Corresponding interaction parameter should satisfy the condition where we had to "rescaled" y components of the magnetic field, α q h 0 µ → h 0 µ + b ⊥ µ , in order to specify x and y components of the transversal magnetization by the same set of parameters θ µ and θ ν . Then, using Equation (18) and noting that to the 1st order in the SO couplingQ +R−Q+ = −Q − , one can find that This expression has been obtained in ref. [56] basically heuristically, by the analogy with isotropic interactions and the expression employing the rotations of the xc field, which will be considered in Section 2.8. Here, we have provided a more rigorous proof of Equation (21). The real space parameters can obtained by the Fourier transform of J q, µν and d z q, µν . For practical purposes, it may be more convenient to calculate in terms of onlyQ ↑↓ q , and then relate it with J q, µν and d z q, µν using the property (13), which yields and Thus, J q, µν is related to the average energy of spin spirals propagating in q and −q, while d z q, µν is related to the energy difference [14].

Sum rule and local exchange-correlation field
The sum rule is obtained from the identity which can be further rearranged aŝ whereb =Ĥ ↑ (k)−Ĥ ↓ (k) is assumed to be local (i.e., site-diagonal and not depending on k). Then, integrating over ε and k, and using the definition 14 for the response tensor, one obtains: This sum rule has very straightforward explanation: q = 0 corresponds to the uniform rotation of the ground-state magnetization, where all spins are rotated in the same direction by the same angle. Therefore, the transversal magnetization is described by the same xc field b as in the ground state. Nevertheless, in the TB representation, such field is not necessary local. For instance, in the local spin-density approximation (LSDA), the splittingĤ ↑ (k)−Ĥ ↓ (k) can have interatomic matrix elements and depend on k. In such a situation, it can be important to reenforce the sum rule, by defining new local xc field as b =Q ↑↓ 0 m, which would yield the given ground-state magnetization m. For instance, this is a simple and transparent alternative to the kernel polynomial method, which was recently proposed to deal with nonlocal matrix elements of the xc field [25]. In fact, if the xc field is nonlocal, the total energy change for the infinitesimal rotations of spins is no longer representable in the form of pairwise interactions.

Right object to rotate: magnetization matrices versus local magnetic moments
So far, we did not properly specify the spin object which should be rotated on the magnetic sites in order to obtain the total energy change (10). All above discussions implied that it is the magnetization matrixm µ , while the spin model is typically formulated in terms of the magnetic moments M µ = Tr L {m µ }. Undoubtedly, the rotation ofm µ , as a whole, by the angle θ µ will rotate M µ by the same angle. However, is this choice unique? Are there other perturbations ofm µ , resulting in the same rotations of M µ but preferably at lower energy cost? Here, we will follow the discussion in ref. [53]. Nevertheless, we would like to note that somewhat similar ideas can be found in the work of Antropov et al [57].
Indeed, for the Hermitian matrixm µ , one can always choose the diagonal representationm µ = diag( . . . , m a µ , . . .) with respect to the orbital indices. In principle, each orbital, a, in such representation can be rotated by its own angle θ a µ . Then, the transversal magnetization in the local coordinate frame, where it is parallel to x, will bê m ⊥ µ = diag( . . . , θ a µ m a µ , . . .). Nevertheless, these angles are subjected to the additional constraint because Tr L {m ⊥ µ } should be equal to θ µ M µ . Importantly, this conditions is softer than the rotation ofm µ as a whole, where all orbitals are rotated by the same θ a µ = θ µ . Therefore, it is reasonable to expect that the energy change will be smaller, so as the exchange parameters.
Mathematically, we have to minimize the energy change (10) with the additional condition a θ a µ − θ µ m a µ = 0 on each site µ: where h 0a µ are the constraining fields acting on m a µ and λ µ are the Lagrange multipliers. Then, minimizing δE with respect to θ a µ , it is straightforward to find that h 0a µ = λ µ . Thus, in order to rotate the spin moments at the minimal energy cost, one have to apply the scalar field h µ , i.e. the same for all orbitals a. Moreover, in this case it is convenient to use the "spherically averaged" version of the linear response, wherem µ is replaced by M µ ,b µ is replaced by B µ = 1 n Tr L {b µ }, and R ↓↑ ab,cd (q) is replaced by The corresponding exchange interaction parameters will be given by and , andR σσ ′ q ≡R σσ ′ (q) is the matrix specified by the atomic indices in the unit cell.
In comparison with Equations (19) and (21), based on rotations of the magnetization matrix, Equations (27) and (28) are expected to be more suitable for the analysis of low-energy excitations, of course, provided that the latter can be described by the spin model (1). In the following, these two methods will be denoted asm and M, after the basic variable describing the infinitesimal rotations of spins.

Rotations of the exchange-correlation field as an alternative perturbation
In this section, we consider the original formulation of the linear response theory, as it was proposed by Liechtenstein et al [10], which is frequently called the magnetic force theorem [29]. However, there are two important points about the work of Liechtenstein et al [10], which should be distinguished from each other [58]: • The general claim that, in SDFT, the energy change caused by the infinitesimal rotations of spins can be related to the KS eigenstates in the ground state is certainly correct and should not be revised. This is what is actually called the "magnetic force theorem" stating that the interaction parameters are the ground state properties and can be found by knowing the electronic structure in the ground state; • Nevertheless, the practical expression (3), which was derived by Liechtenstein et al [10] for the exchange interactions, relies on additional approximations and, in principle, can be improved.
The starting assumption of Liechtenstein et al [10] is that since the rotation of the magnetization results in the rotation of the xc field by the same angle (Section 2.3), it is logical to treat the change of the xc field, δb q =b q −b GS , as a perturbation without the constraining field. Then, we have to consider only δE sp in Equation (9) caused by this δb q . Furthermore, the transversal magnetization, m ⊥ ′ q , which is induced by δb q will generally differ from m ⊥ q because, without the constraining field, m q will tend to relax towards the ground state [28][29][30]. The DM interactions, if any, will tend to additionally rotate m ⊥ ′ q relative to b ⊥ q : m ⊥ ′ q ≈ m ⊥0 ′ q + β q n z × m ⊥0 ′ q . Thus, instead of Equation (10), this method relies on (in the local coordinate frame) arising from the single-particle energies for the perturbations caused by the transversal and longitudinal parts of the xc field. Again, the important point here is that m ⊥ ′ deviates from m ⊥ . Otherwise, the right-hand side of Equation (29) would identically be equal to zero, as in Section 2.3. Then, using the definitions and Alternatively, d z q, µν can be obtained from Equation (20) for h 0 µ = 0 and b ⊥ µ = b µ θ µ . Substituting Equation (14) into Equation (30) and Fourier transforming it to the real space, one obtains the well-known Equation (3). Nevertheless, these are the approximate expressions, which can be formally obtained from the exact ones, Equations (30) and (31), replacing m by b andQ byR, with the additional minus sign. In the following, we will refer to this method as "methodb" or "approximate methodb".
In principle, one can also introduce the "spherically averaged" version of this method replacing b ν by B ν andR by R [29], though it is rarely used. Without the SO coupling,R + =R ↑↓ and corresponding exchange interactionsĴ B q ≡ [J B q, µν ] can be written asĴ B q = − 1 2B (R ↑↓ q +Î −1 )B, whereB = diag( . . . , B ν , . . .) andÎ = diag( . . . , −B ν /M ν , . . .) are the diagonal matrices of, respectively, exchange splittings and effective Stoner parameters. By adapting the same matrix form for the "exact" interactions (27) . .), one can find the following expression, connectingĴ q withĴ B q [29]: Thus,Ĵ q can be indeed replaced byĴ B q at least in two cases: (i) the long wavelength limit q → 0 and (ii) the strong-coupling limitB → ∞. In fact, considering the strongcoupling limit in Equation (30) [59,60], one can derive all known expressions for the double exchange J ij ∼ Ĥ ij [15] [61], etc., where U and V is the on-site and intersite Coulomb repulsion, respectively,Ĥ ij are the transfer integrals (see Appendix A), which are typically associated with the matrix elements of the KS Hamiltonian in LDA, and . . . denotes the expectation value in the ground state. The strong-coupling limit for the DM interaction (31) results in the spin-current model [52,62], which can be viewed as the relativistic counterpart of the double exchange mechanism [56]. The expression for RKKY interactions can be also derived starting from Equation (30), but using slightly different philosophy [9]. In this case, b is the field created by localized spins and acting on conduction electrons. Without b, the conduction bands are non-magnetic (and the tensorR is evaluated in this non-magnetic state [9]).

Elimination of the ligand spins
By knowing J q, µν and d z q, µν , one can, in principle, calculate isotropic and DM interactions operating between all sites in the unit cell. Nevertheless, the magnetism of these sites may have completely different origin. For instance, the transition-metal (T) sites in many oxide materials participate as a source of the magnetism, being primarily responsible for the spontaneous time-reversal symmetry breaking, while the oxygen or any other ligand (L) sites behave as "magnetic slaves": although they can host an appreciable portion of the magnetization, it is solely induced by hybridization with the T sites and strictly follow the change of the magnetization on the T sites. The corresponding energy change for each q can be schematically expressed as whereX AB is the matrix specified by the atomic sites of the types A and B, θ A are the polar angles specifying the rotations of the magnetic moments of the type A in the form of the column vector, and θ † A is the corresponding to it row vector. Then, one can try to eliminate the L degrees of freedom by transferring their effect into the interaction parameters between the T sites. This can be done by employing the ideas of adiabatic spin dynamics [63,64] and assuming that the T spins are sufficiently slow so that the L spins have sufficient time to adjust each change in the system of the T spins. Mathematically, this means that for each θ T , θ L can be found from the condition and The corresponding parameters of isotropic and DM interactions can be obtained from X TT using Equations (23) and (24).
In the method M, the matrix inversion in Equation (35) can be combined with the one of the response matrixQ ↑↓ = R ↑↓ −1 to obtain the following expression forX TT : M T is the diagonal matrix of magnetic moments on the sites T, andÎ L is the diagonal matrix of effective Stoner parameters on the sites L. In this expression, the explicit dependence ofX TT onÎ L is incorporated into ∆X TT , whileX 0 TT formally does not depend onÎ L . The parametersÎ L can play a very important role in the theory of exchange interactions. For instance, according to the GKA rules, they largely contribute to the FM coupling in the systems, where the T-L-T bond angle is close to 90 • [31][32][33][34].
In the linear response theories, these effects are incorporated into ∆X TT [65].
Furthermore, the representation (36) allows us to improve the numerical accuracy. Since in most of the systems the ligand band is filled, the matrix elements ofR ↑↓ associated with the L sites are typically small. Therefore, when we invertR ↑↓ in order to obtainQ ↑↓ , we have to deal with very large numbers even for the T sublattice. Thus, the bare interactionsX TT ∼Q ↑↓ TT are typically large and strongly compensated by the second term in Equation (35) [53,65]. Similar situation occurs in the schemem. On the contrary, the calculation ofX 0 TT and ∆X TT using Equations (37) and (38) involves the inversion of only the T block ofR ↑↓ . This procedure is numerically much more stable rather than the inversion of the whole matrixR ↑↓ in Equation (35). The idea of downfolding, somewhat similar to ours, was previously considered in the literature. For instance, Mryasov et al used this technique in order to eliminate the 5d states of heavy Pt atoms and explain the unusual temperature dependence of the magnetic anisotropy energy in the ordered FePt alloy [66]. A simplified approach in the framework of the schemeb, which did not take into account the effects of ligand-ligand interactions andÎ L , was also considered by Logemann et al [67].

Chromium trihalides
In order to illustrate abilities of considered linear response techniques, we start with the detailed analysis of exchange interactions in chromium trihalides CrX 3 (X= Cr and I). These van der Walls compounds crystallize in the rhombohedral R3 structure, which is built from the honeycomb layers as shown in Figure 2a,b [68,69]. The interactions between the layers are weak, but not negligible. For instance, sizable exchange interactions spread up to 6th coordination sphere, in and between the layers, as shown in Figure 2b.
CrI 3 is the ferromagnet with the Curie temperature T C = 61 K [69-71], while CrCl 3 is a antiferromagnet with the Néel temperature T N = 14.1 K [72,73]. However, the AFM transition in CrCl 3 is believed to be followed by another transition to a pseudo-FM phase at T C ∼ 17 K, where the FM order within the layer coexist with an interlayer disorder [73].
CrI 3 is viewed as the prominent two-dimensional ferromagnet [74], where one of the key ingredients is the strong SO coupling on the I sites [75], which is mainly responsible for the exchange anisotropy and emergence of the long-range FM order at relatively high T C (i.e., contrary to what would be expected from the Mermin-Wagner theorem in the isotropic case [76]). Besides that, CrCl 3 and CrI 3 are regarded as the testbed materials for studying fundamental aspects related to the origin of the ferromagnetism. Namely, why are these materials ferromagnetic? The popular answer is that since the Cr-X-Cr angle is close to 90 • , the interaction is expected to be ferromagnetic due to intraatomic exchange coupling I X on the ligand sites, as prescribed by the GKA rules [31][32][33][34]. However, the GKA rule for the 90 • exchange is not very conclusive, because typically there are several competing mechanisms supporting either ferromagnetism or antiferromagnetism [77]. In fact, Kanamori himself admitted that there are several exceptions from this rule in the 90 • case [34]. Moreover, below we will see that, under certain circumstances, I X can easily become negative and act against the ferromagnetism.
Although CrI 3 is more useful practically, mainly due to the strong SO coupling steaming from the heavy I atoms, CrCl 3 is interesting from the explanatory point of view. Even in LSDA, the electronic structure of CrCl 3 consists of well separated Cr t 2g , Cr e g , and Cl 3p bands, as explained in Figure 2c. 4 Therefore, it is possible to study separately the contributions of each of these bands to the exchange interactions by constructing proper TB models in the basis of Wannier functions [38,78]. Besides that, one can also consider correlated models, both for CrCl 3 and CrI 3 , which include the on-site Coulomb interactions. This procedure is briefly explained in Appendix A. The numerical calculations are performed on the basis of the linear muffin-tin orbital (LMTO) method in the atomic-spheres approximation [79,80], using the mesh of the 10×10×10 points, both for the k-and q-integration.
First, we will review the behavior of interatomic exchange interactions depending on each new ingredient added to the model as well as the type of infinitesimal spin rotations (whether the rotated object isb,m, or M). Then, the detailed comparison with the experimental data will be given in Section 3.6.

3-orbital t 2g model in LSDA
The simplest model for CrCl 3 is the half-filled t 2g model. It contains only occupied ↑-spin t 2g band and unoccupied ↓-spin t 2g band. Since all basis functions of this model are associated with the Cr states, there are no additional complications coming from the ligand states. Even in LSDA, the exchange splitting, B Cr , between the ↑-and ↓-spin t 2g bands is large compared to the bandwidth. Thus, the system should be in the strongcoupling limit. The corresponding parameters of exchange interactions, calculated by rotatingb,m, or M are summarized in Table 1.
In this case we obtain very consistent description as all three methods provide very similar sets of parameters. This is not surprising: the approximate schemeb is justified Table 1. Isotropic exchange interactions in CrCl 3 (in meV): results of t 2g model in LSDA.b,m, and M stand for the methods based on the infinitesimal rotations of, respectively, the xc field, the magnetization matrix, and local spin moments. The notations of J k are explained in Fig. 2.
in the strong coupling limit. Moreover, three t 2g levels are nearly degenerate. Therefore, the asphericity ofm is small and corresponding parameters are practically identical to the ones obtained in the M scheme. However, all interactions are antiferromagnetic, which is quite expected for the half filling [2], but totally contradicts to the experimental situation.

5-orbital Cr 3d model
The next important question is whether the ferromagnetism of CrX 3 can be explained without the ligand states, in the model including both t 2g and e g bands. In LSDA, the ↑-and ↓-spin bands are subjected to the exchange splitting B Cr . At the first sign, there is only a small addition compared to the t 2g model -the unoccupied e g band. However, it changes the story dramatically. First, the crystal-field splitting between t 2g and e g bands is comparable with B Cr . Moreover, since only transitions between occupied ↑bands and unoccupied ↓-bands contribute toR ↑↓ , such exchange interactions do not know about the existence of the ↑-spin e g band. Thus, although B Cr is large (as in the t 2g model), the behavior of the exchange interactions is also controlled by other details of the electronic structure and the system is no longer in the strong coupling limit. Furthermore, the matrixm in the basis of all five 3d orbitals acquires additional degrees of freedom besides M, which is only the spherical part ofm. All these tendencies are reflected in the behavior of interatomic exchange interactions (Table 2). First, the inclusion of the e g band into the model gives rise to the FM interactions. These interactions prevail in the nearest-neighbor (nn) bonds and considerably weaken the AFM interactions in other neighbors, in agreement with the experimental situation. Then, the schemesb,m, and M provide quite a different description. Compared to the methodsm and M, the approximate method b underestimate the FM interaction J 1 . Moreover, the methodm (in comparison with M) has a tendency to overestimate the in-plane interactions J 1 , J 3 , and J 6 , as expected from the analysis in Section 2.7.
Another important question is how to define the xc field in the TB model. The common practice is to useb =Ĥ ↑ −Ĥ ↓ and take only site-diagonal (local) part of this splitting. However, in LSDA, the off-diagonal elements ofĤ can also depend on spin, giving rise to non-local contributions to the xc field [25]. Then, the use of the site- Table 2. Isotropic exchange interactions in CrCl 3 (in meV): results of the Cr 3d model in LSDA.b,m, and M stand for the methods based on the infinitesimal rotations of, respectively, the xc field, the magnetization matrix, and the local spin moments. The xc field was either computed from the on-site spin splitting of the TB Hamiltonian (H) or obtained from the sum rule (sr). T S is corresponding spin transition temperature in RPA (in K). The methodm yields the ferromagnetic ground state, while the methodŝ b and M result in an incommensurate spin spiral propagating both in and between the planes. The notations of J k are explained in Figure 2b.
diagonal part alone will violate the sum rule becauseR + ↑↓ b with such b will no longer reproduce the ground-state magnetization m (see Section 2.6). Another possibility is to reenforce the sum rule by defining the local xc field b =Q ↑↓ 0 m, which would reproduce the ground state magnetization m. In the 3-orbital t 2g model, this results only in minor change of the exchange interactions. However, starting from the 5-orbital model, the difference become significant (see Table 2). Similar analysis can be performed for the correlated model in the Hartree-Fock approximation (see Appendix A). The non-interacting electron part of the model Hamiltonian in this case is evaluated in LDA. The corresponding electronic structure is shown in Figure 2c,d: since in CrCl 3 and CrI 3 the Cr 3d bands are well separated from the ligand ones, such model can be easily constructed for both compounds. The screened Coulomb interactions are evaluated within the constrained random-phase approximation (cRPA) [82], as explained in ref. [78]. The obtained (averaged) parameters of onsite Coulomb repulsion, intraatomic exchange interaction, and nonsphericity (see Appendix A for explanations) are 1.79 (1.15), 0.85 (0.78), and 0.09 (0.07) eV for CrCl 3 (CrI 3 ). The screened U is not particularly large. Moreover, the screening is more efficient in CrI 3 due to proximity of the I 5d and Cr t 2g bands [78]. The corresponding densities of states are shown in Figure 3. As expected [83,84], the Coulomb U increases the band gap in comparison with LSDA. The band gap is smaller in CrI 3 because U is smaller.
Corresponding parameters of exchange interactions are summarized in Tables 3  and 4 for CrCl 3 and CrI 3 , respectively. The approximate methodb systematically underestimates the FM interactions. For instance, all interactions in CrCl 3 are antiferromagnetic, which clearly contradicts to the experimental situation. In CrI 3 , only J 1 is weakly FM, which is not enough to stabilize the FM ground state [85]. The method M systematically improves the situation: the FM interactions clearly prevail and the FM ground state is stabilized both in CrCl 3 and CrI 3 . 5 The corresponding Curie temperature, evaluated in RPA [86] (see Appendix B for details), is also in reasonable agreement with the experimental data. The methodm has a tendency to overestimate the interactions J 1 , J 3 , and J 6 , making the FM ground state in CrI 3 unstable. Table 3. Isotropic exchange interactions in CrCl 3 (in meV): results of correlated Cr 3d model in the Hartree-Fock approximation.b,m, and M stand for the methods based on the infinitesimal rotations of, respectively, the xc field, the magnetization matrix, and the local spin moments. T S is the corresponding spin transition temperature in RPA (in K). The methodsm and M yield the ferromagnetic ground state, while the methodb results in an incommensurate spin-spiral state. The notations of J k are explained in Fig. 2b.  Table 4. The same as Table 3 but for CrI 3 . The method M yields the ferromagnetic ground state, while the methodsb andm result in an incommensurate spin-spiral state.
(A) denotes the same parameters calculated in the antiferromagnetic state. Thus, the minimal model, which can capture the FM character of the exchange interactions in chromium trihalides is the 5-orbital model. Formally, the ferromagnetism can be obtained without invoking the ligand states. Nevertheless, it is crucial to consider the unoccupied e g band. Furthermore, it is crucially important to use the exact technique, formulated in terms of the inverse response function. The approximate methodb, which is linear inR, can lead to an incorrect answer.
Generally, the interatomic exchange interactions, defined via infinitesimal rotations of spins, can depend on the magnetic state. In a number of cases, this dependence can be very strong, reflecting the dependence of the electronic structure on the magnetic state [87,88]. Considering how J k in CrX 3 change depending on the method used for their calculations (for instance, M versusb), it is reasonable to expect that the system is no longer in the strong-coupling limit and, therefore, these parameters can also depend on the magnetic state. However, this appears to be not the case. For instance, the exchange parameters calculated in the AFM state of CrI 3 are practically identical to the ones in the FM state (Table 4). This observation is very important and means, for example, that the parameters derived in the ordered FM state can be also used to evaluate the spin transition temperature, T S . The main reason why theb and M methods yield different J k is related to the fact that the crystal-field splitting, 10Dq, is comparable with the intraatomic exchange splitting. However, this 10Dq does not depend on the magnetic state, and therefore does not contribute to the magnetic-state dependence of J k .

Role of the ligand states
Now we turn to the analysis of the most general model, which explicitly include the contributions of the Cr 3d as well we the ligand Cl 3p or I 5p states.
3.3.1. Cr 3d and X p bands in LSDA We start with the analysis of general TB model, including both Cr 3d and X p bands, in LSDA. First, we note that the Cr 3d and ligand p states are antipolarized. In this case M Cr exceeds the nominal value of 3 µ B . Nevertheless, it is compensated by the negative M X , distributed over three ligand atoms, so that the total moment, M tot = M Cr + 3M X , is integer and equal to 3 µ B , as should be for these FM insulators. This antipolarization is a joint effect of the spin splitting on the Cr atoms and the hybridization between the occupied ligand states and unoccupied Cr e g states. Since the ↑-spin Cr e g states are closer to the ligand band, the hybridization is stronger, resulting in the stronger admixture of the ↑-spin Cr e g states into the occupied ligand band (and transfer of the ↑-spin ligand states into the unoccupied Cr e g band). Therefore, when the ligand band is added to the model, we will have M Cr > 3 µ B but M X < 0, as schematically illustrated in Figure 4. For instance, the LSDA yields M Cr = 3.14 (3.37) µ B and M X = −0.05 (−0.12) µ B for CrCl 3 (CrI 3 ). As expected, the effect is stronger in CrI 3 where the Cr e g states are closer to the I 5p band (Figure 2d). Moreover, the I 5p states are more extended in comparison with the Cl 3p ones, resulting in stronger hybridization.
The new aspect of calculations of the exchange interactions is that now they can also depend on the strength of the Stoner coupling I X , which contributes to the second term of Equation (36). Therefore, the first problem we have to solve is how to properly define I. In fact, there are several possible definitions: whereb µ is associated with the intraatomic spin splitting ofĤ σ andm µ is the ground-state magnetization; ( which is nothing but the spherically averaged version of (I); (III) the same as (I), but takingb µ from the sum rule (25); (IV) the same as (II), but taking B µ from the sum rule M =R ↑↓ B; (V) the same as (I), but takingm µ and M µ from the sum rule (and definingb µ as the intraatomic spin splitting ofĤ σ ); (VI) the same as (II), but taking M µ from the sum rule.
The results are summarized in Table 5. One can see that I Cr only weakly depends on the definition (thought it is different in CrCl 3 and CrI 3 ). However, we do not need this parameter in our calculations. On the other hand, I X is very sensitive to the definition [65]. Apparently, one can discard the definitions V and VI as they yield incorrect M tot violating the fundamental property M tot = 3µ B . 6 Furthermore, it makes sense to enforce the sum rule by using the definitions III and IV instead of I and II. Finally, the definitions II and IV is more appropriate for the spherically averaged method M, while the definitions I and III should be used in the combination with the matrix form of the methodsm andb. While I Cr in CrX 3 is close to atomic values (∼0.7 eV) and can be interpreted as the intraatomic exchange integral responsible for Hund's first rule [89][90][91], I X is definitely not. In most of the cases I X is negative (except scheme III for CrI 3 ). This is another manifestation that the moments M X are merely induced by the hybridization with the Cr 3d states, which acts against B X . According to Equation (38), the negative I X will tend to decrease the FM coupling, contrary to a common believe that the ferromagnetism in CrX 3 is driven by the Hund's rule coupling on the ligand sites, as suggested by the phenomenological GKA rules [31][32][33][34].
The parameters of exchange interactions are summarized in Tables 6 and 7, for CrCl 3 and CrI 3 , respectively. We note the following. In comparison with the 5orbital model (Section 3.2), the parameters J k becomes mostly ferromagnetic, except antiferromagnetic J 6 . The exchange interactions are sensitive to I X . This dependence is illustrated only for the method M, but very similar behavior is observed also for the methodsb andm. The use of the M method in combination with I X obtained in the spherically averaged scheme IV substantially improves the agreement with the experimental data for T C (but not for the spin-wave dispersion, which will be considered in Section 3.6). The exchange parameters obtained in the methodm are unrealistically large [53] and not shown here. Furthermore, in the earlier work [53], the approximate methodb was concluded to provide systematically worse description, especially when considers the contributions of the ligand states [53]. However, the situation crucially depends to the two factors: (i) the proper choice for I X ; (ii) the proper definition of the xc fieldb. For instance, the abilities of this method can be substantially improved by enforcing the sum rule in the definition ofb and I X , and taking into account the aspherical contributions to I X , following the definition III.
The sensitivity of the exchange interactions to the parameter I X can be regarded as the weak point of the downfolding technique. Nevertheless, one can propose an Table 6. Isotropic exchange interactions in CrCl 3 (in meV): results of Cr 3d+Cl 3p model in LSDA.b and M stand for the methods based on the infinitesimal rotations of, respectively, the xc field and the local spin moments. The xc field was either associated with the on-site spin splitting of the TB Hamiltonian (H) or obtained from the sum rule (sr). I Cl is the value of the effective Stoner parameter used in the calculations (the results of the methods III and IV from Table 5, in eV). In the L0 method, the xc field on the ligand sites is set to be zero, while that on the Cr sites is set to reproduce the total magnetization M tot = 3µ B . The T C is corresponding Curie temperature in RPA (in K). The notations of J k are explained in Figure 2b.  Table 7. Isotropic exchange interactions in CrI 3 (in meV): results of the Cr 3d+I 5p model in LSDA.  Tables 6 and 7. In comparison with the M scheme, these parameters somewhat overestimate T C but otherwise fall within the range of typical estimates for J k . We will continue to use this scheme in the future and call it the L0 scheme. It is especially good for semiquantitative estimates, which would illustrate the basic trend in the behavior of interatomic exchange interactions. The 2nd term in Equation (36) also vanishes if Q ↑↓ LL = 0, meaning that there is no ligand-ligand interactions. In this sense, there is certain similarity with the method proposed by Logemann et al [67], but reformulated for the scheme M.

Ligand states in correlated Cr 3d band
In this section, we try to extend the correlated 5-orbital model, by expanding its Wannier basis W into the pseudo-atomic Wannier basis, A, of the more general model, containing Cr 3d as well as the ligand p bands. Namely, we still consider only ten Cr 3d bands (for each spin), but transform |C in Equation (12) from the basis W to the basis A: |C W → |C A =T |C W , whereT is the transformation matrix. Our intension here is to consider explicitly the X p states and all the contributions of these states to the exchange interactions in the 5-orbital model. The magnetic moments in the A basis are M Cr = 2.63 (2.56) µ B and M X = 0.12 (0.15) µ B for CrCl 3 (CrI 3 ). In the 5-orbital model, the only occupied ↑-spin t 2g band provides 3 electrons, which are distributed between one Cr and three ligand atoms. Therefore, the Cr and ligand moments in the t 2g band are ferromagnetically coupled (see also Figure 4), while in order to obtain the opposite polarization of these states, it is essential to consider a more general model, which would explicitly include the ligand p band.
Without the ligand bands, the inversion of the matrixR ↑↓ in the A basis becomes unstable as it contains small matrix elements associated with the ligand sites. Nevertheless, one can still try to estimate J k using the method L0, where it is necessary to invert only the subblock ofR ↑↓ in the subspace of the Cr sites. It yields the magnetic moments: M Cr = 2.64 (2.59) µ B and M X = 0.12 (0.14) µ B for CrCl 3 (CrI 3 ), which are close to the ones reported above.
The corresponding exchange interactions and T C for the 5-orbital model reformulated in the pseudo-atomic basis, are clearly overestimated (see Table 8, the rows denoted 5o). However, this is quite in line with general understanding. The Table 8. Exchange interactions in CrCl 3 and CrI 3 (in meV) obtained in the scheme L0 for the 5-orbital model with the ligand states (5o) and after merging this model with the ligand p bands (5o+p). T C is the corresponding Curie temperature in RPA (in K). The notations of J k are explained in Figure 2b. interactions mediated by the tails of the Wannier functions spreading to the ligand sites and transferring there the FM magnetization can be viewed as an analog of the direct Heisenberg exchange [1]. In order to evaluate these direct exchange contributions numerically, one typically performs the 6-dimensional integration in the real space [35,92,93]. Nevertheless, the downfolding method suggests how these contributions can be naturally evaluated within SDFT. The bare direct exchange integrals are typically large and need to be additionally scaled, in the spirit of cRPA, to account for the screening caused by other bands [35,93]. The same situation is here: the attempt to include the ligand states into 5-orbital model only worsen the description of exchange interactions. The discrepancy can be resolved by considering more general model, which would explicitly include the contributions of the ligand p band.

Merging correlated and ligand bands
The next important step is to merge the correlated Cr 3d bands with the ligand bands in the framework of the Cr 3d + X p model. The correlated 5-orbital model is formulated in the Wannier basis for the Cr 3d bands in LDA. After solving this model in the Hartree-Fock approximation, we want to replace the original Cr 3d bands in the allelectron LDA band structure by these correlated bands in the pseudo-atomic A basis. Since such basis states of the 5-orbital model and the ligand p bands are orthogonal to other states, such merging is not unique and an arbitrary energy shift of the Cr 3d and X p bands relative to each other will not change the magnetization (but will change the exchange interactions). Similar problem arises in the LDA+U method [83], where the relative position of the transition-metal 3d and ligand p states is controlled by the empirical double-counting term [51,84]. Therefore, we have to make some empirical conjecture about this merging and we choose it from the "charge neutrality" condition, by requesting the Fermi energy of the correlated model to coincide with the one in LDA. If the system is insulating, Fermi energy is chosen in the middle of the gap. The electronic structure after such merging is shown in Fig. 5.
We would like to emphasize that this electronic structure is rather artificial and would require different sets of the xc fields for the Cr 3d and X p bands (zero for the X p band, but finite for the Cr 3d one), which is hardly useful for our purposes. Therefore, in order to evaluate the exchange interactions, we again employ the L0 technique, which yields the following magnetic moments: M Cr = 3.35 (3.49) µ B and M X = −0.12 (−0.16) µ B for CrCl 3 (CrI 3 ). Thus, after adding the X p band, the L0 method nicely capture the antipolarization of the Cr 3d and X p states, which is again stronger in CrI 3 . However, it should be understood that this effect is "mimicked" by the specific choice of the xc fields, while from the viewpoint of electronic structure itself, the X p band remains unpolarized and the Cr 3d band is polarized as in Section 3.3.2. In the other words, in the L0 method, the imperfection of the electronic structure is corrected by the physical choice of the xc fields (though the electronic structure itself was obtained not for these fields).
The corresponding exchange parameters are listed in Table 8 (and denoted as 5o+p). There is a systematic improvement in comparison with the 5-orbital model with the ligand states: all exchange interactions become smaller so as T C . However, this improvement is only partial as these parameters are still too large and clearly overestimate the tendency towards the ferromagnetism. This demonstrates the complexity of the problem: the L0 method is only as the starting point of the selfconsistent solution, which should take into account the feedback of correlated Cr 3d bands onto the ligand bands via the xc field and result in magnetic polarization of the ligand bands. However, such field will inevitably mix up the Cr 3d and X p bands, which would require to redefine the Wannier basis for the correlated model and recalculate the model parameters. An alternative solution is LDA+U [83,84], which automatically combines the transition-metal and ligand states takes into account the polarization of the ligand bands. Today there are many applications of this method for calculating the interatomic exchange interactions in various transition-metal oxides and other strongly correlated systems [21][22][23]. However, it should be understood that the LDA+U approach is supplemented with additional approximations of purely empirical character, such as the form of the double counting and the choice of the correlated subspace for the solution of the Hubbard-type model [51].

Dzyaloshinskii-Moriya interactions
In this section, we briefly discuss the behavior of DM interactions in CrI 3 , which are driven by the strong SO coupling of the heavy I atoms (so as other magnetic interactions of the relativistic origin [75]). The DM interactions in CrCl 3 are considerably weaker [44,56] and not significant [72]. In addition to d z , two other components of the DM vector, d x and d y , can be computed by rotating the coordinate frame and applying the same Equation (28) (or similar Equations for theb orm rotations). The symmetry of CrI 3 is such that two Cr sublattices (shown by different colors in Figure 2a,b) are transformed to each other by the spacial inversion. Therefore, all intersublattice interactions will identically vanish [6,7]. On the other hand, the interactions within the sublattices can exist and for each vector R connecting two Cr sites, the interactions in the sublattices I and II are related as: d II (R) = −d I (R). The strongest interaction is the azimuthal angle specifying the direction of the bond k in the xy plane (relative to the bond along a in Figure 2a, corresponding to k = 1) and ψ specifies the direction of the DM vector in the plane relative to this bond. For the R3 symmetry, the parameters d ⊥ , d z , and ψ do not depend on k. They are listed in Table 9.
We note that the schemesb and M provide very consistent description for the DM interactions. The small discrepancy for ψ in the correlated 5-orbital model is probably due to the fact that the angle ψ is ill-defined when d ⊥ is small. There is also surprisingly good agreement between results of the correlated 5-orbital model and the Cr 3d + I 5p model in LSDA. However, such agreement is probably fortuitous because the models are very different. 7 The DM interactions in CrCl 3 are considerably weaker (|d z | ∼ 0.02 meV [56]). However, this is to be expected and the order of magnitude difference of the DM interactions in CrCl 3 and CrI 3 well correlates with the strength of the SO coupling on the ligand sites, which differs by the same order of magnitude: ξ Cl = 95 meV versus ξ Cl = 881 meV.

Comparison with experimental data
The experimental parameters of exchange interactions, derived from inelastic neutron scattering data [70,72,73], are summarized in Table 10. 8 On of the interesting features of the experimental spin-wave dispersion in CrI 3 is a ∼ 4 meV gap at the Dirac (K) point, indicating at the existence of large DM interaction d z between 2nd neighbors in the honeycomb plane [70]. Quite expectably, no such feature was observed in CrCl 3 [72,73], where this DM interaction is small. For the isotropic interactions, the agreement between 7 The correlated 5-orbital model includes orbital polarization caused by on-site Coulomb interactions [54]. On the other hand, the Cr 3d+I 5p model in LSDA includes the additional contributions steaming from the magnetically polarized I 5p band. Apparently, these two effects are comparable with each other, which explain a good agreement in Table 9. 8 In order to be consistent with our definition of the spin model, Equation (1), the experimental parameters have been multiplied by S 2 = (3/2) 2 .  Table 6). The calculations are performed for the hexagonal cell, where six branches of ω(q) correspond to six magnetic Cr sublattices.
theoretical and experimental data depends on the model and approximations employed for treating the Coulomb correlations and the ligand states, which we will discuss below. The theoretical DM interaction in CrI 3 appears to be underestimated by factor two, both in correlated 5-orbital model and LSDA for the Cr 3d+I 5p model. The theoretical spin-wave dispersion in CrCl 3 and CrI 3 , calculated for some representative sets of parameters, is plotted in Figures 6 and 7 Table 7). The calculations are performed for the hexagonal cell, where six branches of ω(q) correspond to six magnetic Cr sublattices.
the FM ground state in both the compounds. In CrCl 3 , the dispersion is overestimated by factor two, while in CrI 3 we obtain an overall fair agreement with the experimental data, except for the band splitting in the point K, which is underestimated by factor two. Nevertheless, the situation changes significantly when we explicitly consider the contributions of the ligand states in the framework of the all-electron Cr 3d+X p model in LSDA (panels c and d). Now, the FM state becomes stable already in the schemeb. The spin-wave dispersion in CrCl 3 remains overestimated by factor two. The agreement with the experimental data is better for CrI 3 , again except for the band splitting in the point K. When we turn to the M scheme, which is expected to be more accurate, the theoretical spin-wave dispersion in CrCl 3 is substantially reduced, so that the lowest experimental band is reproduced pretty well. However, the behavior of theoretical upper bands is not satisfactory. This is related to the fact that the theoretical J 1 = 0.64 meV (Table 6) is underestimated by factor three. A similar situation is realized in CrI 3 , where the theoretical J 1 = 1.32 meV (Table 7) is underestimated by the same amount, which worsens the agreement with the experimental data. The problem is partially related to our choice of I X , which appears to be more negative for the scheme M due to the additional spherical approximation, as discussed in Section 3.3.1. For instance, if one used the same I X as in the schemeb (or the L0 scheme, where we do not need any I X ), one could get a better agreement with the experimental data for the upper bands (but not for the dispersion in the entire energy region).
Regarding the value of d z and the band gap, first we note that there is an alternative mechanism of the gap opening, related to the existence of the long-range isotropic exchange interactions beyond the 6th coordination sphere [94]. These interactions were taken into account also in our calculations. 9 However, their effect is not particularly strong. Therefore, the band gap is mainly controlled by d z and the discrepancy with the experimental data is caused by the underestimation of this interaction in the theoretical calculations, as was also reported by Kvashnin et al [44] and Olsen [95]. In this respect, we note that by explicitly considering the contributions of the ligand states in the correlated 5-orbital model (Section 3.3.2) and using the L0 scheme for the evaluation of the exchange parameters, we can easily get |d z | as large as 1.10 meV, which exceeds the value obtained in the regular correlated 5-orbital model by factor five and the experimental value by factor two [70]. Such enhancement is caused by the additional contribution of the basis functions, which explicitly takes into account the effect of the heavy I atoms. Nevertheless, when we try to combine this effect with another contribution stemming from the I 5p band itself (employing for these purposes the merging of the Cr 3d and I 5p bands, as discussed in Section 3.4), again the framework of the L0 scheme, the DM interaction parameter is reduced till |d z | = 0.16 meV. Such strong reduction is apparently caused by the cancelation of contributions in the Cr 3d and I 5p bands. Since the I 5p states in the Cr 3d and I 5p bands are antipolarized, the fact of the cancelation itself is not surprising. However, such analysis clearly demonstrate the fragility of the situation, where the value of d z appears to depend on the delicate balance of two large contributions arising from the Cr 3d and I 5p bands. The merging of these two bands considered in Section 3.4 was probably too crude to explain the experimental situation.

Brief summary
To summarize this section, we would like to stress again the main points: • The minimal model, which captures the magnetic properties of CrX 3 , is the 5orbital model, constructed for the magnetic Cr 3d bands near the Fermi level. The main advantage of this model is the simplicity: in this case there is simply no contributions associated with the ligand states and all calculations of the exchange interactions become pretty straightforward. Nevertheless, it is absolutely essential to use for these purposes the "exact" approach dealing with rotations of magnetic moments M. The approximate scheme, dealing with rotations of the xc fieldsb, severely underestimate the FM interactions and fails to reproduce the FM ground state. The scheme M improves the situation tremendously; • The FM interactions can be additionally stabilized in the all-electron Cr 3d+X p model. However, the exchange interactions in this case depend on the strength of the effective Stoner coupling I X on the ligand atoms and the proper definition of these Stoner parameters is still not completely resolved problem. As soon as the xc field is corrected to satisfy the sum rules for the given magnetization, the approximateb scheme works reasonably well. The main issue in this context is even not the differences between the schemesb and M, but how to properly define I X within each scheme. As an alternative solution, we have proposed the L0 method, which makes the related to I X contributions inactive. This method can be used for FM insulators and half-metallic materials; • Another open question is how to properly merge the correlated Cr 3d bands and ligand X p bands. The exchange interactions can crucially depend on the details of such merging. The method considered in Section 3.4 is probably only the first step in this direction.

Half-metallic ferromagnets
The half-metallicity is basically the peculiar type of the electronic structure where one spin channel is metallic while another one is semiconducting [96]. Most of such materials are ferro-or ferrimagnets. The fully compensated half-metallic ferrimagnets, where 100% spin polarization of the conduction electrons coexists with zero net magnetization, have been also proposed [97]. In this section, we further explore abilities of the linear response theories for the analysis of interatomic exchange interactions in this type of materials with the emphasis on the Coulomb correlations and the ligand states. We consider two such examples: the canonical CrO 2 [98] and more recent Co 3 Sn 2 S 2 , which has attracted a great deal of attention due to the large anomalous Hall effect and other intriguing properties [99][100][101][102][103][104].

CrO 2
The half-metallic ferromagnetism is not very common in stoichiometric transition-metal oxides. Nevertheless, there are exceptions and CrO 2 is one of them, which is widely considered in various applications related to spintronics and magnetic recording [105]. One of the limitations of CrO 2 from this practical point of view is the relatively low T C ∼ 390 K [105]. CrO 2 crystallizes in the rutile structure (the space group P 4 2 /mnm, Figure 8a) [106]. The exchange interactions remain sizable up to at least the 8th coordination sphere. Moreover, since the rutile structure is nonsymmorphic, there are two types of interactions, J 7 and J 8 , which are denoted by superscripts > and <, as explained in Figure 8b. The calculations are performed for the experimental parameters of the crystal structure [106] on the mesh of the 10×10×16 points both for k and q.
The LDA band structure is featured by three separated bands: O 2p in the occupied part, Cr t 2g near the Fermi level, and Cr e g in the unoccupied part. Therefore, one can consider two types of correlated models: the 3-orbital model for the Cr t 2g bands  and 5-orbital model both for Cr t 2g and Cr e g bands. In the former case, the onsite Coulomb and exchange interactions can be described in terms of two Kanamori parameters [107]: the intraorbital Coulomb interaction U and the exchange interaction J , which can be evaluated within cRPA as 2.84 eV and 0.70 eV, respectively [108]. The third parameter of interorbital Coulomb interaction can be obtained from these two as U ′ = U − 2J [107]. The parameters of 5-orbital model are specified by U = 1.98 eV, J = 0.94 eV, and B = 0.09 eV (see Appendix A). The corresponding densities of states, in the Hartree-Fock approximation, are shown in Figure 9. Alternatively, one Table 11. Parameters of exchange interactions in CrO 2 (in meV), as obtained in correlated 3-and 5-orbital models (respectively, 3o and 5o), and in LSDA for the allelectron Cr 3d+O 2p model using the schemesb and M (the infinitesimal rotations of the xc fields and local spin moments, respectively). In the L0 scheme, xc field was redefined to enforce I O = 0. The notations of exchange parameters are explained in Figure 8b. D xx = D yy and D zz are nonvanishing elements of the spin-stiffness tensor (in meV·Å 2 ). T C is the Curie temperature in RPA (in K). can consider the all-electron Cr 3d+O 2p model in LSDA. The half-metallic character of the electronic structure is evident from the LSDA density of states ( Figure 8c). The Coulomb interactions additionally split the ↑-spin t 2g band, placing the Fermi energy inside a pseudogap, and shift the ↓-spin states to the higher energy region, away from the Fermi level. The parameters of exchange interactions are summarized in Table 11. First, there is a large difference in the parameters obtained in the schemesb and M, especially in the Hartree-Fock approximations for correlated 3-and 5-orbital models, where the more accurate method M strengthens the FM interactions practically in all the bonds. The situation is more modest in LSDA for the all-electron Cr 3d+O 2p model: when going from the methodb to the method M, the nearest FM interactions J 1 and J 2 increase only partially, while other longer-range interactions tend to decrease and become more antiferromagnetic.
In LSDA, the Cr and O atoms are antipolarized due to the joint effect of the hybridization and the intraatomic spin splitting, similar to The spin-wave dispersions calculated with different sets of parameters for the 3-and 5-orbital models is plotted in Figure 10 (see Appendix B for details). The nonvanishing xx and zz elements of the spin-wave stiffness tensorD = [D αβ ], describing the dispersion of the acoustic (A) spin-wave branch near the Γ point, are summarized in Table 11. The experimental estimates for the averaged spin-wave stiffness in CrO 2 vary from 70 to 150 meV·Å 2 [110,111]. The Curie temperature is about 390 K [105]. Then, the LSDA values seem to be overestimated, though the situation is rather subtle. On the one hand, the theoretical values of the spin-wave stiffness are within the experimental scatter. On the other hand, the exchange parameters derived for ground state are not supposed to reproduce T C for this metallic system and more sophisticated methods, considering the temperature dependence of the exchange interactions, may be needed [112][113][114]. Anyways, in the light of well-known limitations of LSDA for the transition-metal oxides [83], such moderate disagreement would not be surprising.
Even more interesting situation is realized in correlated 3-and 5-orbital models. At first glance, the approximate schemeb provides a much better description for the spin-wave stiffness and T C in comparison with the scheme M, which is supposed to be more accurate. However, this is true only in the Hartree-Fock approximation, which has serious limitations for the metallic systems. If one goes beyond the Hartree-Fock approximation, the situation can change dramatically. Particularly, as expected for the half-metallic systems [115], the dynamic correlations lead to a strong redistribution of the electronic states. These effects can be treated in the framework of dynamical meanfield theory (DMFT) [116]. The latter can be regarded as an extension of SDFT in which the KS potential is replaced by the frequency-dependent self-energy Σ ↑,↓ (ω) [117]. This method has been applied to CrO 2 in ref. [108], using the same correlated 3orbital model. The exchange interactions were evaluated basically in theb scheme using Equation (3), by considering the infinitesimal rotations of the frequency-dependent xc field ∆Σ(ω) = Σ ↑ (ω) − Σ ↓ (ω) [118]. The dynamic correlations substantially reduce ∆Σ(ω), resulting in much stronger AFM interactions J 7 and J 8 , which overcome the effect of FM interactions J 1 and J 2 , making the FM state unstable. The corresponding spin-wave dispersion is also plotted in Figure 10, where the negative frequencies ω(q) along the directions Γ-X, Γ-M, and Γ-Z mean that the theoretical ground state should be in one of these points. Therefore, it was concluded in ref. [108] that in order to reproduce the FM ground state of CrO 2 , it is essential to extend the 3-orbital model, by adding new ingredients such as the Cr e g and O 2p bands. For instance, the spin-wave stiffness is systematically larger in the 5-orbital model (Table 11), which makes the FM state more stable. Nevertheless, the present analysis also suggests that the problem may be not in the 3-orbital model itself, but in additional approximations underlying the schemeb for the exchange interactions. The method M systematically improves the stability of the FM ground state in CrO 2 . This tendency may be overestimated in the Hartree-Fock approximation. However, a more rigorous treatment of correlation effects in the framework of DMFT, in the combination with the method M, could possibly bring the situation to a better agreement with the experimental data.
The DM interactions can take place between the atoms in different Cr sublattices. The strongest ones are realized in the 2nd coordination sphere (in the combination with J 2 , as shown in Figure 8b). If ǫ = 1 √  M in both cases). The interesting point is that the DM interactions are expected to be strong (and comparable with the ones in CrI 3 ) even though CrO 2 does not contain heavy elements. This is basically the consequence of two factors: (i) the partially filled t 2g states, opening the room for unquenched orbital magnetization; (ii) half-metallic character of the electronic structure, giving rise to the spin-current mechanism of the DM interactions in the ↑-spin t 2g band [52,62]. Nevertheless, d is much smaller than the isotropic interaction J 2 operating in the same bonds.

Co 3 Sn 2 S 2
Co 3 Sn 2 S 2 is a ferromagnet in which small spontaneous magnetization (about 0.3 µ B per Co atom) coexists with relatively high T C = 177 K [119]. It crystallizes in the rhombohedral R3m structure, hosting the kagome lattice of Co ions (Figure 11a-c) [119]. The S atoms sit on the top of the Co 3 triangles, forming with them the network of alternating tetrahedra. There are two inequivalent Sn sites located in (Sn 1 ) and between (Sn 2 ) the kagome planes. Recently, Co 3 Sn 2 S 2 has attracted a lot of attention as a magnetic Weyl semimetal whose nontrivial topology of the electronic states gives rise to a large anomalous Hall effect [99][100][101]. Furthermore, Co 3 Sn 2 S 2 appears to be halfmetallic, as was demonstrated in experimental [102] and theoretical [103,104] studies. We have chosen Co 3 Sn 2 S 2 to illustrate the complexities one can face in the analysis of interatomic exchange interactions for this type of itinerant electron systems. At the present stage, we have in our disposition only all-electron Co 3d+Sn 5p+S 3p model, which was constructed in ref. [104] within GGA [120], by employing the Vienna ab initio simulation package (VASP) [121] for the electronic structure calculation and the maximal localization technique for the Wannier functions [122]. The corresponding density of states is shown in Figure 11h  However, the attempts to construct a more compact model, which would also include the Coulomb correlations, were so far unsuccessful. Formally, the electronic structure in Figure 11h can be divided into fully occupied 19 bands and remaining 8 bands, which are separated by a band gap. However, the construction of the Wannier basis for the upper eight-band manifold is not straightforward, as these Wannier functions will probably reside not on single atomic sites [123]. A compact model for Co 3 Sn 2 S 2 was constructed in ref. [124], but using empirical arguments. Co 3 Sn 2 S 2 provides an interesting example of itinerant magnetism. According to the GGA calculations [104], finite rotations of spins away from the FM ground state (for instance, by forcing three Co spins to form the umbrella texture) eventually collapse the magnetization, so that the system falls in the nonmagnetic state. 10 A similar behavior was found in fcc Ni [125,126] and several Ru-based oxides [127,128]. Moreover, the size of the magnetic moments is expected to change with the temperature, affecting both electronic structure and interatomic exchange interactions. Therefore, the bilinear Heisenberg model is not defined in the global sense (for instance, for describing T C ). Nevertheless, it can still be defined locally, for the analysis of local stability of the FM state with respect to the infinitesimal rotations of spins. The exchange interactions spread at least up to 9th coordination sphere around each Co site, as explained in Figure 11d-g. Moreover, some of these interactions, formally belonging to the same coordination sphere, can be inequivalent (as J 4 and J ′ 4 in the plane, J 5 and J ′ 5 , and J 9 and J ′ 9 between the planes). The calculations have been performed on the mesh of the 20 × 20 × 20 points both for k and q.
The effective Stoner parameters, calculated by using different definitions are summarized in Table 12. As was discussed in Section 3.3, the definitions II and IV should be used in combination with the spherically averaged scheme M, while the definitions I and III are more suitable for the schemesb andm. One can clearly see that all the parameters in Co 3 Sn 2 S 2 strongly depend on the definition. This holds even for I Co , which can easily change by about 50%. This is clearly in contrast with other considered compounds, where the value of I on the magnetic transition-metal site practically did Table 13. Parameters of exchange interactions in Co 3 Sn 2 S 2 (in meV) calculated using the schemesb,m, and M (the infinitesimal rotations of, respectively, the xc field, magnetization matrix, and local spin moments). In the schemeb, the xc field was either associated with the site-diagonal part of the TB Hamiltonian (H) or derived from the sum rule (sr). The roman number in the parentheses stands for the set of the parameters I Sn and I S from Table 12, which was used in the calculations of J k . In the L0 scheme, xc field was redefined to ensure I = 0 on the sites Sn and S. The notations of J k are explained in Figure 11d  not depend on the definition. Nevertheless, the exchange interactions (36) between the Co sites do not explicitly depend on the choice of I Co . The uncertainty with the choice of the parameters I Sn and I S , which strongly depend on the way how they are defined, posses a more serious problem. The parameters I Sn are mainly negative (except I Sn 1 in the case I) and tend to destabilize the FM interactions. On the other hand, positive I S strengthens the ferromagnetism. However, depending on the definition, I S can change by factor four, and I Sn changes by an order of magnitude. This dependence has a strong impact on the exchange interactions, which are summarized in Table 13. As was pointed out before, the attempts to estimate the Curie temperature entirely from J k , defined via the infinitesimal rotations of spins near the FM ground state, are pretty much meaningless in the case of Co 3 Sn 2 S 2 , where proper theories for T C should consider also the longitudinal change of the magnetization [104]. Nevertheless, these J k can be used to evaluate the spin-wave dispersion and the spin stiffness. The latter has been measured experimentally and the nonvanishing elements of the spin-stiffness tensor are D xx = D yy = 803 ± 46 and D zz = 237 ± 13 meV·Å 2 [129], which can be used for comparison with theoretical data.
First, we note that the results of the schemeb strongly depend on the definition of the xc field: the enforcement of the sum rule forb strengthens the FM interactions, bringing the spin stiffness to a good agreement with the experimental data. In this case, the exchange parameters are given mainly by the bare interactions, corresponding to the first term in Equation (35). The corrections caused by the ligand states are small and do not play a decisive role. An interesting situation is realized in the schemem, based on the rigid rotations of the magnetization matrix: as expected, the individual parameters are overestimated (see Section 2.7). Nevertheless, the interactions are longranged and many of them are antiferromagnetic, resulting in the relatively small spin stiffness. The exchange interactions in the scheme M are very sensitive to I Sn and I S . If one uses the definition II, where the xc field in the expression for I is taken from the site-diagonal part of the TB Hamiltonian, I Sn and I S are relatively small so that the exchange interactions are given basically by the first term of Equation (36). In this case, J 1 -J 3 are weakly antiferromagnetic, while other interactions are ferromagnetic (except J 9 , which is antiferromagnetic in all considered methods). As the result, the FM state is stable, though D xx is somewhat underestimated in comparison with the experiment.
The situation changes dramatically, if one uses the definition IV, where the xc field is taken from the sum rule, which substantially strengthens both I Sn and I S . Apparently, negative I Sn plays a dominant role and is responsible for the AFM character of interactions in some of the bonds, particularly J 1 -J 3 become prominently antiferromagnetic and J 5 changes from strongly ferromagnetic to antiferromagnetic. Then, the FM state become unstable and the spins tend to form the 120 • in the xy plane (at least, on the mean-field level). The corresponding spin stiffness is strongly underestimated compared to the experimental data. This is clearly an artefact of the model analysis, which is probably caused by unrealistic estimate for I Sn . At least, the conclusion does not seem to be consistent with the brute-force GGA calculations, where the 120 • alignment of spins in the xy plane leads to the collapse of the magnetic state [104].
The L0 approach allows us to eliminate the dependence of J k on I Sn and I S , assuming that all magnetic moments are induces by B Co , while B Sn = B S = 0 (so as I Sn and I S ). This results only in a small change of magnetic moments in comparison to the plain GGA values reported above: M Co = 0.36, M Sn 1 = −0.05, M Sn 2 = −0.07, and M S = 0.02 µ B (thus, M tot remains equal to 1 µ B ). The exchange interactions in these case reminiscent the ones obtained in the scheme M using the definition II for I Sn and I S with somewhat stronger tendency towards the ferromagnetism, which is reflected in larger values of D xx and D zz .
In principle, the schemesb and L0 both provide a reasonably good agreement with the experimental data for the spin-stiffness constants. Nevertheless, the behavior of exchange interactions is quite different. In the schemeb, the strongest FM interaction is J 1 . The longer-range interactions operating practically at the same distance in the 4th and 5th coordination spheres, in and between the kagome planes (see Figure 11e,g), are also sizable but smaller than J 1 . In the scheme L0, J 1 is weakly antiferromagnetic, while the strongest FM interactions take place in the 4th and 5th coordination spheres. The corresponding spin-wave dispersions are plotted in Figure 12. In the long wavelength limit q → 0, two methods provide very similar description, while the main difference  [130], whereJ k = 1 3 (J k + 2J ′ k ) stands for the averaged exchange interactions in the 4th and 5th coordination spheres. 11 There is certain similarity with the picture provided by the L0 scheme: at least J 1 is negligibly small, J 2 is antiferromagnetic, and the main FM interactions occurs in the next coordination spheres. Nevertheless, there are also the differences. Especially, experimentalJ 5 appears to be smaller thanJ 4 , while theoretical parameters are at least comparable (in all theoretical approaches summarized in Table 13). According to the theoretical analysis, the values of inequivalent parameters J k and J ′ k can be very different and this difference was not considered in the fitting of the experimental spin-wave spectrum. On the other hand, it is not clear whether it can resolve the discrepancy between theoretical and experimental data. It is also unclear whether the experimental data available only in the vicinity of q = 0 are enough to make a decisive conclusion about large complexity of the exchange interactions in Co 3 Sn 2 S 2 .
The Co 3 Sn 2 S 2 structure has three Co sublattices (which are transformed to each other by the threefold rotations). The spacial inversion transforms each Co atom to the same sublattice. Thus, there can be no DM interactions within the sublattices. Nevertheless, the DM interactions between the sublattices are permitted by the R3m space group. The strongest ones take place between nearest neighbors in the kagome plane, which are explained in Figure 13 (and using for these purposes the same fragment of the crystal structure as in Figure 11d). Using the convention, where the directions of the bonds in the triangles can be transformed to each other by the threefold rotations, the DM vectors can be presented in the form: d = d ⊥ [ǫ × n z ] + d z n z . If the FM moment is parallel to z (the experimental situation), d z does not contribute to the energy change, while d ⊥ is responsible for the canting of spins and tends to form the umbrella texture. If e k = (θ sin πk 3 , θ cos πk 3 , 1 − θ 2 2 ) are the directions of spins in the triangle (k= 1, 2, or 3, as explained in Figure 13), the energy gain due to the DM interaction is δE DM = − √ 3d ⊥ θ (per one Co site). The corresponding energy loss due to the isotropic exchange is δE H = 3 4 J 1,23 θ 2 , where J 1,23 is the effective interaction of the atom in the 1st sublattice with the atoms of the 2nd and 3rd sublattices in all the coordination spheres. The values of the parameters obtained in the scheme L0 are J 1,23 = 6.04 meV, d ⊥ = 0.20 meV, and d z = −0.29 meV. 12 Thus, the angle θ can be estimates as 2 • , which perfectly agrees with θ ∼ 2 • obtained in the brute-force GGA calculations with the SO coupling [104].

Orthorhombic perovskite manganites
Perovskite manganites AMnO 3 (where A is the trivalent rare earth or alkaline earth element) have attracted a great deal of attention. Most of them crystallize in the orthorhombic P bnm structure (Figure 14a). The magnetic Mn 3+ ions in the octahedral environment accommodate four electrons, three of which occupy the t 2g states, while the remaining one resides in the doubly degenerate manifold of the e g states. Therefore, the system is subjected to the Jahn-Teller distortion, resulting in the peculiar orbital ordering (the alternation of occupied e g orbitals, as schematically shown in Figure 14e).
LaMnO 3 is the parent material of colossal magnetoresistive oxides [131] and a popular testbed system for studying abilities of first-principle electronic structure calculations, especially in reproducing the insulating behavior, the cooperative Jahn-Teller distortion, and associated with it A-type AFM ordering, where the FM coupling within the orthorhombic ab plane coexists with weakly AFM coupling between the planes (see Figure 14d) [132][133][134]. The orthorhombic distortion systematically increases   with the decrease of the size of the A 3+ ion in the direction La→Ho. At certain point it makes the A-type AFM phase unstable in the ab plane. For instance, the ground state of TbMnO 3 is a spin spiral with the propagation vector q ≈ (0, 1 4 , 0) [135], whereas HoMnO 3 forms the twofold periodic E-type AFM structure corresponding to q = (0, 1 2 , 0) [136,137]. These magnetic superstructures breaks the inversion symmetry, giving rise to the rich multiferroic activity [138][139][140].
In this section we consider LaMnO 3 and HoMnO 3 as two characteristic example and discuss abilities of the linear-response based techniques for describing the change of the exchange interactions, which would eventually lead to the change of the magnetic structure from A to E. We use the experimental crystal-structure parameters reported in [141] and [136] for LaMnO 3 and HoMnO 3 , respectively. The main parameters are summarized in Table 14. Particularly, the MnO 6 octahedra have two long Mn-O bonds in the ab plane and four short ones, which is the signature of the Jahn-Teller distortion [142]. The Jahn-Teller distortion practically does not change when going from LaMnO 3 to HoMnO 3 . On the other hand, the Mn-O-Mn angles change significantly with much stronger deviation from the ideal cubic value of 180 • in the case of HoMnO 3 . This change is accompanied by some shrinking of the orthorhombic lattice along a and c.
The orthorhombic distortion has a profound effect of the electronic structure in LDA/LSDA (Figure 14f,g), splitting the Mn e g band and opening the band gap in the A-type AFM phase. The latter is the consequence of the Jahn-Teller distortion and quasi-two-dimensional character of the A-type AFM order [143]. The occupied Mn e g band is located around −1 eV and well separated from other bands. The corresponding distribution of the e g electron density has the form of the orbital ordering, which is schematically shown in Figure 14e.
The exchange interactions in LaMnO 3 and HoMnO 3 are rather complex (see Figure 14b,c) [37]. Particularly, in addition to the nn interactions in and between the ab planes (J 1 and J ⊥ 1 , respectively), there are several important longer-range interactions, such as: (i) the next-nn interactions between the planes, J 1 2 and J 2 2 ; (ii) the 2nd neighbor interactions in the plane, J a 2 and J b 2 , operating along orthorhombic axes a and b, respectively; and (iii) the 3rd neighbor interactions in the plane, J 1 3 and J 2 3 . These interactions obey certain symmetry properties. For instance, considering Mn site 1 in Figure 14b as the reference point, J 1 3 and J 2 3 will operate in the bonds ±(a, a, 0) and ±(a, −a, 0), respectively. The parameters J 1 3 and J 2 3 around Mn sites 2, 3, and 4 are obtained by the 180 • rotations of these bonds about a, b, and c combined with the lattice shifts by ( a 2 , b 2 , 0), (0, 0, c 2 ), and ( a 2 , b 2 , c 2 ), respectively. The interactions J ⊥ 1 , J 1 2 and J 2 2 control the AFM coupling between the planes, while the formation of the long-periodic magnetic structures in the plane results from the interplay of J 1 , J a 2 , J b 2 , J 1 3 and J 2 3 . The behavior of J 1 3 and J 2 3 is directly related to the orbital ordering: these 3rd neighbor interactions can be regarded as the super-superexchange, mediated by the states of intermediate Mn sites. If the lobes of the occupied e g orbitals are directed toward each other along the bond, as in the case of J 1 3 , the interaction is strong. If the lobes are parallel to each other and perpendicular to the bond, as in the case of J 2 3 , the interaction is weak. Moreover, there are two crystallographically different types of next-nn bonds between the planes, resulting in slightly different values of the parameters J 1 2 and J 2 2 [53]. From the viewpoint of the electronic structure, one can construct two types of model: the correlated 5-orbital model, where the one-electron part is taken from LDA and combined with the screened Coulomb interactions obtained in cRPA (as explained in Appendix A), and the all electron model within LSDA, which includes Mn 3d as well as O 2p and A 5d bands. The calculations are performed for the A-type AFM phase on the mesh of the 8×8×6 points, both for k and q, unless it is specified otherwise.

Correlated 5-orbital model
First, we investigate abilities of correlated 5-electron model: whether it can reproduce the main tendencies in the behavior of interatomic exchange interactions in LaMnO 3 and HoMnO 3 , stabilizing the A-type AFM state in the former case and making it unstable with respect to the formation of the spin superstructures along the orthorhombic axis b in the latter one. Another important question is how good is the approximate schemê b, which is frequently used for the analysis of interatomic exchange interactions in these orthorhombic manganites [37,132] Thus, the Coulomb U∼2 eV is not particularly strong. 13 This finding is very important as it naturally explains the existence of the long-range exchange interaction, J b 2 and J 1 3 , arising in the higher orders ofĤ/U beyond the conventional superexchange approximation [2], and responsible for the formation of the spin superstructures along b. On the other hand, the relatively small U also means that the strong-coupling limit is hardly satisfied. Therefore, it is reasonable to expect that the approximate schemeb may experience serious limitations for these orthorhombic manganites.
The densities of states obtained in the Hartree-Fock approximation for the A-type AFM state are displayed in Figure 15. LaMnO 3 and HoMnO 3 have similar electronic structure. The only difference is slightly smaller e g bandwidth in the case of HoMnO 3 . Nevertheless, the behavior exchange interactions appears to be very different. These interactions are summarized in Tables 15 and 16 (lines 5o) for LaMnO 3 and HoMnO 3 , respectively.
The schemeb yields strong interlayer coupling J ⊥ = J ⊥ 1 + 2J 1 2 + 2J 2 2 = −10.78 (−10.1) meV for LaMnO 3 (HoMnO 3 ), which is even stronger than the intralayer one J 1 = 2.36 (−6.53) meV. This behavior is inconsistent with the experimental neutron- Table 15. Isotropic exchange interactions in LaMnO 3 (in meV) as obtained in correlated 5-orbital (5o) model in comparison with LSDA for the all-electron Mn 3d+O 2p+La 5d model. The notations of parameters are explained in Figure 14b,c. The AFM interactions J a 2 and J 2 3 are considerably weaker and not shown here. T N is the Néel temperature evaluated in RPA (in K) and q = (0, q b , 0) is the theoretical ground state propagation vector (in units of reciprocal lattice translations, where q = 0 corresponds to the A-type AFM state). method scattering data for LaMnO 3 , indicating that J ⊥ ∼ −4.7 meV is weaker than J 1 = 6.6 meV [145,146]. 14 Furthermore, the longer-range AFM interactions J b 2 and J 1 3 in LaMnO 3 are comparable or even stronger than J 1 , making the A-type AFM structure unstable.
On the contrary, the scheme M systematically improves the description of the interatomic exchange interactions. First, the interlayer coupling becomes considerably weaker: J ⊥ = −2.35 (−4.34) meV for LaMnO 3 (HoMnO 3 ). Then, the itralayer interaction J 1 in LaMnO 3 becomes strongly ferromagnetic, as expected for the "antiferro" orbital ordering [147], 15 and overcomes the AFM long-range interactions J b 2 and J 1 3 . As the result, the A-type AFM phase becomes stable. The theoretical T N = 186 K, evaluated in RPA, is in fair agreement with the experimental value of 140 K [145,146]. In HoMnO 3 , J 1 is significantly reduced due to the additional buckling of the Mn-O-Mn bonds, and becomes comparable with J b 2 and J 1 3 . This makes the A-type AFM phase unstable. Regarding the direction of this instability, it is very important that J a 2 = −0.35 meV is much weaker than J b 2 = −1.53 meV. Therefore, the propagation vector for the expected theoretical ground state is parallel to b, in agreement with the experimental observation. This instability is further enhanced by the AFM interactions J 1 3 . Nevertheless, in order to reproduce the experimental E phase with commensurate q = (0, 1 2 , 0), it is essential to consider other ingredients such as the exchange striction and single-ion anisotropy, which would lock the spin superstructure to the lattice [148][149][150]. In HoMnO 3 , such lock-in transition occurs at 29 K, while the Néel temperature is about 41 K [136,137], which is close to theoretical value of T N = 53 K. Other aspects of stability of the E-type AFM phase will be considered in Section 5.4.

All-electron model in LSDA
The all-electron model in LSDA provides an alternative description for the exchange interactions. We start with the analysis of effective Stoner parameters. As was discussed in Section 3.3.1, amongst several possible definitions of the parameters I, the most relevant seem to be the definitions III and IV, which should be considered in combination with the methodsb and M, respectively. In the A-type AFM phase, the A and apical O atoms, located between the antiferromagnetically coupled layers, are nonmagnetic. Therefore, we set I = 0 for them. The remaining parameters for Mn and planar O atoms are listed in Table 17. As expected, I Mn is practically identical for LaMnO 3  Tables 15 and 16 (lines  LSDA). Unlike for the correlated 5-orbital model, the methodsb and M in all-electron LSDA provide a consistent description. On the one hand, there is an effect of the O 2p band, which is explicitly treated in this model. On the other hand, the underestimation of the FM interactions in the schemeb is partly compensated by I O . Moreover, the correct definition of the xc fieldb using the sum rule is important. For instance, the use ofb defined via site-diagonal elements ofĤ ↑,↓ underestimates the FM interactions and makes the A-type AFM state in LaMnO 3 unstable.
The exchange interactions are generally stronger than in correlated 5-orbital model (except J 1 ), partly due to the fact that in these insulating materials the exchange interactions are expected to increase when the effective Coulomb repulsion decreases, as in the superexchange mechanism [2]. Furthermore, the additional contribution may arise from the oxygen band [113,151]. Anyway, the total interlayer coupling in Figure 16. Form of DM interactions operating between nearest neighbors in the orthorhombic structure, in and between the plane. The numbering of Mn sublattices is the same as in Figure 14. Each DM vector is attached to its Mn-O-Mn bond, starting from the atoms 2 or 4 for the in-plane interactions and atoms 3 or 4 for the out-of-plane interactions.
LaMnO 3 , J ⊥ = −4.82 (−5.41) in the scheme M (b), is consistent with the experimental data [145,146]. J 1 is strongly ferromagnetic, which appears to be sufficient to overcome strong AFM interactions J b 2 and J 1 3 , and stabilize A-type AFM phase. In HoMnO 3 , this J 1 is substantially weaker, making the A-type AFM state unstable with respect to an incommensurate spin structure propagating along b. The corresponding theoretical T N is in fair agreement with experimental data.

Dzyaloshinskii-Moriya interactions
All nn DM interactions in the orthorhombic planes, d ij , are transformed to each other by the symmetry operations of the space group P bnm [132]. The same holds for the interplane interactions d ⊥ ij . Furthermore, since neighboring planes are connected by the mirror reflection, the z (c) component of d ⊥ ij will be equal to zero. Therefore, there will be 5 parameters describing all nn DM interactions: d a , d b , and d c for the in-plane interactions, and d ⊥ a and d ⊥ b for the out-of-plane interactions. The corresponding DM vectors, attached to neighboring Mn-O-Mn bonds, are illustrated in Figure 16. The parameters obtained in the scheme M are summarized in Table 18.
For LaMnO 3 in LSDA we note a good agreement with the results of ref. [132], employing rather different strategy -the mixed perturbation theory, where the rotation Table 18. Parameters of nearest-neighbor DM interactions (in meV) as obtained in correlated 5-orbital (5o) model and all-electron LSDA using the scheme M . The results of ref. [132], employing mixed perturbation theory, are shown for comparison. of the xc field on one Mn sites was combined with the SO coupling on another such site. This means that the 5d states of the heavy A atoms, which are located in the unoccupied part of the spectrum (Figure 14), do not largely contribute to the DM interactions (though these states were taken into account in our calculations) -quite contrary to our expectations based on the previous analysis for CrI 3 (Section 3.5). Partly, this may be due to the fact that the A sites remain nonmagnetic in the A-type AFM state. The parameters obtained in correlated 5-orbital model are generally smaller than in LSDA -similar to what was found for the isotropic interactions. However, this is not surprising and can be again explained by the Coulomb U in the denominator of superexchange interaction [2].
The DM interactions in AMnO 3 mainly contribute to the spin canting. The magnetocrystalline anisotropy tends to align the spins parallel to the b axis [152][153][154]. In LaMnO 3 , they order according to the A-type: e µ = (0, −1, 0) for µ = 1 and 2 and e µ = (0, 1, 0) for µ = 3 and 4 in Figure 16. This perfect AFM alignment will be further deformed by the DM interactions, which additionally rotate the spins along a and c by, respectively, e a = −d c /(J 1 − J 1 2 − J 2 2 ) and e c = d ⊥ a /(J 1 + 2J 1 2 + 2J 2 2 ) [54]. The a components will order according to the G-type, 16 while c components gives rise to the weak ferromagnetism [152,153]. Using the obtained parameters of exchange interactions, |e a | and |e c | can be estimated as 0.019 and 0.110, respectively. Taking into account that in all-electron LSDA M = 3.59 µ B , the weak FM moment can be estimated as 0.4 µ B (being somewhat larger that the experimental 0.18 µ B , derived from magnetization measurements [155]).
An interesting question is whether the multiferroicity associated with the E-type AFM phase in HoMnO 3 can coexist with the weak ferromagnetism. In the E-type AFM structure, each of the J 1 2 and J 2 2 will contribute to two FM bonds and two AFM ones. Therefore, these contributions cancel each other and the spin canting along c will be given by e c = ±d ⊥ a /J 1 (for the sublattices with different spins in the ab plane). Using the parameters for the P bnm structure of HoMnO 3 , |e c | can be estimated as 0.169 and 0.050 in the all-electron LSDA and correlated 5-orbital model, respectively, where the difference mainly comes from d ⊥ a . Nevertheless, this component will order according to the AFM E-type and, therefore, there will be no weak ferromagnetism.

Internal instability of the E phase
An interesting aspect of interatomic exchange interactions defined via infinitesimal rotations of spins near some equilibrium magnetic states is that these exchange interactions depend on the magnetic state in which they are calculated, reflecting the dependence of the electronic structure on the magnetic state. If it happens, the bilinear spin model (1) is ill-defined in the global sense, meaning, for instance, that the same set of the exchange parameters cannot be used for the analysis of the low-temperature spinwave dispersion and the magnetic transition temperature, where the electronic structure is strongly modified by the spin disorder [113,114]. Nevertheless, the model (1) can be defined locally, near each magnetic equilibrium. Then, in principle, one can expect some kine of self-organization phenomenon, when the change of the electronic structure in certain magnetic state itself can stabilize this magnetic state, at least locally. This is what happens at least with some of the exchange interactions in the 50% doped manganites, where the zigzag (CE-type) AFM alignment opens the band gap, induces the orbital ordering, and additionally stabilizes the FM coupling in the zigzag chains [87]. Can we expect a similar behavior in the E phase?
Intuitively, the reason for the formation of this noncentrosymmetric E-type AFM structure can be understood as follow: the AFM interactions J b 2 and J 1 3 tend to align the corner spins, separated by the orthorhombic translations a and b, antiferromagnetically, as shown in Figure 17a. The central Mn atom is located in the inversion center. Then, due to the AFM alignment, the corner atoms should transform to each other by the symmetry operationÎT (where the spacial inversionÎ is combined with the time reversalT ), which is formally compatible with the P bnm space group. However, the same symmetry operation would make the central Mn atom nonmagnetic. This is an eligible scenario from the symmetry point of view, but would lead to gigantic energy loss, Mn ∼ 4 eV, caused by the violation of the first Hund's rule for the atoms, which are potentially expected to be in the high-spin state. Therefore, the more favorable scenario is to keep the Mn atom magnetic, but break the inversion symmetry, which is quite common in multiferroic materials, especially those involving the high-spin ions.
Then, the next question is whether the E-type AFM order is compatible with the change of the exchange interactions induced by this order.
Below we consider results of correlated 5-orbital model (in the scheme M) for HoMnO 3 . A qualitatively similar behavior has been found in LSDA for the all-electron model. The E-type AFM order results in the additional narrowing of the t 2g and e g bands (Figure 17), which is consistent with the decrease of the number of FM bonds around each Mn site (2 instead of 4 in the A-type AFM phase). The inversion symmetry breaking makes the nn interactions inequivalent, where J ↑↑ 1 in the FM bond is generally different from J ↑↓ 1 in the AFM one. In order to stabilize the E state, one would need at least J ↑↑ 1 > J ↑↓ 1 (the FM coupling is stronger in the FM bond). However, we have found the opposite tendency: J ↑↑ 1 = 2.99 meV and J ↑↓ 1 = 4.00 meV, meaning that the E state corresponds to the energy maximum and any small rotations of spins will slide the system away from this equilibrium to a new state. Furthermore, there is an intrinsic mechanism inducing these rotations of spins. Indeed, the E-type AFM order, producing some changes in the electronic structure, can induce the electric polarization even in the centrosymmetric P bnm structure [148]. Then, it is reasonable to expect that the same changes in the electronic structure may lead to the appearance of new DM interactions, which would otherwise be forbidden in the P bnm structure. Particularly, we have found finite interactions d b 2 = (±0.01, 0, ±0.02) meV and d 1 3 = (±0.01, ±0.01, ±0.02) meV (being analogs of isotropic J b 2 and J 1 3 , respectively, where the ± signs depend on the origin and direction of the bond). Although these interactions are not particularly strong, they are sufficient to shift the system of spins from the extremum (maximum) point, so that it will start to relax to a new magnetic equilibrium.
Thus, the E-type AFM state in HoMnO 3 is intrinsically unstable and cannot be stabilized by purely electronic mechanisms. For these purposes, it is essential to consider the exchange striction or single-ion anisotropy (or both of them) [148][149][150].

Brief summary
• The minimal correlated 5-orbital model provides a consistent description for interatomic exchange interactions in LaMnO 3 and HoMnO 3 , explaining stability of the A-type AFM phase in the former material and the tendency towards formation of spin superstructures along the orthorhombic axis b in the case of HoMnO 3 . The driving force behind this change of the magnetic structure is the buckling of the Mn-O-Mn bonds, which is stronger in HoMnO 3 . The use of the scheme M appears to be crucial within this 5-orbital model, while the approximate schemeb strongly underestimates the FM interactions.
• The all-electron model, explicitly treating the Mn 3d, O 2p, and A 5d bands in LSDA provides an alternative description for the exchange interactions. Both 5-orbital model and all-electron LSDA capture the experimental situation pretty well. Why do we have such consistent description? Apparently, this is due to the cancellation of many contributions. For instance, the polarization of the oxygen states in LSDA strengthens the FM interactions [37]. In addition to them, there are FM superexchange interactions, which are controlled by the ratio of the on-site exchange and Coulomb repulsion, J/U, and expected to be stronger for smaller U [147]. However, these effects are compensated by stronger AFM interactions, also expected in the superexchange theory for smaller U [2]. Furthermore, the AFM interactions are additionally stabilized by correlations effects beyond the Hartree-Fock approximation [37].
• Unlike in the correlated 5-orbital model, the schemesb and M provide a consistent description within all-electron LSDA. Nevertheless, the special attention should be paid to the definition of the xc field and the choice of the parameters I O . The incorrect choice of these parameters can worsens the description.
• The inversion symmetry breaking, caused by the E-type AFM order in HoMnO 3 , gives rise to new DM interactions, operating across the inversion centers in the P bnm structure, which, in the combination with the magnetic-state dependence of the isotropic interactions, act to destabilize this E-type AFM state.

Symmetric anisotropic exchange interactions
The spin-spiral concept, which assumes that the perturbation of the magnetic ground state can be described in the form of an incommensurate spin spiral, is applicable only for the calculations of isotropic and DM interactions [14]. All these calculations are based on the generalized Bloch theorem, which allows us to deal with this spinspiral periodicity [41]. Unfortunately, the theorem cannot be used for calculations of the exchange anisotropy or any other symmetric anisotropic interaction, emerging in the second order of the SO coupling and involving spin diagonal as well as off-diagonal elements of this coupling. As was pointed out in Section 2.2, the DM interactions support the spin-spiral propagation, while the symmetric anisotropic interactions act against it. Nevertheless, one can still consider the perturbation caused by the infinitesimal rotations of the xc field and evaluate the 3×3 exchange tensorĴ ij separately for each magnetic bond. Then, 1 2 (Ĵ ij −Ĵ ji ), which has only three inequivalent matrix elements, can be related to the DM vector, while the traceless part of 1 2 (Ĵ ij +Ĵ ji ) gives rise to the exchange anisotropy. This method was successfully implemented in several computational packages, which are actively used for calculations of isotropic as well anisotropic exchange interactions [27,[44][45][46][47]. Nevertheless, one should remember that this method has limitations inherent to the schemeb, which is an approximation. Unfortunately, only the schemeb can be easily reformulated in the real space, where the exchange interactions can be calculated separately for each bond. Due to the additional inversion of the response matrix in the scheme M, the exchange parameters in the bond will depend on the matrix elements of the response function in other bonds.

Dynamic electron correlations
As was already pointed out in Section 4.1, another interesting direction is to go beyond the conventional SDFT by incorporating the effects of dynamic electron correlations. The latter are typically evaluated in the framework of DMFT [116], which provides a formal extension of the KS equations, where the local potential is replaced by also local, but frequency-dependent self-energy [117]. Then, the exchange interactions can be still evaluated using the schemeb, but with the frequency-dependent xc field [118]. This approach is especially important for metallic systems, where the static Hartree-Fock approximation is clearly insufficient and, as long as the Coulomb interactions are taken into account, they should be treated on a more rigorous footing. For instance, the dynamic correlations have a profound effect on the electronic structure of half-metallic compounds [115], which is reflected in the behavior of exchange interactions, as was demonstrated for CrO 2 [108,156]. However, as the most interesting applications of this method are beyond the strong-coupling limit, the schemeb can have serious limitations. Therefore, it is very important to reformulate the interatomic exchange interactions in terms of the inverse response function, as it is done in the scheme M. In this respect, Equation (10) seems to be very general and could be a good starting point for such extension.
Another advantage of this DMFT based approach is that it provides a natural extension for the analysis of temperature dependence of the exchange interactions [112].

Summary and conclusions
The linear response theory becomes a powerful tool for calculations of interatomic exchange interactions in various substances. By treating infinitesimal rotations of spins as a perturbation, it allows us to present the total energy change caused by these rotations in the form of pairwise interactions. The main purpose of this topical review was to clarify basic principles of this technique and provide a transparent explanation to more recent developments and controversies related to its practical realization.
The first group of questions is related to the validity of the magnetic force theorem for the infinitesimal rotations of spins [29,58]. The magnetic force theorem is certainly valid and the total energy change underlying calculations of the exchange interactions can be replaced by the change of the single-particle energies, which seems to be a general fundamental property. However, it would be a mistake to make the equality between the magnetic force theorem and Equation (3). Such attempts are obviously misleading as the correct use of the magnetic force theorem for the exchange interactions should also include the contributions of the external magnetic field, which is needed to control the direction of the magnetization. Equation (3) does not take into account such contribution. This is an approximation leading to the linear dependence of the exchange interactions on the response tensor, which can be justified only in the limits of long wavelength and strong coupling. In fact, the exact expression (10) for the energy change is extremely simple. We only need to know how to find the external magnetic field and for these purposes we need the linear response theory. The total energy change (and interatomic exchange interactions) in this case will be proportional to the inverse response tensor. The procedure is applicable for calculations of the isotropic exchange as well as the DM interactions. In the latter case, it requires the additional constraining field in order to compensate the rotations caused by the DM interactions and we have shown how this field can be found using the linear response theory.
Besides fundamental aspects of the magnetic force theorem, there is a number of more practically oriented questions. The first one is which object is more suitable for the description of infinitesimal rotations of spins. In this respect, we have considered rigid rotations of the magnetization matrix and local magnetic moments (the schemesm and M, respectively). Since the local magnetic moment is nothing but the spherical part of the magnetization matrix and only this spherical part is subjected to the constraint in the scheme M (while other components of the magnetization matrix are allowed to relax) such perturbations are expected to cost less energy and, therefore, should be more suitable for the analysis of low-energy excitations.
Another important question is what to do with the ligand states. In most of the applications, these states play a role of effective medium participating in the electron transfer between magnetic transition-metal sites. The exchange interactions are typically calculated only between the transition-metal sites, while the contributions of the ligand sites, which can carry an appreciable portion of the magnetization due to the hybridization with the transition-metal sites, are ignores. In this respect, we have proposed the downfolding method [53], which allows us to eliminate the ligand states, by transferring their effect to the exchange interactions between the transition-metal sites. This can be achieved by employing the adiabaticity concept and assuming that the fast ligand states instantaneously follow the slow change of the magnetization on the transition-metal sites.
An important aspect of the downfolding method is that it can naturally incorporate the dependence of the exchange interactions on the strength of the Stoner coupling I L on the ligand sites, which is regarded to be the key ingredient of phenomenological GKA rules for the 90 • -exchange and in certain cases primarily responsible for the FM character of this exchange. Nevertheless, the covalent mixing can easily make the effective coupling I L negative. Such situation is realized, for instance, in CrCl 3 , CrI 3 , and CrO 2 , where the magnetic polarization of the ligand states is antiparallel to the one of the transitionmetal states. In such case, the negative I L acts against the FM coupling and the ferromagnetism is stabilized by other mechanisms involving the Cr e g states.
For most of the applications, we have considered two pictures, which provide a supplementary description for the exchange interactions. The first one is the correlated minimal model constructed for the magnetic, mainly transition-metal, bands near the Fermi level. The second one is based on all-electron LSDA, which take into consideration both magnetic transition-metal and ligand states. Once the Coulomb interactions are added to the TB model, they should be treated rigorously: if static Hartree-Fock approximation is applicable for the analysis of exchange interactions in insulating systems with lifted orbital degeneracy, the situation in metallic systems can be very different. In certain cases, the use of the Hartree-Fock approximation can lead to a misleading answer, as was demonstrated for half-metallic CrO 2 [108]. However, this is not the only source of the error for the exchange interactions. Another problem is related to the fact that most of the real materials are far from the strong-coupling limit and the commonly used Equation (3) for the exchange interactions is no longer justified. A more rigorous approach is to evaluate the exchange interactions via the inverse response (the scheme M), which is certainly more difficult from the computational point of view. However, such scheme tremendously improves the description of interatomic exchange interactions in the correlated model, as was demonstrated for insulating CrX 3 (X= Cl, I) and AMnO 3 (A= La, Ho) in the framework of the Hartree-Fock approximation. The application of correlated models for the exchange interactions in metallic systems should consider two aspects: (i) The method should be beyond the Hartree-Fock approximation. The commonly used alternative is DMFT [116,117]; (ii) The calculations of the exchange interactions themselves should be based on the inverse response, as in the scheme M.
The exchange interactions in all-electron LSDA (GGA) obey quite different principles. Certainly, this is a simple approximation derived in the limit of homogeneous electron gas. Nevertheless, it satisfies certain fundamental sum rules and on many occasions provides an insightful view on the properties of systems far beyond this limit [157]. Thus, it can be regarded as a supplementary view in addition to the correlated model, though not necessarily the perfect one. For the magnetic systems, LSDA and strong-coupling limit are basically incompatible with each other (or this limit can be strongly underestimated in LSDA). Therefore, there is absolutely no guarantee that the approximate schemeb can properly capture the behavior of interatomic exchange interactions and a more rigorous scheme M formulated in terms of the inverse response looks more preferable. Nevertheless, the improvement expected in the scheme M is somewhat diminished by the effects of other parameters controlling the properties of the exchange interactions in the all-electron case, particularly the choice of the xc field b and effective Stoner coupling I L on the ligand sites. By taking the xc field from the sum rule, b =Q + 0 m, one can substantially improve the description in the schemeb (even despite approximate form in terms of the response function). However, the choice of the parameters I L , which can be strongly depend on the definition, poses a more serious problem. We have considered several such definitions. In a number of cases, they provide a consistent description for the interatomic exchange interactions, but not always. It seems that there is still an ambiguity with the choice of I L and the isotropic exchange interactions depend on such ambiguity. On the other hand, the dependence of DM interactions on I L is considerably weaker. As an alternative approach, we have proposed the scheme L0, which assumes that the magnetic moments on the ligand sites are induced solely by the hybridization with the magnetic transition-metal sites, while all Stoner parameters I L can be set to zero. Such scheme is also formulated in terms of the linear response and most suitable for the analysis interatomic exchange interactions in FM insulators and half-metals.
The merging of correlated model for the transition-metal bands with LDA electronic structure for the ligand bands is another largely unresolved problem. Although such merging is expected to improve the description of the exchange interactions, by combining the effects of Coulomb correlations in transition-metal bands with the explicit treatment of the ligand states, the strategy suffers from many ambiguities, as was demonstrated for CrCl 3 and CrI 3 . We are always forced to compromise between genuine physical effect and intrinsic error of such merging caused by additional approximations and assumptions, which are inevitable in this case. On the other hand, on many occasions such merging is not really needed as the physically meaningful picture for the exchange interactions can be obtained in correlated 5-orbital model constructed only for the magnetic 3d bands. As a supplementary step, one can always consider the picture of all-electron LSDA (GGA), which gives an idea about the role of the ligand states.
Amongst other interesting properties, CrI 3 has attracted a considerable attention due to the strong DM interaction induced by the large SO coupling on the heavy I atoms [70]. On the other hand, CrI 3 is a closed shell material, where orbital degrees of freedom are quenched by the large crystal-field splitting between the occupied t 2g and unoccupied e g states. Therefore, from the practical point of view, comparable or even stronger DM interactions can be expected in the open shell materials even without heavy 5p elements, as was demonstrated, for instance, for CrO 2 , LaMnO 3 , and HoMnO 3 .
corresponding KS Hamiltonian in the Wannier basis. In LSDA, these parameters depend on spin. Furthermore,Ĥ σ = [H σ ia,jb ] can include spin-diagonal part of the SO coupling, which is needed to calculate DM interactions [14,56]. In this case,Ĥ σ also depends on spin, even in the non-magnetic LDA, where it holds:Ĥ ↓ =Ĥ ↑ * . For practical purposes we use mainly the LMTO method [79,80]. Then, the Wannier functions can be generated using the projector-operator technique [38,78]. The spin-diagonal part of the SO coupling is typically added to the site-diagonal part of the LMTO Hamiltonian asĤ ↑,↓ ii →Ĥ ↑,↓ ii ± 1 2 ξ iL z i . After that, the TB Hamiltonian is constructed in the Wannier basis for the target bands. Thus, even though the target bands are formally of the transition-metal 3d type, the corresponding Wannier functions and the TB Hamiltonian include some contributions of the SO coupling of the heavy ligand atoms (if any), which are admixed into the target bands via the hybridization effects.
The parameters of screened on-site Coulomb interactions,Û = [U abcd i ], are evaluated in cRPA starting with the LDA band structure [82]. The basic idea of the constraint in this case is to get rid of nonphysical metallic screening in LDA emerging from the bands near the Fermi level. For the d electrons, the matrixÛ 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), where U is responsible for the charge stability of given electronic configuration, while J and B are responsible for the first and second Hund's rules, respectively [78]. Strictly speaking, the parametrization in terms of U, J, and B is valid only for the isolated atoms in the spherical environment [159]. It is used only for the explanatory purposes, while all numerical calculations are performed with the matrices of Coulomb interactionsÛ extracted from cRPA without additional fitting.
After the construction, the model is solved in the mean-field Hartree-Fock approximation, where the second term in Equation (A.1) is replaced by and the potentialV σ i = [V σ i,ab ] is found self-consistently. Then, the obtained electronic structure is used for calculating the exchange interactions. The xc field in this case is defined asb i =V ↑ i −V ↓ i . When the SO interaction is added toĤ σ , the potentialV σ i is recalculated self-consistently to include the effects or the SO coupling. After that the obtained electronic structure andV σ i are used to calculate the DM interactions. The Hartree-Fock approximation is believed to be a good starting points for the analysis of magnetic properties of insulating materials, where the orbital degeneracy is lifted by the lattice distortions. Other details can be found in ref. [78].
Appendix B. Random-phase approximation for the magnetic transition temperature Here, we generalize the RPA expression for the critical transition temperature, T S , of compounds with multiple magnetic sublattices to the case, where the ground state is a noncollinear spin spiral with the propagation vector q. Moreover, the sublattices are allowed to acquire the additional phases γ q, µ , describing relative rotations of the magnetization in different sublattices relative to each other. In all other respects, we follow the derivation considered for the collinear case by Rusz et al [160].
Our first goal is to find the ground state, corresponding to the minimum of the energy where c q, µ = e iγq, µ and J q, µν is the Fourier transform of J ij between the sublattices µ and ν. Equation (B.1) can be generalized to include the DM interactions d z ij by replacing J q, µν with J q, µν − id z q, µν . For each q, we minimize E(q) with respect to γ q, µ using the gradient descent method. Then, we pick up q corresponding to the global minimum of E(q) and for each k redefine J k, µν as J k, µν → c * q, µ J k, µν c q, ν with the coefficient c q, ν determined for the ground state. This is nothing but the transformation to the new local coordinate frame corresponding to the energy minimum.
Then, the spin-wave frequencies can be associated with the eigenvalues of the positive-defined matrixΩ k =Ŝ −1N k , where N k, µν = δ µν ν ′ J q, µν ′ − J k, µν , (B.2) S =M /2, andM is the diagonal matrix of magnetic moments. In RPA, the corresponding magnetic transition temperature T µ S for the sublattice µ is given by [160]: (Ω BZ being the volume of the first Brillouin zone). For the inequivalent sublattices, this T µ S depends on the ratio of the magnetic moments, M µ /M ν , which can also depend on the temperature. In order to find true transition temperature, these rations should be adjusted to make the same T µ S = T ν S ≡ T S for all the sublattices [160]. The spectral theorem is employed in order to calculateN −1 k . Then, a small imaginary part, iδ with δ = 0.01 meV (corresponding to δ/k B ∼ 0.1 K) is added to the eigenvalues ofN k for k close to q. Finally, the k-space integration is replaced by the summation on a very dense mesh of k-points (for instance, a typical mesh used for CrCl 3 and CrI 3 is 258×258×258).