# Fileset

[Solovyev_2024_J._Phys.__Condens._Matter_36_223001.pdf](https://mdr.nims.go.jp/filesets/355e1d31-e744-41b2-a2ad-153eae417728/download)

## Creator

[I V Solovyev](https://orcid.org/0000-0002-2010-9877)

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Linear response theories for interatomic exchange interactions](https://mdr.nims.go.jp/datasets/4cbdbbc4-8b9b-4450-bb9d-4d5a0ee883df)

## Fulltext

Linear response theories for interatomic exchange interactionsJournal of Physics: CondensedMatter     TOPICAL REVIEW • OPEN ACCESSLinear response theories for interatomic exchangeinteractionsTo cite this article: I V Solovyev 2024 J. Phys.: Condens. Matter 36 223001 View the article online for updates and enhancements.You may also likeC-type antiferromagnetic structure oftopological semimetal CaMnSb$_2$Bo Li, Xu-Tao Zeng, Qianhui Xu et al.-Observation of $\Omega(2012) \to\Xi(1530)\bar{K}$ and measurement of theeffective couplings of $\Omega(2012)$ to$\Xi(1530)\bar{K}$ and $\Xi\bar{K}$Chengping Shen-Sodium-based di-chalcogenide: Apromising material for tandem solar cellsDanilo Gómez-Ríos, Santiago Perez-Walton, Francisco López-Giraldo et al.-This content was downloaded from IP address 144.213.253.16 on 07/03/2024 at 23:43https://doi.org/10.1088/1361-648X/ad215a/article/10.1088/0256-307X/41/3/037104/article/10.1088/0256-307X/41/3/037104/article/10.1088/1674-1137/aca0de/article/10.1088/1674-1137/aca0de/article/10.1088/1674-1137/aca0de/article/10.1088/1674-1137/aca0de/article/10.1088/2516-1075/ad2f5b/article/10.1088/2516-1075/ad2f5bJournal of Physics: Condensed MatterJ. Phys.: Condens. Matter 36 (2024) 223001 (36pp) https://doi.org/10.1088/1361-648X/ad215aTopical ReviewLinear response theories for interatomicexchange interactionsI V SolovyevResearch Center for Materials Nanoarchitectonics (MANA), National Institute for Materials Science(NIMS), 1-1 Namiki, Tsukuba, Ibaraki 305-0044, JapanE-mail: SOLOVYEV.Igor@nims.go.jpReceived 8 August 2023, revised 14 December 2023Accepted for publication 22 January 2024Published 7 March 2024AbstractThe linear response is a perturbation theory establishing the relationship between given physicalvariable and the external field inducing this variable. A well-known example of the linearresponse theory in magnetism is the susceptibility relating the magnetization with the magneticfield. In 1987, Liechtenstein et al came up with the idea to formulate the problem of interatomicexchange interactions, which would describe the energy change caused by the infinitesimalrotations of spins, in terms of this susceptibility. The formulation appears to be very genericand, 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 theelectronic structure obtained, for instance, in the first-principles calculations based on thedensity functional theory. The purpose of this review is to elaborate basic ideas of the linearresponse theories for the exchange interactions as well as more recent developments. The specialattention is paid to the approximations underlying the original method of Liechtenstein et al incomparison with its more recent and more rigorous extensions, the roles of the on-site Coulombinteractions and the ligand states, and calculations of antisymmetric Dzyaloshinskii–Moriyainteractions, which can be performed alongside with the isotropic exchange, within onecomputational scheme. The abilities of the linear response theories as well as many theoreticalnuances, which may arise in the analysis of interatomic exchange interactions, are illustrated onmagnetic van der Walls materials CrX3 (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.Keywords: electronic structure, linear response theory, exchange interactions,Dzyaloshinskii–Moriya interactions, ligand states,transition-metal oxides and related compoundsOriginal Content from this work may be used under theterms of the Creative Commons Attribution 4.0 licence. Anyfurther distribution of this work must maintain attribution to the author(s) andthe title of the work, journal citation and DOI.1 © 2024 The Author(s). Published by IOP Publishing Ltdhttps://doi.org/10.1088/1361-648X/ad215ahttps://orcid.org/0000-0002-2010-9877mailto:SOLOVYEV.Igor@nims.go.jphttp://crossmark.crossref.org/dialog/?doi=10.1088/1361-648X/ad215a&domain=pdf&date_stamp=2024-3-7https://creativecommons.org/licenses/by/4.0/J. Phys.: Condens. Matter 36 (2024) 223001 Topical Review1. IntroductionOn many occasions our image of magnetism rests on thepicture interacting spins (ei and ej) attached to the atomic sites(i and j). If the system is isotropic, such interactions have aform of the scalar products ei · ej [1–5]. The interactions areferromagnetic (FM) if they force ei · ej > 0 and antiferromag-netic (AFM) if ei · ej < 0. If i and j are no longer connectedby the spacial inversion, the spins tend to align neither fer-romagnetically nor antiferromagnetically.1 The correspond-ing interaction, which is called Dzyaloshinskii–Moriya (DM)interaction, is given by the cross product [ei× ej] and drivenby the relativistic spin–orbit (SO) coupling [6, 7]. Thus, itis always nice to have a transparent toy model, which wouldexplain that certainmaterial has a particular magnetic structurebecause some interactions are strong or weak, FM or AFM,etc. The experimental inelastic neutron scattering data are typ-ically fitted to extract parameters of such physically meaning-ful model. In theory, the spin model can be constructed byaveraging the energies of interatomic interactions over non-magnetic degrees of freedom [1–5, 7]. In this article we willexplain how the interatomic exchange interactions can be gen-erally derived starting from the electronic structure obtained inthe first-principles calculations.Let us consider the simplest possible spin model with theenergyE=−12∑i ̸=j(Jijei · ej− dij · [ei× ej]) , (1)where Jij is the isotropic exchange, dij = (dxij,dyij,dzij) is DMvector, and the spin moments ei are normalized to the unity:|ei|= 1. Although we will be primarily interested in the beha-vior of isotropic interactions, it appears to be possible to con-sider Jij in the combination with dij within one computa-tional scheme. The reason will become clear in a moment. Ourgoal is to find parameters of this model using the informa-tion about the electronic structure. Of course, the model (1)is an approximation as there is no reason why the energyof a general magnetic system should have such a simpleform and be described exclusively by the bilinear interac-tions. There are only few microscopic mechanisms, which areconsistent with the form of equation (1). These are the dir-ect 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 energy, related to the antisym-metry of fermionic wave functions. In the last two examples,this is the consequence of the 2nd order perturbation the-ory with respect to, respectively, transfer integrals and intra-atomic exchange interactions, that couples localized core spinsto the outer conduction electrons.1 The collinear FM or AFM alignment is the consequence of the spacial inver-sion in the bond. In the former case, the spacial inversion should be among thesymmetry operations. In the latter case, it is combined with the time reversal.Therefore, if the inversion symmetry is broken, neither FM nor AFM align-ment satisfies the symmetry properties.Nevertheless, there is one more, very special case, wherethe magnetic energy can be also described by equation (1).These are the infinitesimal rotations of spins near the equilib-rium, as was realized by Liechtenstein et al [10, 11]. Indeed,considering rotations ei = (θ cosqRi,θ sinqRi,1− θ22 ) near theground state eGS = (0,0,1) (Ri being the position of the sitei, q being the spin-spiral propagation vector), one can evalu-ate the energy change (per one unit cell) caused by interac-tions between the transversal (xy) components of spins. Forthe model (1), this energy change is given byδEq =−12(Jq− idzq)θ2, (2)where Xq =∑jX0j exp(−iq ·Rj) is the Fourier image of Xij.Then, the basic idea is to extract the same energy changefrom the electronic structure calculations, typically withinspin-density functional theory (SDFT) [12–14], and map it onequation (2). This should give us the parameters of exchangeinteractions Jq and idzq (or Jij and dzij after the Fourier trans-form to the real space). Since the rotations of ei are chosenin the form of the conical spin spiral (which is compatiblewith the DM interactions), Jij can be considered in the com-bination with dzij in the one computational scheme, where theenergy change is uniquely specified by q [15]. Basically, thisis a perturbation theory, which can be formulated in terms ofthe response function (or the susceptibility).One of the most attractive points of the infinitesimalrotations of spins is that the bilinear form of equation (1)remains valid irrespectively on which microscopic mechan-ism stands behind the interaction parameters. It can be thesuperexchange [2], RKKY [3–5], double exchange [16] orany other mechanism, provided that the rotations are small.Even biquadratic exchange [17] for small θ can be reformu-lated in the bilinear form (1). Without SO coupling, the energychange caused by the infinitesimal rotations of spins can bealways described by the Heisenberg model and in this sensethe method is very universal.Another important point of the work of Liechtensteinet al [11] is that they have proposed a practical scheme forcalculating the exchange parameters and proved for these pur-poses the magnetic force theorem [18–20], which justifies theuse of the single-particle energies, obtained from the Kohn–Sham (KS) equations in SDFT [13, 14], for evaluating theenergy change caused by the infinitesimal rotations of spins.The theorem greatly simplifies the calculations and improvesthe numerical accuracy.The basic variable of SDFT is the magnetization dens-ity m̂. In the ground state, m̂ is controlled by the exchange–correlation (xc) field b̂. Therefore, instead of rotating m̂,Liechtenstein et al [11] have proposed to rotate b̂, assumingthat m̂ will automatically rotate by the same angle. This leadsto the commonly used expression for Jij:Jij =12πImˆ εF−∞dεTrL{Ĝ↑ij (ε) b̂j Ĝ↓ji (ε) b̂i}, (3)which is nothing but the 2nd order perturbation theoryfor the single-particle energy, formulated in terms of the2J. Phys.: Condens. Matter 36 (2024) 223001 Topical Reviewsingle-particle Green’s functions Ĝσij with the spins σ = ↑ or ↓(εF being the Fermi energy, TrL stands for the trace over orbitalindices, and Im denotes the imaginary part). Although this res-ult was anticipated by the previous works on the RKKY inter-actions [5, 8], the exchange interactions in the paramagneticmedium [21, 22], as well as general theories of the itinerantmagnetism [23, 24], the expression (3) can be relatively easilycombined with the first-principles electronic structure calcu-lations for the ground state, in the framework of SDFT or itsrefinements. Today, it is known for almost 40 years and wassuccessfully applied for the analysis of interatomic magneticinteractions in various substances [25–32].Nevertheless, there are also open questions. Particularly,several authors have raised doubts that the true energy changecaused by the infinitesimal rotations of spins can be describedby rotating only b̂ and suggested that it should include theadditional contribution steaming from the external magneticfield, which is needed to control the direction of the mag-netization [33–35]. Thus, equation (3) may be incomplete.Presumably, the most persuasive arguments were given in2003 by Bruno [34], who proposed how equation (3) should becorrected. Surprisingly, however, that even 20 years later afterthis publication there is no systematic analysis of the prob-lem: equation (3) is widely used, but little is known how goodit is. The problem is complicated by rather common misunder-standing putting the equality between equation (3) and morefundamental magnetic force theorem.Another question 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 withthe 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 magneticsite, one can define the magnetization matrix and rotate it. Forinstance, equation (3) implies such matrix form. Nevertheless,which construction is more suitable, based on rotations of thescalar moments or the magnetization matrices, is absolutelyunclear.Then, what shall we do with the ligand sites, whichcan hardly be the source of the magnetism, but frequentlycarry an appreciable magnetization due to the hybridizationwith the magnetic transition-metal sites? The contributionsof such ligand states are typically ignored, and the interac-tions (3) are computed only between the transition-metal siteswithout the justification. On the other hand, there are well-known Goodenough–Kanamori–Anderson (GKA) rules [36–39], which state, among others, that in certain circumstancesthe exchange interactions between the transition-metal sitescan be controlled by the effective Stoner coupling on the inter-mediate ligand sites. Such coupling is typically added empiric-ally to correct Jij given by equation (3) between the transition-metal sites [40–42]. However, if the theory is general enough,it should include all such contributions automatically.The aim of this article is to review some basic ideasas well as more recent developments related to the use ofthe linear response methods for the analysis of interatomicexchange interactions and give clear answers to all above ques-tions. The general theory is discussed in section 2. Then,sections 3–5 deal with practical examples for several typesof compounds, where our main goal is to explain the theor-etical nuances, which may arise in various parts of calcula-tions of the interatomic exchange interactions: the applicab-ility of equation (3) and its refinements, the role of on-siteCoulomb correlations, the contributions of the ligand sites,the merging of correlated and uncorrelated bands, etc. Shortsection 6 outline other developments beyond the main scopesof this review. The article is summarized in section 7. Twoappendices deal with the construction of the tight-binding(TB) Hamiltonians using the band structure obtained in first-principles calculations and the evaluation of magnetic trans-ition temperature for the spin model in the random phaseapproximation (RPA).2. Infinitesimal spin rotations and exchangeinteractions2.1. Basic idea, notations, and conventionsIn the magnetic equilibrium, the 1st derivative of the totalenergy with respect to a small change of the magnetizationm(r) vanishes and the energy changes is described by the2nd derivative. In a general sense, the magnetic force theoremstates that not only the total energy but also its 2nd derivat-ive with respect to the infinitesimal rotations of the magnet-ization is the ground-state property as it can be expressed viathe eigenvalues and eigenfunctions of the ground state. Thisstatement can be traced back to fundamentals of the quantummechanics, where the system can be measured only via per-turbations. Therefore, it is logical that the 2nd derivative can beconnected to the properties of the ground state. The responsefunction (or the susceptibility) is the useful tool, which estab-lishes such connection by means of the perturbation theory.In practical terms, we will deal mainly with SDFT [12–14], where the ground-state magnetization density and thetotal energy are described with the help of the single-particlespin-dependent KS Hamiltonian Hσ(r) with some local self-consistent potential incorporating all effects of exchange andcorrelations [13, 14]. This locality implies that the change ofthe potential in certain point r depends only on the changeof magnetization in the same point r. As a consequence, theenergy change caused by the infinitesimal rotations of themagnetization can be presented in the form of pairwise inter-actions. Without SO coupling it corresponds to the isotropicbilinear Heisenberg model. This is a general property of the2nd order perturbation theory with the local potentials.From the viewpoint of analysis and interpretation, it ismore convenient to adopt the TB representation, which dealswith the atomically resolved properties emerging from thesolution of some lattice model. Another advantage of theTB representation is the on-site Coulomb interactions, whichcan be easily incorporated into the model. The purpose ofthese interactions is to correct limitations of the local dens-ity approximation (LDA) or the generalized gradient approx-imation (GGA), which are derived in the limit of homogen-eous electron gas and typically used to describe the effectsof exchange and correlations in SDFT. Mathematically, this3J. Phys.: Condens. Matter 36 (2024) 223001 Topical Reviewcan be done by constructing the orthonormal basis of local-ized Wannier functions centered on the atomic sites [43]. TheKS Hamiltonian in this basis is the matrix specified by thelattice (i, j) and orbital (a, b) indices, Ĥσ ≡ [Hσia,jb], so thatall other matrices can be obtained from Ĥσ . For instance, theGreen function is Ĝσ(ε) = [ε− Ĥσ]−1 and the magnitude ofthe magnetization is given by m̂=Θ(εF − Ĥ↑)−Θ(εF − Ĥ↓),in terms of the Heaviside function Θ. We assume that the xcpotential in this TB representation remains local in the sensethat on each atomic site i it depends only on the magnetizationm̂i on the same site. It is a common practice to relate the mag-netic part of such potential, the so-called xc field b̂, with thesite-diagonal elements of the TB Hamiltonian: b̂i = Ĥ↑ii − Ĥ↓ii.However, such b̂ is ill-defined because it ignores the non-localcontributions steaming from the off-diagonal part of Ĥσij [30].A more consistent definition of the local b̂ in terms of theresponse function and the ground-state magnetization will begiven in section 2.6.Other conventions can be formulated as follows:• For periodic systems, [Ĥσia,jb] can be Fourier transformedto [Hσµa,νb(k)], with µ and ν denoting the atomic posi-tions within the primitive cell. Furthermore, the analysisthroughout this paper assumes the use of the periodic gaugeHσµa,νb(k+G) = Hσµa,νb(k) for any reciprocal lattice trans-lation G [44].• The n× nmatrix âµ, specified by the n orbital indices, can beviewed as the column vector a⃗µ of the length n2. The scalarproduct a⃗†µ · b⃗µ is the shorthand notation for TrL{âµb̂µ}(where a⃗†µ is the row vector corresponding to the columnvector a⃗µ). In the case of Cartesian vectors (such as the mag-netization matrix or the magnetic field interacting with themagnetization), the notation a⃗†µ · b⃗µ stands for the regularscalar product with the summation over the orbital indices.The notation a⃗† · b⃗ implies the summation over the atomicindices as well.• The n× n×m×m tensorA= [Aab,cd], with first two orbit-als (ab) residing on the site µ and last two orbitals (cd) resid-ing on the site ν, can be viewed as the n2 ×m2 matrix Âµν .The construction Âµν b⃗ν implies the summation over the twoorbital indices on the site ν.• The 2nd derivative is a local probe and the interaction para-meter depends on the point in which it is calculated. Forinstance, considering the simplest interaction energy E=−Jcosφ between two spins in the bond, the 2nd deriv-ative near FM (φ= 0) and AFM (φ = π) configurationsof spins will be, respectively, J and −J. Nevertheless,as it is typically done, we will additionally change thesign of interaction parameters for the antiferromagneticallycoupled bonds, thus adopting the universal definition whereJ> 0 and < 0 stands for the FM and AFM interactions,respectively.2.2. Spin spirals, SO coupling, and DM interactionsThe SO interaction is known to consist of the spin-diagonal,ξ2 L̂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 vectorof angular momenta, and σ̂ = (σ̂x, σ̂y, σ̂z) is that of the Paulimatrices. The antisymmetric DM interaction d= (dx,dy,dz)emerges in the 1st order of ξ and generally one should be ableto calculate all three vector projections onto x, y, and z (unlessthey are related by the symmetry properties). Nevertheless,by proper rotations of the coordinate frame, which transformxyz to zyx and zxy, dx and dy can be viewed as dz in the newcoordinate frame. Therefore, we need the numerical proced-ure only for calculating dz. The important point in this respectwas realized by Sandratskii [15], who suggested that in orderto calculate dz, it is sufficient to consider only spin-diagonalpart of the SO coupling. His idea was based on a simple obser-vation that DM interactions give rise to spiral magnetic struc-tures [45], which can be regarded as ‘eigenstates’ of the spinmodel (1). Therefore, the energies of the model (1) can beuniquely specified by the vectors q, describing propagationof the spin spiral. Then, the same should hold for the elec-tronic model, which is used for the mapping onto the spin one,and the spin spirals should be amount possible magnetic solu-tions of such model. In practical terms, this means that theelectronic states should obey the generalized Bloch theorem,which combines translations with the SU(2) rotations of spinsin the spiral texture [46]. Nevertheless, this theorem can beapplied only if the spin is the good quantum number so thatthe Hamiltonian Ĥσ remains diagonal with respect to the spinindices σ. Therefore, such Ĥσ can include the diagonal part ofthe SO coupling, but not the off-diagonal one.The method is suitable for the DM interactions, but not forthe magnetic anisotropy, which emerges in the 2nd order ofthe SO coupling and typically include both diagonal and off-diagonal contributions. This is again in line with the idea of thespin-spiral approach: the DM interactions give rise to the spinspirals, while the magnetic anisotropy acts against them, bydeforming the spin spirals and locking them to the crystallo-graphic lattice [47, 48]. The alternative to the spin-spiral tech-nique is to work in the real space separately for each magneticbond [49–52]. Such methods, which have certain limitations,will be briefly considered in section 6.1.2.3. General expression for the energy changeAs was already pointed out before, our basic idea is to ‘excite’the spin spiral, rotating the ground-state magnetization m̂GS =(0,0, m̂) asm̂q,i =(θ cosqRi,θ sinqRi,1−θ22)m̂, (4)(see figure 1) and evaluate the interactions between thetransversal ‘fluctuations’ of the magnetization m̂⊥q,i =(cosqRi,sinqRi,0)θm̂ for small θ. In order to induce suchm̂⊥q,i, we have to apply the external fieldĥq,i = (cos(qRi+αq) ,sin(qRi+αq) ,0) ĥ (5)in the direction perpendicular to the ground state magnetiza-tion. If the SO coupling is included, ĥq,i is not necessarily par-allel to m̂⊥q,i as the latter can experience the effect of the DM4J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewFigure 1. (a) Conical spin spiral, which is assumed in calculations of interatomic exchange interactions: q is the propagation vector, h0 isthe constraining field inducing the transversal magnetization m⊥i without the spin–orbit coupling, and mi is the rotated magnetization on thesite i. Reprinted figure with permission from [53], Copyright (2023) by the American Physical Society. (b) Configuration of constrainingfield in the xy plane: h0 is required to induce given transversal magnetization m⊥, while the additional perpendicular field, αnz× h⃗0, isrequired to compensate the additional rotation of m⊥ caused by the DM interaction dz.interaction dz, which tend to additionally rotate the magnetiz-ation in the xy plane. Therefore, the phases αq are needed tocompensate the effect of DM interactions (see figure 1). Forsmall αq, ĥq,i can be written asĥq,i ≈ ĥ0q,i+αqnz× ĥ0q,i, (6)where ĥ0q,i = (cosqRi,sinqRi,0) ĥ and nz = (0,0,1).The corresponding energy can be evaluated in the frame-work of constrained SDFT [34, 54] as:E [m⃗] = T [m⃗] + Exc [m⃗] +12h⃗†q ·(m⃗− m⃗⊥q), (7)where T and Exc are, respectively, the kinetic and xc energies,while the last term controls the size of the transversal magnet-ization m⃗⊥q . For simplicity, we drop here all irrelevant depend-encies of E on the charge density.Then, the kinetic energy can be expressed as the sum of theoccupied KS single-particle energies, Esp, calculated for theexternal field h⃗q and the xc fieldb⃗q = 2∂Exc [m⃗]∂m⃗∣∣∣m⃗=m⃗q,minus the interaction energy of m⃗q with h⃗q and b⃗q [13, 14],yieldingE [m⃗q] = Esp(h⃗q+ b⃗q)− 12(h⃗q+ b⃗q)†· m⃗q+ Exc [m⃗q] . (8)The important property of the xc energy in this respect,being the consequence of fundamental gauge invariance ofthe density functional theory [55], is that rotations of thespin magnetization m̂GS → m̂q,i do not change Exc: Exc[m̂q,i] =Exc[m̂GS] [56–58]. This is a general property of SDFT,which becomes especially transparent in the local spin-density approximation (LSDA), based on the picture of homo-geneous electron gas. In this case, Exc in each point rdepends only on the magnitude of the magnetization, ELSDAxc ≡ELSDAxc [|m(r)|] [59, 60], and therefore does not change underrotations of m(r). Since Exc[m̂q,i] = Exc[m̂GS], the rotation ofthe magnetization will rotate the xc field by the same angle:b̂q,i =(θ cosqRi,θ sinqRi,1−θ22)b̂. (9)Therefore, b⃗†q · m⃗q does not change either and equation (8) willlead to the following energy change:δEq = δEsp(h⃗q+ b⃗q)− 12h⃗†q · m⃗q.Then, δEsp can be evaluated by treating h⃗q+ b⃗q as a perturba-tion, to the 2nd order in h⃗q + b⃗⊥q and the 1st order in the lon-gitudinal change of the xc field, − 12 b̂θ2. The details are elab-orated in [61], leading to the simple but general expression:δEq =−14h⃗0†q · m⃗⊥q . (10)In fact, this result is well anticipated. On the one hand, δEqshould be proportional to h⃗q as there would be no energychange without the external field. On the other hand, thereonly possible interaction of h⃗q with the constrained magnetiz-ation m⃗⊥q is the scalar product given by equation (10). It maylook incomplete because δEq does not seem to know anythingabout αq and the DM interactions. Nevertheless, all necessaryinformation is in equation (10) and in section 2.5 we will showhow it should be used to derive practical expressions for theisotropic exchange and DM interactions.In addition to the rotations given by (4), the magnetiza-tion can experience the longitudinal change, which is caused5J. Phys.: Condens. Matter 36 (2024) 223001 Topical Reviewby these rotations. It will affect m̂, resulting in an additionalchange of each of the terms in equation (8). Nevertheless, thesecontributions can be shown to cancel out in the lowest orderof θ [11, 57].2.4. Response tensorThe response theory is basically the perturbation theory relat-ing the small change of the potential v⃗with the induced densityn⃗: n⃗= R̂v⃗. Spin-dependent v⃗ can be generally specified by fourelements:v⃗=(v⃗↑↑ v⃗↑↓v⃗↓↑ v⃗↓↓).Then, each v⃗σσ′induces the corresponding change n⃗σσ ′:n⃗σσ ′= R̂σσ ′v⃗σσ′, (11)where the rank-4 tensor R̂σσ ′can be found in terms of the1st-order perturbation theory for the wave functions [62]. Inour case, the perturbation is h⃗q+ b⃗⊥q and our goal is to eval-uate m⃗⊥q . Then, it is convenient to use the local coordinateframe where ĥi = (1,αq,0)ĥ0, which is obtained by rotatingĥq,i about z by the angles −qRi, and employ the generalizedBloch theorem, combining lattice translations with the SU(2)rotations of spins [46]. This will lead to the additional shiftof the k-mesh for the states with σ =↑ relative to those withσ =↓. Moreover, since ĥi = (1,αq,0)ĥ0 corresponds tov⃗h =12(0 1− iαq1+ iαq 0)h⃗0,we have to consider only R̂↑↓ and R̂↓↑. Then, the perturbationtheory yieldsR↑↓ab,cd (q) =∑mlkf↑mk− f↓lk+qε↑mk− ε↓lk+q(Ca↑mk)∗Cb↓lk+q(Cc↓lk+q)∗Cd↑mk,(12)where εσmk and |Cσlk⟩= [Caσlk ] are, respectively, the eigenval-ues and eigenvectors of Ĥσ, and fσmk ≡Θ(εF − εσmk) is theFermi distribution function. The orbital indices in each of thepairs ab and cd belong to the same atomic sites in the unitcell. Furthermore, R↓↑ can be obtained from R↑↓ using thepropertyR↑↓ab,cd (q) =[R↓↑ba,dc (−q)]∗. (13)If Ĥσ remains invariant under the time reversal (e.g. withoutSO interaction), equation (13) is reduced to R↑↓ab,cd(q) =R↓↑ba,dc(q).22 Without SO interaction, R↑↓ab,cd and R↓↑ab,cd are considered only in the com-bination with the symmetric matrices m̂ and b̂. Therefore, one can writeR↑↓ab,cd =R↓↑ab,cd.R↑↓ab,cd(q) can be also related to the Green functionĜσ(ε,k) = [ε− Ĥσ(k)]−1 as [63]R↑↓ab,cd (q) =− 1πBZ∑kImˆ εF−∞dε{G↑da (ε,k)G↓bc (ε,k+ q)}(14)with the summation running over the 1st Brillouin zone (BZ).Then, using the definition (11), one can find:n⃗↑↓ =12R̂↑↓q(h⃗0 + b⃗⊥ − iαqh⃗0)andn⃗↓↑ =12R̂↓↑q(h⃗0 + b⃗⊥ + iαqh⃗0),where R̂σσ ′q ≡ R̂σσ ′(q). In the local coordinate frame, thesen⃗↑↓ and n⃗↓↑ should give us themagnetization m⃗⊥ = n⃗↑↓ + n⃗↓↑along x:m⃗⊥ = R̂+q(h⃗0 + b⃗⊥)− iαqR̂−q h⃗0, (15)where R̂±q = 12 (R̂↑↓q ±R̂↓↑q ). Another equation,iR̂−q(h⃗0 + b⃗⊥)+αqR̂+q h⃗0 = 0, (16)requires that the perpendicular to it magnetization along y,i (⃗n↑↓ − n⃗↓↑), should vanish (so as the y component of thexc field) according to our constraint conditions. These are theequations for h⃗0 and αq for given m⃗⊥ = θm⃗. Their meaning isvery straightforward. For instance, in equation (16), the iso-tropic part of the magnetization αqR̂+q h⃗0, which is induced byαqh⃗0 along y, is compensated by the one, which is induced dueto the DM interaction by the field h⃗0 + b⃗⊥ acting in the perpen-dicular direction x. The same is with equation (15), where m⃗⊥has two components: the isotropic one, induced by h⃗0 + b⃗⊥along x, and the one caused by the DM interaction, transfer-ring the effect of the magnetic field αqh⃗0, applied along y, tothemagnetization along x. This explains how one can naturallyseparate the contributions of the isotropic and DM interactionsin equation (10).2.5. Exchange interactionsThe next step is the mapping of the total energy change (10)onto the spin model:δEq =−12∑µν(Jq,µν − idzq,µν)θµθν , (17)where we explicitly consider the possibility of having sev-eral magnetic sublattices. In the local coordinate frame,equation (10) can be rearranged asδEq =−14(h⃗0 + b⃗⊥)†· m⃗⊥ +14b⃗⊥† · m⃗⊥,6J. Phys.: Condens. Matter 36 (2024) 223001 Topical Reviewwhere we have added and subtracted the xc field b⃗⊥. Ourstrategy is to start with the expression for m⃗⊥ without the SOcoupling, which is given by the 1st term in equation (15), andthen consider the corrections arising in the 1st order of the SOcoupling, which are given by the 2nd term.3 Then, noting thath⃗0 + b⃗⊥ = Q̂+q m⃗⊥, (18)where Q̂σσ ′q = [R̂σσ ′q ]−1, Q̂±q = 12 (Q̂↑↓q ±Q̂↓↑q ), andwithout SO coupling Q̂+q = [R̂+q ]−1, one immediatelyfinds the following expression for the isotropic exchangeinteractions:Jq,µν =12(m⃗†µ · Q̂+q,µνm⃗ν − b⃗†µ · m⃗µδµν). (19)Considering the 2nd term in equation (15), the constructioni (⃗h0 + b⃗⊥)† · R̂−q αqh⃗0 describes the interaction between x andy components of the magnetic field caused by the DM interac-tions. Corresponding interaction parameter should satisfy theconditiondzq,µν θµθν =i2(h⃗0µ + b⃗⊥µ)†· R̂−q,µν(h⃗0ν + b⃗⊥ν), (20)where we had to ‘rescale’ y components of the magnetic field,αqh⃗0µ → h⃗0µ + b⃗⊥µ , in order to specify x and y components ofthe transversal magnetization by the same set of parametersθµ and θν . Then, using equation (18) and noting that to the1st order in the SO coupling Q̂+R̂−Q̂+ =−Q̂−, one can findthatdzq,µν =− i2m⃗†µ · Q̂−q,µνm⃗ν . (21)This expressionwas obtained in [53] basically heuristically, bythe analogy with isotropic interactions and similar expressionformulated in terms of the xc fields, which will be consideredin section 2.8. Here, we have provided a more rigorous proofof equation (21). The real space parameters can be obtainedby the Fourier transform of Jq,µν and dzq,µν .Thus, the exchange interactions are proportional to theinverse response function. For the isotropic exchange, this isbasically the result of Bruno [34]. For theHubbardmodel, sim-ilar relationship has been established by Szczech et al [64].For practical purposes, it may be more convenient to calcu-lateXq,µν =12(m⃗†µ · Q̂↑↓q,µνm⃗ν − b⃗†µ · m⃗µδµν), (22)in terms of only Q̂↑↓q , and then relate it with Jq,µν and dzq,µνusing the property (13), which yieldsJq,µν =12(Xq,µν +X∗−q,µν)(23)3 Since we neglect spin-off-diagonal elements of the SO coupling (seesection 2.2), higher-order corrections are meaningless.anddzq,µν =− i2(Xq,µν −X∗−q,µν). (24)Thus, Jq,µν is related to the average energy of spin spiralspropagating in q and −q, while dzq,µν is related to the energydifference [15].2.6. Sum rule and local xc fieldThe sum rule is obtained from the identity[Ĝ↓ (ε,k)]−1−[Ĝ↑ (ε,k)]−1= b̂,which can be further rearranged asĜ↑ (ε,k)− Ĝ↓ (ε,k) = Ĝ↓ (ε,k) b̂ Ĝ↑ (ε,k)= Ĝ↑ (ε,k) b̂ Ĝ↓ (ε,k) ,where b̂= Ĥ↑(k)− Ĥ↓(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 can find:m⃗= R̂↑↓0 b⃗. (25)This sum rule has very straightforward meaning: q= 0 cor-responds to the uniform rotation of the ground-state magnet-ization, where all spins are rotated in the same direction bythe same angle. Therefore, the transversal magnetization isdescribed by the same xc field b⃗ as in the ground state (withoutany constraining fields).Nevertheless, in the TB representation, such xc field is notnecessary local. For instance, in LSDA, the splitting Ĥ↑(k)−Ĥ↓(k) can have interatomic matrix elements and depend onk. In such a situation, it can be important to reenforce the sumrule, by defining new local xc field as b⃗= Q̂↑↓0 m⃗, which wouldyield the given ground-state magnetization m⃗. For instance,this is a simple and transparent alternative to the kernel poly-nomial method, which was recently proposed to deal with non-local matrix elements of the xc field [30]. In fact, if the xc fieldis nonlocal, the total energy change for the infinitesimal rota-tions of spins is no longer representable in the form of pairwiseinteractions.2.7. Right object to rotate: magnetization matrices versuslocal magnetic momentsSo far, we did not properly specify the spin object whichshould be rotated on the magnetic sites in order to obtain thetotal energy change (10). All above discussions implied thatit is the magnetization matrix m̂µ, while the spin model istypically formulated in terms of the magnetic momentsMµ =TrL{m̂µ}. Undoubtedly, the rotation of m̂µ, as a whole, by theangle θµ will rotate Mµ by the same angle. However, is thischoice unique? Are there other perturbations of m̂µ, resultingin the same rotations of Mµ but preferably at lower energy7J. Phys.: Condens. Matter 36 (2024) 223001 Topical Reviewcost? Here, wewill follow the discussion in [61]. Nevertheless,wewould like to note that somewhat similar ideas can be foundin the work of Antropov et al [65].Indeed, for the Hermitian matrix m̂µ, one can alwayschoose the diagonal representation m̂µ = diag( . . . , maµ, . . .)with respect to the orbital indices. In principle, each orbital a insuch 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 be m̂⊥µ = diag( . . . , θaµmaµ, . . .).Nevertheless, these angles are subjected to the additionalconstraint because TrL{m̂⊥µ } should be equal to θµMµ.Importantly, this condition is softer than the rotation of m̂µ asa whole, where all orbitals are rotated by the same θaµ = θµ.Therefore, it is reasonable to expect that the energy changewill be smaller, so as the exchange parameters.Mathematically, we have to minimize the energychange (10) with the additional condition∑a(θaµ − θµ)maµ =0 on each site µ:δE =−14∑µa{θaµmaµh0aµ −(θaµ − θµ)maµλµ}, (26)where h0aµ are the constraining fields acting on maµ and λµ arethe Lagrange multipliers. Then, minimizing δE with respect toθaµ, it is straightforward to find that h0aµ = λµ. Thus, in order torotate the spin moments at the minimal energy cost, one haveto apply the scalar field hµ, i.e. the same for all orbitals a.Moreover, in this case it is convenient to use the ‘sphericallyaveraged’ version of the linear response, where m̂µ is replacedby Mµ, b̂µ is replaced by Bµ = 1nTrL{b̂µ}, and R↑↓ab,cd(q) isreplaced byR↑↓µν (q) =∑a∈µ,c∈νR↑↓aa,cc (q) .The corresponding exchange interaction parameters will begiven byJq,µν =12(MµQ+q,µνMν −BµMµδµν)(27)anddzq,µν =− i2MµQ−q,µνMν , (28)where Q̂±q = 12 (Q̂↑↓q ± Q̂↓↑q ), Q̂σσ ′q = [R̂σσ ′q ]−1, and R̂σσ ′q ≡R̂σσ ′(q) is the matrix specified by the atomic indices in theunit cell.In comparison with equations (19) and (21), based on rota-tions of the magnetization matrix, equations (27) and (28) areexpected to be more suitable for the analysis of low-energyexcitations, of course, provided that the latter can be describedby the spin model (1). In the following, these twomethods willbe denoted as m̂ andM, after the basic variable describing theinfinitesimal rotations of spins.2.8. Rotations of the xc field as an alternative perturbationIn this section, we consider the original formulation of thelinear response theory, as it was proposed by Liechtensteinet al [11], which is frequently called the magnetic force the-orem [34]. However, there are two important points about thework of Liechtenstein et al [11], which should be distinguishedfrom each other [66]:• The general claim that, in SDFT, the energy change causedby the infinitesimal rotations of spins can be related to theKS eigenstates in the ground state is certainly correct andshould not be revised. This is what is actually called the‘magnetic force theorem’ stating that the interaction para-meters are the ground state properties and can be found byknowing the electronic structure in the ground state;• Nevertheless, the practical expression (3), which wasderived by Liechtenstein et al [11] for the exchange interac-tions, relies on additional approximations and, in principle,can be improved.The starting assumption of Liechtenstein et al [11] is that sincethe rotation of the magnetization results in the rotation of thexc field by the same angle (section 2.3), it is logical to treatthe change of the xc field, δb̂q = b̂q− b̂GS, as a perturbationwithout the constraining field. Then, we have to consider onlyδEsp in equation (9) caused by this δb̂q. Furthermore, the trans-versal magnetization, m⃗⊥ ′q , which is induced by δb̂q will gen-erally differ from m⃗⊥q because, without the constraining field,m⃗q will tend to relax toward the ground state [33–35]. TheDM interactions, if any, will tend to additionally rotate m⃗⊥ ′qrelative to b⃗⊥q : m⃗⊥ ′q ≈ m⃗⊥0 ′q +βqnz× m⃗⊥0 ′q . Thus, instead ofequation (10), this method relies on (in the local coordinateframe)δEq ≈14b⃗⊥† · m⃗⊥ ′− 14b⃗† · m⃗θ2, (29)arising from the single-particle energies for the perturbationscaused 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 identic-ally be equal to zero, as was discussed in section 2.3. Then,using the definitions m⃗⊥0 ′= R̂+q b⃗⊥ and βqm⃗⊥0 ′= iR̂−q b⃗⊥,one can find thatJq,µν =−12(b⃗†µ · R̂+q,µν b⃗ν − b⃗†µ · m⃗µδµν), (30)anddzq,µν =i2b⃗†µ · R̂−q,µν b⃗ν . (31)Alternatively, dzq,µν can be obtained from equation (20)for h⃗0µ = 0 and b⃗⊥µ = b⃗µθµ. Substituting equation (14) intoequation (30) and Fourier transforming it to the real space, oneobtains the well-known equation (3). Nevertheless, these arethe approximate expressions, which can be formally obtainedfrom the exact ones, equations (30) and (31), replacing m⃗ by8J. Phys.: Condens. Matter 36 (2024) 223001 Topical Reviewb⃗ and Q̂ by R̂, with the additional minus sign. In the follow-ing, we will refer to this method as ‘method b̂’ or ‘approximatemethod b̂’.In principle, one can also introduce the ‘spherically aver-aged’ version of this method (the so-called method B) repla-cing b⃗ν by Bν and R̂ by R [34], though it is rarelyused. Without the SO coupling, R̂+ = R̂↑↓ and correspond-ing exchange interactions ĴBq ≡ [JBq,µν ] can be written as ĴBq =− 12 B̂( R̂↑↓q + Î−1 )B̂, where B̂= diag( . . . , Bν , . . .) and Î =diag( . . . ,−Bν/Mν , . . .) are the diagonal matrices of, respect-ively, exchange splittings and effective Stoner parameters. Byadapting the samematrix form for the ‘exact’ interactions (27),Ĵq = 12M̂( [R̂↑↓q ]−1 + Î )M̂ with M̂= diag( . . . ,Mν , . . .), onecan find the following expression, connecting Ĵq with ĴBq [34]:Ĵq = ĴBq(1− 2B̂−1M̂−1ĴBq)−1. (32)Thus, Ĵq can be indeed replaced by ĴBq at least in two cases:(i) the long wavelength limit q→ 0 and (ii) the strong-coupling limit B̂→∞. Therefore, the spin-wave stiffness inthe limit q→ 0 is expected to be the same in both meth-ods. Nevertheless, this statement should not be exaggeratedbecause equation (32) holds only in the spherical case, wherethe xc field and the magnetization on each magnetic site aregiven by the scalar parameters Bν andMν . In the matrix case,the simple relationship (32) is no longer valid [61]. That is whyeven the spin-wave stiffness in themethods b̂ andM can be dif-ferent. Furthermore, we will see that there is indeed a numberexamples, where the approximate method b̂ fails to reproducethe correct magnetic ground state, while the method M dra-matically improves the description.Of course, it is reasonable to ask what are the right objectsto rotate in this case: whether they should be the wholematrices b̂ν or only the spherical parts of these matricesBν? If in the case of the magnetization, the answer can befound by minimizing the energy change (10) (see section 2.7),equation (10) is not applicable for rotations of the xc field.Therefore, the answer is open. However, historically most ofthe applications deal with the rotations of the matrices b̂ν .Considering the strong-coupling limit in equation (30) [67,68], one can derive all known expressions for the doubleexchange Jij ∼ ⟨Ĥij⟩ [16], superexchange Jij ∼−⟨Ĥ2ij⟩/U [2],superexchange with the interatomic Coulomb repulsion Jij ∼−⟨Ĥ2ij⟩/(U−V) [69], etc where U and V is the on-site andintersite Coulomb repulsion, respectively, Ĥij are the transferintegrals (see appendix A), which are typically associated withthe matrix elements of the KS Hamiltonian in LDA (GGA),and ⟨. . .⟩ denotes the expectation value in the ground state.The strong-coupling limit for the DM interaction (31) resultsin the spin-current model [58, 70], which can be viewed asthe relativistic counterpart of the double exchange mechan-ism [53]. The expression for RKKY interactions can be alsoderived starting from equation (30), but using slightly differentphilosophy [9]. In this case, b⃗ is the field created by localizedcore spins and acting on outer conduction electrons. Withoutb⃗, the conduction bands are non-magnetic (and the tensor R̂ isevaluated in this non-magnetic state [9]).2.9. Relationship to the spin-wave spectraIn the previous sections, we have considered how the spinmodel can be generally derived from the electronic one usingthe concept of infinitesimal rotations of spins. The paramet-ers of such spin model are expressed in terms of the staticspin susceptibility (or the response function). On the otherhand, the spin-wave dispersion, ωq, which is the experiment-ally measurable quantity, can be derived in the framework ofRPA from the poles of the dynamic spin susceptibility [71–73]. However, this ωq does not necessary coincide with theone of the spin model with the parameters derived from thestatic spin susceptibility [63, 64]. In this respect, Katsnelsonand Lichtenstein [63] have argued that although the methodproposed by Bruno [34] is more consistent with the staticresponse formulation, the method b̂ should more suitable forthe analysis of the spin-wave spectra. Here, we will brieflyconsider this problem. For simplicity, we assume that there isonly one magnetic site in the unit cell and drop all matrix nota-tions. Then, in the spherical case, the dynamic response func-tion is given by R̃↑↓q (ω) = R↑↓q (ω)[1+ IR↑↓q (ω)]−1, whereR↑↓q (ω) is obtained from equation (12) replacing in the denom-inator (ε↑mk− ε↓lk+q) by (ω+ ε↑mk− ε↓lk+q). Therefore, one hasto solve the equation1+ IR↑↓q (ωq) = 0, (33)which can be equivalently rearranged as: 2M−1Jq(ωq) = 0,where Jq(ω) is given by equation (27) with Q↑↓q (ω) instead ofQ↑↓q =Q+q . Thus, the problem is that the spin-wave energiesare given by the zeros of Jq(ω) and do not necessary coincidewith ωMq = 2M−1Jq(0), expected from the solution of the spinmodel.In the limit ω≪ ε↓lk+q− ε↑mk (which takes place, forinstance, for insulating and half-metallic materials, wherethe occupied and unoccupied states with opposite projectionsof spins are separated by an energy gap), one can use thelinearization R↑↓q (ω)≈ R↑↓q (0)+ωṘ↑↓q (0), where Ṙ↑↓q (0) =∂∂ωR↑↓q (ω)|ω=0, and find the following expression: ωq =ωBq I−1[Ṙ↑↓q (0)]−1B−1, where ωBq = 2M−1JBq (0) is the spin-wave dispersion calculated with the parameters of the schemeB. This example clearly shows that ωBq should be addition-ally renormalized, though this renormalization is generallydifferent from the one given by equation (32), connectingthe parameters of the methods M and B. Ṙ↑↓q (0) depends onthe details of the electronic structure. In practical terms, itcan be calculated replacing (ε↑mk− ε↓lk+q) by −(ε↑mk− ε↓lk+q)2in equation (12). Then, in certain circumstances, the methodM can be a good starting point for the analysis of the spin-wave dispersion. For instance, if the ↑-spin (↓-spin) statesare fully occupied (empty) and B is large compared to theband dispersion, it is straightforward to obtain that Ṙ↑↓q (0)≈−B−1R↑↓q (0) and ωq ≈ ωMq .Thus, it would be fair to conclude that the analysis of thespin-wave dispersion requires the additional renormalizationof the parameters derived from the static spin susceptibil-ity [63, 64]. This conclusion applies to allmethods (M, B, and9J. Phys.: Condens. Matter 36 (2024) 223001 Topical Reviewb̂). Therefore, this is an open question which method servesbetter for the analysis of the spin-wave dispersion. It is cer-tainly true that, in LSDA, the method b̂ better reproduces theexperimental spin-wave dispersion in the canonical case ofbcc-Fe and fcc-Ni [61, 63]. However, this conclusion does notseem to be general and for other materials the comparison canbe less favorable.Finally, we note that equation (33) can be further rearrangedas R↑↓q (0)(B+ hω) =M, where hω = [R↑↓q (0)]−1[R↑↓q (ω)−R↑↓q (0)]B≈ ω[R↑↓q (0)]−1Ṙ↑↓q (0)B has a meaning of the con-straining field, which is needed to correct the effect of thexc field B in order to reproduce the ground-state magneticmoment for an arbitrary q. In this sense, there is an analogybetween the search of the poles of the dynamic susceptibilityand the constrained SDFT considered in section 2.3.2.10. Elimination of the ligand spinsBy knowing Jq,µν and dzq,µν , one can, in principle, calculateisotropic and DM interactions operating between all sites inthe unit cell. Nevertheless, themagnetization at these sites mayhave completely different origin. For instance, the transition-metal (T) sites in many oxide materials participate as a sourceof the magnetism, being primarily responsible for the spon-taneous time-reversal symmetry breaking, while the oxygenor any other ligand (L) sites behave as ‘magnetic slaves’:although they can host an appreciable portion of the magnetiz-ation, it is solely induced by hybridization with the T sites andstrictly follow the change of the magnetization on the T sites.The corresponding energy change for each q can be schemat-ically expressed asδE =−12(θ⃗ †T · X̂TTθ⃗T + θ⃗ †T · X̂TLθ⃗L + θ⃗ †L · X̂LTθ⃗T + θ⃗ †L · X̂LLθ⃗L),(34)where X̂AB is the matrix specified by the atomic sites of thetypes A or B, θ⃗A are the polar angles specifying the rotations ofmagneticmoments of the typeA in the form of the column vec-tor, and θ⃗ †A is the corresponding to it row vector. Then, one cantry to eliminate the L degrees of freedom by transferring theireffect into the interaction parameters between the T sites. Thiscan be done by employing the ideas of adiabatic spin dynam-ics [74, 75] and assuming that the T spins are sufficiently slowso that the L spins have sufficient time to adjust each changein the system of the T spins. Mathematically, this means thatfor each θ⃗T, θ⃗L can be found from the condition ∂∂θ⃗ †LδE = 0,yieldingθ⃗L =−[X̂LL]−1X̂LTθ⃗T (35)and δE =− 12 θ⃗†T ·ˆ̃XTTθ⃗T withˆ̃XTT = X̂TT − X̂TL[X̂LL]−1X̂LT. (36)The corresponding parameters of isotropic and DM interac-tions can be obtained from ˆ̃XTT using equations (23) and (24).In the method M, the matrix inversion in equation (36)can be combined with the one of the response matrix Q̂↑↓ =[R̂↑↓]−1 to obtain the following expression for ˆ̃XTT:ˆ̃XTT =ˆ̃X0TT +∆ˆ̃XTT, (37)whereˆ̃X0TT = M̂T[R̂↑↓TT]−1M̂T, (38)and∆ˆ̃XTT = M̂T[R̂↑↓TT]−1R̂↑↓TLQ̂↑↓LL(Q̂↑↓LL + ÎL)−1× ÎLR̂↑↓LT[R̂↑↓TT]−1M̂T. (39)Here, M̂T is the diagonal matrix of magnetic moments on thesites T, and ÎL is the diagonal matrix of effective Stoner para-meters on the sites L. In this expression, the explicit depend-ence of ˆ̃XTT on ÎL is incorporated into∆ˆ̃XTT, while ˆ̃X0TT form-ally does not depend on ÎL. The parameters ÎL can play avery important role in the theory of exchange interactions. Forinstance, according to GKA rules, they largely contribute tothe FM coupling in the systems, where the T–L–T bond angleis close to 90◦ [36–39]. In the linear response theories, theseeffects are incorporated into ∆ˆ̃XTT [76].Furthermore, the representation (37) allows us to improvethe numerical accuracy. Since in most of the systems the lig-and band is filled, the matrix elements of R̂↑↓ associated withthe L sites are typically small. Therefore, when we invert R̂↑↓,we have to deal with very large numbers even for the T sub-lattice. This is the reason why the bare interactions X̂TT∼Q̂↑↓TTare typically large and strongly compensated by the secondterm in equation (36) [61, 76]. Similar situation occurs in thescheme m̂. On the contrary, the calculation of ˆ̃X0TT and ∆ˆ̃XTTusing equations (38) and (39) involves the inversion of onlythe T block of R̂↑↓. This procedure is numerically much morestable rather than the inversion of the whole matrix R̂↑↓ inequation (36).The idea of downfolding somewhat similar to ours was pre-viously considered by Mryasov et al in order to eliminate the5d states of the heavy Pt atoms and explain the unusual tem-perature dependence of the magnetic anisotropy energy in theordered FePt alloy [77]. A simplified approach in the frame-work of the scheme b̂, which did not take into account theeffects of ligand–ligand interactions and ÎL, was also con-sidered by Logemann et al [78].3. Chromium trihalidesIn order to illustrate abilities of considered linear responsetechniques, we start with the detailed analysis of exchangeinteractions in chromium trihalides CrX3 (X=Cl and I). Thesevan der Walls compounds crystallize in the rhombohedral R3structure, which is built from the honeycomb layers as shownin figures 2(a) and (b) [79, 80]. The interactions between10J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewFigure 2. (a) Top view on the CrCl3 (CrI3) layer. The hexagonal unit cell is denoted by the broken line. Reprinted figure with permissionfrom [61], Copyright (2021) by the American Physical Society. (b) Stacking of adjacent honeycomb layers with the notation of mainexchange interactions. Two Cr atoms in the primitive rhombohedral unit cell are denoted by different colors. Reprinted figure withpermission from [61], Copyright (2021) by the American Physical Society. (c) and (d) Densities of states (DOS) of CrCl3 and CrI3 in LDA(top) and LSDA for the ferromagnetic state (bottom). Shaded areas show partial contributions of the Cr 3d states. The Fermi level is at zeroenergy (the middle of the band gap in the insulating phase).the layers are weak, but not negligible. For instance, sizableexchange interactions spread up to 6th coordination sphere, inand between the layers, as shown in figure 2(b).CrI3 is the ferromagnet with the Curie temperature TC =61K [80–82], while CrCl3 is a antiferromagnet with the Néeltemperature TN = 14.1K [83, 84]. The AFM transition inCrCl3 is followed by another transition to a pseudo-FM phasewith TC ∼ 17K. In both cases, the magnetic moments tend toorder ferromagnetically within the honeycomb layers. BelowTN, the interlayer coupling is weakly AFM, while in the tem-perature interval TN < T< TC the magnetic behavior of CrCl3is explained by the interlayer disorder [84].CrI3 is viewed as the prominent two-dimensional ferro-magnet [85], where one of the key ingredients is the strongSO coupling steaming from the heavy I atoms [86], whichis mainly responsible for the exchange anisotropy and emer-gence of the long-range FM order at relatively high TC(i.e. contrary to what would be expected from the Mermin–Wagner theorem in the isotropic case [87]). Besides that, CrCl3and CrI3 are regarded as the testbed materials for studyingfundamental aspects related to the origin of the ferromag-netism. Namely, why are these materials FM? The popularanswer is that since the Cr–X–Cr angle is close to 90◦, theinteraction is expected to be FM due to intra-atomic exchangecoupling IX on the ligand sites, as prescribed by the GKArules [36–39]. However, the GKA rule for the 90◦ exchangeis not very conclusive, because typically there are severalcompeting mechanisms supporting either ferromagnetism orantiferromagnetism [88]. In fact, Kanamori himself admittedthat there are several exceptions from this rule in the 90◦case [39]. Moreover, below we will see that, under certain cir-cumstances, IX can easily become negative and act against theferromagnetism.Although CrI3 is more useful practically, CrCl3 is inter-esting from the explanatory point of view. Even in LSDA,the electronic structure of CrCl3 consists of well separatedCr t2g, Cr eg, and Cl 3p bands, as explained in figure 2(c).4Therefore, it is possible to study separately the contributions ofeach of these bands to the exchange interactions by construct-ing proper TB models in the basis of Wannier functions [43,89]. Besides that, one can also consider correlated models,both for CrCl3 and CrI3, which explicitly consider the on-siteCoulomb interactions. This procedure is briefly explained inappendix A. The numerical calculations are performed on thebasis of the linear muffin-tin orbital (LMTO) method in theatomic-spheres approximation [90, 91], using the mesh of the10× 10× 10 points, both for the k- and q-integration.First, we will review the behavior of interatomic exchangeinteractions depending on each new ingredient added to themodel as well as the type of infinitesimal spin rotations(whether the rotated object is b̂, m̂, or M). Having estim-ated the parameters of interatomic exchange interactions, onecan find the magnetic ground state and evaluate the magnetictransition temperature as explained in appendix B. Then, thedetailed comparison with the experimental data will be givenin section 3.6.3.1. Three-orbital t2g model in LSDAThe simplest model for CrCl3 is the half-filled t2g model. Itcontains only occupied ↑-spin t2g band and unoccupied ↓-spint2g band. Since all basis functions of this model are associ-ated with the Cr states, there are no additional complicationscoming from the ligand states. Even in LSDA, the exchangesplitting BCr between the ↑- and ↓-spin t2g bands is large4 In the octahedral environment, the Cr 3d states are split into triply-degenerate t2g levels and doubly-degenerate eg levels. In many transition-metal oxides or related materials, the strong crystal-field splitting, 10Dq, iscaused by the hybridization between 3d and ligand p states [92].11J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewTable 1. Isotropic exchange interactions in CrCl3 (in meV): resultsof the t2g model in LSDA. b̂, m̂, and M stand for the methods basedon the infinitesimal rotations of, respectively, the xc field, themagnetization matrix, and the local spin moments. The notations ofJk are explained in figure 2.Method J1 J2 J3 J4 J5 J6b̂ −6.44 −0.50 −0.69 −0.17 −0.18 −0.63m̂ −6.32 −0.49 −0.68 −0.17 −0.18 −0.62M −6.30 −0.49 −0.69 −0.17 −0.18 −0.63compared to the bandwidth. Thus, the system should be inthe strong-coupling limit. The corresponding parameters ofexchange interactions, calculated by rotating b̂, m̂, or M aresummarized in table 1.In this casewe obtain very consistent description as all threemethods provide very similar sets of parameters. This is notsurprising: the approximate scheme b̂ is justified in the strongcoupling limit. Moreover, three t2g levels are nearly degener-ate. Therefore, the asphericity of m̂ is small and correspondingparameters are practically identical to the ones obtained in theM scheme. However, all interactions are AFM, which is quiteexpected for the half filling [2], but totally contradicts to theexperimental situation.3.2. Five-orbital Cr 3d modelThe next important question is whether the ferromagnetismof CrX3 can be explained without the ligand states, in themodel including both t2g and eg bands. In LSDA, the ↑- and↓-spin bands are subjected to the exchange splitting BCr. Atthe first sign, there is only a small addition in comparisonwith the t2g model—the unoccupied eg bands. However, itchanges the story dramatically. First, the crystal-field splittingbetween t2g and eg bands is comparable with BCr. Moreover,since only the transitions between the occupied ↑-bands andunoccupied ↓-bands contribute to R̂↑↓, such exchange inter-actions do not know anything about the existence of the ↑-spin eg band. Thus, although BCr is large (as in the t2g model),the behavior of the exchange interactions is also controlledby other details of the electronic structure and the system isno longer in the strong coupling limit. Furthermore, the mat-rix m̂ in the basis of all five 3d orbitals acquires additionaldegrees of freedom besidesM, which is only the spherical partof m̂.All these tendencies are reflected in the behavior ofinteratomic exchange interactions (table 2). First, the inclu-sion of the eg band into the model gives rise to the FM interac-tions. These interactions prevail in the nearest-neighbor (nn)bonds and considerably weaken the AFM interactions in otherneighbors, in agreement with the experimental situation. Then,the schemes b̂, m̂, andM provide quite a different description.Compared to the methods m̂ and M, the approximate methodb̂ underestimates the FM interaction J1. Moreover, the methodm̂ (in comparison with M) has a tendency to overestimate thein-plane interactions J1, J3, and J6, as expected from the ana-lysis in section 2.7.Another important question is how to define the xc fieldin the TB model. The common practice is to use b̂= Ĥ↑ −Ĥ↓ and take only site-diagonal (local) part of this splitting.However, in LSDA, the off-diagonal elements of Ĥ can alsodepend on spin, giving rise to non-local contributions to thexc field [30]. Then, the use of the site-diagonal part alone willviolate the sum rule because R̂↑↓0 b⃗ with such b⃗ will no longerreproduce the ground-state magnetization m⃗ (see section 2.6).Another possibility is to reenforce the sum rule by defining thelocal xc field b⃗= Q̂↑↓0 m⃗, which would reproduce the groundstate magnetization m⃗. In the three-orbital (3o) t2g model, thisresults only in minor change of the exchange interactions.However, starting from the five-orbital (5o) model, the differ-ence becomes more significant (see table 2).Similar analysis can be performed for the correlated modelin the Hartree–Fock approximation (see appendix A). Thenon-interacting electron part of the model Hamiltonian inthis case is evaluated in LDA. The corresponding electronicstructure is shown in figures 2(c) and (d): since in CrCl3and CrI3 the Cr 3d bands are well separated from the lig-and ones, such model can be easily constructed for both com-pounds. The screened Coulomb interactions are evaluatedwithin constrained RPA (cRPA) [93], as explained in [89].The obtained (averaged) parameters of on-site Coulomb repul-sion U, intra-atomic exchange interaction J, and nonspheri-city B (see appendix A for explanations) are U= 1.79 (1.15),J= 0.85 (0.78), and B= 0.09 (0.07) eV for CrCl3 (CrI3). Thescreened U is not particularly large. Moreover, the screeningis more efficient in CrI3 due to proximity of the I 5d and Crt2g bands [89]. The corresponding densities of states (DOSs)are shown in figure 3. As expected [94, 95], the Coulomb Ufurther increases the band gap in comparison with LSDA. Theband gap is smaller in CrI3 because U is smaller.Corresponding parameters of exchange interactions aresummarized in tables 3 and 4 for CrCl3 and CrI3, respect-ively. The approximate method b̂ systematically underestim-ates the FM interactions. For instance, all interactions in CrCl3are AFM, which clearly contradicts to the experimental situ-ation. In CrI3, only J1 is weakly FM, which is not enough tostabilize the FM ground state [96]. The methodM systematic-ally improves the situation: the FM interactions clearly prevailand the FM ground state is stabilized both in CrCl3 and CrI3.5The corresponding Curie temperature, evaluated in RPA [97](see appendix B for details), is also in reasonable agreementwith the experimental data. The method m̂ has a tendency tooverestimate the interactions J1, J3, and J6, making the FMground state in CrI3 unstable.Thus, the minimal model, which can capture the FM char-acter of the exchange interactions in the chromium trihalidesis the 5o model. Formally, the ferromagnetism can be obtainedwithout invoking the ligand states. Nevertheless, it is crucial toconsider the unoccupied eg bands. Furthermore, it is cruciallyimportant to use the exact technique, formulated in terms of the5 The situation with the AFM interlayer coupling in the case of CrCl3 is ratherfragile: the interaction J2 is indeed weakly AFM. However, in theM method,it is counterbalanced by two weakly FM interaction J4 and J5.12J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewTable 2. Isotropic exchange interactions in CrCl3 (in meV): results of the Cr3d model in LSDA. b̂, m̂, and M stand for the methods basedon the infinitesimal rotations of, respectively, the xc field, the magnetization matrix, and the local spin moments. The xc field was eithercomputed from the on-site spin splitting of the TB Hamiltonian (H) or obtained from the sum rule (sr). TS is corresponding spin transitiontemperature in RPA (in K). The method m̂ yields the ferromagnetic ground state, while the methods b̂ and M result in an incommensuratespin-spiral structure propagating both in and between the planes. The notations of Jk are explained in figure 2(b).Method J1 J2 J3 J4 J5 J6 TSb̂ (H) 2.39 −0.27 −0.45 −0.03 −0.04 −0.50 8b̂ (sr) 2.27 −0.16 −0.37 0.01 0 −0.46 5m̂ 11.94 −0.01 −0.69 0.08 0.09 −0.69 30M 4.68 −0.21 −0.40 0 −0.01 −0.48 5Figure 3. Densities of states (DOS) for the FM state of five-orbital Cr 3d model in (a) LSDA for CrCl3, (b) Hartree–Fock approximation forCrCl3, and (c) Hartree–Fock approximation for CrI3. The zero energy is in the middle of the band gap.Table 3. Isotropic exchange interactions in CrCl3 (in meV): results of the 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 spinmoments. TS is the corresponding spin transition temperature in RPA (in K). The methods m̂ and M yield the ferromagnetic ground state,while the method b̂ results in an incommensurate spin-spiral state. The notations of Jk are explained in figure 2(b).Method J1 J2 J3 J4 J5 J6 TSb̂ −1.69 −0.13 −0.28 −0.02 −0.03 −0.18 6m̂ 6.66 −0.06 −0.58 0.02 0.02 −0.35 11M 3.76 −0.02 −0.13 0.05 0.05 −0.14 27Table 4. The same as table 3 but for CrI3. The method M yields the ferromagnetic ground state, while the methods b̂ and m̂ result in anincommensurate spin-spiral state. (A) denotes the same parameters calculated in the antiferromagnetic state.Method J1 J2 J3 J4 J5 J6 TSb̂ 0.76 −0.30 −0.39 0.04 0 −0.46 10m̂ 5.87 −0.10 −0.21 0.33 0.16 −0.71 14M 4.58 −0.08 −0.06 0.37 0.25 −0.37 51M (A) 4.52 −0.07 −0.06 0.36 0.25 −0.41 49inverse response function. The approximate method b̂, whichis linear in R̂, can lead to an incorrect magnetic ground state(see tables 3 and 4).Generally, the interatomic exchange interactions, definedvia infinitesimal rotations of spins, can depend on the mag-netic state. In a number of cases, this dependence can be verystrong, reflecting the dependence of the electronic structureon the magnetic arrangement [98, 99]. Considering how Jk inCrX3 change depending on the method used for their calcu-lations (for instance, M versus b̂), it is reasonable to expectthat the system is no longer in the strong-coupling limit and,therefore, these parameters can also depend on the magneticstate. However, this appears to be not the case. For instance,the exchange parameters calculated in the AFM state of CrI3are 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 alsoused to evaluate the spin transition temperature, TS. The mainreason why the b̂ and M methods yield different Jk is relatedto the fact that the crystal-field splitting, 10Dq, is comparable13J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewFigure 4. Schematic view on the distribution of magnetic momentsin polarized transition-metal (T) and ligand (L) bands: M0T and M0Ldenote the spin moments on the T and L sites in the T bands, andδMT and δML denote those in the L band. The total moments areMT(L) =M0T(L) + δMT(L). In the T band, the hybridization betweenT and L sites induces M0L in the same direction as M0T. Then, δMT inthe L band emerges as the joint effect of hybridization andintraatomic exchange interactions (denoted by dashed lines). Sincethe L bands is fully occupied, δML =−δMT, resulting in theantipolarization of M0L and δML (and also ML and MT).with the intra-atomic exchange splitting. However, this 10Dqdoes not depend on the magnetic state, and therefore does notcontribute to the magnetic-state dependence of Jk.3.3. Role of the ligand statesNow we turn to the analysis of the most general model, whichexplicitly includes the contributions of the Cr 3d as well wethe ligand Cl 3p or I 5p states.3.3.1. Cr 3d and Xp bands in LSDA. We start with the ana-lysis of general TBmodel, including both Cr 3d andX p bands,in LSDA. First, we note that the Cr 3d and ligand p states areantipolarized. In this case MCr exceeds the nominal value of3µB (where µB denotes the Bohr magneton). Nevertheless, itis compensated by the negative MX on the ligand atoms, sothat the total moment,Mtot =MCr + 3MX, is integer and equalto 3µB. This antipolarization is a joint effect of the spin split-ting on the Cr atoms and the hybridization between the occu-pied ligand states and unoccupied Cr eg states. Since the ↑-spin Cr eg states are closer to the ligand band, the hybridiza-tion is stronger, resulting in the stronger admixture of the ↑-spin Cr eg states into the occupied ligand band (and transferof the ↑-spin ligand states into the unoccupied Cr eg band).Therefore, when the ligand band is added to the model, wehave MCr > 3µB but MX < 0, as schematically illustrated infigure 4. For instance, the LSDA yields MCr = 3.14 (3.37) µBand MX =−0.05 (−0.12) µB for CrCl3 (CrI3). As expected,the effect is stronger in CrI3 where the Cr eg states are closerto the I 5p band (see figure 2(d)). Moreover, the I 5p states aremore extended in comparison with the Cl 3p ones, resulting instronger hybridization.Table 5. The effective Stoner parameters in CrCl3 and CrI3 (in eV),obtained using six possible definitions (as explained in the text).DefinitionCrCl3 CrI3ICr ICl ICr III 0.82 −3.08 0.76 −0.55II 0.81 −0.75 0.75 −0.13III 0.83 −0.28 0.77 0.31IV 0.83 −3.32 0.77 −0.46V 0.83 −2.50 0.77 −0.51VI 0.83 −0.64 0.77 −0.13The new aspect of calculations of the exchange interac-tions is that now they can also depend on the strength of theStoner coupling IX, which contributes to the second term ofequation (37). Therefore, the first problem we have to solveis how to properly define I. In fact, there are several possibledefinitions:(I) Iµ =−TrL{b̂µm̂µ}/(Mµ)2, where b̂µ is associated withthe intra-atomic spin splitting of Ĥσ and m̂µ is theground-state magnetization;(II) Iµ =−Bµ/Mµ, where Bµ = 1nTrL{b̂µ} and Mµ =TrL{m̂µ}, which is nothing but the spherically averagedversion of (I);(III) the same as (I), but taking b̂µ 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 taking m̂µ andMµ from the sum rule(and defining b̂µ as the intra-atomic spin splitting of Ĥσ);(VI) the same as (II), but taking Mµ from the sum rule.The results are summarized in table 5. One can see that ICronly weakly depends on the definition (thought it is differentin CrCl3 and CrI3). However, we do not need this parameter inour calculations. On the other hand, IX is very sensitive to thedefinition [76]. Apparently, one can discard the definitions Vand VI as they yield incorrect Mtot violating the fundamentalproperty Mtot = 3µB.6 Furthermore, it makes sense to enforcethe sum rule by using the definitions III and IV instead of I andII. Finally, the definitions II and IV are more appropriate forthe spherically averaged methodM, while the definitions I andIII should be used in the combination with the matrix form ofthe methods m̂ and b̂.While ICr in CrX3 is close to atomic values (∼0.7 eV)and can be interpreted as the intra-atomic exchange integ-ral responsible for Hund’s first rule [100–102], IX is defin-itely not. In most of the cases IX is negative (except schemeIII for CrI3). This is another manifestation of the fact thatthe moments MX are merely induced by the hybridizationwith the Cr 3d states, which acts against BX . According toequation (39), the negative IX will tend to decrease the FMcoupling, contrary to rather common believe that the ferro-magnetism in CrX3 is driven by the Hund’s rule coupling on6 Mtot = 2.93 (2.93) and 2.91 (2.93) µB for CrCl3 (CrI3) in the schemes V andVI, respectively.14J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewTable 6. Isotropic exchange interactions in CrCl3 (in meV): results of the Cr3d+Cl3p model in LSDA. b̂ and M stand for the methodsbased on the infinitesimal rotations of, respectively, the xc field and the local spin moments. The xc field was either associated with theon-site spin splitting of the TB Hamiltonian (H) or obtained from the sum rule (sr). ICl is the value of the effective Stoner parameter used inthe 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 Mtot = 3µB (further details are given in the text). TC is correspondingCurie temperature in RPA (in K). The notations of Jk are explained in figure 2(b).Method ICl J1 J2 J3 J4 J5 J6 TCb̂ (H) −0.28 −0.51 0.26 0.34 0.16 0.15 −0.31 21b̂ (sr) −0.28 4.21 0.28 0.26 0.13 0.13 −0.34 49M −0.28 4.36 0.28 0.25 0.16 0.16 −0.32 52M −3.32 0.64 0.25 0.53 0.08 0.10 −0.42 20L0 0 4.73 0.29 0.25 0.17 0.17 −0.32 56Table 7. Isotropic exchange interactions in CrI3 (in meV): results of the Cr3d+ I5p model in LSDA (see table 6 for the notations).Method II J1 J2 J3 J4 J5 J6 TCb̂ (H) 0.31 3.93 0.77 0.75 0.87 0.58 −0.44 97b̂ (sr) 0.31 3.40 0.76 0.71 0.85 0.56 −0.46 88M 0.31 4.29 0.87 0.89 1.03 0.67 −0.41 113M −0.46 1.32 0.75 0.77 0.82 0.54 −0.51 62L0 0 3.21 0.84 0.86 0.96 0.63 −0.46 96the ligand sites, as suggested by the phenomenological GKArules [36–39].The parameters of exchange interactions are summarizedin tables 6 and 7, for CrCl3 and CrI3, respectively. We notethe following. In comparison with the 5o model (section 3.2),the parameters Jk becomes mostly FM, except AFM J6. Theexchange interactions are sensitive to IX. This dependenceis illustrated only for the method M, but very similar beha-vior is observed also for the methods b̂ and m̂. The use of theM method in the combination with IX obtained in the spher-ically averaged scheme IV substantially improves the agree-ment with the experimental data for TC (but not for the spin-wave dispersion, which will be considered in section 3.6). Theexchange parameters obtained in the method m̂ are unreal-istically large [61] and not shown here. Furthermore, in theearlier work [61], the approximate method b̂ was concludedto provide systematically worse description, especially whenconsiders contributions of the ligand states. However, the situ-ation crucially depends to two factors: (i) the proper choicefor IX; (ii) the proper definition of the xc field b̂. For instance,abilities of thismethod can be substantially improved by enfor-cing the sum rule in the definition of b̂ and IX, and taking intoaccount the aspherical contributions to IX, following the defin-ition III.Recently, Ke and Katsnelson [103] reported the exchangeparameters for CrI3, which were also derived using the inverseresponse function (i.e. similar to our method M). In SDFT,they have found (using our notations, in meV): J1 = 6.58,J2 = 0.38, J3 = 1.14, J4 = 0.62, J5 = 0.94, and J6 =−0.14.These parameters are in reasonable agreement with our data,perhaps except J1, which is systematically smaller in ourcase. Nevertheless, our experience shows that in order toobtain reasonable parameters Jk in CrI3 it is absolutely cru-cial to consider the contributions of the ligand states (using thedownfolding technique described in section 2.10). The bareexchange interactions between the Cr 3d states are stronglyAFM and fails to reproduce the correct FM ground state [61].Unfortunately, it is not clear how this problem was tackledby Ke and Katsnelson [103]. Furthermore, our parameters aresensitive to the choice of II, as seen in table 7.The sensitivity of the exchange interactions to the para-meter IX can be regarded as the weak point of the downfold-ing technique. Nevertheless, one can propose an alternativeoption, which can be viewed as an extension of the methodM for the FM insulators and half-metals. In the future, we willcall it the L0 scheme. The crucial observation is that MX isinduced by the hybridization between the Cr 3d and ligandp states (which is totally in line with the general idea of thedownfolding technique, considered in section 2.10). Then, itis reasonable to enforce BX = 0 and find the remaining para-meter of the xc field, BCr, from the equation M⃗= R̂↑↓B⃗. ThisBCr induces the magnetic moments on the Cr sites as well asthe ligand sites. Then, the parameter BCr can be chosen toreproduce Mtot = 3µB. For CrCl3 (CrI3) in LSDA, such pro-cedure yieldsMCr = 3.18 (3.40) µB andMX =−0.06 (−0.13)µB, which are pretty close to the values of spin magneticmoments in the LSDA ground state. The good aspect of thisapproximation is that IX = 0 and the exchange interactions aresolely determined by the 1st term of equation (37) (but withthe redefinedMCr for BX = 0). The corresponding values of Jkand TC are also listed in tables 6 and 7. In comparison with theM scheme, these parameters somewhat overestimate TC butotherwise fall within the range of typical estimates for Jk. Wewill continue to use this L0 scheme in the future. It is espe-cially good for semiquantitative estimates, which would illus-trate the basic trend in the behavior of interatomic exchangeinteractions. The 2nd term in equation (37) also vanishes ifQ̂↑↓LL = 0, meaning that there is no ligand–ligand interactions.15J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewIn this sense, there is certain similarity with the downfoldingmethod proposed by Logemann et al [78], but reformulatedfor the scheme M.3.3.2. Ligand states in correlated Cr 3d band. In thissection, we try to extend the correlated 5omodel, by expandingits Wannier basis W into the pseudo-atomic Wannier basis, A,of the more general model, containing Cr 3d as well as the lig-and p bands. Namely, we still consider only ten Cr 3d bands(for each spin), but transform |C⟩ in equation (12) from thebasis W to the basis A: |CW⟩ → |CA⟩= T̂|CW⟩, where T̂ is thetransformation matrix. Our intension here is to consider expli-citly theX p states and all the contributions of these states to theexchange interactions in the 5o model. The magnetic momentsin the A basis areMCr = 2.63 (2.56) µB andMX = 0.12 (0.15)µB for CrCl3 (CrI3). In the 5o model, the only occupied ↑-spin t2g band provides three electrons, which are distributedbetween one Cr and three ligand atoms. Therefore, the Cr andligand moments in the t2g band are ferromagnetically coupled(see also figure 4), while in order to obtain the opposite polar-ization of these states, it is essential to consider a more generalmodel, which would explicitly include the ligand p band.Without the ligand bands, the inversion of the matrix R̂↑↓in the A basis becomes unstable as it contains small matrixelements associated with the ligand sites. While for calcula-tions of the exchange interactions for the given IX such inver-sion can be largely avoided using equations (38) and (39), theinverse matrix is still needed to evaluate IX using the sumrule. Nevertheless, one can still try to estimate Jk using themethod L0, where it is sufficient to invert only the subblockof R̂↑↓ in the subspace of the Cr sites. It yields the magneticmoments: MCr = 2.64 (2.59) µB and MX = 0.12 (0.14) µB forCrCl3 (CrI3), which are close to the ones reported above.The corresponding exchange interactions and TC for the5o model reformulated in the pseudo-atomic basis, are clearlyoverestimated (see table 8, the rows denoted 5o). However, thisis quite in line with general understanding. The interactionsmediated by the tails of theWannier functions spreading to theligand sites and transferring there the FMmagnetization can beviewed as an analog of the direct Heisenberg exchange [1]. Inorder to evaluate these direct exchange contributions numer-ically, one typically performs the six-dimensional integrationin the real space [40, 104, 105]. Nevertheless, the downfold-ing method suggests how these contributions can be natur-ally evaluated within SDFT. The bare direct exchange integ-rals are typically large and need to be additionally scaled, inthe spirit of cRPA, to account for the screening caused byother bands [40, 105]. The same situation is here: the attemptto include the ligand states into 5o model only worsen thedescription of exchange interactions. The discrepancy can beresolved by considering more general model, which wouldexplicitly include the contributions of the ligand p band.3.4. Merging correlated and ligand bandsThe next important step is to merge the correlated Cr 3dbands with the ligand bands in the framework of the CrTable 8. Exchange interactions in CrCl3 and CrI3 (in meV)obtained in the scheme L0 for the 5o model with the ligand states(5o) and after merging this model with the ligand p bands (5o+ p).TC is the corresponding Curie temperature in RPA (in K). Thenotations of Jk are explained in figure 2(b).System J1 J2 J3 J4 J5 J6 TCCrCl3 (5o) 14.27 0.85 1.65 0.42 0.45 0.21 242CrCl3 (5o+ p) 10.81 0.28 0.35 0.18 0.15 −0.27 102CrI3 (5o) 14.03 2.37 3.63 1.89 1.59 0.76 484CrI3 (5o+ p) 6.87 0.72 0.91 1.14 0.63 −0.31 1193d+X p model. The correlated 5o model is formulated inthe Wannier basis for the Cr 3d bands in LDA. After solv-ing this model in the Hartree–Fock approximation, we wantto replace the original Cr 3d bands in the all-electron LDAband structure by these correlated bands in the pseudo-atomicA basis. Since such basis states of the 5o model and the lig-and p bands are orthogonal to other states, such merging isnot unique and an arbitrary energy shift of the Cr 3d and Xp bands relative to each other will not change the magnet-ization (but will change the exchange interactions). Similarproblem arises in the LDA+U method [94], where the rel-ative position of the transition-metal 3d and ligand p statesis controlled by the empirical double-counting term [57, 95].Therefore, we have to make some empirical conjecture aboutthis merging and we choose it from the ‘charge neutral-ity’ condition, by requesting the Fermi energy of the correl-ated model to coincide with the one in LDA. If the systemis insulating, Fermi energy is chosen in the middle of thegap. The electronic structure after such merging is shown infigure 5.We would like to emphasize that this electronic structureis rather artificial and would require different sets of the xcfields 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 ourpurposes. Therefore, in order to evaluate the exchange inter-actions, we again employ the L0 technique. It yields the fol-lowing magnetic moments: MCr = 3.35 (3.49) µB and MX =−0.12 (−0.16) µB for CrCl3 (CrI3). Thus, after adding the X pband, the L0 method nicely capture the antipolarization of theCr 3d and X p states, which is again stronger in CrI3. However,it should be understood that this effect is ‘mimicked’ by thespecific choice of the xc fields, while from the viewpoint ofelectronic structure itself, the X p band remains unpolarizedand the Cr 3d band is polarized as in section 3.3.2. In theother words, in the L0 method, the imperfection of the elec-tronic structure is corrected by the physical choice of the xcfields (though the electronic structure itself was obtained notfor these fields).The corresponding exchange parameters are listed in table 8(and denoted as 5o+p). We note a systematic improvement incomparison with the 5o model with the ligand states (denotedas 5o): all exchange interactions become smaller so as TC.However, this improvement is only partial as these parametersare still too large and clearly overestimate the tendency towardthe ferromagnetism. This demonstrates the complexity of the16J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewFigure 5. Electronic structure obtained by merging the correlated Cr 3d bands, treated in the model Hartree–Fock approximation, anduncorrelated ligand p bands in LDA for CrCl3 (left) and CrI3 (right). The zero energy is in the middle of the band gap.problem: the L0method is only as the starting point of the self-consistent solution, which should take into account the feed-back of correlated Cr 3d bands onto the ligand bands via thexc field and results in the magnetic polarization of the ligandbands. However, such field will inevitably mix up the Cr 3dandX p bands. After that, one has to redefine theWannier basisfor the correlated model and recalculate the model parameters.An alternative solution is LDA+U [94, 95], which considersboth transition-metal and ligand states, and takes into accountthe polarization of the ligand bands. Today there are manyapplications of this method for calculating the interatomicexchange interactions in various transition-metal oxides andother strongly correlated systems [26–28]. However, it shouldbe understood that the LDA+U approach is supplementedwith additional approximations of purely empirical character,such as the form of the double counting as well as the choiceof the Coulomb interaction parameters and the correlated sub-space for the solution of the Hubbard-type model [57].3.5. DM interactionsIn this section, we briefly discuss the behavior of DM inter-actions in CrI3, which are driven by the strong SO couplingof the heavy I atoms (so as other magnetic interactions of therelativistic origin [86]). The DM interactions in CrCl3 are con-siderably weaker [49, 53] and not significant [83]. In addi-tion to dz, two other components of the DM vector, dx anddy, can be computed by rotating the coordinate frame andapplying the same equation (28) (or similar equations for theb̂ or m̂ rotations). The symmetry of CrI3 is such that twoCr sublattices (shown by different colors in figures 2(a) and(b)) are transformed to each other by the spacial inversion.Therefore, all intersublattice interactions will identically van-ish [6, 7]. On the other hand, the interactions within the sub-lattices can exist and for each vector R connecting two Crsites, the interactions in the sublattices I and II are relatedas: dII(R) =−dI(R). The strongest interaction is d3, whichoccurs in the 2nd coordination sphere of the honeycomb plane(together with the isotropic interaction J3). In total, there aresix bonds n connecting the central atom with the atoms inTable 9. Parameters of Dzyaloshinskii–Moriya interactions in CrI3(dxy and dz are in meV and ψ is in degrees) in correlated five-orbitalmodel (5o) and LSDA for the all-electron Cr3d+ I5p model.Method Model |dxy| ψ |dz|b̂ 5o 0.01 39 0.22M 5o 0.02 19 0.20b̂ LSDA 0.08 −2 0.25M LSDA 0.08 0 0.28the 2nd coordination sphere. The corresponding DM vectorsare d3,n = (dxy cos(ϕn+ψ ),dxy sin(ϕn+ψ ),(−1)ndz), whereϕn = (2n+ 1)π6 is the azimuthal angle specifying the direc-tion of the bond n in the xy plane (relative to the bond along ain figure 2(a), corresponding to n= 1) and ψ specifies the dir-ection of the DM vector in the plane relative to this bond. Forthe R3 symmetry, the parameters dxy, dz, and ψ do not dependon n. They are listed in table 9.We note that the schemes b̂ and M provide very consist-ent description for the DM interactions. The small discrep-ancy for ψ in the correlated 5o model is probably due to thefact that the angle ψ is ill-defined when dxy is small. Thereis also surprisingly good agreement between results of thecorrelated 5o model and the Cr 3d+ I5p model in LSDA.However, such agreement is probably fortuitous because themodels are very different7. The DM interactions in CrCl3 areconsiderably weaker (|dz| ∼ 0.02meV [53]). However, this isto be expected and the order of magnitude difference of theDM interactions in CrCl3 and CrI3 well correlates with thestrength of the SO coupling on the ligand sites, which dif-fers by the same order of magnitude: ξCl = 95meV versusξI = 881meV.7 The correlated 5o model includes orbital polarization caused by on-siteCoulomb interactions [62]. On the other hand, the Cr3d+ I5pmodel in LSDAincludes the additional contributions steaming from the magnetically polar-ized I 5p band. Apparently, these two effects are comparable with each other,which explain a good agreement in table 9.17J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewTable 10. Experimental parameters of exchange interactions inCrCl3 and CrI3.Compound Reference J1 J2 J3 J6 dzCrCl3 [83] 2.14 −0.02 0.05 −0.11 0.03CrCl3 [84] 2.12 0.08 −0.16CrI3 [81] 4.52 1.33 0.36 −0.18 0.503.6. Comparison with experimental dataThe experimental parameters of exchange interactions,derived from the inelastic neutron scattering data [81, 83, 84],are summarized in table 10.8 One of the interesting featuresof the experimental spin-wave dispersion in CrI3 is a ∼4meVgap at the Dirac (K) point, indicating at the existence of largeDM interaction dz between 2nd neighbors in the honeycombplane [81]. Quite expectably, no such feature was observedin CrCl3 [83, 84], where this DM interaction is small. For theisotropic interactions, the agreement between theoretical andexperimental data depends on the model and approximationsemployed for treating the Coulomb correlations and the lig-and states, which we will discuss below. The theoretical DMinteraction in CrI3 appears to be underestimated by factor two,both in the correlated 5o model and LSDA for the Cr3d+ I5pmodel.The theoretical spin-wave dispersion in CrCl3 and CrI3,calculated for some representative sets of parameters, is plot-ted in figures 6 and 7, respectively, in comparison with theexperimental data. First, let us consider results of correlated5o model (shown in the panels a and b of these figures),where there are no additional complications caused by the lig-and states. Within this 5o model, the approximate scheme b̂severely underestimates the FM interactions, making the FMground state unstable both in CrCl3 and CrI3. The situation isdramatically improved when we use the ‘exact’ scheme M: itstrengthens the FM interactions and stabilizes the FM groundstate in both the compounds. In CrCl3, the dispersion is over-estimated by factor two, while in CrI3 we obtain an overallfair agreement with the experimental data, except for the bandsplitting in the point K, which is underestimated by factortwo. Nevertheless, the situation changes significantly whenwe explicitly consider the contributions of the ligand states inthe framework of the all-electron Cr3d+Xp model in LSDA(panels c and d). Now, the FM state becomes stable already inthe scheme b̂. The spin-wave dispersion in CrCl3 remains over-estimated by factor two. The agreement with the experimentaldata is better for CrI3, again except for the band splitting inthe point K. When we turn to the M scheme, which is expec-ted to be more accurate, the theoretical spin-wave dispersionin CrCl3 is substantially reduced, so that the lowest experi-mental band is reproduced pretty well. However, the behaviorof theoretical upper bands is not satisfactory. This is related tothe fact that the theoretical J1 = 0.64meV (table 6) is under-estimated by factor three. A similar situation occurs in CrI3,8 In order to be consistent with our definition of the spin model, equation (1),the experimental parameters have been multiplied by S2 = (3/2)2.where the theoretical J1 = 1.32meV (table 7) is underestim-ated by the same amount, which worsens the agreement withthe experimental data. The problem is partially related to ourchoice of IX, which appears to bemore negative for the schemeM due to the additional spherical approximation, as discussedin section 3.3.1. For instance, if one uses the same IX as in thescheme b̂ (or the L0 scheme, where we do not need any IX),one could get a better agreement with the experimental data forthe upper bands (but not for the dispersion in the entire energyregion).Regarding the value of dz and the band gap, first we notethat there is an alternative mechanism of opening the bandgap, related to long-range isotropic exchange interactions bey-ond the 6th coordination sphere [103]. These interactions weretaken into account also in our calculations.9 However, theireffect is not particularly strong. Therefore, the band gap ismainly controlled by dz and the discrepancy with the exper-imental data is caused by the underestimation of this inter-action in the theoretical calculations, as was also reported byKvashnin et al [49] and Olsen [106]. In this respect, we notethat by explicitly considering the contributions of the ligandstates in the correlated 5o model (section 3.3.2) and using theL0 scheme for the evaluation of the exchange parameters, wecan easily get |dz| as large as 1.10meV, which exceeds thevalue obtained in the regular correlated 5omodel by factor fiveand the experimental value by factor two [81]. Such enhance-ment is caused by the additional contribution of the basis func-tions, which explicitly takes into account the effect of theheavy I atoms. Nevertheless, when we try to combine thiseffect with another contribution stemming from the I 5p banditself (employing for these purposes the merging of the Cr 3dand I 5p bands, as discussed in section 3.4), again in the frame-work of the L0 scheme, the DM interaction is reduced till|dz|= 0.16meV. Such strong reduction is apparently causedby the cancellation of contributions in the Cr 3d and I 5p bands.Since the I 5p states in the Cr 3d and I 5p bands are anti-polarized, the fact of the cancellation itself is not surprising.However, such analysis clearly demonstrates the fragility ofthe situation, where the value of dz appears to depend on thedelicate balance of two large contributions arising from theCr 3d and I 5p bands. The merging of these two bands con-sidered in section 3.4 was probably too crude to explain theexperimental situation.3.7. Brief summaryTo summarize this section, we would like to stress again themain points:• The minimal model, which captures the magnetic proper-ties of CrX3, is the 5o model, constructed for the magneticCr 3d bands near the Fermi level. The main advantage ofthis model is the simplicity: in this case there is simply no9 We have considered all the interactions within the coordination sphere ofabout 16 (20)Å in CrCl3 (CrI3). Figure 2(b) shows only representative para-meters up to 6th coordination sphere.18J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewFigure 6. Theoretical and experimental spin-wave dispersion for CrCl3: results of correlated five-orbital model with the isotropicparameters (J) obtained in the schemes b̂ (a) andM (b); results of all-electron Cr3d+Cl3p model in LSDA with the parameters obtained inthe schemes b̂ (c) and M (d) (data rows 2 and 4 in table 6). The calculations are performed for the hexagonal cell, where six branches ofω(q) correspond to six magnetic Cr sublattices.Figure 7. Theoretical and experimental spin-wave dispersion for CrI3: results of correlated five-orbital model with the isotropic only (J) aswell as both isotropic and DM (J+ dz) parameters obtained in the schemes b̂ (a) and M (b); results of all-electron Cr Cr3d+ I5p model inLSDA with the parameters obtained in the schemes b̂ (c) and M (d) (data rows 2 and 4 in table 7). The calculations are performed for thehexagonal cell, where six branches of ω(q) correspond to six magnetic Cr sublattices.19J. Phys.: Condens. Matter 36 (2024) 223001 Topical Reviewcontributions associated with the ligand states and all calcu-lations of the exchange interactions become pretty straight-forward. Nevertheless, it is absolutely essential to use forthese purposes the ‘exact’ approach dealing with rotationsof magnetic momentsM. The approximate scheme, dealingwith rotations of the xc fields b̂, severely underestimate theFM 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 Cr3d+Xp model. However, the exchange inter-actions in this case depend on the strength of the effect-ive Stoner coupling IX on the ligand atoms and the properdefinition of these Stoner parameters is still not completelyresolved problem. As soon as the xc field is corrected to sat-isfy the sum rules for the given magnetization, the approx-imate b̂ scheme works reasonably well. The main issue inthis context is even not the differences between the schemesb̂ andM, but how to properly define IX within each scheme.As an alternative solution, we have proposed the L0 method,which makes the related to IX contributions inactive. Thismethod can be used for FM insulators and half-metallicmaterials.• Another open question is how to properly merge the correl-ated Cr 3d bands and ligand X p bands. The exchange inter-actions can crucially depend on details of such merging. Themethod considered in section 3.4 is probably only the firststep in this direction.4. Half-metallic ferromagnetsThe half-metallicity is basically the peculiar type of the elec-tronic structure where one spin channel is metallic whileanother one is semiconducting [107]. Most of such materi-als are ferro- or ferrimagnets. The fully compensated half-metallic ferrimagnets, where 100% spin polarization of theconduction electrons coexists with zero net magnetization,have been also proposed [108]. In this section, we furtherexplore abilities of the linear response theories for the ana-lysis of interatomic exchange interactions in this type ofmater-ials with the emphasis on the Coulomb correlations and theligand states. We consider two such examples: the canonicalCrO2 [109] and more recent Co3Sn2S2, which has attracted agreat deal of attention due to the large anomalous Hall effectand other intriguing properties [110–115].Although the electronic structure is metallic, the responsetensor R̂↑↓ involves the transitions between occupied andempty states with opposite projections of spins. For the half-metallic compounds, these transitions are gaped. In such situ-ation, fine details of the electronic structure near the Fermilevel do not play a primary role in the behavior of interatomicexchange interactions, as it could be expected, for instance, inthe RKKY theory for the regular metals [8, 9].We start our analysis of interatomic exchange interactionswith either LSDA or GGA picture. It should be noted that thealternative point of view on the electronic structure of CrO2and other rutile oxides is based on the LDA+U concept [116,117]. Nevertheless, the situation is probably disputable asFigure 8. (a) Fragment of the crystal structure of CrO2, illustratingthe arrangement of the CrO6 octahedra. Reprinted figure withpermission from [61], Copyright (2021) by the American PhysicalSociety. (b) The lattice of Cr atoms with the notations of exchangeinteractions. Reprinted figure with permission from [61], Copyright(2021) by the American Physical Society. (c) Densities of states(DOS) for CrO2 in LDA (top) and LSDA for the ferromagnetic state(bottom). Shaded areas show partial contributions of the Cr 3dstates. The Fermi level is at zero energy.many experimental data of CrO2 are described reasonablywell already on the level of LSDA with no clear indicationsof strong correlations of the Hubbard type [118]. Moreover,the dynamic correlations are known to play a very importantrole in the half-metallic materials [119] and can substantiallyrevise the picture based on the static LDA+U approach [120].We will illustrate this idea by considering the behavior ofinteratomic exchange interactions in the case of CrO2.4.1. CrO2The half-metallic ferromagnetism is not very common in stoi-chiometric transition-metal oxides. Nevertheless, there areexceptions and CrO2 is one of them, which is widely con-sidered in various applications related to spintronics andmagnetic recording [121]. One of the limitations of CrO2from this practical point of view is the relatively lowTC ∼ 390K [121].CrO2 crystallizes in the rutile structure (the space groupP42/mnm, figure 8(a)) [122]. The exchange interactionsremain sizable up to at least the 8th coordination sphere.Moreover, since the rutile structure is nonsymmorphic, thereare two types of interactions, J7 and J8, which are denoted bysuperscripts > and <, as explained in figure 8(b). The calcu-lations are performed for the experimental parameters of thecrystal structure [122] on the mesh of the 10× 10× 16 pointsboth for k and q.The LDA band structure is featured by three separatedbands: O2p in the occupied part, Cr t2g near the Fermi level,and Creg in the unoccupied part. Therefore, one can considertwo types of correlated models: the 3o model for the Cr t2gbands and 5o model both for Cr t2g and Cr eg bands. In theformer case, the on-site Coulomb and exchange interactionscan be described in terms of two Kanamori parameters [123]:20J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewFigure 9. Densities of states (DOS) for CrO2 in the FM state: (a)Hartree–Fock approximation for the three-orbital model and (b) thesame for the five-orbital model. The Fermi level is at zero energy.the intraorbital Coulomb interaction U and the exchange inter-action J , which can be evaluated within cRPA as 2.84 eV and0.70 eV, respectively [120]. The third parameter of interorbitalCoulomb interaction can be obtained from these two as U ′ =U − 2J [123]. The parameters of 5o model are specified byU= 1.98 eV, J= 0.94 eV, and B= 0.09 eV (see appendix A).The corresponding DOSs, in the Hartree–Fock approximation,are shown in figure 9. Alternatively, one can consider the all-electron Cr3d+O2pmodel in LSDA. The half-metallic char-acter of the electronic structure is evident from the LSDADOSs (figure 8(c)). The Coulomb interactions additionallysplit the ↑-spin t2g band, placing the Fermi energy inside apseudogap,10 and shift the ↓-spin states to the higher energyregion, away from the Fermi level.The parameters of exchange interactions are summarizedin table 11. First, there is a large difference in the parametersobtained in the schemes b̂ and M, especially in the Hartree–Fock approximations for correlated 3o and 5o models, wherethe more accurate method M strengthens the FM interactionspractically in all the bonds. The situation is more modest inLSDA for the all-electron Cr3d+O2p model: when goingfrom the method b̂ to the method M, the nearest FM inter-actions J1 and J2 increase only partially, while other longer-range interactions tend to decrease and become more AFM.In LSDA, the Cr and O atoms are antipolarized due tothe joint effect of the hybridization and the intraatomic spinsplitting, similar to CrX3. The corresponding spin magneticmoments areMCr = 2.15µB andMO =−0.08µB. The effectiveStoner parameters were calculated using the definitions III andIV, which should be used in the combination with the meth-ods b̂ and M, respectively (see section 3.3.1). As expected,10 In the rutile structure, three t2g orbitals on the Cr sites belong to three dif-ferent irreducible one-dimensional representations of the point groupmmm=D2h, meaning that the local (site-diagonal) part of the model Hamiltonian willbe diagonal with respect to the t2g orbital indices. In such situation, the on-siteCoulomb repulsion will tend to occupy two t2g levels with lower crystal-fieldenergies and empty the remaining one. On the other hand, there is a stronghopping operating between ‘occupied’ and ‘empty’ orbitals of the differentCr sites, which leads to the formation of the pseudogap instead of the realgap. Further details can be found in [120].ICr = 0.87 eV practically does not depend on the definition.On the other hand, the definitions III and IV yield different val-ues of IO:−0.78 eV and−0.48 eV, respectively. Nevertheless,the parameters are pretty small and the exchange interactionsin CrO2 appears to be less sensitive to the ligand states (atleast, in comparison with CrX3) [61]. For instance, the schemeL0, where we enforce IO = 0, provides basically the sameset of parameters as the M scheme with IO =−0.48 eV (seetable 11).The spin-wave dispersions calculated with different sets ofparameters for the 3o and 5o models is plotted in figure 10 (seeappendix B for details). The nonvanishing xx and zz elementsof the spin-wave stiffness tensor D̂= [Dαβ ],ωA (q)≈∑α,βDαβqαqβ ,describing the dispersion of the acoustic (A) spin-wave branchnear the Γ point, are summarized in table 11. The experi-mental estimates for the averaged spin-wave stiffness in CrO2vary from 70 to 150meV·Å2 [125, 126]. The Curie temperat-ure is about 390K [121]. Then, the LSDA values of TC seemto be overestimated, though the situation is rather subtle. Onthe one hand, the theoretical values of the spin-wave stiff-ness are within the experimental scatter. On the other hand,the exchange parameters derived for ground state are not sup-posed to reproduce TC for this metallic system and more soph-isticated methods, considering the temperature dependenceof the exchange interactions, may be needed [22, 127, 128].Anyways, in the light of well-known limitations of LSDA forthe transition-metal oxides [94], such moderate disagreementwould not be surprising.Even more interesting situation is realized in correlated3o and 5o models. At first glance, the approximate schemeb̂ provides a much better description for the spin-wave stiff-ness and TC in comparison with the scheme M, which is sup-posed to be more accurate. However, this is true only in theHartree–Fock approximation, which has serious limitationsfor the metallic systems. If one goes beyond the Hartree–Fock approximation, the situation can change dramatically.Particularly, as expected for the half-metallic systems [119],the dynamic correlations lead to a strong redistribution of theelectronic states. These effects can be treated in the frame-work of dynamical mean-field theory (DMFT) [129]. Thelatter can be regarded as an extension of SDFT in whichthe KS potential is replaced by the frequency-dependentself-energy Σ↑,↓(ω) [130]. This method has been appliedto CrO2 in [120], using the same correlated 3o model. Theexchange interactions were evaluated basically in the b̂ schemeusing equation (3), by considering the infinitesimal rota-tions of the frequency-dependent xc field ∆Σ(ω) = Σ↑(ω)−Σ↓(ω) [131]. The dynamic correlations substantially reduce∆Σ(ω), resulting in much stronger AFM interactions J7 andJ8, which overcome the effect of FM interactions J1 and J2,making the FM state unstable. The corresponding spin-wavedispersion is also plotted in figure 10, where the negative fre-quencies ω(q) along the directions Γ–X, Γ–M, and Γ–Mmeanthat the theoretical ground state should be in one of these21J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewTable 11. Parameters of exchange interactions in CrO2 (in meV), as obtained in correlated three- and five-orbital models (respectively, 3oand 5o), and in LSDA for the all-electron Cr3d+O2p model using the schemes b̂ andM (corresponding to the infinitesimal rotations of thexc fields and local spin moments, respectively). In the L0 scheme, xc field was redefined to enforce IO = 0. The notations of exchangeparameters are explained in figure 8(b). Dxx = Dyy and Dzz are nonvanishing elements of the spin-stiffness tensor (in meV·Å2). TC is theCurie temperature in RPA (in K).3o 5o LSDAb̂ M b̂ M b̂ M L0J1 9.18 15.76 5.24 17.47 30.87 36.86 37.63J2 10.94 23.15 11.47 31.17 21.39 24.80 25.32J3 1.12 1.55 0.33 0.44 3.04 2.54 2.57J4 0.84 1.51 0.71 1.54 1.49 1.22 1.20J5 −0.34 −0.38 −0.12 0.04 −0.79 −1.72 −1.83J6 −1.92 −1.78 −1.75 −1.73 −3.58 −4.96 −5.14J>7 −3.41 −3.36 −2.50 −2.79 −6.25 −8.39 −8.61J<7 −1.09 −1.77 −0.38 −0.67 −2.09 −3.00 −3.11J>8 −0.02 0.19 0.10 0.62 −1.50 −1.98 −2.08J<8 −0.52 −0.64 −0.39 −0.30 −0.53 −0.66 −0.73Dxx 37 258 113 532 123 113 109Dzz 11 168 39 362 161 159 152TC 310 1026 430 1563 868 880 875Figure 10. Theoretical spin-wave dispersion for CrO2 with theexchange parameters obtained in the Hartree–Fock approximationfor correlated three-orbital (3o) and five-orbital (5o) models, byrotating either the xc field (b̂) or the local magnetic moments (M).Σ̂(ω) denotes results of [120] (the extension of the scheme b̂ forcorrelated three-orbital model, where the scalar xc field wasreplaced by the frequency-dependent self-energy evaluated withinthe dynamical mean-field theory). The notations of thehigh-symmetry points of the Brillouin zone are taken from [124].points. Therefore, it was concluded in [120] that in order toreproduce the FMground state of CrO2, it is essential to extendthe 3o model, by adding new ingredients such as the Cr egand O 2p bands. For instance, the spin-wave stiffness is sys-tematically larger in the 5o model (table 11), which makesthe FM state more stable. Nevertheless, the present analysisalso suggests that the problem may be not in the 3o modelitself, but in additional approximations underlying the schemeb̂ for the exchange interactions. The methodM systematicallyimproves the stability of the FM ground state in CrO2. Thistendency may be overestimated in the Hartree–Fock approx-imation. However, a more rigorous treatment of correlationeffects in the framework of DMFT, in the combination withthe method M, could possibly bring the situation to a betteragreement with the experimental data.The DM interactions can take place between atoms in dif-ferent Cr sublattices. The strongest ones are realized in the2nd coordination sphere (in the combination with J2, as shownin figure 8(b). If ϵ= 1√2a2+c2(±a,±a,±c) are the directionsof such bonds, the corresponding to them DM vectors havethe following form d= dϵz[ϵ×nz] (i.e. the vectors lie in thexy plane, are perpendicular to the bonds, and have oppos-ite signs for the bonds in the directions ±z). The parameterd can be estimated as 0.85 and 1.38meV, thus yielding d=(±0.23,±0.23,0) and (±0.37,±0.37,0)meV in the Hartree–Fock approximation for correlated 5o model and all-electronLSDA, respectively (employing the scheme M in both cases).The interesting point is that the DM interactions are expec-ted to be strong (and comparable with the ones in CrI3) eventhough CrO2 does not contain heavy elements. This is basic-ally the consequence of two factors: (i) the partial filling of thet2g states, opening the room for unquenched orbital magnetiz-ation; (ii) half-metallic character of the electronic structure,giving rise to the spin-current mechanism of the DM interac-tions in the ↑-spin t2g band [58, 70]. Nevertheless, d is muchsmaller than the isotropic interaction J2 operating in the samebonds.4.2. Co3Sn2S2Co3Sn2S2 is a ferromagnet in which small spontaneous mag-netization (about 0.3µB per Co atom) coexists with relat-ively high TC = 177K [132]. It crystallizes in the rhombo-hedral R3m structure, hosting the kagome lattice of Co ions(figures 11(a)–(c)) [132]. The S atoms sit on the top of the22J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewFigure 11. (a)–(c) Fragments of the crystal structure of Co3Sn2S2. (a) Top view on the kagome Cr layer surrounded by Sn and S atoms. (b)Side view on the same layer. (c) Relative arrangement of adjacent kagome layers. (d)–(g) Main exchange interactions. Reprinted figure withpermission from [115], Copyright (2022) by the American Physical Society. (d) and (e) Interactions operating in the kagome plane. (f) and(g) Interactions between adjacent planes (top view). The Co atoms located in different planes are denoted by different colors, the same as in(c). The coordination spheres of Co atoms around the origin are denoted by dotted circles. (d) and (f) Interactions, which are the same for allthe bonds in the given coordination sphere. (e) and (g) Interactions, which are characterized by two different values for two types of theinequivalent bonds in the same coordination sphere. The behavior of Jk around two other Co sites in the primitive cell are obtained by thethreefold rotations. (h) GGA density of states (DOS) for Co3Sn2S2 in the ferromagnetic state. Shaded areas show partial contributions of theCo 3d states. The Fermi level is at zero energy.Co3 triangles, forming with them the network of alternatingtetrahedra. There are two inequivalent Sn sites located in (Sn2)and between (Sn1) the kagome planes. Recently, Co3Sn2S2has attracted a lot of attention as a magnetic Weyl semi-metal whose nontrivial topology of the electronic states givesrise to a large anomalous Hall effect [110–112]. Furthermore,Co3Sn2S2 appears to be half-metallic, as was demonstrated inexperimental [113] and theoretical [114, 115] studies.We have chosen Co3Sn2S2 to illustrate the complexitiesone can face in the analysis of interatomic exchange inter-actions for this type of itinerant electron systems. At thepresent stage, we have in our disposition only all-electronCo3d+Sn5p+S3p model, which was constructed in [115]within GGA [133], by employing the Vienna ab initio sim-ulation package [134] for the electronic structure calculationand the maximal localization technique for the Wannier func-tions [135]. The corresponding DOS is shown in figure 11(h):the Fermi level crosses the ↑-spin band and falls inside thegap for the ↓-spin channel. The magnetic moments areMCo =0.33,MSn1 =−0.03,MSn2 =−0.04, andMS = 0.04µB, so thatthe total moment per one unit cell is 1µB. It is believed tobe still an open question how good is GGA for Co3Sn2S2.However, the attempts to construct a more compact model,which would also include the Coulomb correlations, wereso far unsuccessful. Formally, the electronic structure infigure 11(h) can be divided into fully occupied 19 bandsand remaining 8 bands, which are separated by a band gap.However, the construction of the Wannier basis for the uppereight-band manifold is not straightforward, as these Wannierfunctions will probably reside not on single atomic sites [136].A compact model for Co3Sn2S2 was constructed in [137], butusing empirical arguments.Co3Sn2S2 provides an interesting example of itinerant mag-netism. According to the GGA calculations [115], finite rota-tions of spins away from the FM ground state (for instance,by forcing three Co spins to form the umbrella texture) even-tually leads to the collapse of magnetization, so that the sys-tem falls in the nonmagnetic state.11 A similar behavior wasfound in fcc Ni [138, 139] and several Ru-based oxides [140,141]. Moreover, the size of the magnetic moments is expectedto shrink due the temperature disorder, affecting both elec-tronic structure and interatomic exchange interactions [115].Therefore, the bilinear Heisenberg model cannot be defined inthe global sense (for instance, for describing simultaneouslylow-temperature spin-wave dispersion and TC). Nevertheless,it can still be defined locally, for the analysis of local sta-bility of the FM state with respect to the infinitesimal rota-tions of spins. The exchange interactions spread at least up to9th coordination sphere around each Co site, as explained infigures 11(d)–(g). Moreover, some of these interactions, form-ally belonging to the same coordination sphere, can be inequi-valent (as J4 and J ′4 in the plane, J5 and J ′5, and J9 and J ′9between the planes). The calculations have been performed onthe mesh of the 20× 20× 20 points both for k and q.The effective Stoner parameters, calculated by using differ-ent definitions are summarized in table 12. As was discussedin section 3.3, the definitions II and IV should be used in com-bination with the spherically averaged scheme M, while thedefinitions I and III are more suitable for the matrix schemes b̂and m̂. One can clearly see that all the parameters in Co3Sn2S2strongly depend on the definition. This holds even for ICo,which can easily change by about 50%. This is clearly incontrast with other considered compounds, where the valueof I on the magnetic transition-metal site practically did notdepend on the definition. Nevertheless, the exchange interac-tions (37) between the Co sites do not explicitly depend on the11 In the umbrella texture, three Co spins are rotated from the FM axis z by thesame angle such that their projections in the xy plane form the 120◦ texture.23J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewTable 12. The effective Stoner parameters (in eV) obtained usingdefinitions I–IV (as explained in section 3.3.1). Sn2 and Sn1 arelocated, respectively, in and between kagome planes.Definition ICo ISn1 ISn2 ISI 1.30 0.17 −0.10 1.22II 0.98 −0.06 −0.20 0.88III 1.55 −0.65 −0.65 2.84IV 1.45 −3.44 −2.53 3.85Table 13. Parameters of exchange interactions in Co3Sn2S2 (inmeV) calculated using the schemes b̂, m̂, and M (the infinitesimalrotations of, respectively, the xc field, magnetization matrix, andlocal spin moments). In the scheme b̂, the xc field was eitherassociated with the site-diagonal part of the TB Hamiltonian (H) orderived from the sum rule (sr). The roman number in theparentheses stands for the set of the parameters ISn and IS fromtable 12, which was used in the calculations of Jk. In the L0 scheme,xc field was redefined to ensure I = 0 on the sites Sn and S. Thenotations of Jk are explained in figures 11(d)–(g). Dxx = Dyy and Dzzare nonvanishing elements of the spin-stiffness tensor (in meV·Å2).b̂ (H,I) b̂ (sr,III) m̂ (III) M (II) M (IV) L0J1 1.60 2.54 7.71 −0.08 −1.51 −0.09J2 0.05 0.08 −0.14 −0.06 −0.25 −0.06J3 0.09 0.14 0.22 −0.06 −1.11 −0.07J4 0.18 0.22 −0.01 0.34 0.35 0.41J ′4 0.54 0.77 0.77 0.61 0.52 0.73J5 0.42 0.58 0.88 0.83 −0.23 0.98J ′5 0.65 0.90 1.56 0.95 0.92 1.13J6 0.29 0.41 0.64 0.44 0.51 0.53J7 0.10 0.14 0.08 0.12 0.18 0.15J8 0.04 0.06 −0.13 0.10 0.07 0.11J9 −0.13 −0.19 −0.66 −0.27 −0.25 −0.31J ′9 0.07 0.10 0.03 0.02 0.16 0.03Dxx 521 727 383 490 36 582Dzz 246 317 73 310 48 374choice of ICo. The uncertainty with the choice of the paramet-ers ISn and IS, which strongly depend on the way how theyare defined, posses a more serious problem. The parametersISn are mainly negative (except ISn1 in the case I) and tend todestabilize the FM interactions. On the other hand, positive ISstrengthens the ferromagnetism. However, depending on thedefinition, IS can change by factor four, and ISn changes byan order ofmagnitude. This dependence has a strong impact onthe exchange interactions, which are summarized in table 13.As was pointed out before, the attempts to estimate theCurie temperature entirely from Jk, defined via the infinites-imal rotations of spins near the FM ground state, are prettymuch meaningless in the case of Co3Sn2S2, where proper the-ories for TC should consider also the longitudinal change ofthe magnetization [115]. Nevertheless, these Jk can be used toevaluate the spin-wave dispersion and the spin stiffness. Thelatter has been measured experimentally and the nonvanishingelements of the spin-stiffness tensor areDxx = Dyy = 803± 46and Dzz = 237± 13meV·Å2 [142], which can be used forcomparison with theoretical data.First, we note that the results of the scheme b̂ stronglydepend on the definition of the xc field: the enforcement ofthe sum rule for b̂ strengthens the FM interactions, bring-ing the spin stiffness to a good agreement with the experi-mental data. In this case, the exchange parameters are givenmainly by the bare interactions, corresponding to the first termin equation (36). The corrections caused by the ligand statesare small and do not play a decisive role. An interesting situ-ation is realized in the scheme m̂, based on the rigid rota-tions of the magnetization matrix: as expected, the individualparameters are overestimated (see section 2.7). Nevertheless,the interactions are long-ranged and many of them are AFM,resulting in the relatively small spin stiffness. The exchangeinteractions in the scheme M are very sensitive to ISn andIS. If one uses the definition II, where the xc field in theexpression for I is taken from the site-diagonal part of theTB Hamiltonian, ISn and IS are relatively small so that theexchange interactions are given basically by the first term ofequation (37). In this case, J1–J3 are weakly AFM, while otherinteractions are FM (except J9, which is AFM in all con-sidered methods). As the result, the FM state is stable, thoughDxx is somewhat underestimated in comparison with theexperiment.The situation changes dramatically if one uses the definitionIV, where the xc field is taken from the sum rule, which sub-stantially strengthens both ISn and IS. Apparently, negativeISn plays a dominant role and is responsible for the AFM char-acter of interactions in some of the bonds: particularly, J1–J3become prominently AFM and J5 changes from strongly FMto AFM. Then, the FM state become unstable and the spinstend to form the 120◦ in the xy plane (at least, on the mean-field level). The corresponding spin stiffness is strongly under-estimated compared to the experimental data. This is clearlyan artifact of the model analysis, which is probably caused byunrealistic estimate for ISn. At least, the conclusion does notseem to be consistent with the brute-force GGA calculations,where the 120◦ alignment of spins in the xy plane leads to thecollapse of the magnetic state [115].The L0 approach allows us to eliminate the dependence ofJk on ISn and IS, assuming that all magnetic moments areinduces by BCo, while BSn = BS = 0 (so as ISn and IS). Thisresults only in a small change of magnetic moments in com-parison to the plain GGA values reported above: MCo = 0.36,MSn1 =−0.05, MSn2 =−0.07, and MS = 0.02µB (thus, Mtotremains equal to 1µB). The exchange interactions in these casereminiscent the ones obtained in the schemeM using the defin-ition II for ISn and IS with somewhat stronger tendency towardthe ferromagnetism, which is reflected in larger values of Dxxand Dzz.In principle, the schemes b̂ and L0 both provide a reason-ably good agreement with the experimental data for the spin-stiffness constants. Nevertheless, the behavior of exchangeinteractions is quite different. In the scheme b̂, the strongestFM interaction is J1. The longer-range interactions operatingpractically at the same distance in the 4th and 5th coordinationspheres, in and between the kagome planes (see figures 11(e)and (g)), are also sizable but smaller than J1. In the schemeL0, J1 is weakly AFM, while the strongest FM interactions24J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewFigure 12. Theoretical spin-wave dispersion for Co3Sn2S2 with theexchange parameters derived in the schemes b̂ and L0. Thenotations are taken from [143] for the hexagonal unit cell. Reprintedfigure with permission from [115], Copyright (2022) by theAmerican Physical Society.take place in the 4th and 5th coordination spheres. The corres-ponding spin-wave dispersions are plotted in figure 12. In thelong wavelength limit q→ 0, two methods provide very sim-ilar description, while the main difference is in the positionof optical branches in the higher-energy region. The resultsof recent inelastic neutron scattering data (which were never-theless limited by the behavior of the acoustic branch in thevicinity of q= 0) where interpreted in terms of four interac-tions (in our notations): J2 =−0.32, J3 = 0.64, J̄4 = 1.66, andJ̄5 = 0.09meV [143], where J̄k = 13 (Jk+ 2J ′k) stands for theaveraged exchange interactions in the 4th and 5th coordin-ation spheres.12 There is certain similarity with the pictureprovided by the L0 scheme: at least J1 is negligibly small,J2 is AFM, and the main FM interactions occurs in the nextcoordination spheres. Nevertheless, there are also the differ-ences. Especially, experimental J̄5 appears to be smaller thanJ̄4, while theoretical parameters are at least comparable (seetable 13). According to the theoretical analysis, the values ofinequivalent parameters Jk and J ′k can be very different andthis difference was not considered in the fitting of the experi-mental spin-wave spectrum. On the other hand, it is not clearwhether it can resolve the discrepancy between theoretical andexperimental data. It is also unclear whether the experimentaldata available only in the vicinity of q= 0 are enough to makea decisive conclusion about the complexity of the exchangeinteractions in Co3Sn2S2.The Co3Sn2S2 structure has three Co sublattices (which aretransformed to each other by the threefold rotations). The spa-cial 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 arepermitted by the R3m space group. The strongest ones takeplace between nearest neighbors in the kagome plane. Theyare explained in figure 13. Using the convention, where the12 In order to be consistent with our definition of the spin model, equation (1),the exchange parameters reported in [143] were additionally multiplied by2S2 = 118.Figure 13. Dzyaloshinskii–Moriya interactions between nearestneighbors in Co3Sn2S2. The Co atoms belonging to one of the threesublattices are denoted by numbers. Blue vectors show thedirections of the bonds. Green vectors are projections of the DMvectors onto the kagome plane for each such bond. The dzcomponent is perpendicular to this plane and has the same directionfor all such bonds.directions of the bonds (ϵ) in the triangles can be transformedto each other by the threefold rotations, the DM vectors canbe presented in the form: d= dxy[ϵ×nz] + dznz. If the FMmoment is parallel to z (the experimental situation), dz doesnot contribute to the energy change, while dxy is responsiblefor the canting of spins and tends to form the umbrella texture.If ek = (θ sin π k3 ,θ cosπ k3 ,1−θ22 ) are the directions of spinsin the triangle (k= 1, 2, or 3, as explained in figure 13), theenergy gain due to the DM interaction is δEDM =−√3dxyθ(per one Co site). The corresponding energy loss due to theisotropic exchange is δEH = 34J1,23θ2, where J1,23 is the effect-ive interaction of the atom in the 1st sublattice with the atomsof the 2nd and 3rd sublattices in all the coordination spheres.The values of the parameters obtained in the scheme L0are J1,23 = 6.04meV, dxy = 0.20meV, and dz =−0.29meV.13Thus, the angle θ = 2dxy√3J1,23can be estimates as 2◦, which per-fectly agrees with θ ∼ 2◦ obtained in the brute-force GGA cal-culations with the SO coupling [115].5. Orthorhombic perovskite manganitesPerovskite manganites AMnO3 (where A is the trivalent rareearth or alkaline earth element) have attracted a great dealof attention. Most of them crystallize in the orthorhombicPbnm structure (figure 14(a)). The magnetic Mn3+ ions in theoctahedral environment accommodate four electrons, three ofwhich occupy the t2g states, while the remaining one residesin the doubly degenerate manifold of the eg states. Therefore,the system is subjected to the Jahn–Teller distortion, resultingin the peculiar orbital ordering (the alternation of occupied egorbitals, as schematically shown in figure 14(e)).13 The DM parameters are not sensitive to the definition and very similar val-ues are obtained in the schemeM, taking ISn and IS from the set IV. However,the same procedure would give us J1,23 =−4.42meV, meaning that the FMstate is unstable, as was already explained in the main text.25J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewFigure 14. (a) Fragment of distorted perovskite structure, illustrating the arrangement of the MnO6 octahedra. Reprinted figure withpermission from [61], Copyright (2021) by the American Physical Society. (b) Main exchange interactions in the orthorhombic cell.Reprinted figure with permission from [61], Copyright (2021) by the American Physical Society. (c) Main exchange interactions in theorthorhombic ab plane. Reprinted figure with permission from [61], Copyright (2021) by the American Physical Society. (d) AFMalignment of the types A and E. (e) Schematic view on the orbital ordering underlaying the behavior of exchange interactions in the abplane. Atoms forming four Mn sublattices are denoted by numbers. Reproduced from [42]. CC BY 4.0. (f) and (g) Densities of states (DOS)for LaMnO3 and HoMnO3 in LDA (top) and LSDA for the A-type AFM state (bottom). Shaded areas show partial contributions of the Mn3d states (from two magnetic sublattices with the spins ↑ and ↓ in the case of A-type AFM order). The Fermi level is at zero energy (themiddle of the band gap in the insulating phase).Table 14. Selected parameters of the orthorhombic structure for LaMnO3 [154] and HoMnO3 [149]: the orthorhombic lattice parameters a,b, and c; the Mn–O bondlengths (dMn–O) in the ab plane (the first two values) and along c (the third value); and the Mn–O–Mn angles(∠Mn–O–Mn) in the ab plane (first value) and between the planes (second value).Compound a (Å) b (Å) c (Å) dMn–O (Å) ∠Mn–O–Mn (◦)LaMnO3 5.532 5.742 7.668 1.906, 2.118, 1.959 154, 157HoMnO3 5.257 5.835 7.361 1.905, 2.222, 1.943 144, 142LaMnO3 is the parent material of colossal magnetores-istive oxides [144] and a popular testbed system for study-ing abilities of first-principle electronic structure calcula-tions, especially in reproducing the insulating behavior, thecooperative Jahn–Teller distortion, and associated with itA-type AFM ordering, where the FM coupling within theorthorhombic ab plane coexists with weakly AFM coup-ling between the planes (see figure 14(d)) [145–147]. Theorthorhombic distortion systematically increases with thedecrease of the size of the A3+ ions in the direction La→Ho.At certain point it makes the A-type AFM phase unstable inthe ab plane. For instance, the ground state of TbMnO3 isbelieved to be a spin spiral with the propagation vector q≈(0, 14 ,0) [148], whereas HoMnO3 forms the twofold periodicE-type AFM structure corresponding to q= (0, 12 ,0) [149,150]. These magnetic superstructures breaks the inversionsymmetry, giving rise to rich multiferroic activity [151–153].In this section we consider LaMnO3 and HoMnO3 astwo characteristic example and discuss abilities of the linear-response based techniques for describing the change of theexchange interactions, which would eventually lead to thechange of the magnetic structure from A to E. We use theexperimental crystal-structure parameters reported in [154]and [149] for LaMnO3 and HoMnO3, respectively. The mainparameters are summarized in table 14. Particularly, theMnO6octahedra have two long Mn–O bonds in the ab plane andfour short ones, which is the signature of the Jahn–Teller dis-tortion [155]. The Jahn-Teller distortion practically does notchange when going from LaMnO3 to HoMnO3. On the otherhand, the Mn–O–Mn angles change significantly with muchstronger deviation from the ideal cubic value of 180◦ in thecase of HoMnO3. This change is accompanied by some shrink-ing of the orthorhombic lattice along a and c.The orthorhombic distortion has a profound effect of theelectronic structure in LDA/LSDA (figures 14(f) and (g)),splitting the Mn eg 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 theA-type AFM order [156]. The occupied Mn eg band is loc-ated around −1 eV and well separated from other bands. Thecorresponding distribution of the eg electron density has theform of the orbital ordering, which is schematically shown infigure 14(e).The exchange interactions in LaMnO3 and HoMnO3 arerather complex (see figures 14(b) and (c)) [42]. Particularly, inaddition to the nn interactions in and between the ab planes (J∥1and J⊥1 , respectively), there are several important longer-rangeinteractions, such as: (i) the next-nn interactions between theplanes, J12 and J22; (ii) the 2nd neighbor interactions in theplane, Ja2 and Jb2, operating along orthorhombic axes a andb, respectively; and (iii) the 3rd neighbor interactions in theplane, J13 and J23. These interactions obey certain symmetry26https://creativecommons.org/licenses/by/4.0/J. Phys.: Condens. Matter 36 (2024) 223001 Topical Reviewproperties. For instance, considering Mn site 1 in figure 14(b)as the reference point, J13 and J23 will operate in the bonds±(a,a,0) and±(a,−a,0), respectively. The parameters J13 andJ23 around Mn sites 2, 3, and 4 are obtained by the 180◦ rota-tions of these bonds about a, b, and c combined with the lat-tice shifts by ( a2 ,b2 ,0), (0,0,c2 ), and (a2 ,b2 ,c2 ), respectively. Theinteractions J⊥1 , J12 and J22 control the AFM coupling betweenthe planes, while the formation of the long-periodic magneticstructures in the plane results from the interplay of J∥1 , Ja2, Jb2,J13 and J23. The behavior of J13 and J23 is directly related tothe orbital ordering: these 3rd neighbor interactions can beregarded as the super-superexchange ones, mediated by thestates of intermediate Mn sites. If the lobes of the occupiedeg orbitals are directed toward each other along the bond, as inthe case of J13, the interaction is strong. If the lobes are parallelto each other and perpendicular to the bond, as in the case ofJ23, the interaction is weak. Moreover, two crystallographicallydifferent types of next-nn bonds between the planes result inslightly different values of the parameters J12 and J22 [61].From the viewpoint of the electronic structure, one can con-struct two types of model: the correlated 5o model, wherethe one-electron part is taken from LDA and combined withthe screened Coulomb interactions obtained in cRPA (asexplained in appendix A), and the all electron model withinLSDA, which includes Mn 3d as well as O 2p and A 5d bands.The calculations are performed for the A-type AFM phase onthe mesh of the 8× 8× 6 points, both for k and q, unless it isspecified otherwise.5.1. Correlated 5o modelFirst, we investigate abilities of correlated 5o model: whetherit can reproduce the main tendencies in the behavior ofinteratomic exchange interactions in LaMnO3 and HoMnO3,stabilizing the A-type AFM state in the former case andmaking it unstable with respect to the formation of the spinsuperstructures along the orthorhombic axis b in the latterone. Another important question is how good is the approx-imate scheme b̂, which is frequently used for the analysisof interatomic exchange interactions in these orthorhombicmanganites [42, 145], in comparison with the more accuratescheme M. The averaged on-site parameters of the Coulombrepulsion, intraatomic exchange interaction, and nonspheri-city (see appendix A for definitions) calculated in cRPA are,respectively, 2.15 (2.16), 0.85 (0.85), and 0.09 (0.09) eV forLaMnO3 (HoMnO3) [42]. Thus, the Coulomb U∼ 2 eV is notparticularly strong.14 This finding is very important as it natur-ally explains the existence of the long-range exchange interac-tion, Jb2 and J13, arising in higher orders of Ĥ/U beyond the con-ventional superexchange approximation [2], and responsiblefor the formation of the spin superstructures along b. On the14 This is qualitatively consistent with old estimates, based in the constrainedLSDA [157]. The key point is that Mn eg electrons are efficiently screened byother Mn 3d electrons from other bands.Figure 15. Densities of states (DOS) for the A-type AFM order inLaMnO3 and HoMnO3 as obtained in the Hartree–Fockapproximation for correlated five-orbital model. Zero energy is inthe middle of the band gap.other hand, the relatively small U also means that the strong-coupling limit is hardly satisfied. Therefore, it is reasonable toexpect that the approximate scheme b̂ may experience seriouslimitations for these orthorhombic manganites.The DOSs obtained in the Hartree–Fock approximation forthe A-type AFM order are displayed in figure 15. LaMnO3 andHoMnO3 have similar electronic structure. The only differ-ence is slightly smaller eg bandwidth in the case of HoMnO3.Nevertheless, the behavior exchange interactions appears to bevery different. These interactions are summarized in tables 15and 16 (lines 5o) for LaMnO3 and HoMnO3, respectively.The scheme b̂ yields strong interlayer coupling J⊥ = J⊥1 +2J12 + 2J22 =−10.78 (−10.1)meV for LaMnO3 (HoMnO3),which is even stronger than the intralayer one J∥1 = 2.36(−6.53)meV. This behavior is inconsistent with the experi-mental neutron-scattering data for LaMnO3, indicating thatJ⊥ ∼−4.7meV is weaker than J∥1 = 6.6meV [158, 159].15Furthermore, the longer-range AFM interactions Jb2 and J13 inLaMnO3 are comparable or even stronger than J∥1 , making theexperimental A-type AFM structure unstable.On the contrary, the scheme M systematically improvesthe description of the interatomic exchange interactions.First, the interlayer coupling becomes considerably weaker:J⊥ =−2.35 (−4.34)meV for LaMnO3 (HoMnO3). Then, theitralayer interaction J∥1 in LaMnO3 becomes strongly FM, asexpected for the ‘antiferro’ orbital ordering [160],16 and over-comes the AFM long-range interactions Jb2 and J13. As the res-ult, the experimental A-type AFM phase becomes stable. Thetheoretical TN = 186K, evaluated in RPA, is in fair agreementwith the experimental value of 140K [158, 159]. In HoMnO3,J∥1 is significantly reduced due to the additional buckling ofthe Mn–O–Mn bonds to become comparable with Jb2 and J13.This makes the A-type AFM phase unstable. Regarding thedirection of this instability, it is very important that Ja2 =−0.35meV is much weaker than Jb2 =−1.53meV. Therefore,the propagation vector for the expected theoretical ground15 In order to be consistent with our definition of the spin model, equation (1),the experimental parameters have been multiplied by 2S2 = 8.16 The alternation of the 3x2–r2 and 3y2–r2 orbitals in the ab plane (seefigure 14(e)).27J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewTable 15. Isotropic exchange interactions in LaMnO3 (in meV) as obtained in the correlated five-orbital (5o) model in comparison withLSDA for the all-electron Mn3d+O2p+La5d model. The notations of parameters are explained in figures 14(b) and (c). The AFMinteractions Ja2 and J23 are considerably weaker and not shown here. TN is the Néel temperature evaluated in RPA (in K) and q= (0,qb,0) isthe theoretical ground state propagation vector (in units of reciprocal lattice translations, where q= 0 corresponds to the A-type AFM state).Method J∥1 J⊥1 J12 J22 Jb2 J13 TN qbb̂ in 5o 2.36 −6.64 −0.97 −1.10 −1.00 −3.27 80 0.32M in 5o 17.23 3.51 −1.42 −1.51 −1.13 −3.56 186 0b̂ in LSDA 15.10 7.37 −3.14 −3.25 −1.78 −9.46 131 0M in LSDA 17.09 8.20 −3.24 −3.27 −1.83 −9.21 185 0Table 16. The same as table 15 but for HoMnO3. The LSDA results are obtained in the Mn3d+O2p+Ho5d model.Method J∥1 J⊥1 J12 J22 Jb2 J13 TN qbb̂ in 5o −6.53 −7.58 −0.54 −0.72 −1.44 −2.22 65 0.22M in 5o 4.02 −1.44 −0.63 −0.82 −1.53 −2.36 53 0.26b̂ in LSDA 2.23 −1.64 −1.72 −1.85 −2.57 −6.90 79 0.39M in LSDA 4.15 −1.03 −1.71 −1.81 −2.46 −6.50 79 0.34state is parallel to b, in agreement with the experimental obser-vation. This instability is further enhanced by the AFM inter-actions J13. Nevertheless, in order to reproduce the experi-mental E phase with commensurate q= (0, 12 ,0), it is essentialto consider other ingredients such as the exchange strictionand single-ion anisotropy, which would lock the spin super-structure to the lattice [161–163]. In HoMnO3, such lock-in transition occurs at 29K, while the Néel temperature isabout 41K [149, 150], which is close to theoretical value ofTN = 53K. Other aspects of stability of the E-type AFM phasewill be considered in section 5.4.5.2. All-electron model in LSDAThe all-electron model in LSDA provides an alternativedescription for the exchange interactions. We start with theanalysis of effective Stoner parameters. As was discussed insection 3.3.1, among several possible definitions of the para-meters I, the most relevant seem to be III and IV, whichshould be considered in combination with the methods b̂ andM, respectively. In the A-type AFM phase, the A and apicalO atoms, located between the antiferromagnetically coupledlayers, are nonmagnetic. Therefore, we set I = 0 for them.The remaining parameters for Mn and planar O atoms are lis-ted in table 17. IMn is practically identical for LaMnO3 andHoMnO3, and does not depend on the definition. IO appearsto be more sensitive to the environment and the definition.Nevertheless, the change of IO is rather modest. For bothdefinitions, IO is large and positive, that will additionallystrengthens the FM interactions.The exchange interactions are summarized in tables 15and 16 (lines LSDA). Unlike for the correlated 5o model, themethods b̂ and M in the all-electron LSDA provide a consist-ent description. On the one hand, there is an effect of the O2p band, which is explicitly treated in this model. On the otherhand, the underestimation of the FM interactions in the schemeb̂ is partly compensated by IO. Moreover, the correct definitionof the xc field b̂ using the sum rule is important. For instance,Table 17. The effective Stoner parameters (in eV) for Mn and planarO atoms in LaMnO3 and HoMnO3, obtained using definitions IIIand IV (as explained in section 3.3.1). The A and apical O atoms arenonmagnetic in the A-type AFM phase and not considered here.DefinitionLaMnO3 HoMnO3IMn IO IMn IOIII 0.98 3.71 0.97 4.20IV 0.98 4.35 0.97 4.23the use of b̂ defined via site-diagonal elements of Ĥ↑,↓ under-estimates the FM interactions and makes the experimental A-type AFM state unstable in LaMnO3.The exchange interactions are generally stronger than incorrelated 5o model (except J∥1 ), partly due to the fact thatin these insulating materials the exchange interactions areexpected to increase when the effective Coulomb repulsiondecreases [2]. Furthermore, the oxygen band can also contrib-ute to the exchange interactions [22, 164]. Anyway, the totalinterlayer coupling in LaMnO3, J⊥ =−4.82 (−5.41) in thescheme M (b̂), is consistent with the experimental data [158,159]. J∥1 is strongly FM, which appears to be sufficient to over-come the strong AFM interactions Jb2 and J13, and stabilizethe experimental A-type AFM phase. In HoMnO3, this J∥1 issubstantially weaker, making the A-type AFM phase unstablewith respect to an incommensurate spin structure propagatingalong b. The theoretical TN is in fair agreement with experi-mental data.5.3. DM interactionsAll nn DM interactions in the orthorhombic planes, d∥ij , aretransformed to each other by the symmetry operations ofthe space group Pbnm [145]. The same holds for the inter-plane interactions d⊥ij . Furthermore, since neighboring planesare connected by the mirror reflection, the z (c) componentof d⊥ij is equal to zero. Therefore, there are five parameters28J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewFigure 16. Form of DM interactions operating between nearestneighbors in the orthorhombic structure, in and between the plane.The numbering of Mn sublattices is the same as in figure 14. EachDM vector is attached to its Mn–O–Mn bond, starting from theatoms 2 or 4 for the in-plane interactions and atoms 3 or 4 for theout-of-plane interactions. Reprinted figure with permission from[145], Copyright (1996) by the American Physical Society.Table 18. Parameters of nearest-neighbor DM interactions (in meV)as obtained in the correlated five-orbital (5o) model and all-electronLSDA using the scheme M. The results of [145], employing mixedperturbation theory, are shown for comparison.Compound Model d∥a d∥b d∥c d⊥a d⊥bLaMnO3 [145] 0.44 0.33 0.53 0.45 0.71LaMnO3 LSDA 0.46 0.35 0.63 0.53 0.60HoMnO3 LSDA 0.54 0.38 0.73 0.70 0.47HoMnO3 5o 0.23 0.09 0.28 0.20 0.11describing 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 inter-actions. The corresponding DM vectors, attached to neighbor-ing Mn–O–Mn bonds, are shown in figure 16. The values ofthe parameters, obtained in the scheme M, are summarized intable 18.For LaMnO3 in LSDA we note a good agreement withthe results of [145], employing the mixed perturbation theory,where the rotation of the xc field on one Mn sites was com-bined with the SO coupling on another such site. This meansthat the 5d states of the heavy A atoms, which are located farin the unoccupied part of the spectrum (see figure 14), do notstrongly contribute to the DM interactions. Partly, this may bedue to the fact that the A sites remain nonmagnetic in the A-type AFM phase.The parameters obtained in the correlated 5omodel are gen-erally smaller than in LSDA—similar to what was found forthe isotropic interactions. However, this is not surprising andcan be again explained by the Coulomb U in the denominatorof the superexchange interactions [2].The DM interactions in AMnO3 mainly contribute to thespin canting. The magnetocrystalline anisotropy tends to alignthe spins parallel to the b axis [165–167]. In LaMnO3, theyorder according to the A-type: eµ = (0,−1,0) for µ= 1 and2 and eµ = (0,1,0) for µ= 3 and 4 in figure 16. This per-fect AFM alignment will be further deformed by the DMinteractions, which additionally rotate the spins along a and cby, respectively, ea =−d∥c /(J∥1 − J12 − J22) and ec = d⊥a /(J∥1 +2J12 + 2J22) [62]. The a components will order according to theG-type,17 while c components gives rise to the weak ferromag-netism [165, 166]. Using the obtained parameters of exchangeinteractions, |ea| and |ec| can be estimated as 0.019 and 0.110,respectively. Taking into account that in all-electron LSDAM= 3.59µB, the weak FMmoment can be estimated as 0.4µBper Mn atom (being somewhat larger that the experimental0.18µB, derived from magnetization measurements [168]).An interesting question is whether the multiferroicity asso-ciated with the E-type AFM phase in HoMnO3 can coexistwith the weak ferromagnetism. In the E-type AFM structure,each of the J12 and J22 will contribute to two FM bonds and twoAFM ones. Therefore, these contributions cancel each otherand the spin canting along c will be given by ec =±d⊥a /J∥1(for the sublattices with different spins in the ab plane). Usingthe parameters for the Pbnm structure of HoMnO3, |ec| canbe estimated as 0.169 and 0.050 in the all-electron LSDA andcorrelated 5o model, respectively, where the difference mainlycomes from d⊥a . Nevertheless, this component will repeat thepattern of the E-type AFM phase in each ab plane and, there-fore, there will be no weak ferromagnetism.5.4. Internal instability of the E phaseAn interesting aspect of interatomic exchange interactionsdefined via infinitesimal rotations of spins near some equi-librium magnetic states is that these exchange interactionsdepend on the magnetic state in which they are calculated,reflecting the dependence of the electronic structure on themagnetic state. If it happens, the bilinear spin model (1) isill-defined in the global sense, meaning, for instance, thatthe same set of the exchange parameters cannot be usedfor the analysis of the low-temperature spin-wave dispersionand the magnetic transition temperature, where the electronicstructure is strongly modified by the spin disorder [22, 128].Nevertheless, the model (1) can be defined locally, near eachmagnetic equilibrium. Then, in principle, one can expect somekind of self-organization phenomenon, when the change ofthe electronic structure in certain magnetic state itself can sta-bilize this magnetic state, at least locally. This is what hap-pens at least with some of the exchange interactions in the50% doped manganites, where the zigzag (CE-type) AFMalignment opens the band gap, induces the orbital ordering,17 The AFM coupling between all nearest neighbors, in and between theplanes.29J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewFigure 17. (Left) Inequivalent nearest-neighbor interactions in theab plane of HoMnO3 due to the E-type AFM order. (Right)Densities of states (DOS) for the E-type AFM state in HoMnO3 asobtained in the Hartree–Fock approximation for the correlatedfive-orbital model. Zero energy is in the middle of the band gap.and additionally stabilizes the FM coupling in the zigzagchains [98]. Can we expect a similar behavior in the E phase?Intuitively, the reason for the formation of this non-centrosymmetric E-type AFM structure can be understood asfollows: theAFM interactions Jb2 and J13 tend to align the cornerspins, separated by the orthorhombic translations a and b, anti-ferromagnetically, as shown in figure 17(a). The central Mnatom is located in the inversion center. Then, due to the AFMalignment, the corner atoms should transform to each otherby the symmetry operation ÎT̂ (where the spacial inversion Î iscombined with the time reversal T̂), which is formally compat-ible with the Pbnm space group. However, the same symmetryoperation would make the central Mn atom nonmagnetic. Thisis an eligible scenario from the symmetry point of view, butwould lead to gigantic energy loss, 14IMnM2Mn ∼ 4 eV, causedby the violation of the first Hund’s rule for the atoms, whichare potentially expected to be in the high-spin state. Therefore,the more favorable scenario is to keep the Mn atom mag-netic, but break the inversion symmetry, which is quite com-mon in multiferroic materials, especially those with high-spinions.Then, the next question is whether the E-type AFM orderis compatible with the change of the exchange interactionsinduced by this order.Below we consider results of correlated 5o model (in theschemeM) for HoMnO3. A qualitatively similar behavior hasbeen found in LSDA for the all-electron model. The E-typeAFM order results in the additional narrowing of the t2g andeg bands (figure 17), which is consistent with the decrease ofthe number of FM bonds around each Mn site (2 instead of4 in the A-type AFM phase). The inversion symmetry break-ing makes the nn interactions inequivalent, where J↑↑1 in theFM bond is generally different from J↑↓1 in the AFM one. Inorder 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.99meV andJ↑↓1 = 4.00meV, meaning that the E state corresponds tothe energy maximum and any small rotations of spins willslide the system away from this equilibrium towards a newmagnetic state. Furthermore, there is an intrinsic mechanisminducing these rotations of spins. Indeed, the E-type AFMorder, producing some changes in the electronic structure,can induce the electric polarization even in the centrosym-metric Pbnm structure [161]. Then, it is reasonable to expectthat the same changes in the electronic structure may leadto the appearance of new DM interactions, which wouldotherwise be forbidden in the Pbnm structure. Particularly,we have found finite interactions db2 = (±0.01,0,±0.02)meVand d13 = (±0.01,±0.01,±0.02)meV (being analogs of iso-tropic Jb2 and J13, respectively, where the ± signs depend onthe origin and direction of the bond). Although these interac-tions are not particularly strong, they are sufficient to shift thesystem of spins away from the extremum (maximum) point,so that it will start to relax to a new magnetic equilibrium.Thus, the E-type AFM state in HoMnO3 is intrinsicallyunstable and cannot be stabilized by purely electronic mech-anisms. For these purposes, it is essential to consider theexchange striction or the single-ion anisotropy (or both ofthem) [161–163].5.5. Brief summary• The minimal correlated 5o model provides a consistentdescription for the interatomic exchange interactions inLaMnO3 and HoMnO3, explaining stability of the A-type AFM phase in the former material and the tend-ency toward formation of the spin superstructures alongthe orthorhombic axis b in the case of HoMnO3. The driv-ing force behind this change of the magnetic structure isthe buckling of the Mn–O–Mn bonds, which is strongerin HoMnO3. The use of the scheme M appears to be cru-cial within this 5o model, while the approximate schemeb̂ strongly underestimates the FM interactions and fails toexplain the stability of the A-type AFM phase in LaMnO3.• The all-electron model, explicitly treating the Mn 3d, O2p, and A 5d bands in LSDA provides an alternativedescription for the exchange interactions. Both 5o modeland all-electron LSDA capture the experimental situationpretty well. Why do we have such consistent descrip-tion? Apparently, this is due to the cancellation of manycontributions. For instance, the polarization of the oxy-gen states in LSDA strengthens the FM interactions [42].In addition to them, there are FM superexchange inter-actions, which are controlled by the ratio of the on-siteexchange and Coulomb repulsion, J/U, and expected to bestronger for smaller U [160]. However, these effects arecompensated by stronger AFM interactions, also expectedin the superexchange theory for smallerU [2]. Furthermore,the AFM interactions are additionally stabilized by correla-tions effects beyond the Hartree–Fock approximation [42].• Unlike in the correlated 5o model, the schemes b̂ and Mprovide a consistent description within all-electron LSDA.Nevertheless, the special attention should be paid to thedefinition of the xc field and the choice of the parametersIO. The incorrect choice of these parameters can worsensthe description.30J. Phys.: Condens. Matter 36 (2024) 223001 Topical Review• The inversion symmetry breaking, caused by the E-typeAFM order in HoMnO3, gives rise to the new DM inter-actions, operating across the inversion centers in the Pbnmstructure, which, in the combination with the magnetic-statedependence of the isotropic interactions, act against thisE-type AFM state. Thus, the latter state can be stabilizedonly by extrinsic mechanisms such as the exchange strictionand/or the single-ion anisotropy.6. Other developments6.1. Symmetric anisotropic exchange interactionsThe spin-spiral concept, which assumes that the perturbationof the magnetic ground state can be described in the form of anincommensurate spin spiral, is applicable only for calculationsof isotropic and DM interactions [15]. All these calculationsare based on the generalized Bloch theorem, which allows usto deal with this spin-spiral periodicity [46]. Unfortunately,this theorem is no longer applicable for calculations of theexchange anisotropy or any other symmetric anisotropic inter-action, emerging in the second order of the SO coupling andinvolving spin diagonal as well as off-diagonal elements of thiscoupling. As was pointed out in section 2.2, the DM interac-tions support the spin-spiral propagation, while the symmet-ric anisotropic interactions act against it. Nevertheless, onecan still consider the perturbations caused by the infinites-imal rotations of the xc field and evaluate the 3× 3 exchangetensor Ĵij in the real space, separately for each magnetic bond.Then, 12 (Ĵij− Ĵji), which has only three inequivalent matrixelements, can be related to the DM vector, while the trace-less part of 12 (Ĵij+ Ĵji) gives rise to the exchange anisotropy.This method was successfully implemented in several com-putational packages, which are actively used for calculationsof isotropic as well anisotropic exchange interactions [32, 49–52]. Nevertheless, one should remember that this method haslimitations inherent to the scheme b̂, which is an approxima-tion. Unfortunately, only the scheme b̂ can be easily reformu-lated in the real space so that the exchange interactions canbe calculated separately for each bond. Due to the additionalinversion of the responsematrix in the schemeM, the exchangeparameters in the bond will depend on the matrix elements ofthe response function in other bonds.6.2. Dynamic electron correlationsAs was already pointed out in section 4.1, another interest-ing direction is to go beyond the conventional SDFT by incor-porating the effects of dynamic electron correlations. The lat-ter are typically evaluated in the framework of DMFT [129],which provides a formal extension of the KS equations, wherethe local potential is replaced by also local, but frequency-dependent self-energy [130]. Then, the exchange interac-tions can be still evaluated using the scheme b̂, but withthe frequency-dependent xc field [131]. This approach isespecially important for metallic systems, where the staticHartree–Fock approximation is clearly insufficient and, aslong as the Coulomb interactions are taken into account, theyshould be treated on a more rigorous footing. For instance, thedynamic correlations have a profound effect on the electronicstructure of half-metallic compounds [119], which is reflectedin the behavior of exchange interactions, as was demonstratedfor CrO2 [120, 169]. However, as the most interesting applic-ations of this method are beyond the strong-coupling limit,the scheme b̂ can have serious limitations. Therefore, it canbe important to reformulate the interatomic exchange interac-tions in terms of the inverse response function, as it is done inthe scheme M. In this respect, equation (10) seems to be verygeneral and could be a good starting point for such extension.Another advantage of this DMFT based approach is that itprovides a natural extension for the analysis of temperaturedependence of the exchange interactions [127].7. Summary and conclusionsThe linear response theory becomes a powerful tool for cal-culations of interatomic exchange interactions in various sub-stances. By treating infinitesimal rotations of spins as a per-turbation, it allows us to present the total energy change causedby these rotations in the form of pairwise interactions. Themain purpose of this topical review was to clarify basic prin-ciples of this technique and provide a transparent explanationto more recent developments and controversies related to itspractical realization.The first group of questions is related to the validity ofthe magnetic force theorem for the infinitesimal rotations ofspins [34, 66]. The magnetic force theorem is certainly validand the total energy change underlying calculations of theexchange interactions can be replaced by the change of thesingle-particle energies, which seems to be a general funda-mental property. However, it would be a mistake to make theequality between the magnetic force theorem and equation (3).Such attempts are obviously misleading as the correct useof the magnetic force theorem for the exchange interactionsshould also include contributions of the external magneticfield, which is needed to control the direction of themagnetiza-tion. Equation (3) does not take into account such contribution.This is an approximation leading to the linear dependence ofthe exchange interactions on the response tensor, which canbe justified only in the long-wavelength and strong-couplinglimits. In fact, the exact expression (10) for the energy changeis extremely simple. We only need to know how to find theexternal magnetic field and this can be done by using the linearresponse theory. The total energy change (and the interatomicexchange interactions) in this case will be proportional to theinverse response tensor. The procedure is applicable for cal-culations of the isotropic exchange as well as the DM interac-tions. In the latter case, it requires the additional constrainingfield in order to compensate rotations caused by the DM inter-actions and we have shown how this field can be found usingthe linear response theory.31J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewBesides fundamental aspects of the magnetic force the-orem, there is a number of more practically oriented ques-tions. The first one is which object is more suitable for thedescription of infinitesimal rotations of spins. In this respect,we have considered rigid rotations of the magnetization matrixand local magnetic moments (the schemes m̂ and M, respect-ively). Since the local magnetic moment is nothing but thespherical part of the magnetization matrix and only this spher-ical part is subjected to the constraint conditions in the schemeM (while other components of the magnetization matrix areallowed to relax) such perturbations are expected to cost lessenergy and, therefore, should be more suitable for the analysisof low-energy excitations.Another important question is what to do with the ligandstates. In most of the applications, these states play a role ofeffectivemedium participating in the electron transfer betweenmagnetic transition-metal sites. The exchange interactions aretypically calculated only between the transition-metal sites,while the contributions of the ligand sites, which can carry anappreciable portion of the magnetization due to the hybridiz-ation with transition-metal sites, are ignored. In this respect,we have proposed the downfolding method [61], which allowsus to eliminate the ligand states, by transferring their effect tothe exchange interactions between the transition-metal sites.This can be achieved by employing the adiabaticity conceptand assuming that the fast ligand states instantaneously followthe slow change of the magnetization on the transition-metalsites.An important aspect of the downfolding method is thatit can naturally incorporate the dependence of the exchangeinteractions on the strength of the Stoner coupling IL on theligand sites, which is regarded to be the key ingredient of phe-nomenological GKA rules for the 90◦-exchange and in cer-tain cases primarily responsible for the FM character of thisexchange. Nevertheless, the covalent mixing can easily makethe effective coupling IL negative. Such situation is realized,for instance, in CrCl3, CrI3, and CrO2, where the magneticpolarization of the ligand states is antiparallel to the one ofthe transition-metal states. In such case, the negative IL actsagainst the FM coupling and the ferromagnetism is stabilizedby other mechanisms involving the unoccupied Cr eg states.For most of the applications, we have considered two pic-tures, which provide a supplementary to each other descrip-tion for the exchange interactions. The first one is the cor-related minimal model constructed for the magnetic, mainlytransition-metal, bands near the Fermi level. The second oneis based on the all-electron LSDA, which takes into considera-tion both magnetic transition-metal and ligand states. Once theCoulomb interactions are added to the TB model, they shouldbe treated rigorously: if the static Hartree–Fock approxim-ation is applicable for the analysis of exchange interactionsin insulating systems with lifted orbital degeneracy, the situ-ation in metallic systems can be very different. In certaincases, the use of the Hartree–Fock approximation can leadto a misleading answer, as was demonstrated for the half-metallic CrO2 [120]. However, this is not the only source ofthe error for the exchange interactions. Another problem isrelated to the fact that most of the real materials are far fromthe strong-coupling limit and the commonly used scheme b̂,which leads to equation (3) for the exchange interactions, isno longer justified. A more rigorous approach is to evalu-ate the exchange interactions via the inverse response (thescheme M), which is certainly more difficult from the com-putational point of view. However, such scheme tremend-ously improves the description of interatomic exchange inter-actions in the correlated model, as was demonstrated for insu-lating CrX3 (X= Cl, I) and AMnO3 (A= La, Ho) in theframework of Hartree–Fock approximation. The applicationof correlated models for the exchange interactions in metallicsystems should consider two aspects: (i) the method shouldbe beyond the Hartree–Fock approximation. The commonlyused alternative is DMFT [129, 130]; (ii) the calculations ofthe exchange interactions themselves should be based on theinverse response, as in the scheme M.The exchange interactions in all-electron LSDA (GGA)obey quite different principles. Certainly, this is a simpleapproximation derived in the limit of homogeneous electrongas. Nevertheless, it satisfies certain fundamental sum rulesand on many occasions provides an insightful view on theproperties of systems far beyond this limit [170]. Thus, it canbe regarded as an alternative view on the problem of exchangeinteractions, though not necessarily the perfect one. For themagnetic systems, LSDA and strong-coupling limit are basic-ally incompatible with each other (or this limit can be stronglyunderestimated in LSDA). Therefore, there is absolutely noguarantee that the approximate scheme b̂ can properly capturethe behavior of interatomic exchange interactions and a morerigorous scheme M looks more preferable. Nevertheless, theimprovement expected in the scheme M is somewhat dimin-ished by the effects of other parameters controlling the proper-ties of the exchange interactions in the all-electron case, partic-ularly the choice of the xc field b⃗ and the effective Stoner coup-ling IL on the ligand sites. By taking the xc field from the sumrule, b⃗= Q̂+0 m⃗, one can substantially improve the descrip-tion in the framework of the method b̂. However, the choiceof the parameters IL, which can be strongly depend on thedefinition, poses a more serious problem. We have consideredseveral such definitions. In a number of cases, they providea consistent description for the interatomic exchange interac-tions, but not always. It seems that there is still an ambiguitywith the choice of IL and the isotropic exchange interactionsdepend on such ambiguity. On the other hand, the depend-ence of DM interactions on IL is considerably weaker. As analternative approach, we have proposed the scheme L0, whichassumes that the magnetic moments on the ligand sites areinduced solely by the hybridization with the transition-metalsites, while all Stoner parameters IL can be set to zero. Suchscheme is also formulated in terms of the linear response andmost suitable for the analysis interatomic exchange interac-tions in FM insulators and half-metals.The merging of correlated model for the transition-metalbands with LDA electronic structure for the ligand bands isanother largely unresolved problem. Although such mergingis expected to improve the description of the exchange inter-32J. Phys.: Condens. Matter 36 (2024) 223001 Topical Reviewactions, by combining the effects of Coulomb correlationsin transition-metal bands with the explicit treatment of theligand states, the strategy suffers from many ambiguities, aswas demonstrated for CrCl3 and CrI3. We are always forcedto compromise between genuine physical effect and intrinsicerror of such merging caused by additional approximationsand assumptions, which are inevitable in this case. On the otherhand, on many occasions such merging is not really neededas the physically meaningful picture for the exchange interac-tions can be obtained in correlated 5o model constructed onlyfor the magnetic 3d bands. As a supplementary step, one canalways consider the all-electron LSDA (GGA), which gives anidea about the role played by the ligand states.Among other interesting properties, CrI3 has attracteda considerable attention due to the strong DM interactioninduced by large SO coupling of the heavy I atoms [81].On the other hand, CrI3 is a closed shell material, whereorbital degrees of freedom are quenched by the large crystal-field splitting between occupied t2g and unoccupied eg states.Therefore, from the practical point of view, comparable oreven stronger DM interactions can be expected in the openshell materials evenwithout heavy 5p elements, as was demon-strated, for instance, for CrO2, LaMnO3, and HoMnO3.Data availability statementThe data cannot be made publicly available upon publicationbecause they are not available in a format that is sufficientlyaccessible or reusable by other researchers. The data that sup-port the findings of this study are available upon reasonablerequest from the author.AcknowledgmentsI am grateful to Mikhail Katsnelson, Alexander Liechtensteinand Vladimir Antropov for valuable comments. The TBHamiltonian for Co3Sn2S2 in section 4.2 was constructed bySergey Nikolaev [115]. MANA is supported byWorld PremierInternational Research Center Initiative, MEXT, Japan.Appendix A. Construction and solution of effectiveHubbard-type modelsThe Hubbard-type model,Ĥ=∑ij∑σ∑abHσia,jbĉσ†i a ĉσj b+12∑i∑σσ ′∑abcdUabcdi ĉ†σi a ĉ†σ ′i c ĉσi b ĉσ ′i d , (A.1)is constructed in the Wannier basis for some particular groupbands (the so-called target bands) [43, 89, 171]. The oper-ator ĉ†σi a (ĉσi a ) stands for the creation (annihilation) of an elec-tron with the spin σ in the Wannier orbital a on the site i. Itis assumed that LDA (or LSDA) is the good starting pointfor the noninteracting one-electron part of the model so thatthe parameters Hσia,jb can be associated with the matrix ele-ments of corresponding KS Hamiltonian in the Wannier basis.In LSDA, these parameters depend on spin. Furthermore,Ĥσ = [Hσia,jb] can include spin-diagonal part of the SO coup-ling, which is needed to calculate DM interactions [15, 53].In this case, Ĥσ also depends on spin, even in the non-magnetic LDA, where it holds: Ĥ↓ = Ĥ↑∗. For practical pur-poses we use mainly the LMTO method [90, 91]. Then,the Wannier functions can be generated using the projector-operator technique [43, 89]. The spin-diagonal part of the SOcoupling is typically added to the site-diagonal part of theLMTO Hamiltonian as Ĥ↑,↓ii → Ĥ↑,↓ii ± 12ξiL̂zi . After that, theTBHamiltonian is constructed in theWannier basis for the tar-get bands. Thus, even though the target bands are formally ofthe transition-metal 3d type, the corresponding Wannier func-tions and the TB Hamiltonian include some contributions ofthe SO coupling of the heavy ligand atoms (if any), which areadmixed into the target bands via the hybridization effects.The parameters of screened on-site Coulomb interactions,Û= [Uabcdi ], are evaluated in cRPA starting with the LDA bandstructure [93]. The basic idea of the constraint in this case isto get rid of nonphysical metallic screening in LDA emergingfrom the bands near the Fermi level. For the d electrons, thematrix Û can be fitted in terms of the on-site Coulomb repul-sion U= F 0, the intraatomic exchange interaction J= (F 2 +F 4)/14, and ‘nonsphericity’ B= (9F 2 − 5F 4)/441 (F0, F2,and F4 being radial Slater’s integrals), where U is respons-ible for the charge stability of given electronic configuration,while J and B are responsible for the first and second Hund’srules, respectively [89]. Strictly speaking, the parametrizationin terms ofU, J, andB is valid only for the isolated atoms in thespherical environment [172]. It is used only for the explanatorypurposes, while all numerical calculations are performed withthe matrices of Coulomb interactions Û extracted from cRPAwithout additional fitting. For the t2g electrons alone, it is con-venient to use the Kanamori parametrization in terms of theintraorbital Coulomb repulsionU and the exchange interactionJ [123]. For the d electrons, these parameters can be relatedto the above U and J as U ≈ U+ 8J/7 and J ≈ 0.77J [89].Moreover, for the t2g model,U can be reduced due to additionalchannels of screening by the eg electrons, which are typicallyconsidered in the t2g model, but not in the more general oneconstructed for all d (i.e. t2g + eg) electrons [89].After the construction, the model is solved in the mean-field Hartree–Fock approximation, where the second term inequation (A.1) is replaced by∑i∑σ∑abVσi,abĉσ†i a ĉσi b (A.2)and the potential V̂σi = [Vσi,ab] is found self-consistently. Then,the obtained electronic structure is used to calculate theexchange interactions. The xc field in this case is defined asb̂i = V̂↑i − V̂↓i . When the SO interaction is added to Ĥσ, thepotential V̂σi is recalculated self-consistently to include theeffects or the SO coupling. After that the obtained electronicstructure and V̂σi are used to calculate theDM interactions. The33J. Phys.: Condens. Matter 36 (2024) 223001 Topical ReviewHartree–Fock approximation is believed to be a good start-ing points for the analysis of magnetic properties of insulatingmaterials, where the orbital degeneracy is lifted by the latticedistortions. Other details can be found in [89].Appendix B. Magnetic ground state andrandom-phase approximation for the magnetictransition temperatureHere, we generalize the RPA expression for the critical trans-ition temperature, TS, of compounds with multiple magneticsublattices to the case, where the ground state is a noncollin-ear spin spiral with the propagation vector q. Moreover, thesublattices are allowed to acquire the additional phases γq,µ,describing relative rotations of the magnetization in differentsublattices relative to each other. In all other respects, we fol-low the derivation considered for the collinear case by Ruszet al [173].Our first goal is to find the ground state, corresponding tothe minimum of the energyE(q) =−12∑µνc∗q,µJq,µν cq,ν , (B.1)where cq,µ = eiγq,µ and Jq,µν is the Fourier transform ofJij between the sublattices µ and ν. Equation (B.1) can begeneralized to include the DM interactions dzij by replacingJq,µν with Jq,µν − idzq,µν . For each q, we minimize E(q) withrespect to γq,µ using the gradient descent method. Then, wepick up q corresponding to the global minimum of E(q) andfor each k redefine Jk,µν as Jk,µν → c∗q,µJk,µνcq,ν with thecoefficient cq,ν determined for the ground state. This is noth-ing but the transformation to the new local coordinate framecorresponding to the energy minimum.Then, the spin-wave energies can be associated withthe eigenvalues of the positive-defined matrix Ω̂k = Ŝ−1N̂k,whereNk,µν = δµν∑ν ′Jq,µν ′ − Jk,µν , (B.2)Ŝ= M̂/2, and M̂ is the diagonal matrix of magneticmoments.18 In RPA, the corresponding magnetic transitiontemperature TµS for the sublattice µ is given by [173]:TµS =1+ 1/Sµ3kB(1ΩBZˆdk[N̂−1k]µµ)−1(B.3)(ΩBZ being the volume of the first BZ). For the inequival-ent sublattices, this TµS depends on the ratio of the magneticmoments, Mµ/Mν , which can also depend on the temperat-ure. In order to find true transition temperature, these rationsshould be adjusted to make the same TµS = TνS ≡ TS for all thesublattices [173].18 Here, we use the exchange parameters derived from the static responsefunction and do not consider any renormalization effects [63, 64].The spectral theorem is employed in order to calculate N̂−1k .Then, a small imaginary part, iδ with δ= 0.01meV (corres-ponding to δ/kB ∼ 0.1K) is added to the eigenvalues of N̂kfor k close to q. Finally, the k-space integration is replaced bythe summation on a very dense mesh of k-points (for instance,a typical mesh used for CrCl3 and CrI3 is 258× 258× 258).ORCID iDI V Solovyev https://orcid.org/0000-0002-2010-9877References[1] Heisenberg W 1928 Z. Phys. 49 619[2] Anderson P W 1959 Phys. Rev. 115 2[3] Ruderman M A and Kittel C 1954 Phys. Rev. 96 99[4] Kasuya T 1956 Prog. Theor. Phys. 16 45[5] Yosida K 1957 Phys. Rev. 106 893[6] Dzyaloshinsky I 1958 J. Chem. Phys. Solids 4 241[7] Moriya T 1960 Phys. Rev. 120 91[8] Roth L M, Zeiger H J and Kaplan T A 1966 Phys. Rev.149 519[9] Bruno P and Chappert C 1992 Phys. Rev. B 46 261[10] Liechtenstein A I, Katsnelson M I and Gubanov V A 1984 J.Phys. F: Met. Phys. 14 L125[11] Liechtenstein A I, Katsnelson M I, Antropov V P andGubanov V A 1987 J. Magn. Magn. Mater. 67 65[12] Hohenberg P and Kohn W 1964 Phys. Rev. 136 B864[13] Kohn W and Sham L J 1965 Phys. Rev. 140 A1133[14] von Barth U and Hedin L 1972 J. Phys. C: Solid State Phys.5 1629[15] Sandratskii L M 2017 Phys. Rev. B 96 024450[16] de Gennes P-G 1960 Phys. Rev. 118 141[17] Nagaev E L 1982 Sov. Phys. Usp. 25 31[18] Mackintosh A K and Andersen O K 1975 Electrons at theFermi Surface ed M Springford (Cambridge UniversityPress)[19] Heine V 1980 Solid State Physics vol 35, ed H Ehrenreich,F Seitz and D Turnbull (Academic)[20] Oswald A, Zeller R, Braspenning P J and Dederichs P H1985 J. Phys. F 15 193[21] Oguchi T, Terakura K and Hamada N 1983 J. Phys. F: Met.Phys. 13 145[22] Oguchi T, Terakura K and Williams A R 1983 Phys. Rev. B28 6443[23] Liu S H 1977 Phys. Rev. B 15 4281[24] Prange R E and Korenman V 1979 Phys. Rev. B 19 4691[25] Solovyev I V 2003Magnetic Interactions in Transition-MetalOxides (Recent Research Developments in Magnetism &Magnetic Materials vol 1) (Transworld ResearchNetwork) p 253[26] Kvashnin Y O, Grånäs O, Di Marco I, Katsnelson M I,Lichtenstein A I and Eriksson O 2015 Phys. Rev. B91 125133[27] Korotin Dm M, Mazurenko V V, Anisimov V I andStreltsov S V 2015 Phys. Rev. B 91 224405[28] Yoon H, Kim T J, Sim J H, Jang S W, Ozaki T and Han M J2018 Phys. Rev. B 97 125132[29] Matsumoto M and Akai H 2020 Phys. Rev. B 101 144402[30] Nomoto T, Koretsune T and Arita R 2020 Phys. Rev. B102 014444[31] Grytsiuk S, Hanke J-P, Hoffmann M, Bouaziz J,Gomonay O, Bihlmayer G, Lounis S, Mokrousov Y andBlügel S 2020 Nat. Commun. 11 511[32] He X, Helbig N, Verstraete M J and Bousquet E 2021Comput. Phys. Commun. B 264 10793834https://orcid.org/0000-0002-2010-9877https://orcid.org/0000-0002-2010-9877https://doi.org/10.1007/BF01328601https://doi.org/10.1007/BF01328601https://doi.org/10.1103/PhysRev.115.2https://doi.org/10.1103/PhysRev.115.2https://doi.org/10.1103/PhysRev.96.99https://doi.org/10.1103/PhysRev.96.99https://doi.org/10.1143/PTP.16.45https://doi.org/10.1143/PTP.16.45https://doi.org/10.1103/PhysRev.106.893https://doi.org/10.1103/PhysRev.106.893https://doi.org/10.1016/0022-3697(58)90076-3https://doi.org/10.1016/0022-3697(58)90076-3https://doi.org/10.1103/PhysRev.120.91https://doi.org/10.1103/PhysRev.120.91https://doi.org/10.1103/PhysRev.149.519https://doi.org/10.1103/PhysRev.149.519https://doi.org/10.1103/PhysRevB.46.261https://doi.org/10.1103/PhysRevB.46.261https://doi.org/10.1088/0305-4608/14/7/007https://doi.org/10.1088/0305-4608/14/7/007https://doi.org/10.1016/0304-8853(87)90721-9https://doi.org/10.1016/0304-8853(87)90721-9https://doi.org/10.1103/PhysRev.136.B864https://doi.org/10.1103/PhysRev.136.B864https://doi.org/10.1103/PhysRev.140.A1133https://doi.org/10.1103/PhysRev.140.A1133https://doi.org/10.1088/0022-3719/5/13/012https://doi.org/10.1088/0022-3719/5/13/012https://doi.org/10.1103/PhysRevB.96.024450https://doi.org/10.1103/PhysRevB.96.024450https://doi.org/10.1103/PhysRev.118.141https://doi.org/10.1103/PhysRev.118.141https://doi.org/10.1070/PU1982v025n01ABEH004495https://doi.org/10.1070/PU1982v025n01ABEH004495https://doi.org/10.1088/0305-4608/15/1/021https://doi.org/10.1088/0305-4608/15/1/021https://doi.org/10.1088/0305-4608/13/1/018https://doi.org/10.1088/0305-4608/13/1/018https://doi.org/10.1103/PhysRevB.28.6443https://doi.org/10.1103/PhysRevB.28.6443https://doi.org/10.1103/PhysRevB.15.4281https://doi.org/10.1103/PhysRevB.15.4281https://doi.org/10.1103/PhysRevB.19.4691https://doi.org/10.1103/PhysRevB.19.4691https://doi.org/10.1103/PhysRevB.91.125133https://doi.org/10.1103/PhysRevB.91.125133https://doi.org/10.1103/PhysRevB.91.224405https://doi.org/10.1103/PhysRevB.91.224405https://doi.org/10.1103/PhysRevB.97.125132https://doi.org/10.1103/PhysRevB.97.125132https://doi.org/10.1103/PhysRevB.101.144402https://doi.org/10.1103/PhysRevB.101.144402https://doi.org/10.1103/PhysRevB.102.014444https://doi.org/10.1103/PhysRevB.102.014444https://doi.org/10.1038/s41467-019-14030-3https://doi.org/10.1038/s41467-019-14030-3https://doi.org/10.1016/j.cpc.2021.107938https://doi.org/10.1016/j.cpc.2021.107938J. Phys.: Condens. Matter 36 (2024) 223001 Topical Review[33] Stocks G M, Ujfalussy B, Wang X, Nicholson D M C,Shelton W A, Wang Y, Canning A and Gyorffy B L 1998Phil. Mag. 78 665[34] Bruno P 2003 Phys. Rev. Lett. 90 087205[35] Streib S, Borisov V, Pereiro M, Bergman A, Sjöqvist E,Delin A, Eriksson O and Thonig D 2020 Phys. Rev. B102 214407[36] Anderson P W 1950 Phys. Rev. 79 350[37] Goodenough J B 1955 Phys. Rev. 100 564[38] Goodenough J B 1958 J. Phys. Chem. Solids 6 287[39] Kanamori J 1959 J. Phys. Chem. Solids 10 87[40] Mazurenko V V, Skornyakov S L, Kozhevnikov A V , Mila Fand Anisimov V I 2007 Phys. Rev. B 75 224408[41] Streltsov S V and Khomskii D I 2008 Phys. Rev. B 77 064405[42] Solovyev I 2009 J. Phys. Soc. Japan 78 054710[43] Marzari N, Mostofi A A, Yates J R, Souza I and Vanderbilt D2012 Rev. Mod. Phys. 84 1419[44] King-Smith R D and Vanderbilt D 1993 Phys. Rev. B47 1651[45] Dzyaloshinskii I 1964 Sov. Phys. JETP 19 960 (available at:http://jetp.ras.ru/cgi-bin/dn/e_019_04_0960.pdf)[46] Sandratskii L M 1998 Adv. Phys. 47 91[47] Koehler W C, Cable J W, Wilkinson M K and Wollan E O1966 Phys. Rev. 151 414[48] Solovyev I V 2011 Phys. Rev. B 83 054404[49] Kvashnin Y O, Bergman A, Lichtenstein A I andKatsnelson M I 2020 Phys. Rev. B 102 115162[50] Ebert H and Mankovsky S 2009 Phys. Rev. B 79 045209[51] Mahfouzi F and Kioussis N 2021 Phys. Rev. B 103 094410[52] Lange H, Mankovsky S, Polesya S, Weißenhofer M,Nowak U and Ebert H 2023 Phys. Rev. B 107 115176[53] Solovyev I V 2023 Phys. Rev. B 107 054442[54] Dederichs P H, Blügel S, Zeller R and Akai H 1984 Phys.Rev. Lett. 53 2512[55] Vignale G and Rasolt M 1987 Phys. Rev. Lett. 59 2360[56] Vignale G and Rasolt M 1988 Phys. Rev. B 37 10685[57] Solovyev I V and Terakura K 1998 Phys. Rev. B 58 15496[58] Katsnelson M I, Kvashnin Y O, Mazurenko V V andLichtenstein A I 2010 Phys. Rev. B 82 100403(R)[59] Kübler J, Höck K-H, Sticht J and Williams A R 1988 J. Phys.F 18 469[60] Eich F G and Gross E K U 2013 Phys. Rev. Lett. 111 156401[61] Solovyev I V 2021 Phys. Rev. B 103 104428[62] Solovyev I V 2014 Phys. Rev. B 90 024417[63] Katsnelson M I and Lichtenstein A I 2004 J. Phys.: Condens.Matter 16 7439[64] Szczech Y H, Tusch M A and Logan D E 1995 Phys. Rev.Lett. 74 2804[65] Antropov V P, van Schilfgaarde M, Brink S and Xu J L 2006J. Appl. Phys. 99 08F507[66] Antropov V P 2004 No new renormalized magnetic forcetheorem (arXiv:cond-mat/0407739 [cond-mat.mtrl-sci])[67] Solovyev I V and Terakura K 1999 Phys. Rev. Lett. 82 2959[68] Solovyev I V 1999 Phys. Rev. B 60 8550[69] van den Brink J, Meinders M B J, Lorenzana J, Eder R andSawatzky G A 1995 Phys. Rev. Lett. 75 4658[70] Kikuchi T, Koretsune T, Arita R and Tatara G 2016 Phys.Rev. Lett. 116 247201[71] Cooke J F 1973 Phys. Rev. B 7 1108[72] Callaway J, Wang C S and Laurent D G 1981 Phys. Rev. B24 6491[73] Savrasov S Y 1998 Phys. Rev. Lett. 81 2570[74] Antropov V P, Katsnelson M I, Harmon B N,van Schilfgaarde M and Kusnezov D 1996 Phys. Rev. B54 1019[75] Halilov S V, Eschrig H, Perlov A Y and Oppeneer P M 1998Phys. Rev. B 58 293[76] Solovyev I V, Ushakov A V and Streltsov S V 2022 Phys.Rev. B 106 L180401[77] Mryasov O N, Nowak U, Guslienko K Y and Chantrell R W2005 Europhys. Lett. 69 805[78] Logemann R, Rudenko A N, Katsnelson M I and Kirilyuk A2017 J. Phys.: Condens. Matter 29 335801[79] Morosin B and Narath A 1964 J. Chem. Phys. 40 1958[80] McGuire M A, Dixit H, Cooper V R and Sales B C 2015Chem. Mater. 27 612[81] Chen L, Chung J H, Gao B, Chen T, Stone M B,Kolesnikov A I, Huang Q and Dai P 2018 Phys. Rev. X8 041028[82] McGuire M A 2017 Crystals 7 121[83] Chen L, Stone M B, Kolesnikov A I, Winn B, Shon W, Dai Pand Chung J H 2022 2D Mater. 9 015006[84] Schneeloch J A, Tao Y, Cheng Y, Daemen L, Xu G, Zhang Qand Louca D 2022 npj Quantum Mater. 7 66[85] Huang B et al 2017 Nature 546 270[86] Lado J L and Fernández-Rossier J 2017 2D Mater. 4 035002[87] Mermin N D and Wagner H 1966 Phys. Rev. Lett. 17 1133[88] Chaloupka J, Jackeli G and Khaliullin G 2013 Phys. Rev.Lett. 110 097204[89] Solovyev I V 2008 J. Phys.: Condens. Matter 20 293201[90] Andersen O K 1975 Phys. Rev. B 12 3060[91] Gunnarsson O, Jepsen O and Andersen O K 1983 Phys. Rev.B 27 7144[92] Kanamori J 1957 Prog. Theor. Phys. 17 177[93] Aryasetiawan F, Imada M, Georges A, Kotliar G, Biermann Sand Lichtenstein A I 2004 Phys. Rev. B 70 195104[94] Anisimov V I, Zaanen J and Andersen O K 1991 Phys. Rev.B 44 943[95] Solovyev I V, Dederichs P H and Anisimov V I 1994 Phys.Rev. B 50 16861[96] Besbes O, Nikolaev S, Meskini N and Solovyev I 2019 Phys.Rev. B 99 104432[97] Solovyev I V and Streltsov S V 2019 Phys. Rev. Mater.3 114402[98] Solovyev I V and Terakura K 1999 Phys. Rev. Lett. 83 2825[99] Solovyev I V and Terakura K 2019 Phys. Rev. B63 174425[100] Gunnarsson O 1976 J. Phys. F: Met. Phys. 6 587[101] Janak J F 1977 Phys. Rev. B 16 255[102] Brooks M S S and Johansson B 1983 J. Phys. F: Met. Phys.13 L197[103] Ke L and Katsnelson M I 2021 npj Comput. Mater. 7 4[104] Ku W, Rosner H, Pickett W E and Scalettar R T 2002 Phys.Rev. Lett. 89 167204[105] Badrtdinov D I, Nikolaev S A, Katsnelson M I andMazurenko V V 2016 Phys. Rev. B 94 224418[106] Olsen T 2021 Phys. Rev. Lett. 127 166402[107] de Groot R A, Mueller F M, van Engen P G andBuschow K H J 1983 Phys. Rev. Lett. 50 2024[108] van Leuken H and de Groot R A 1995 Phys. Rev. Lett.74 1171[109] Schwarz K 1986 J. Phys. F: Met. Phys. 16 L211[110] Liu E et al 2018 Nat. Phys. 14 1125[111] Liu D F et al 2019 Science 365 1282[112] Minami S, Ishii F, Hirayama M, Nomoto T, Koretsune T andArita R 2020 Phys. Rev. B 102 205128[113] Jiao L, Xu Q, Cheon Y, Sun Y, Felser C, Liu E and Wirth S2019 Phys. Rev. B 99 245158[114] Yanagi Y, Ikeda J, Fujiwara K, Nomura K, Tsukazaki A andSuzuki M-T 2021 Phys. Rev. B 103 205112[115] Solovyev I V, Nikolaev S A, Ushakov A V, Irkhin V Y,Tanaka A and Streltsov S V 2022 Phys. Rev. B 105 014415[116] Korotin M A, Anisimov V I, Khomskii D I andSawatzky G A 1998 Phys. Rev. Lett. 80 4305[117] Maurya V, Sharma G and Joshi K B 2021 Phys. Scr.96 055807[118] Mazin I I, Singh D and Ambrosch-Draxl C 1999 Phys. Rev. B59 41135https://doi.org/10.1080/13642819808206775https://doi.org/10.1080/13642819808206775https://doi.org/10.1103/PhysRevLett.90.087205https://doi.org/10.1103/PhysRevLett.90.087205https://doi.org/10.1103/PhysRevB.102.214407https://doi.org/10.1103/PhysRevB.102.214407https://doi.org/10.1103/PhysRev.79.350https://doi.org/10.1103/PhysRev.79.350https://doi.org/10.1103/PhysRev.100.564https://doi.org/10.1103/PhysRev.100.564https://doi.org/10.1016/0022-3697(58)90107-0https://doi.org/10.1016/0022-3697(58)90107-0https://doi.org/10.1016/0022-3697(59)90061-7https://doi.org/10.1016/0022-3697(59)90061-7https://doi.org/10.1103/PhysRevB.75.224408https://doi.org/10.1103/PhysRevB.75.224408https://doi.org/10.1103/PhysRevB.77.064405https://doi.org/10.1103/PhysRevB.77.064405https://doi.org/10.1143/JPSJ.78.054710https://doi.org/10.1143/JPSJ.78.054710https://doi.org/10.1103/RevModPhys.84.1419https://doi.org/10.1103/RevModPhys.84.1419https://doi.org/10.1103/PhysRevB.47.1651https://doi.org/10.1103/PhysRevB.47.1651http://jetp.ras.ru/cgi-bin/dn/e_019_04_0960.pdfhttps://doi.org/10.1080/000187398243573https://doi.org/10.1080/000187398243573https://doi.org/10.1103/PhysRev.151.414https://doi.org/10.1103/PhysRev.151.414https://doi.org/10.1103/PhysRevB.83.054404https://doi.org/10.1103/PhysRevB.83.054404https://doi.org/10.1103/PhysRevB.102.115162https://doi.org/10.1103/PhysRevB.102.115162https://doi.org/10.1103/PhysRevB.79.045209https://doi.org/10.1103/PhysRevB.79.045209https://doi.org/10.1103/PhysRevB.103.094410https://doi.org/10.1103/PhysRevB.103.094410https://doi.org/10.1103/PhysRevB.107.115176https://doi.org/10.1103/PhysRevB.107.115176https://doi.org/10.1103/PhysRevB.107.054442https://doi.org/10.1103/PhysRevB.107.054442https://doi.org/10.1103/PhysRevLett.53.2512https://doi.org/10.1103/PhysRevLett.53.2512https://doi.org/10.1103/PhysRevLett.59.2360https://doi.org/10.1103/PhysRevLett.59.2360https://doi.org/10.1103/PhysRevB.37.10685https://doi.org/10.1103/PhysRevB.37.10685https://doi.org/10.1103/PhysRevB.58.15496https://doi.org/10.1103/PhysRevB.58.15496https://doi.org/10.1103/PhysRevB.82.100403https://doi.org/10.1103/PhysRevB.82.100403https://doi.org/10.1088/0305-4608/18/3/018https://doi.org/10.1088/0305-4608/18/3/018https://doi.org/10.1103/PhysRevLett.111.156401https://doi.org/10.1103/PhysRevLett.111.156401https://doi.org/10.1103/PhysRevB.103.104428https://doi.org/10.1103/PhysRevB.103.104428https://doi.org/10.1103/PhysRevB.90.024417https://doi.org/10.1103/PhysRevB.90.024417https://doi.org/10.1088/0953-8984/16/41/023https://doi.org/10.1088/0953-8984/16/41/023https://doi.org/10.1103/PhysRevLett.74.2804https://doi.org/10.1103/PhysRevLett.74.2804https://doi.org/10.1063/1.2176392https://doi.org/10.1063/1.2176392https://arxiv.org/abs/cond-mat/040773https://doi.org/10.1103/PhysRevLett.82.2959https://doi.org/10.1103/PhysRevLett.82.2959https://doi.org/10.1103/PhysRevB.60.8550https://doi.org/10.1103/PhysRevB.60.8550https://doi.org/10.1103/PhysRevLett.75.4658https://doi.org/10.1103/PhysRevLett.75.4658https://doi.org/10.1103/PhysRevLett.116.247201https://doi.org/10.1103/PhysRevLett.116.247201https://doi.org/10.1103/PhysRevB.7.1108https://doi.org/10.1103/PhysRevB.7.1108https://doi.org/10.1103/PhysRevB.24.6491https://doi.org/10.1103/PhysRevB.24.6491https://doi.org/10.1103/PhysRevLett.81.2570https://doi.org/10.1103/PhysRevLett.81.2570https://doi.org/10.1103/PhysRevB.54.1019https://doi.org/10.1103/PhysRevB.54.1019https://doi.org/10.1103/PhysRevB.58.293https://doi.org/10.1103/PhysRevB.58.293https://doi.org/10.1103/PhysRevB.106.L180401https://doi.org/10.1103/PhysRevB.106.L180401https://doi.org/10.1209/epl/i2004-10404-2https://doi.org/10.1209/epl/i2004-10404-2https://doi.org/10.1088/1361-648X/aa7b00https://doi.org/10.1088/1361-648X/aa7b00https://doi.org/10.1063/1.1725428https://doi.org/10.1063/1.1725428https://doi.org/10.1021/cm504242thttps://doi.org/10.1021/cm504242thttps://doi.org/10.1103/PhysRevX.8.041028https://doi.org/10.1103/PhysRevX.8.041028https://doi.org/10.3390/cryst7050121https://doi.org/10.3390/cryst7050121https://doi.org/10.1088/2053-1583/ac2e7ahttps://doi.org/10.1088/2053-1583/ac2e7ahttps://doi.org/10.1038/s41535-022-00473-3https://doi.org/10.1038/s41535-022-00473-3https://doi.org/10.1038/nature22391https://doi.org/10.1038/nature22391https://doi.org/10.1088/2053-1583/aa75edhttps://doi.org/10.1088/2053-1583/aa75edhttps://doi.org/10.1103/PhysRevLett.17.1133https://doi.org/10.1103/PhysRevLett.17.1133https://doi.org/10.1103/PhysRevLett.110.097204https://doi.org/10.1103/PhysRevLett.110.097204https://doi.org/10.1088/0953-8984/20/29/293201https://doi.org/10.1088/0953-8984/20/29/293201https://doi.org/10.1103/PhysRevB.12.3060https://doi.org/10.1103/PhysRevB.12.3060https://doi.org/10.1103/PhysRevB.27.7144https://doi.org/10.1103/PhysRevB.27.7144https://doi.org/10.1143/PTP.17.177https://doi.org/10.1143/PTP.17.177https://doi.org/10.1103/PhysRevB.70.195104https://doi.org/10.1103/PhysRevB.70.195104https://doi.org/10.1103/PhysRevB.44.943https://doi.org/10.1103/PhysRevB.44.943https://doi.org/10.1103/PhysRevB.50.16861https://doi.org/10.1103/PhysRevB.50.16861https://doi.org/10.1103/PhysRevB.99.104432https://doi.org/10.1103/PhysRevB.99.104432https://doi.org/10.1103/PhysRevMaterials.3.114402https://doi.org/10.1103/PhysRevMaterials.3.114402https://doi.org/10.1103/PhysRevLett.83.2825https://doi.org/10.1103/PhysRevLett.83.2825https://doi.org/10.1103/PhysRevB.63.174425https://doi.org/10.1103/PhysRevB.63.174425https://doi.org/10.1088/0305-4608/6/4/018https://doi.org/10.1088/0305-4608/6/4/018https://doi.org/10.1103/PhysRevB.16.255https://doi.org/10.1103/PhysRevB.16.255https://doi.org/10.1088/0305-4608/13/10/003https://doi.org/10.1088/0305-4608/13/10/003https://doi.org/10.1038/s41524-020-00469-2https://doi.org/10.1038/s41524-020-00469-2https://doi.org/10.1103/PhysRevLett.89.167204https://doi.org/10.1103/PhysRevLett.89.167204https://doi.org/10.1103/PhysRevB.94.224418https://doi.org/10.1103/PhysRevB.94.224418https://doi.org/10.1103/PhysRevLett.127.166402https://doi.org/10.1103/PhysRevLett.127.166402https://doi.org/10.1103/PhysRevLett.50.2024https://doi.org/10.1103/PhysRevLett.50.2024https://doi.org/10.1103/PhysRevLett.74.1171https://doi.org/10.1103/PhysRevLett.74.1171https://doi.org/10.1088/0305-4608/16/9/002https://doi.org/10.1088/0305-4608/16/9/002https://doi.org/10.1038/s41567-018-0234-5https://doi.org/10.1038/s41567-018-0234-5https://doi.org/10.1126/science.aav2873https://doi.org/10.1126/science.aav2873https://doi.org/10.1103/PhysRevB.102.205128https://doi.org/10.1103/PhysRevB.102.205128https://doi.org/10.1103/PhysRevB.99.245158https://doi.org/10.1103/PhysRevB.99.245158https://doi.org/10.1103/PhysRevB.103.205112https://doi.org/10.1103/PhysRevB.103.205112https://doi.org/10.1103/PhysRevB.105.014415https://doi.org/10.1103/PhysRevB.105.014415https://doi.org/10.1103/PhysRevLett.80.4305https://doi.org/10.1103/PhysRevLett.80.4305https://doi.org/10.1088/1402-4896/abec00https://doi.org/10.1088/1402-4896/abec00https://doi.org/10.1103/PhysRevB.59.411https://doi.org/10.1103/PhysRevB.59.411J. Phys.: Condens. Matter 36 (2024) 223001 Topical Review[119] Katsnelson M I, Irkhin V Y, Chioncel L, Lichtenstein A I andde Groot R A 2008 Rev. Mod. Phys. 80 315[120] Solovyev I V, Kashin I V and Mazurenko V V 2015 Phys.Rev. B 92 144407[121] Skomski R 2008 Simple Models of Magnetism (OxfordUniversity Press)[122] Porta P, Marezio M, Remeika J P and Dernier P D 1972Mater. Res. Bull. 7 157[123] Kanamori J 1963 Prog. Theor. Phys. 30 275[124] Bradley C J and Cracknell A P 1972 The MathematicalTheory of Symmetry in Solids (Clarendon)[125] Sims H, Oset S J, Butler W H, MacLaren J M andMarsman M 2010 Phys. Rev. B 81 224436[126] Beliayev E Y, Horielyi V A and Kolesnichenko Y A 2021Low Temp. Phys. 47 355[127] Katanin A A, Belozerov A S, Lichtenstein A I andKatsnelson M I 2023 Phys. Rev. B 107 235118[128] Staunton J, Gyorffy B L, Pindor A J, Stocks G M andWinter H 1984 J. Magn. Magn. Mater. 45 14[129] Georges A, Kotliar G, Krauth W and Rozenberg M J 1996Rev. Mod. Phys. 68 13[130] Kotliar G, Savrasov S Y, Haule K, Oudovenko V S,Parcollet O and Marianetti C A 2006 Rev. Mod. Phys.78 865[131] Katsnelson M I and Lichtenstein A I 2000 Phys. Rev. B61 8906[132] Vaqueiro P and Sobany G G 2009 Solid State Sci. 11 513[133] Perdew J P, Burke K and Ernzerhof M 1996 Phys. Rev. Lett.77 3865[134] Kresse G and Hafner J 1993 Phys. Rev. B 47 558[135] Mostofi A A, Yates J R, Pizzi G, Lee Y S, Souza I,Vanderbilt D and Marzari N 2014 Comput. Phys.Commun. 185 2309[136] Rossi A et al 2021 Phys. Rev. B 104 155115[137] Ozawa A and Nomura K 2019 J. Phys. Soc. Japan 88 123703[138] Turzhevskii S A, Lichtenstein A I and Katsnelson M I 1990Sov. Phys. Solid State 32 1138[139] Singer R, Fähnle M and Bihlmayer G 2005 Phys. Rev. B71 214435[140] Streltsov S, Mazin I I and Foyevtsova K 2015 Phys. Rev. B92 134408[141] Schnelle W et al 2021 Phys. Rev. B 103 214413[142] Liu C et al 2021 Sci. China Phys. Mech. Astron. 64 217062[143] Zhang Q et al 2021 Phys. Rev. Lett. 127 117201[144] Tokura Y (ed) 2000 Colossal Magnetoresistive Oxides(Gordon and Breach Science publishers)[145] Solovyev I, Hamada N and Terakura K 1996 Phys. Rev. Lett.76 4825[146] Sawada H, Morikawa Y, Terakura K and Hamada N 1997Phys. Rev. B 56 12154[147] He J and Franchini C 2012 Phys. Rev. B 86 235117[148] Kimura T, Goto T, Shintani H, Ishizaka K, Arima T andTokura Y 2003 Nature 426 55[149] Muñoz A, Casáis M T, Alonso J A, Martínez-Lope M J,Martínez J L and Fernández-Díaz M T 2001 Inorg. Chem.40 1020[150] Ishiwata S, Kaneko Y, Tokunaga Y, Taguchi Y, Arima T Hand Tokura Y 2010 Phys. Rev. B 81 100411(R)[151] Cheong S-W and Mostovoy M 2007 Nat. Mater. 6 13[152] Khomskii D 2009 Physics 2 20[153] Tokura Y and Seki S 2010 Adv. Mater 22 1554[154] Elemans J B A A, van Laar B, van der Veen K R andLoopstra B O 1971 J. Solid State Chem. 3 238[155] Kanamori J 1960 J. Appl. Phys. 31 S14[156] Gor’kov L P and Kresin V Z 1998 JETP Lett. 67 986[157] Solovyev I, Hamada N and Terakura K 1996 Phys. Rev. B53 7158[158] Hirota K, Kaneko N, Nishizawa A and Endoh Y 1996 J.Phys. Soc. Japan 65 3736[159] Moussa F, Hennion M, Rodriguez-Carvajal J, Moudden H,Pinsard L and Revcolevschi A 1996 Phys. Rev. B54 15149[160] Kugel K I and Khomskii D I 1982 Sov. Phys.-Usp. 25 231[161] Picozzi S, Yamauchi K, Sanyal B, Sergienko I A andDagotto E 2007 Phys. Rev. Lett. 99 227201[162] Okuyama D et al 2011 Phys. Rev. B 84 054440[163] Solovyev I V, Valentyuk M V and Mazurenko V V 2012Phys. Rev. B 86 144406[164] Zaanen J and Sawatzky G A 1987 Can. J. Phys. 65 1262[165] Bozorth R M 1958 Phys. Rev. Lett. 1 362[166] Treves D 1962 Phys. Rev. 125 1843[167] Matsumoto G 1970 J. Phys. Soc. Japan 29 606[168] Skumryev V, Ott F, Coey J M D, Anane A, Renard J-P,Pinsard-Gaudart L and Revcolevschi A 1999 Eur. Phys. J.B 11 401[169] Solovyev I V, Kashin I V and Mazurenko V V 2016 J. Phys.:Condens. Matter 28 216001[170] Gunnarsson O and Lundqvist B I 1976 Phys. Rev. B13 4274[171] Imada M and Miyake T 2010 J. Phys. Soc. Japan 79 112001[172] Slater J C 1929 Phys. Rev. 34 1293[173] Rusz J, Turek I and Divǐs M 2005 Phys. Rev. B 71 17440836https://doi.org/10.1103/RevModPhys.80.315https://doi.org/10.1103/RevModPhys.80.315https://doi.org/10.1103/PhysRevB.92.144407https://doi.org/10.1103/PhysRevB.92.144407https://doi.org/10.1016/0025-5408(72)90272-3https://doi.org/10.1016/0025-5408(72)90272-3https://doi.org/10.1143/PTP.30.275https://doi.org/10.1143/PTP.30.275https://doi.org/10.1103/PhysRevB.81.224436https://doi.org/10.1103/PhysRevB.81.224436https://doi.org/10.1063/10.0004228https://doi.org/10.1063/10.0004228https://doi.org/10.1103/PhysRevB.107.235118https://doi.org/10.1103/PhysRevB.107.235118https://doi.org/10.1016/0304-8853(84)90367-6https://doi.org/10.1016/0304-8853(84)90367-6https://doi.org/10.1103/RevModPhys.68.13https://doi.org/10.1103/RevModPhys.68.13https://doi.org/10.1103/RevModPhys.78.865https://doi.org/10.1103/RevModPhys.78.865https://doi.org/10.1103/PhysRevB.61.8906https://doi.org/10.1103/PhysRevB.61.8906https://doi.org/10.1016/j.solidstatesciences.2008.06.017https://doi.org/10.1016/j.solidstatesciences.2008.06.017https://doi.org/10.1103/PhysRevLett.77.3865https://doi.org/10.1103/PhysRevLett.77.3865https://doi.org/10.1103/PhysRevB.47.558https://doi.org/10.1103/PhysRevB.47.558https://doi.org/10.1016/j.cpc.2014.05.003https://doi.org/10.1016/j.cpc.2014.05.003https://doi.org/10.1103/PhysRevB.104.155115https://doi.org/10.1103/PhysRevB.104.155115https://doi.org/10.7566/JPSJ.88.123703https://doi.org/10.7566/JPSJ.88.123703https://doi.org/10.1103/PhysRevB.71.214435https://doi.org/10.1103/PhysRevB.71.214435https://doi.org/10.1103/PhysRevB.92.134408https://doi.org/10.1103/PhysRevB.92.134408https://doi.org/10.1103/PhysRevB.103.214413https://doi.org/10.1103/PhysRevB.103.214413https://doi.org/10.1007/s11433-020-1597-6https://doi.org/10.1007/s11433-020-1597-6https://doi.org/10.1103/PhysRevLett.127.117201https://doi.org/10.1103/PhysRevLett.127.117201https://doi.org/10.1103/PhysRevLett.76.4825https://doi.org/10.1103/PhysRevLett.76.4825https://doi.org/10.1103/PhysRevB.56.12154https://doi.org/10.1103/PhysRevB.56.12154https://doi.org/10.1103/PhysRevB.86.235117https://doi.org/10.1103/PhysRevB.86.235117https://doi.org/10.1038/nature02018https://doi.org/10.1038/nature02018https://doi.org/10.1021/ic0011009https://doi.org/10.1021/ic0011009https://doi.org/10.1103/PhysRevB.81.100411https://doi.org/10.1103/PhysRevB.81.100411https://doi.org/10.1038/nmat1804https://doi.org/10.1038/nmat1804https://doi.org/10.1103/Physics.2.20https://doi.org/10.1103/Physics.2.20https://doi.org/10.1002/adma.200901961https://doi.org/10.1002/adma.200901961https://doi.org/10.1016/0022-4596(71)90034-Xhttps://doi.org/10.1016/0022-4596(71)90034-Xhttps://doi.org/10.1063/1.1984590https://doi.org/10.1063/1.1984590https://doi.org/10.1007/BF03178245https://doi.org/10.1007/BF03178245https://doi.org/10.1103/PhysRevB.53.7158https://doi.org/10.1103/PhysRevB.53.7158https://doi.org/10.1143/JPSJ.65.3736https://doi.org/10.1143/JPSJ.65.3736https://doi.org/10.1103/PhysRevB.54.15149https://doi.org/10.1103/PhysRevB.54.15149https://doi.org/10.1070/PU1982v025n04ABEH004537https://doi.org/10.1070/PU1982v025n04ABEH004537https://doi.org/10.1103/PhysRevLett.99.227201https://doi.org/10.1103/PhysRevLett.99.227201https://doi.org/10.1103/PhysRevB.84.054440https://doi.org/10.1103/PhysRevB.84.054440https://doi.org/10.1103/PhysRevB.86.144406https://doi.org/10.1103/PhysRevB.86.144406https://doi.org/10.1139/p87-201https://doi.org/10.1139/p87-201https://doi.org/10.1103/PhysRevLett.1.362https://doi.org/10.1103/PhysRevLett.1.362https://doi.org/10.1103/PhysRev.125.1843https://doi.org/10.1103/PhysRev.125.1843https://doi.org/10.1143/JPSJ.29.606https://doi.org/10.1143/JPSJ.29.606https://doi.org/10.1007/BF03219176https://doi.org/10.1007/BF03219176https://doi.org/10.1088/0953-8984/28/21/216001https://doi.org/10.1088/0953-8984/28/21/216001https://doi.org/10.1103/PhysRevB.13.4274https://doi.org/10.1103/PhysRevB.13.4274https://doi.org/10.1143/JPSJ.79.112001https://doi.org/10.1143/JPSJ.79.112001https://doi.org/10.1103/PhysRev.34.1293https://doi.org/10.1103/PhysRev.34.1293https://doi.org/10.1103/PhysRevB.71.174408https://doi.org/10.1103/PhysRevB.71.174408 Linear response theories for interatomic exchange interactions 1. Introduction 2. Infinitesimal spin rotations and exchange interactions 2.1. Basic idea, notations, and conventions 2.2. Spin spirals, SO coupling, and DM interactions 2.3. General expression for the energy change 2.4. Response tensor 2.5. Exchange interactions 2.6. Sum rule and local xc field 2.7. Right object to rotate: magnetization matrices versus local magnetic moments 2.8. Rotations of the xc field as an alternative perturbation 2.9. Relationship to the spin-wave spectra 2.10. Elimination of the ligand spins 3. Chromium trihalides 3.1. Three-orbital t2g model in LSDA 3.2. Five-orbital Cr 3d model 3.3. Role of the ligand states 3.3.1. Cr 3d and Xp bands in LSDA. 3.3.2. Ligand states in correlated Cr 3d band. 3.4. Merging correlated and ligand bands 3.5. DM interactions 3.6. Comparison with experimental data 3.7. Brief summary 4. Half-metallic ferromagnets 4.1. CrO2 4.2. Co3Sn2S2 5. Orthorhombic perovskite manganites 5.1. Correlated 5o model 5.2. All-electron model in LSDA 5.3. DM interactions 5.4. Internal instability of the E phase 5.5. Brief summary 6. Other developments 6.1. Symmetric anisotropic exchange interactions 6.2. Dynamic electron correlations 7. Summary and conclusions Appendix A. Construction and solution of effective Hubbard-type models Appendix B. Magnetic ground state and random-phase approximation for the magnetic transition temperature References