Linear response theories for interatomic exchange interactions

In the magnetic equilibrium, the 1st derivative of the total energy with respect to a small change of the magnetization $\boldsymbol{m}(\boldsymbol{r})$ vanishes and the energy changes is described by the 2nd derivative. In a general sense, the magnetic force theorem states that not only the total energy but also its 2nd derivative with respect to the infinitesimal rotations of the magnetization is the ground-state property as it can be expressed via the eigenvalues and eigenfunctions of the ground state. This statement can be traced back to fundamentals of the quantum mechanics, where the system can be measured only via perturbations. Therefore, it is logical that the 2nd derivative can be connected to the properties of the ground state. The response function (or the susceptibility) is the useful tool, which establishes such connection by means of the perturbation theory.

In practical terms, we will deal mainly with SDFT [12–14], where the ground-state magnetization density and the total energy are described with the help of the single-particle spin-dependent KS Hamiltonian $H^{\sigma}(\boldsymbol{r})$ with some local self-consistent potential incorporating all effects of exchange and correlations [13, 14]. This locality implies that the change of the potential in certain point r depends only on the change of magnetization in the same point r . As a consequence, the energy change caused by the infinitesimal rotations of the magnetization can be presented in the form of pairwise interactions. Without SO coupling it corresponds to the isotropic bilinear Heisenberg model. This is a general property of the 2nd order perturbation theory with the local potentials.

From the viewpoint of analysis and interpretation, it is more convenient to adopt the TB representation, which deals with the atomically resolved properties emerging from the solution of some lattice model. Another advantage of the TB representation is the on-site Coulomb interactions, which can be easily incorporated into the model. The purpose of these interactions is to correct limitations of the local density approximation (LDA) or the generalized gradient approximation (GGA), which are derived in the limit of homogeneous electron gas and typically used to describe the effects of exchange and correlations in SDFT. Mathematically, this can be done by constructing the orthonormal basis of localized Wannier functions centered on the atomic sites [43]. The KS Hamiltonian in this basis is the matrix specified by the lattice (i, j) and orbital (a, b) indices, $\hat{H}^{\sigma} \equiv [H_{ia,jb}^{\sigma}]$, so that all other matrices can be obtained from $\hat{H}^{\sigma}$. For instance, the Green function is $\hat{G}^{\sigma}(\varepsilon) = [ \varepsilon - \hat{H}^{\sigma} ]^{-1}$ and the magnitude of the magnetization is given by $\hat{m} = \Theta(\varepsilon_{\mathrm{F}}-\hat{H}^{\uparrow}) - \Theta(\varepsilon_{\mathrm{F}}-\hat{H}^{\downarrow})$, in terms of the Heaviside function Θ. We assume that the xc potential in this TB representation remains local in the sense that on each atomic site i it depends only on the magnetization $\hat{m}_{i}$ on the same site. It is a common practice to relate the magnetic part of such potential, the so-called xc field $\hat{b}$, with the site-diagonal elements of the TB Hamiltonian: $\hat{b}_{i} = \hat{H}^{\uparrow}_{ii}-\hat{H}^{\downarrow}_{ii}$. However, such $\hat{b}$ is ill-defined because it ignores the non-local contributions steaming from the off-diagonal part of $\hat{H}^{\sigma}_{ij}$ [30]. A more consistent definition of the local $\hat{b}$ in terms of the response function and the ground-state magnetization will be given in section 2.6.

Other conventions can be formulated as follows:

• For periodic systems, $[\hat{H}_{ia,jb}^{\sigma}]$ can be Fourier transformed to $[H_{\mu a, \nu b}^{\sigma}(\boldsymbol{k})]$, with µ and ν denoting the atomic positions within the primitive cell. Furthermore, the analysis throughout this paper assumes the use of the periodic gauge $H_{\mu a, \nu b}^{\sigma}(\boldsymbol{k}+\boldsymbol{G}) = H_{\mu a, \nu b}^{\sigma}(\boldsymbol{k})$ for any reciprocal lattice translation G [44].
• The n × n matrix $\hat{a}_{\mu}$, specified by the n orbital indices, can be viewed as the column vector $\vec{a}_{\mu}$ of the length n². The scalar product $\vec{a}_{\mu}^{\, \dagger} \cdot \vec{b}_{\mu}^{\phantom{\dagger}}$ is the shorthand notation for ${\mathrm{Tr}}_{L} \{\hat{a}_{\mu}\hat{b}_{\mu} \}$ (where $\vec{a}_{\mu}^{\, \dagger}$ is the row vector corresponding to the column vector $\vec{a}_{\mu}^{\, \phantom{\dagger}}$). In the case of Cartesian vectors (such as the magnetization matrix or the magnetic field interacting with the magnetization), the notation $\vec{\boldsymbol{a}}_{\mu}^{\, \dagger} \cdot \vec{\boldsymbol{b}}_{\mu}^{\phantom{\dagger}}$ stands for the regular scalar product with the summation over the orbital indices. The notation $\vec{\boldsymbol{a}}_{\phantom{\mu}}^{\, \dagger} \cdot \vec{\boldsymbol{b}}_{\phantom{\mu}}$ implies the summation over the atomic indices as well.
• The $n\times n \times m \times m$ tensor $\mathcal{A} = [\mathcal{A}_{ab,cd}]$, with first two orbitals (ab) residing on the site µ and last two orbitals (cd) residing on the site ν, can be viewed as the $n^{2} \times m^{2}$ matrix $\hat{\mathcal{A}}_{\mu \nu}$. The construction $\hat{\mathcal{A}}_{\mu \nu} \vec{b}_{\nu}$ implies the summation over the two orbital indices on the site ν.
• The 2nd derivative is a local probe and the interaction parameter depends on the point in which it is calculated. For instance, considering the simplest interaction energy $E = -J\text{cos} \varphi$ between two spins in the bond, the 2nd derivative near FM (ϕ = 0) and AFM ($\varphi = \pi$) configurations of spins will be, respectively, J and −J. Nevertheless, as it is typically done, we will additionally change the sign of interaction parameters for the antiferromagnetically coupled bonds, thus adopting the universal definition where J > 0 and $\lt0$ stands for the FM and AFM interactions, respectively.

The SO interaction is known to consist of the spin-diagonal, $\frac{\xi}{2}\hat{L}^{z}\hat{\sigma}^{z}$, as well as off-diagonal, $\frac{\xi}{2}(\hat{L}^{x}\hat{\sigma}^{x} + \hat{L}^{y}\hat{\sigma}^{y})$, parts, where ξ is the SO coupling parameter, $\hat{\boldsymbol{L}} = (\hat{L}^{x},\hat{L}^{y},\hat{L}^{z})$ is the vector of angular momenta, and $\hat{\boldsymbol{\sigma}} = (\hat{\sigma}^{x},\hat{\sigma}^{y},\hat{\sigma}^{z})$ is that of the Pauli matrices. The antisymmetric DM interaction $\boldsymbol{d} = (d^{x},d^{y},d^{z})$ emerges in the 1st order of ξ and generally one should be able to calculate all three vector projections onto x, y, and z (unless they are related by the symmetry properties). Nevertheless, by proper rotations of the coordinate frame, which transform xyz to $\overline{z}yx$ and zxy, d^x and d^y can be viewed as d^z in the new coordinate frame. Therefore, we need the numerical procedure only for calculating d^z . The important point in this respect was realized by Sandratskii [15], who suggested that in order to calculate d^z , it is sufficient to consider only spin-diagonal part of the SO coupling. His idea was based on a simple observation that DM interactions give rise to spiral magnetic structures [45], which can be regarded as 'eigenstates' of the spin model (1). Therefore, the energies of the model (1) can be uniquely specified by the vectors q , describing propagation of the spin spiral. Then, the same should hold for the electronic model, which is used for the mapping onto the spin one, and the spin spirals should be amount possible magnetic solutions of such model. In practical terms, this means that the electronic states should obey the generalized Bloch theorem, which combines translations with the SU(2) rotations of spins in the spiral texture [46]. Nevertheless, this theorem can be applied only if the spin is the good quantum number so that the Hamiltonian $\hat{H}^{\sigma}$ remains diagonal with respect to the spin indices σ. Therefore, such $\hat{H}^{\sigma}$ can include the diagonal part of the SO coupling, but not the off-diagonal one.

The method is suitable for the DM interactions, but not for the magnetic anisotropy, which emerges in the 2nd order of the SO coupling and typically include both diagonal and off-diagonal contributions. This is again in line with the idea of the spin-spiral approach: the DM interactions give rise to the spin spirals, while the magnetic anisotropy acts against them, by deforming the spin spirals and locking them to the crystallographic lattice [47, 48]. The alternative to the spin-spiral technique is to work in the real space separately for each magnetic bond [49–52]. Such methods, which have certain limitations, will be briefly considered in section 6.1.

As was already pointed out before, our basic idea is to 'excite' the spin spiral, rotating the ground-state magnetization $\hat{\boldsymbol{m}}_{\mathrm{GS}} = (0,0,\hat{m})$ as

$$\begin{align} \hat{\boldsymbol{m}}_{\boldsymbol{q},i} = \left( \theta \text{cos} \boldsymbol{q} \boldsymbol{R}_{i}, \theta \text{sin} \boldsymbol{q} \boldsymbol{R}_{i}, 1-\frac{\theta^2}{2} \right) \, \hat{m}, \end{align} \tag{ 4 }$$

(see figure 1) and evaluate the interactions between the transversal 'fluctuations' of the magnetization $\hat{\boldsymbol{m}}_{\boldsymbol{q},i}^{\perp} = ( \text{cos} \boldsymbol{q} \boldsymbol{R}_{i}, \text{sin} \boldsymbol{q} \boldsymbol{R}_{i}, 0 ) \, \theta \hat{m}$ for small θ. In order to induce such $\hat{\boldsymbol{m}}_{\boldsymbol{q},i}^{\perp}$, we have to apply the external field

$$\begin{align} \hat{\boldsymbol{h}}_{\boldsymbol{q},i} = \left( \text{cos} \left(\boldsymbol{q} \boldsymbol{R}_{i} + \alpha_{\boldsymbol{q}} \right), \text{sin} \left(\boldsymbol{q} \boldsymbol{R}_{i} + \alpha_{\boldsymbol{q}} \right), 0 \right) \hat{h} \end{align} \tag{ 5 }$$

Zoom In Zoom Out Reset image size

in the direction perpendicular to the ground state magnetization. If the SO coupling is included, $\hat{\boldsymbol{h}}_{\boldsymbol{q},i}$ is not necessarily parallel to $\hat{\boldsymbol{m}}_{\boldsymbol{q},i}^{\perp}$ as the latter can experience the effect of the DM interaction d^z , which tend to additionally rotate the magnetization in the xy plane. Therefore, the phases $\alpha_{\boldsymbol{q}}$ are needed to compensate the effect of DM interactions (see figure 1). For small $\alpha_{\boldsymbol{q}}$, $\hat{\boldsymbol{h}}_{\boldsymbol{q},i}$ can be written as

$$\begin{align} \hat{\boldsymbol{h}}_{\boldsymbol{q},i} \approx \hat{\boldsymbol{h}}_{\boldsymbol{q},i}^{0} + \alpha_{\boldsymbol{q}} \boldsymbol{n}^{z} \times \hat{\boldsymbol{h}}_{\boldsymbol{q},i}^{0}, \end{align} \tag{ 6 }$$

where $\hat{\boldsymbol{h}}_{\boldsymbol{q},i}^{0} = \left( \text{cos} \boldsymbol{q} \boldsymbol{R}_{i} , \text{sin} \boldsymbol{q} \boldsymbol{R}_{i}, 0 \right) \hat{h}$ and $\boldsymbol{n}^{z} = (0,0,1)$.

The corresponding energy can be evaluated in the framework of constrained SDFT [34, 54] as:

$$\begin{align} {\cal E}\left[\vec{\boldsymbol{m}}\right] = {\cal T}\left[\vec{\boldsymbol{m}}\right] + {\cal E}_{\mathrm{xc}}\left[\vec{\boldsymbol{m}}\right] + \frac{1}{2} \vec{\boldsymbol{h}}_{\boldsymbol{q}}^{\dagger} \cdot \left( \vec{\boldsymbol{m}} - \vec{\boldsymbol{m}}_{\boldsymbol{q}}^{\perp} \right), \end{align} \tag{ 7 }$$

where ${\cal T}$ and ${\cal E}_{\mathrm{xc}}$ are, respectively, the kinetic and xc energies, while the last term controls the size of the transversal magnetization $\vec{\boldsymbol{m}}_{\boldsymbol{q}}^{\perp}$. For simplicity, we drop here all irrelevant dependencies of ${\cal E}$ on the charge density.

Then, the kinetic energy can be expressed as the sum of the occupied KS single-particle energies, ${\cal E}_{\mathrm{sp}}$, calculated for the external field $\vec{\boldsymbol{h}}_{\boldsymbol{q}}$ and the xc field

$$\begin{align*} \vec{\boldsymbol{b}}_{\boldsymbol{q}} = 2 \frac{\partial {\cal E}_{\mathrm{xc}}\left[\vec{\boldsymbol{m}}\right]}{\partial \vec{\boldsymbol{m}}} \Big\vert_{\vec{\boldsymbol{m}} = \vec{\boldsymbol{m}}_{\boldsymbol{q}}}, \end{align*}$$

minus the interaction energy of $\vec{\boldsymbol{m}}_{\boldsymbol{q}}$ with $\vec{\boldsymbol{h}}_{\boldsymbol{q}}$ and $\vec{\boldsymbol{b}}_{\boldsymbol{q}}$ [13, 14], yielding

$$\begin{align} {\cal E}\left[\vec{\boldsymbol{m}}_{\boldsymbol{q}}\right] = {\cal E}_{\mathrm{sp}} \left( \vec{\boldsymbol{h}}_{\boldsymbol{q}} + \vec{\boldsymbol{b}}_{\boldsymbol{q}} \right) - \frac{1}{2} \left(\vec{\boldsymbol{h}}_{\boldsymbol{q}} + \vec{\boldsymbol{b}}_{\boldsymbol{q}}\right)^{\dagger} \cdot \vec{\boldsymbol{m}}_{\boldsymbol{q}} + {\cal E}_{\mathrm{xc}}\left[\vec{\boldsymbol{m}}_{\boldsymbol{q}}\right]. \end{align} \tag{ 8 }$$

The important property of the xc energy in this respect, being the consequence of fundamental gauge invariance of the density functional theory [55], is that rotations of the spin magnetization $\hat{\boldsymbol{m}}_{\mathrm{GS}} \rightarrow \hat{\boldsymbol{m}}_{\boldsymbol{q},i}$ do not change ${\cal E}_{\mathrm{xc}}$: ${\cal E}_{\mathrm{xc}}[ \hat{\boldsymbol{m}}_{\boldsymbol{q},i}] = {\cal E}_{\mathrm{xc}}[ \hat{\boldsymbol{m}}_{\mathrm{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 homogeneous electron gas. In this case, ${\cal E}_{\mathrm{xc}}$ in each point r depends only on the magnitude of the magnetization, ${\cal E}_{\mathrm{xc}}^{\mathrm{LSDA}} \equiv {\cal E}_{\mathrm{xc}}^{\mathrm{LSDA}}[|\boldsymbol{m}(\boldsymbol{r})|]$ [59, 60], and therefore does not change under rotations of $\boldsymbol{m}(\boldsymbol{r})$. Since ${\cal E}_{\mathrm{xc}}[ \hat{\boldsymbol{m}}_{\boldsymbol{q},i}] = {\cal E}_{\mathrm{xc}}[ \hat{\boldsymbol{m}}_{\mathrm{GS}} ]$, the rotation of the magnetization will rotate the xc field by the same angle:

$$\begin{align} \hat{\boldsymbol{b}}_{\boldsymbol{q},i} = \left( \theta \text{cos} \boldsymbol{q} \boldsymbol{R}_{i}, \theta \text{sin} \boldsymbol{q} \boldsymbol{R}_{i}, 1-\frac{\theta^2}{2} \right) \, \hat{b}. \end{align} \tag{ 9 }$$

Therefore, $\vec{\boldsymbol{b}}_{\boldsymbol{q}}^{\dagger} \cdot \vec{\boldsymbol{m}}_{\boldsymbol{q}}$ does not change either and equation (8) will lead to the following energy change:

$$\begin{align*} \delta{\cal E}_{\boldsymbol{q}} = \delta{\cal E}_{\mathrm{sp}} \left( \vec{\boldsymbol{h}}_{\boldsymbol{q}} + \vec{\boldsymbol{b}}_{\boldsymbol{q}} \right) - \frac{1}{2} \vec{\boldsymbol{h}}_{\boldsymbol{q}}^{\dagger} \cdot \vec{\boldsymbol{m}}_{\boldsymbol{q}}. \end{align*}$$

Then, $\delta {\cal E}_{\mathrm{sp}}$ can be evaluated by treating $\vec{\boldsymbol{h}}_{\boldsymbol{q}} + \vec{\boldsymbol{b}}_{\boldsymbol{q}}$ as a perturbation, to the 2nd order in $\vec{\boldsymbol{h}}_{\boldsymbol{q}}^{\phantom{\perp}} + \vec{\boldsymbol{b}}_{\boldsymbol{q}}^{\perp}$ and the 1st order in the longitudinal change of the xc field, $-\frac{1}{2} \hat{b} \theta^2$. The details are elaborated in [61], leading to the simple but general expression:

$$\begin{align} \delta {\cal E}_{\boldsymbol{q}} = -\frac{1}{4} \, \vec{\boldsymbol{h}}_{\boldsymbol{q}}^{\,0 \dagger} \cdot \vec{\boldsymbol{m}}_{\boldsymbol{q}}^{\perp }. \end{align} \tag{ 10 }$$

In fact, this result is well anticipated. On the one hand, $\delta {\cal E}_{\boldsymbol{q}}$ should be proportional to $\vec{\boldsymbol{h}}_{\boldsymbol{q}}$ as there would be no energy change without the external field. On the other hand, there only possible interaction of $\vec{\boldsymbol{h}}_{\boldsymbol{q}}$ with the constrained magnetization $\vec{\boldsymbol{m}}_{\boldsymbol{q}}^{\perp }$ is the scalar product given by equation (10). It may look incomplete because $\delta {\cal E}_{\boldsymbol{q}}$ does not seem to know anything about $\alpha_{\boldsymbol{q}}$ and the DM interactions. Nevertheless, all necessary information is in equation (10) and in section 2.5 we will show how it should be used to derive practical expressions for the isotropic exchange and DM interactions.

In addition to the rotations given by (4), the magnetization can experience the longitudinal change, which is caused by these rotations. It will affect $\hat{m}$, resulting in an additional change of each of the terms in equation (8). Nevertheless, these contributions can be shown to cancel out in the lowest order of θ [11, 57].

The response theory is basically the perturbation theory relating the small change of the potential $\vec{v}$ with the induced density $\vec{n}$: $\vec{n} = \hat{\mathcal{R}} \vec{v}$. Spin-dependent $\vec{v}$ can be generally specified by four elements:

$$\begin{align*} \vec{v} = \left( \begin{array}{cc} \vec{v}^{\, \uparrow \uparrow} & \vec{v}^{\, \uparrow \downarrow} \\ \vec{v}^{\, \downarrow \uparrow} & \vec{v}^{\, \downarrow \downarrow} \end{array} \right). \end{align*}$$

Then, each $\vec{v}^{\, \sigma \sigma^{\prime}}$ induces the corresponding change $\vec{n}^{\, \sigma \sigma^{\prime}}$:

$$\begin{align} \vec{n}^{\, \sigma \sigma^{\prime}} = \hat{\mathcal{R}}^{\sigma \sigma^{\prime}} \vec{v}^{\, \sigma \sigma^{\prime}}, \end{align} \tag{ 11 }$$

where the rank-4 tensor $\hat{\mathcal{R}}^{\sigma \sigma^{\prime}}$ can be found in terms of the 1st-order perturbation theory for the wave functions [62]. In our case, the perturbation is $\vec{\boldsymbol{h}}_{\boldsymbol{q}} + \vec{\boldsymbol{b}}_{\boldsymbol{q}}^{\perp}$ and our goal is to evaluate $\vec{\boldsymbol{m}}_{\boldsymbol{q}}^{\perp}$. Then, it is convenient to use the local coordinate frame where $\hat{\boldsymbol{h}}_{i} = ( 1, \alpha_{\boldsymbol{q}}, 0 ) \hat{h}^{0}$, which is obtained by rotating $\hat{\boldsymbol{h}}_{\boldsymbol{q},i}$ about z by the angles $- \boldsymbol{q} \boldsymbol{R}_{i}$, and employ the generalized Bloch theorem, combining lattice translations with the SU(2) rotations of spins [46]. This will lead to the additional shift of the k -mesh for the states with $\sigma = \uparrow$ relative to those with $\sigma = \downarrow$. Moreover, since $\hat{\boldsymbol{h}}_{i} = ( 1, \alpha_{\boldsymbol{q}}, 0 )\hat{h}^{0}$ corresponds to

$$\begin{align*} \vec{v}_{h} = \frac{1}{2} \left( \begin{array}{cc} 0 & 1-i\alpha_{\boldsymbol{q}} \\ 1+i\alpha_{\boldsymbol{q}} & 0 \end{array} \right)\vec{h}^{0}, \end{align*}$$

we have to consider only $\hat{\mathcal{R}}^{\uparrow \downarrow}$ and $\hat{\mathcal{R}}^{\downarrow \uparrow}$. Then, the perturbation theory yields

$$\begin{align} {\cal R}_{ab,cd}^{\uparrow \downarrow} \left(\boldsymbol{q}\right) = \sum_{ml \boldsymbol{k}} \frac{f_{m \boldsymbol{k}}^{\uparrow} - f_{l \boldsymbol{k}+\boldsymbol{q}}^{\downarrow}}{\varepsilon_{m \boldsymbol{k}}^{\uparrow} - \varepsilon_{l \boldsymbol{k}+\boldsymbol{q}}^{\downarrow}} \left(C_{m \boldsymbol{k}}^{a \, \uparrow}\right)^{*}C_{l \boldsymbol{k}+\boldsymbol{q}}^{b \, \downarrow} \left(C_{l \boldsymbol{k}+\boldsymbol{q}}^{c \, \downarrow}\right)^{*}C_{m \boldsymbol{k}}^{d \, \uparrow}, \end{align} \tag{ 12 }$$

where $\varepsilon_{m \boldsymbol{k}}^{\sigma}$ and $| C_{l \boldsymbol{k}}^{\sigma} \rangle = [ C_{l \boldsymbol{k}}^{a\sigma}]$ are, respectively, the eigenvalues and eigenvectors of $\hat{H}^{\sigma}$, and $f_{m \boldsymbol{k}}^{\sigma} \equiv \Theta(\varepsilon_{\mathrm{F}} - \varepsilon_{m \boldsymbol{k}}^{\sigma})$ is the Fermi distribution function. The orbital indices in each of the pairs ab and cd belong to the same atomic sites in the unit cell. Furthermore, ${\cal R}^{\downarrow \uparrow}$ can be obtained from ${\cal R}^{\uparrow \downarrow}$ using the property

$$\begin{align} {\cal R}_{ab,cd}^{\uparrow \downarrow} \left(\boldsymbol{q}\right) = \left[ {\cal R}_{ba,dc}^{\downarrow \uparrow} \left(-\boldsymbol{q}\right) \right]^{*}. \end{align} \tag{ 13 }$$

If $\hat{H}^{\sigma}$ remains invariant under the time reversal (e.g. without SO interaction), equation (13) is reduced to ${\cal R}_{ab,cd}^{\uparrow \downarrow} (\boldsymbol{q}) = {\cal R}_{ba,dc}^{\downarrow \uparrow} (\boldsymbol{q})$. ²

${\cal R}_{ab,cd}^{\uparrow \downarrow} (\boldsymbol{q})$ can be also related to the Green function $\hat{G}^{\sigma}(\varepsilon,\boldsymbol{k}) = [ \varepsilon - \hat{H}^{\sigma}(\boldsymbol{k}) ]^{-1}$ as [63]

$$\begin{align} \mathcal{R}^{\uparrow \downarrow}_{ab , cd} \left(\boldsymbol{q}\right) = -\frac{1}{\pi} \sum_{\boldsymbol{k}}^{\mathrm{BZ}} {\mathrm{Im}} \int_{- \infty}^{\varepsilon_\mathrm{F}} \mathrm{d} \varepsilon \left\{G_{da}^{\uparrow}\left(\varepsilon,\boldsymbol{k}\right) G_{bc}^{\downarrow}\left(\varepsilon,\boldsymbol{k}+\boldsymbol{q}\right) \right\} \end{align} \tag{ 14 }$$

with the summation running over the 1st Brillouin zone (BZ).

Then, using the definition (11), one can find:

$$\begin{align*} \vec{n}^{\, \uparrow \downarrow} = \frac{1}{2} \hat{\mathcal{R}}^{\uparrow \downarrow}_{\boldsymbol{q}} \left( \vec{h}^{0} + \vec{b}^{\perp} -i\alpha_{\boldsymbol{q}} \vec{h}^{0} \right) \end{align*}$$

and

$$\begin{align*} \vec{n}^{\, \downarrow \uparrow} = \frac{1}{2} \hat{\mathcal{R}}^{\downarrow \uparrow}_{\boldsymbol{q}} \left( \vec{h}^{0} + \vec{b}^{\perp} +i\alpha_{\boldsymbol{q}} \vec{h}^{0} \right), \end{align*}$$

where $\hat{\mathcal{R}}^{\sigma \sigma^{\prime}}_{\boldsymbol{q}} \equiv \hat{\mathcal{R}}^{\sigma \sigma^{\prime}}({\boldsymbol{q}})$. In the local coordinate frame, these $\vec{n}^{\, \uparrow \downarrow}$ and $\vec{n}^{\, \downarrow \uparrow}$ should give us the magnetization $\vec{m}^{\perp} = \vec{n}^{\, \uparrow \downarrow} + \vec{n}^{\, \downarrow \uparrow}$ along x:

$$\begin{align} \vec{m}^{\perp} = \hat{\mathcal{R}}^{+}_{\boldsymbol{q}} \left(\vec{h}^{0} + \vec{b}^{\perp}\right) - i \alpha_{\boldsymbol{q}} \hat{\mathcal{R}}^{-}_{\boldsymbol{q}} \vec{h}^{0}, \end{align} \tag{ 15 }$$

where $\hat{\mathcal{R}}^{\pm}_{\boldsymbol{q}} = \frac{1}{2} ( \hat{\mathcal{R}}^{\uparrow \downarrow}_{\boldsymbol{q}} \pm \hat{\mathcal{R}}^{\downarrow \uparrow}_{\boldsymbol{q}} )$. Another equation,

$$\begin{align} i \hat{\mathcal{R}}^{-}_{\boldsymbol{q}} \left(\vec{h}^{0} + \vec{b}^{\perp}\right)+ \alpha_{\boldsymbol{q}} \hat{\mathcal{R}}^{+}_{\boldsymbol{q}} \vec{h}^{0} = 0, \end{align} \tag{ 16 }$$

requires that the perpendicular to it magnetization along y, $i (\vec{n}^{\, \uparrow \downarrow}-\vec{n}^{\, \downarrow \uparrow} )$, should vanish (so as the y component of the xc field) according to our constraint conditions. These are the equations for $\vec{h}^{0}$ and $\alpha_{\boldsymbol{q}}$ for given $\vec{m}^{\perp} = \theta \vec{m}$. Their meaning is very straightforward. For instance, in equation (16), the isotropic part of the magnetization $\alpha_{\boldsymbol{q}} \hat{\mathcal{R}}^{+}_{\boldsymbol{q}} \vec{h}^{0}$, which is induced by $\alpha_{\boldsymbol{q}} \vec{h}^{0}$ along y, is compensated by the one, which is induced due to the DM interaction by the field $\vec{h}^{0}+\vec{b}^{\perp}$ acting in the perpendicular direction x. The same is with equation (15), where $\vec{m}^{\perp}$ has two components: the isotropic one, induced by $\vec{h}^{0}+\vec{b}^{\perp}$ along x, and the one caused by the DM interaction, transferring the effect of the magnetic field $\alpha_{\boldsymbol{q}} \vec{h}^{0}$, applied along y, to the magnetization along x. This explains how one can naturally separate the contributions of the isotropic and DM interactions in equation (10).

The next step is the mapping of the total energy change (10) onto the spin model:

$$\begin{align} \delta {\cal E}_{\boldsymbol{q}} = - \frac{1}{2} \sum_{\mu \nu} \left( J_{\boldsymbol{q}, \, \mu \nu} - id_{\boldsymbol{q}, \, \mu \nu}^{z} \right) \theta_{\mu} \theta_{\nu}, \end{align} \tag{ 17 }$$

where we explicitly consider the possibility of having several magnetic sublattices. In the local coordinate frame, equation (10) can be rearranged as

$$\begin{align*} \delta {\cal E}_{\boldsymbol{q}} = -\frac{1}{4} \, \left(\vec{h}^{\,0} + \vec{b}^{\perp }\right)^{\dagger} \cdot \vec{m}^{\perp } + \frac{1}{4}\vec{b}^{\perp \dagger} \cdot \vec{m}^{\perp }, \end{align*}$$

where we have added and subtracted the xc field $\vec{b}^{\perp }$. Our strategy is to start with the expression for $\vec{m}^{\perp}$ without the SO coupling, which is given by the 1st term in equation (15), and then consider the corrections arising in the 1st order of the SO coupling, which are given by the 2nd term. ³ Then, noting that

$$\begin{align} \vec{h}^{0} + \vec{b}^{\perp} = \hat{\mathcal{Q}}^{+}_{\boldsymbol{q}} \vec{m}^{\perp}, \end{align} \tag{ 18 }$$

where $\hat{\mathcal{Q}}^{\sigma \sigma^{\prime}}_{\boldsymbol{q}} = [ \hat{\mathcal{R}}^{\sigma \sigma^{\prime}}_{\boldsymbol{q}} ]^{-1}$, $\hat{\mathcal{Q}}^{\pm}_{\boldsymbol{q}} = \frac{1}{2} ( \hat{\mathcal{Q}}^{\uparrow \downarrow}_{\boldsymbol{q}} \pm \hat{\mathcal{Q}}^{\downarrow \uparrow}_{\boldsymbol{q}} )$, and without SO coupling $\hat{\mathcal{Q}}^{+}_{\boldsymbol{q}} = [\hat{\mathcal{R}}^{+}_{\boldsymbol{q}}]^{-1}$, one immediately finds the following expression for the isotropic exchange interactions:

$$\begin{align} J_{\boldsymbol{q}, \, \mu \nu} = \frac{1}{2} \left( \vec{m}^{\, \dagger}_{\mu} \cdot \hat{\mathcal{Q}}^{+}_{\boldsymbol{q}, \, \mu \nu} \vec{m}^{\, \phantom{\dagger}}_{\nu} - \vec{b}^{\, \dagger}_{\mu} \cdot \vec{m}^{\, \phantom{\dagger}}_{\mu} \delta_{\mu \nu} \right). \end{align} \tag{ 19 }$$

Considering the 2nd term in equation (15), the construction $i (\vec{h}^{\,0}+\vec{b}^{\perp})^{\dagger} \cdot \hat{\mathcal{R}}^{-}_{\boldsymbol{q}} \alpha_{\boldsymbol{q}} \vec{h}^{0}$ describes the interaction between x and y components of the magnetic field caused by the DM interactions. Corresponding interaction parameter should satisfy the condition

$$\begin{align} d_{\boldsymbol{q}, \, \mu \nu}^{z} \, \theta_{\mu} \theta_{\nu} = \frac{i}{2} \left(\vec{h}^{\,0}_{\mu}+\vec{b}^{\perp}_{\mu}\right)^{\dagger} \cdot \hat{\mathcal{R}}^{-}_{\boldsymbol{q}, \, \mu \nu} \left(\vec{h}^{\,0}_{\nu}+\vec{b}^{\perp}_{\nu}\right), \end{align} \tag{ 20 }$$

where we had to 'rescale' y components of the magnetic field, $\alpha_{\boldsymbol{q}} \vec{h}^{\, 0}_{\mu} \to \vec{h}^{\,0}_{\mu}+\vec{b}^{\perp}_{\mu}$, in order to specify x and y components of the transversal magnetization by the same set of parameters θ_µ and θ_ν . Then, using equation (18) and noting that to the 1st order in the SO coupling $\hat{\mathcal{Q}}^{+} \hat{\mathcal{R}}^{-} \hat{\mathcal{Q}}^{+} = -\hat{\mathcal{Q}}^{-}$, one can find that

$$\begin{align} d_{\boldsymbol{q}, \, \mu \nu}^{z} = -\frac{i}{2} \vec{m}^{\, \dagger}_{\mu} \cdot \hat{\mathcal{Q}}^{-}_{\boldsymbol{q}, \, \mu \nu} \vec{m}^{\, \phantom{\dagger}}_{\nu}. \end{align} \tag{ 21 }$$

This expression was obtained in [53] basically heuristically, by the analogy with isotropic interactions and similar expression formulated in terms of the xc fields, which will be considered in section 2.8. Here, we have provided a more rigorous proof of equation (21). The real space parameters can be obtained by the Fourier transform of $J_{\boldsymbol{q}, \, \mu \nu}$ and $d_{\boldsymbol{q}, \, \mu \nu}^{z}$.

Thus, the exchange interactions are proportional to the inverse response function. For the isotropic exchange, this is basically the result of Bruno [34]. For the Hubbard model, similar relationship has been established by Szczech et al [64].

For practical purposes, it may be more convenient to calculate

$$\begin{align} X_{\boldsymbol{q}, \, \mu \nu} = \frac{1}{2} \left( \vec{m}^{\, \dagger}_{\mu} \cdot \hat{\mathcal{Q}}^{\uparrow \downarrow}_{\boldsymbol{q}, \, \mu \nu} \vec{m}^{\, \phantom{\dagger}}_{\nu} - \vec{b}^{\, \dagger}_{\mu} \cdot \vec{m}^{\, \phantom{\dagger}}_{\mu} \delta_{\mu \nu} \right), \end{align} \tag{ 22 }$$

in terms of only $\hat{\mathcal{Q}}^{\uparrow \downarrow}_{\boldsymbol{q}}$, and then relate it with $J_{\boldsymbol{q}, \, \mu \nu}$ and $d_{\boldsymbol{q}, \, \mu \nu}^{z}$ using the property (13), which yields

$$\begin{align} J_{\boldsymbol{q}, \, \mu \nu} = \frac{1}{2} \left( X_{\boldsymbol{q}, \, \mu \nu} + X_{-\boldsymbol{q}, \, \mu \nu}^{*} \right) \end{align} \tag{ 23 }$$

and

$$\begin{align} d_{\boldsymbol{q}, \, \mu \nu}^{z} = -\frac{i}{2} \left( X_{\boldsymbol{q}, \, \mu \nu} - X_{-\boldsymbol{q}, \, \mu \nu}^{*} \right). \end{align} \tag{ 24 }$$

Thus, $J_{\boldsymbol{q}, \, \mu \nu}$ is related to the average energy of spin spirals propagating in q and $-\boldsymbol{q}$, while $d_{\boldsymbol{q}, \, \mu \nu}^{z}$ is related to the energy difference [15].

The sum rule is obtained from the identity

$$\begin{align*} \left[ \hat{G}^{\downarrow}\left(\varepsilon,\boldsymbol{k}\right) \right]^{-1} - \left[ \hat{G}^{\uparrow}\left(\varepsilon,\boldsymbol{k}\right) \right]^{-1} = \hat{b}, \end{align*}$$

which can be further rearranged as

$$\begin{align*} \hat{G}^{\uparrow}\left(\varepsilon,\boldsymbol{k}\right) - \hat{G}^{\downarrow}\left(\varepsilon,\boldsymbol{k}\right) & = \hat{G}^{\downarrow}\left(\varepsilon,\boldsymbol{k}\right) \, \hat{b} \, \hat{G}^{\uparrow}\left(\varepsilon,\boldsymbol{k}\right)\nonumber\\ & = \hat{G}^{\uparrow}\left(\varepsilon,\boldsymbol{k}\right) \, \hat{b} \, \hat{G}^{\downarrow}\left(\varepsilon,\boldsymbol{k}\right), \end{align*}$$

where $\hat{b} = \hat{H}^{\uparrow}(\boldsymbol{k})-\hat{H}^{\downarrow}(\boldsymbol{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:

$$\begin{align} \vec{m} = \hat{\mathcal{R}}^{\uparrow \downarrow}_{0} \vec{b}. \end{align} \tag{ 25 }$$

This sum rule has very straightforward meaning: $\boldsymbol{q} = 0$ corresponds to the uniform rotation of the ground-state magnetization, where all spins are rotated in the same direction by the same angle. Therefore, the transversal magnetization is described by the same xc field $\vec{b}$ as in the ground state (without any constraining fields).

Nevertheless, in the TB representation, such xc field is not necessary local. For instance, in LSDA, the splitting $\hat{H}^{\uparrow}(\boldsymbol{k})-\hat{H}^{\downarrow}(\boldsymbol{k})$ can have interatomic matrix elements and depend on k . In such a situation, it can be important to reenforce the sum rule, by defining new local xc field as $\vec{b} = \hat{\mathcal{Q}}^{\uparrow \downarrow}_{0} \vec{m}$, which would yield the given ground-state magnetization $\vec{m}$. For instance, this is a simple and transparent alternative to the kernel polynomial method, which was recently proposed to deal with nonlocal matrix elements of the xc field [30]. In fact, if the xc field is nonlocal, the total energy change for the infinitesimal rotations of spins is no longer representable in the form of pairwise interactions.

So far, we did not properly specify the spin object which should be rotated on the magnetic sites in order to obtain the total energy change (10). All above discussions implied that it is the magnetization matrix $\hat{m}_{\mu}$, while the spin model is typically formulated in terms of the magnetic moments $M_{\mu} = {\mathrm{Tr}}_{L} \{\hat{m}_{\mu} \}$. Undoubtedly, the rotation of $\hat{m}_{\mu}$, as a whole, by the angle θ_µ will rotate M_µ by the same angle. However, is this choice unique? Are there other perturbations of $\hat{m}_{\mu}$, resulting in the same rotations of M_µ but preferably at lower energy cost? Here, we will follow the discussion in [61]. Nevertheless, we would like to note that somewhat similar ideas can be found in the work of Antropov et al [65].

Indeed, for the Hermitian matrix $\hat{m}_{\mu}$, one can always choose the diagonal representation $\hat{m}_{\mu} = {\mathrm{diag}} ( \, \dots, \, m^{a}_{\mu}, \, \dots )$ with respect to the orbital indices. In principle, each orbital a in such representation can be rotated by its own angle $\theta_{\mu}^{a}$. Then, the transversal magnetization in the local coordinate frame, where it is parallel to x, will be $\hat{m}_{\mu}^{\perp} = {\mathrm{diag}} ( \, \dots, \, \theta_{\mu}^{a} m^{a}_{\mu}, \, \dots )$. Nevertheless, these angles are subjected to the additional constraint because ${\mathrm{Tr}}_{L} \{\hat{m}_{\mu}^{\perp} \}$ should be equal to $\theta_{\mu} M_{\mu}$. Importantly, this condition is softer than the rotation of $\hat{m}_{\mu}$ as a whole, where all orbitals are rotated by the same $\theta_{\mu}^{a} = \theta_{\mu}$. Therefore, it is reasonable to expect that the energy change will be smaller, so as the exchange parameters.

Mathematically, we have to minimize the energy change (10) with the additional condition $\sum_{a} \left( \theta^{a}_{\mu} - \theta_{\mu}^{\phantom{a}} \right) m^{a}_{\mu} = 0$ on each site µ:

$$\begin{align} \delta {\cal E} = - \frac{1}{4} \sum_{\mu a} \left\{\theta^{a}_{\mu} m^{a}_{\mu} h^{0 a}_{\mu} - \left( \theta^{a}_{\mu} - \theta_{\mu}^{\phantom{a}} \right) m^{a}_{\mu} \lambda_{\mu} \right\}, \end{align} \tag{ 26 }$$

where $h^{0 a}_{\mu}$ are the constraining fields acting on $m^{a}_{\mu}$ and λ_µ are the Lagrange multipliers. Then, minimizing $\delta {\cal E}$ with respect to $\theta^{a}_{\mu}$, it is straightforward to find that $h^{0 a}_{\mu} = \lambda_{\mu}$. Thus, in order to rotate the spin moments at the minimal energy cost, one have to apply the scalar field h_µ , i.e. the same for all orbitals a. Moreover, in this case it is convenient to use the 'spherically averaged' version of the linear response, where $\hat{m}_{\mu}$ is replaced by M_µ , $\hat{b}_{\mu}$ is replaced by $B_{\mu} = \frac{1}{n}{\mathrm{Tr}}_{L} \{\hat{b}_{\mu} \}$, and ${\cal R}_{ab,cd}^{\uparrow \downarrow}(\boldsymbol{q})$ is replaced by

$$\begin{align*} \mathbb{R}^{\uparrow \downarrow}_{\mu \nu}\left(\boldsymbol{q}\right) = \sum_{a \in \mu, c \in \nu} {\cal R}_{aa,cc}^{\uparrow \downarrow} \left(\boldsymbol{q}\right). \end{align*}$$

The corresponding exchange interaction parameters will be given by

$$\begin{align} J_{\boldsymbol{q}, \, \mu \nu} = \frac{1}{2} \left( M_{\mu} \unicode{x211A}^{+}_{\boldsymbol{q}, \, \mu \nu} M_{\nu} - B_{\mu} M_{\mu} \delta_{\mu \nu} \right) \end{align} \tag{ 27 }$$

and

$$\begin{align} d_{\boldsymbol{q}, \, \mu \nu}^{z} = -\frac{i}{2} M_{\mu} \unicode{x211A}^{-}_{\boldsymbol{q}, \, \mu \nu} M_{\nu}, \end{align} \tag{ 28 }$$

where $\hat{\unicode{x211A}}^{\pm}_{\boldsymbol{q}} = \frac{1}{2} ( \hat{\unicode{x211A}}^{\uparrow \downarrow}_{\boldsymbol{q}} \pm \hat{\unicode{x211A}}^{\downarrow \uparrow}_{\boldsymbol{q}} )$, $\hat{\unicode{x211A}}^{\sigma \sigma^{\prime}}_{\boldsymbol{q}} = [ \hat{\mathbb{R}}^{\sigma \sigma^{\prime}}_{\boldsymbol{q}} ]^{-1}$, and $\hat{\mathbb{R}}^{\sigma \sigma^{\prime}}_{\boldsymbol{q}} \equiv \hat{\mathbb{R}}^{\sigma \sigma^{\prime}}({\boldsymbol{q}})$ is the matrix specified by the atomic indices in the unit cell.

In comparison with equations (19) and (21), based on rotations of the magnetization matrix, equations (27) and (28) are expected to be more suitable for the analysis of low-energy excitations, of course, provided that the latter can be described by the spin model (1). In the following, these two methods will be denoted as $\hat{m}$ and M, after the basic variable describing the infinitesimal rotations of spins.

In this section, we consider the original formulation of the linear response theory, as it was proposed by Liechtenstein et al [11], which is frequently called the magnetic force theorem [34]. However, there are two important points about the work of Liechtenstein et al [11], which should be distinguished from each other [66]:

• The general claim that, in SDFT, the energy change caused by the infinitesimal rotations of spins can be related to the KS eigenstates in the ground state is certainly correct and should not be revised. This is what is actually called the 'magnetic force theorem' stating that the interaction parameters are the ground state properties and can be found by knowing the electronic structure in the ground state;
• Nevertheless, the practical expression (3), which was derived by Liechtenstein et al [11] for the exchange interactions, relies on additional approximations and, in principle, can be improved.

The starting assumption of Liechtenstein et al [11] is that since the rotation of the magnetization results in the rotation of the xc field by the same angle (section 2.3), it is logical to treat the change of the xc field, $\delta \hat{\boldsymbol{b}}_{\boldsymbol{q}} = \hat{\boldsymbol{b}}_{\boldsymbol{q}} - \hat{\boldsymbol{b}}_{\mathrm{GS}}$, as a perturbation without the constraining field. Then, we have to consider only $\delta {\cal E}_{\mathrm{sp}}$ in equation (9) caused by this $\delta \hat{\boldsymbol{b}}_{\boldsymbol{q}}$. Furthermore, the transversal magnetization, $\vec{\boldsymbol{m}}_{\boldsymbol{q}}^{\perp ^{\prime}}$, which is induced by $\delta \hat{\boldsymbol{b}}_{\boldsymbol{q}}$ will generally differ from $\vec{\boldsymbol{m}}_{\boldsymbol{q}}^{\perp}$ because, without the constraining field, $\vec{\boldsymbol{m}}_{\boldsymbol{q}}$ will tend to relax toward the ground state [33–35]. The DM interactions, if any, will tend to additionally rotate $\vec{\boldsymbol{m}}_{\boldsymbol{q}}^{\perp^{^{\prime}}}$ relative to $\vec{\boldsymbol{b}}_{\boldsymbol{q}}^{\perp}$: $\vec{\boldsymbol{m}}_{\boldsymbol{q}}^{\perp^{^{\prime}}} \approx \vec{\boldsymbol{m}}_{\boldsymbol{q}}^{\perp 0^{^{\prime}}} + \beta_{\boldsymbol{q}} \boldsymbol{n}^{z} \times \vec{\boldsymbol{m}}_{\boldsymbol{q}}^{\perp 0^{^{\prime}}}$. Thus, instead of equation (10), this method relies on (in the local coordinate frame)

$$\begin{align} \delta {\cal E}_{\boldsymbol{q}} \approx \frac{1}{4} \, \vec{b}^{\perp \dagger} \cdot \vec{m}^{\perp^{^{\prime}}} - \frac{1}{4} \, \vec{b}^{\, \dagger} \cdot \vec{m} \, \theta^2, \end{align} \tag{ 29 }$$

arising from the single-particle energies for the perturbations caused by the transversal and longitudinal parts of the xc field. Again, the important point here is that $\vec{m}^{\perp^{^{\prime}}}$ deviates from $\vec{m}^{\perp}$. Otherwise, the right-hand side of equation (29) would identically be equal to zero, as was discussed in section 2.3. Then, using the definitions $\vec{m}^{\perp 0^{^{\prime}}} = \hat{\mathcal{R}}^{+}_{\boldsymbol{q}} \vec{b}^{\perp}$ and $\beta_{\boldsymbol{q}} \vec{m}^{\perp 0^{^{\prime}}} = i\hat{\mathcal{R}}^{-}_{\boldsymbol{q}} \vec{b}^{\perp}$, one can find that

$$\begin{align} J_{\boldsymbol{q}, \, \mu \nu} = -\frac{1}{2} \left( \vec{b}^{\, \dagger}_{\mu} \cdot \hat{\mathcal{R}}^{+}_{\boldsymbol{q}, \, \mu \nu} \vec{b}^{\, \phantom{\dagger}}_{\nu} - \vec{b}^{\, \dagger}_{\mu} \cdot \vec{m}^{\, \phantom{\dagger}}_{\mu} \delta_{\mu \nu} \right), \end{align} \tag{ 30 }$$

and

$$\begin{align} d_{\boldsymbol{q}, \, \mu \nu}^{z} = \frac{i}{2} \vec{b}^{\, \dagger}_{\mu} \cdot \hat{\mathcal{R}}^{-}_{\boldsymbol{q}, \, \mu \nu} \vec{b}^{\, \phantom{\dagger}}_{\nu}. \end{align} \tag{ 31 }$$

Alternatively, $d_{\boldsymbol{q}, \, \mu \nu}^{z}$ can be obtained from equation (20) for $\vec{h}^{\,0}_{\mu} = 0$ and $\vec{b}^{\perp}_{\mu} = \vec{b}_{\mu} \theta_{\mu}$. Substituting equation (14) into equation (30) and Fourier transforming it to the real space, one obtains the well-known equation (3). Nevertheless, these are the approximate expressions, which can be formally obtained from the exact ones, equations (30) and (31), replacing $\vec{m}$ by $\vec{b}$ and $\hat{\mathcal{Q}}$ by $\hat{\mathcal{R}}$, with the additional minus sign. In the following, we will refer to this method as 'method $\hat{b}$' or 'approximate method $\hat{b}$'.

In principle, one can also introduce the 'spherically averaged' version of this method (the so-called method B) replacing $\vec{b}_{\nu}$ by B_ν and $\hat{\mathcal{R}}$ by $\mathbb{R}$ [34], though it is rarely used. Without the SO coupling, $\hat{\mathbb{R}}^{+} = \hat{\mathbb{R}}^{\uparrow \downarrow}$ and corresponding exchange interactions $\hat{J}_{\boldsymbol{q}}^{B} \equiv [ J_{\boldsymbol{q}, \, \mu \nu}^{B} ]$ can be written as $\hat{J}_{\boldsymbol{q}}^{B} = - \frac{1}{2} \hat{B} ( \, \hat{\mathbb{R}}_{\boldsymbol{q}}^{\uparrow \downarrow} + \hat{\cal I}^{-1} \, ) \hat{B}$, where $\hat{B} = {\mathrm{diag}} ( \, \dots, \, B_{\nu}, \, \dots )$ and $\hat{\cal I} = {\mathrm{diag}} ( \, \dots, \, -B_{\nu}/M_{\nu}, \, \dots )$ are the diagonal matrices of, respectively, exchange splittings and effective Stoner parameters. By adapting the same matrix form for the 'exact' interactions (27), $\hat{J}_{\boldsymbol{q}} = \frac{1}{2} \hat{M} ( \, [\hat{\mathbb{R}}_{\boldsymbol{q}}^{\uparrow \downarrow} ]^{-1} + \hat{\cal I} \, ) \hat{M}$ with $\hat{M} = {\mathrm{diag}} ( \, \dots, \, M_{\nu}, \, \dots )$, one can find the following expression, connecting $\hat{J}_{\boldsymbol{q}}^{\phantom{B}}$ with $\hat{J}_{\boldsymbol{q}}^{B}$ [34]:

$$\begin{align} \hat{J}_{\boldsymbol{q}}^{\phantom{B}} = \hat{J}_{\boldsymbol{q}}^{B} \left( 1 - 2 \hat{B}^{-1} \hat{M}^{-1} \hat{J}_{\boldsymbol{q}}^{B} \right)^{-1}. \end{align} \tag{ 32 }$$

Thus, $\hat{J}_{\boldsymbol{q}}^{\phantom{B}}$ can be indeed replaced by $\hat{J}_{\boldsymbol{q}}^{B}$ at least in two cases: (i) the long wavelength limit $\boldsymbol{q} \to 0$ and (ii) the strong-coupling limit $\hat{B} \to \infty$. Therefore, the spin-wave stiffness in the limit $\boldsymbol{q} \to 0$ is expected to be the same in both methods. Nevertheless, this statement should not be exaggerated because equation (32) holds only in the spherical case, where the xc field and the magnetization on each magnetic site are given by the scalar parameters B_ν and M_ν . In the matrix case, the simple relationship (32) is no longer valid [61]. That is why even the spin-wave stiffness in the methods $\hat{b}$ and M can be different. Furthermore, we will see that there is indeed a number examples, where the approximate method $\hat{b}$ fails to reproduce the correct magnetic ground state, while the method M dramatically improves the description.

Of course, it is reasonable to ask what are the right objects to rotate in this case: whether they should be the whole matrices $\hat{b}_{\nu}$ or only the spherical parts of these matrices B_ν ? If in the case of the magnetization, the answer can be found 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 of the applications deal with the rotations of the matrices $\hat{b}_{\nu}$.

Considering the strong-coupling limit in equation (30) [67, 68], one can derive all known expressions for the double exchange $J_{ij} \sim \langle \hat{H}_{ij} \rangle$ [16], superexchange $J_{ij} \sim -\langle \hat{H}_{ij}^{2} \rangle/U$ [2], superexchange with the interatomic Coulomb repulsion $J_{ij} \sim -\langle \hat{H}_{ij}^{2} \rangle/(U-V)$ [69], etc where U and V is the on-site and intersite Coulomb repulsion, respectively, $\hat{H}_{ij}$ are the transfer integrals (see appendix A), which are typically associated with the matrix elements of the KS Hamiltonian in LDA (GGA), and $\langle \dots \rangle$ denotes the expectation value in the ground state. The strong-coupling limit for the DM interaction (31) results in the spin-current model [58, 70], which can be viewed as the relativistic counterpart of the double exchange mechanism [53]. The expression for RKKY interactions can be also derived starting from equation (30), but using slightly different philosophy [9]. In this case, $\vec{b}$ is the field created by localized core spins and acting on outer conduction electrons. Without $\vec{b}$, the conduction bands are non-magnetic (and the tensor $\hat{\mathcal{R}}$ is evaluated in this non-magnetic state [9]).

In the previous sections, we have considered how the spin model can be generally derived from the electronic one using the concept of infinitesimal rotations of spins. The parameters of such spin model are expressed in terms of the static spin susceptibility (or the response function). On the other hand, the spin-wave dispersion, $\omega_{\boldsymbol{q}}$, which is the experimentally measurable quantity, can be derived in the framework of RPA from the poles of the dynamic spin susceptibility [71–73]. However, this $\omega_{\boldsymbol{q}}$ does not necessary coincide with the one of the spin model with the parameters derived from the static spin susceptibility [63, 64]. In this respect, Katsnelson and Lichtenstein [63] have argued that although the method proposed by Bruno [34] is more consistent with the static response formulation, the method $\hat{b}$ should more suitable for the analysis of the spin-wave spectra. Here, we will briefly consider this problem. For simplicity, we assume that there is only one magnetic site in the unit cell and drop all matrix notations. Then, in the spherical case, the dynamic response function is given by $\tilde{\mathbb{R}}_{\boldsymbol{q}}^{\uparrow \downarrow}(\omega) = \mathbb{R}_{\boldsymbol{q}}^{\uparrow \downarrow}(\omega)[ 1 + {\cal I} \mathbb{R}_{\boldsymbol{q}}^{\uparrow \downarrow}(\omega)]^{-1}$, where $\mathbb{R}_{\boldsymbol{q}}^{\uparrow \downarrow}(\omega)$ is obtained from equation (12) replacing in the denominator $(\varepsilon_{m \boldsymbol{k}}^{\uparrow}-\varepsilon_{l \boldsymbol{k}+\boldsymbol{q}}^{\downarrow})$ by $(\omega+\varepsilon_{m \boldsymbol{k}}^{\uparrow}-\varepsilon_{l \boldsymbol{k}+\boldsymbol{q}}^{\downarrow})$. Therefore, one has to solve the equation

$$\begin{align} 1+ {\cal I} \mathbb{R}_{\boldsymbol{q}}^{\uparrow \downarrow}\left(\omega_{\boldsymbol{q}}\right) = 0, \end{align} \tag{ 33 }$$

which can be equivalently rearranged as: $2M^{-1}J_{\boldsymbol{q}}(\omega_{\boldsymbol{q}}) = 0$, where $J_{\boldsymbol{q}}(\omega)$ is given by equation (27) with $\unicode{x211A}^{\uparrow \downarrow}_{\boldsymbol{q}}(\omega)$ instead of $\unicode{x211A}^{\uparrow \downarrow}_{\boldsymbol{q}} = \unicode{x211A}^{+}_{\boldsymbol{q}}$. Thus, the problem is that the spin-wave energies are given by the zeros of $J_{\boldsymbol{q}}(\omega)$ and do not necessary coincide with $\omega_{\boldsymbol{q}}^{M} = 2M^{-1}J_{\boldsymbol{q}}(0)$, expected from the solution of the spin model.

In the limit $\omega \ll \varepsilon_{l \boldsymbol{k}+\boldsymbol{q}}^{\downarrow}-\varepsilon_{m \boldsymbol{k}}^{\uparrow}$ (which takes place, for instance, for insulating and half-metallic materials, where the occupied and unoccupied states with opposite projections of spins are separated by an energy gap), one can use the linearization $\mathbb{R}_{\boldsymbol{q}}^{\uparrow \downarrow}(\omega) \approx \mathbb{R}_{\boldsymbol{q}}^{\uparrow \downarrow}(0) + \omega \dot{\mathbb{R}}_{\boldsymbol{q}}^{\uparrow \downarrow}(0)$, where $\dot{\mathbb{R}}_{\boldsymbol{q}}^{\uparrow \downarrow}(0) = \frac{\partial}{\partial \omega} \mathbb{R}_{\boldsymbol{q}}^{\uparrow \downarrow}(\omega)\vert_{\omega = 0}$, and find the following expression: $\omega_{\boldsymbol{q}} = \omega_{\boldsymbol{q}}^{B} {\cal I}^{-1} [\dot{\mathbb{R}}_{\boldsymbol{q}}^{\uparrow \downarrow}(0)]^{-1}B^{-1}$, where $\omega_{\boldsymbol{q}}^{B} = 2M^{-1}J^B_{\boldsymbol{q}}(0)$ is the spin-wave dispersion calculated with the parameters of the scheme B. This example clearly shows that $\omega_{\boldsymbol{q}}^{B}$ should be additionally renormalized, though this renormalization is generally different from the one given by equation (32), connecting the parameters of the methods M and B. $\dot{\mathbb{R}}_{\boldsymbol{q}}^{\uparrow \downarrow}(0)$ depends on the details of the electronic structure. In practical terms, it can be calculated replacing $(\varepsilon_{m \boldsymbol{k}}^{\uparrow}-\varepsilon_{l \boldsymbol{k}+\boldsymbol{q}}^{\downarrow})$ by $-(\varepsilon_{m \boldsymbol{k}}^{\uparrow}-\varepsilon_{l \boldsymbol{k}+\boldsymbol{q}}^{\downarrow})^2$ in equation (12). Then, in certain circumstances, the method M can be a good starting point for the analysis of the spin-wave dispersion. For instance, if the $\uparrow$-spin ($\downarrow$-spin) states are fully occupied (empty) and B is large compared to the band dispersion, it is straightforward to obtain that $\dot{\mathbb{R}}_{\boldsymbol{q}}^{\uparrow \downarrow}(0) \approx -B^{-1} \mathbb{R}_{\boldsymbol{q}}^{\uparrow \downarrow}(0)$ and $\omega_{\boldsymbol{q}} \approx \omega_{\boldsymbol{q}}^{M}$.

Thus, it would be fair to conclude that the analysis of the spin-wave dispersion requires the additional renormalization of the parameters derived from the static spin susceptibility [63, 64]. This conclusion applies to all methods (M, B, and $\hat{b}$). Therefore, this is an open question which method serves better for the analysis of the spin-wave dispersion. It is certainly true that, in LSDA, the method $\hat{b}$ better reproduces the experimental spin-wave dispersion in the canonical case of bcc-Fe and fcc-Ni [61, 63]. However, this conclusion does not seem to be general and for other materials the comparison can be less favorable.

Finally, we note that equation (33) can be further rearranged as $\mathbb{R}_{\boldsymbol{q}}^{\uparrow \downarrow}(0)(B+h_{\omega}) = M$, where $h_{\omega} = [\mathbb{R}_{\boldsymbol{q}}^{\uparrow \downarrow}(0)]^{-1} [\mathbb{R}_{\boldsymbol{q}}^{\uparrow \downarrow}(\omega) - \mathbb{R}_{\boldsymbol{q}}^{\uparrow \downarrow}(0)] B \approx \omega [\mathbb{R}_{\boldsymbol{q}}^{\uparrow \downarrow}(0)]^{-1}\dot{\mathbb{R}}_{\boldsymbol{q}}^{\uparrow \downarrow}(0) B$ has a meaning of the constraining field, which is needed to correct the effect of the xc field B in order to reproduce the ground-state magnetic moment for an arbitrary q . In this sense, there is an analogy between the search of the poles of the dynamic susceptibility and the constrained SDFT considered in section 2.3.

By knowing $J_{\boldsymbol{q}, \, \mu \nu}$ and $d_{\boldsymbol{q}, \, \mu \nu}^{z}$, one can, in principle, calculate isotropic and DM interactions operating between all sites in the unit cell. Nevertheless, the magnetization at these sites may have completely different origin. For instance, the transition-metal (${\mathrm{T}}$) sites in many oxide materials participate as a source of the magnetism, being primarily responsible for the spontaneous time-reversal symmetry breaking, while the oxygen or any other ligand (${\mathrm{L}}$) sites behave as 'magnetic slaves': although they can host an appreciable portion of the magnetization, it is solely induced by hybridization with the ${\mathrm{T}}$ sites and strictly follow the change of the magnetization on the ${\mathrm{T}}$ sites. The corresponding energy change for each q can be schematically expressed as

$$\begin{align} \delta {\cal E} = - \frac{1}{2} \left( \vec{\theta}_{\mathrm{T}}^{\, \dagger} \cdot \hat{X}_{\mathrm{TT}}^{\phantom{T}} \vec{\theta}_{\mathrm{T}}^{\phantom{\dagger}} + \vec{\theta}_{\mathrm{T}}^{\, \dagger} \cdot \hat{X}_{\mathrm{TL}}^{\phantom{T}} \vec{\theta}_{\mathrm{L}}^{\phantom{\dagger}} + \vec{\theta}_{\mathrm{L}}^{\, \dagger} \cdot \hat{X}_{\mathrm{LT}}^{\phantom{T}} \vec{\theta}_{\mathrm{T}}^{\phantom{\dagger}} + \vec{\theta}_{\mathrm{L}}^{\, \dagger} \cdot \hat{X}_{\mathrm{LL}}^{\phantom{T}} \vec{\theta}_{\mathrm{L}}^{\phantom{\dagger}} \right), \end{align} \tag{ 34 }$$

where $\hat{X}_{\mathrm{AB}}^{\phantom{T}}$ is the matrix specified by the atomic sites of the types ${\mathrm{A}}$ or ${\mathrm{B}}$, $\vec{\theta}_{\mathrm{A}}^{\phantom{\dagger}}$ are the polar angles specifying the rotations of magnetic moments of the type ${\mathrm{A}}$ in the form of the column vector, and $\vec{\theta}_{\mathrm{A}}^{\, \dagger}$ is the corresponding to it row vector. Then, one can try to eliminate the ${\mathrm{L}}$ degrees of freedom by transferring their effect into the interaction parameters between the ${\mathrm{T}}$ sites. This can be done by employing the ideas of adiabatic spin dynamics [74, 75] and assuming that the ${\mathrm{T}}$ spins are sufficiently slow so that the ${\mathrm{L}}$ spins have sufficient time to adjust each change in the system of the ${\mathrm{T}}$ spins. Mathematically, this means that for each $\vec{\theta}_{\mathrm{T}}^{\phantom{\dagger}}$, $\vec{\theta}_{\mathrm{L}}^{\phantom{\dagger}}$ can be found from the condition $\frac{\partial}{\partial \vec{\theta}_{\mathrm{L}}^{\, \dagger}} \delta {\cal E} = 0$, yielding

$$\begin{align} \vec{\theta}_{\mathrm{L}} = - \left[ \hat{X}_{\mathrm{LL}} \right]^{-1} \hat{X}_{\mathrm{LT}} \vec{\theta}_{\mathrm{T}} \end{align} \tag{ 35 }$$

and $\delta {\cal E} = - \frac{1}{2} \vec{\theta}_{\mathrm{T}}^{\, \dagger} \cdot \hat{\tilde{X}}_{\mathrm{TT}}^{\phantom{T}} \vec{\theta}_{\mathrm{T}}^{\phantom{\dagger}}$ with

$$\begin{align} \hat{\tilde{X}}_{\mathrm{TT}} = \hat{X}_{\mathrm{TT}} - \hat{X}_{\mathrm{TL}} \left[ \hat{X}_{\mathrm{LL}} \right]^{-1} \hat{X}_{\mathrm{LT}}. \end{align} \tag{ 36 }$$

The corresponding parameters of isotropic and DM interactions can be obtained from $\hat{\tilde{X}}_{\mathrm{TT}}$ 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 $\hat{\unicode{x211A}}^{\uparrow \downarrow} = [ \hat{\mathbb{R}}^{\uparrow \downarrow}]^{-1}$ to obtain the following expression for $\hat{\tilde{X}}_{\mathrm{TT}}$:

$$\begin{align} \hat{\tilde{X}}_{\mathrm{TT}} = \hat{\tilde{X}}_{\mathrm{TT}}^{0} + \Delta \hat{\tilde{X}}_{\mathrm{TT}}, \end{align} \tag{ 37 }$$

where

$$\begin{align} \hat{\tilde{X}}_{\mathrm{TT}}^{0} = \hat{M}_{\mathrm{T}}^{\phantom{\uparrow}} \left[ \hat{\mathbb{R}}^{\uparrow \downarrow}_{\mathrm{TT}} \right]^{-1} \hat{M}_{\mathrm{T}}^{\phantom{\uparrow}}, \end{align} \tag{ 38 }$$

and

$$\begin{align} \Delta \hat{\tilde{X}}_{\mathrm{TT}} & = \hat{M}_{\mathrm{T}}^{\phantom{\uparrow}} \left[ \hat{\mathbb{R}}^{\uparrow \downarrow}_{\mathrm{TT}} \right]^{-1} \hat{\mathbb{R}}^{\uparrow \downarrow}_{\mathrm{TL}} \hat{\unicode{x211A}}^{\uparrow \downarrow}_{\mathrm{LL}} \left( \hat{\unicode{x211A}}^{\uparrow \downarrow}_{\mathrm{LL}} + \hat{\cal I}_{\mathrm{L}}^{\phantom{\uparrow}} \right)^{-1}\nonumber\\ & \quad \times \hat{\cal I}_{\mathrm{L}}^{\phantom{\uparrow}} \hat{\mathbb{R}}^{\uparrow \downarrow}_{\mathrm{L}T} \left[ \hat{\mathbb{R}}^{\uparrow \downarrow}_{\mathrm{TT}} \right]^{-1} \hat{M}_{\mathrm{T}}^{\phantom{\uparrow}}. \end{align} \tag{ 39 }$$

Here, $\hat{M}_{\mathrm{T}}$ is the diagonal matrix of magnetic moments on the sites ${\mathrm{T}}$, and $\hat{\cal I}_{\mathrm{L}}$ is the diagonal matrix of effective Stoner parameters on the sites ${\mathrm{L}}$. In this expression, the explicit dependence of $\hat{\tilde{X}}_{\mathrm{TT}}$ on $\hat{\cal I}_{\mathrm{L}}$ is incorporated into $\Delta \hat{\tilde{X}}_{\mathrm{TT}}$, while $\hat{\tilde{X}}_{\mathrm{TT}}^{0}$ formally does not depend on $\hat{\cal I}_{\mathrm{L}}$. The parameters $\hat{\cal I}_{\mathrm{L}}$ can play a very important role in the theory of exchange interactions. For instance, according to GKA rules, they largely contribute to the FM coupling in the systems, where the ${\mathrm{T}}$–${\mathrm{L}}$–${\mathrm{T}}$ bond angle is close to 90° [36–39]. In the linear response theories, these effects are incorporated into $\Delta \hat{\tilde{X}}_{\mathrm{TT}}$ [76].

Furthermore, the representation (37) allows us to improve the numerical accuracy. Since in most of the systems the ligand band is filled, the matrix elements of $\hat{\mathbb{R}}^{\uparrow \downarrow}$ associated with the ${\mathrm{L}}$ sites are typically small. Therefore, when we invert $\hat{\mathbb{R}}^{\uparrow \downarrow}$, we have to deal with very large numbers even for the ${\mathrm{T}}$ sublattice. This is the reason why the bare interactions $\hat{X}_{\mathrm{TT}}^{\phantom{\uparrow}}$ $\sim$ $\hat{\unicode{x211A}}^{\uparrow \downarrow}_{\mathrm{TT}}$ are typically large and strongly compensated by the second term in equation (36) [61, 76]. Similar situation occurs in the scheme $\hat{m}$. On the contrary, the calculation of $\hat{\tilde{X}}_{\mathrm{TT}}^{0}$ and $\Delta \hat{\tilde{X}}_{\mathrm{TT}}$ using equations (38) and (39) involves the inversion of only the ${\mathrm{T}}$ block of $\hat{\mathbb{R}}^{\uparrow \downarrow}$. This procedure is numerically much more stable rather than the inversion of the whole matrix $\hat{\mathbb{R}}^{\uparrow \downarrow}$ in equation (36).

The idea of downfolding somewhat similar to ours was previously considered by Mryasov et al in order to eliminate the 5d states of the heavy Pt atoms and explain the unusual temperature dependence of the magnetic anisotropy energy in the ordered FePt alloy [77]. A simplified approach in the framework of the scheme $\hat{b}$, which did not take into account the effects of ligand–ligand interactions and $\hat{\cal I}_{\mathrm{L}}$, was also considered by Logemann et al [78].

The simplest model for CrCl₃ is the half-filled $t_\mathrm{2g}$ model. It contains only occupied $\uparrow$-spin $t_\mathrm{2g}$ band and unoccupied $\downarrow$-spin $t_\mathrm{2g}$ band. Since all basis functions of this model are associated with the Cr states, there are no additional complications coming from the ligand states. Even in LSDA, the exchange splitting $B_{\mathrm{Cr}}$ between the $\uparrow$- and $\downarrow$-spin $t_\mathrm{2g}$ bands is large compared to the bandwidth. Thus, the system should be in the strong-coupling limit. The corresponding parameters of exchange interactions, calculated by rotating $\hat{b}$, $\hat{m}$, or M are summarized in table 1.

Table 1. Isotropic exchange interactions in CrCl₃ (in meV): results of the $t_\mathrm{2g}$ model in LSDA. $\hat{b}$, $\hat{m}$, and M stand for the methods based on the infinitesimal rotations of, respectively, the xc field, the magnetization matrix, and the local spin moments. The notations of J_k are explained in figure 2.

Method	J1	J2	J3	J4	J5	J6
$\hat{b}$	−6.44	−0.50	−0.69	−0.17	−0.18	−0.63
$\hat{m}$	−6.32	−0.49	−0.68	−0.17	−0.18	−0.62
M	−6.30	−0.49	−0.69	−0.17	−0.18	−0.63

In this case we obtain very consistent description as all three methods provide very similar sets of parameters. This is not surprising: the approximate scheme $\hat{b}$ is justified in the strong coupling limit. Moreover, three $t_\mathrm{2g}$ levels are nearly degenerate. Therefore, the asphericity of $\hat{m}$ is small and corresponding parameters are practically identical to the ones obtained in the M scheme. However, all interactions are AFM, which is quite expected for the half filling [2], but totally contradicts to the experimental situation.

The next important question is whether the ferromagnetism of CrX₃ can be explained without the ligand states, in the model including both $t_\mathrm{2g}$ and $e_\mathrm{g}$ bands. In LSDA, the $\uparrow$- and $\downarrow$-spin bands are subjected to the exchange splitting $B_{\mathrm{Cr}}$. At the first sign, there is only a small addition in comparison with the $t_\mathrm{2g}$ model—the unoccupied $e_\mathrm{g}$ bands. However, it changes the story dramatically. First, the crystal-field splitting between $t_\mathrm{2g}$ and $e_\mathrm{g}$ bands is comparable with $B_{\mathrm{Cr}}$. Moreover, since only the transitions between the occupied $\uparrow$-bands and unoccupied $\downarrow$-bands contribute to $\hat{\mathcal{R}}^{\uparrow \downarrow}$, such exchange interactions do not know anything about the existence of the $\uparrow$-spin $e_\mathrm{g}$ band. Thus, although $B_{\mathrm{Cr}}$ is large (as in the $t_\mathrm{2g}$ model), the behavior of the exchange interactions is also controlled by other details of the electronic structure and the system is no longer in the strong coupling limit. Furthermore, the matrix $\hat{m}$ in the basis of all five 3d orbitals acquires additional degrees of freedom besides M, which is only the spherical part of $\hat{m}$.

All these tendencies are reflected in the behavior of interatomic exchange interactions (table 2). First, the inclusion of the $e_\mathrm{g}$ band into the model gives rise to the FM interactions. These interactions prevail in the nearest-neighbor (nn) bonds and considerably weaken the AFM interactions in other neighbors, in agreement with the experimental situation. Then, the schemes $\hat{b}$, $\hat{m}$, and M provide quite a different description. Compared to the methods $\hat{m}$ and M, the approximate method $\hat{b}$ underestimates the FM interaction J₁. Moreover, the method $\hat{m}$ (in comparison with M) has a tendency to overestimate the in-plane interactions J₁, J₃, and J₆, as expected from the analysis in section 2.7.

Table 2. Isotropic exchange interactions in CrCl₃ (in meV): results of the ${\mathrm{Cr}} \, 3d$ model in LSDA. $\hat{b}$, $\hat{m}$, and M stand for the methods based on the infinitesimal rotations of, respectively, the xc field, the magnetization matrix, and the local spin moments. The xc field was either computed from the on-site spin splitting of the TB Hamiltonian (H) or obtained from the sum rule (sr). $T_{\mathrm{S}}$ is corresponding spin transition temperature in RPA (in K). The method $\hat{m}$ yields the ferromagnetic ground state, while the methods $\hat{b}$ and M result in an incommensurate spin-spiral structure propagating both in and between the planes. The notations of J_k are explained in figure 2(b).

Method	J1	J2	J3	J4	J5	J6	$T_{\mathrm{S}}$
$\hat{b}$ (H)	$\phantom{1}2.39$	−0.27	−0.45	−0.03	−0.04	−0.50	8
$\hat{b}$ (sr)	$\phantom{1}2.27$	−0.16	−0.37	$\phantom{-}0.01$	0	−0.46	5
$\hat{m}$	11.94	−0.01	−0.69	$\phantom{-}0.08$	$\phantom{-}0.09$	−0.69	30
M	$\phantom{1}4.68$	−0.21	−0.40	0	−0.01	−0.48	5

Another important question is how to define the xc field in the TB model. The common practice is to use $\hat{b} = \hat{H}^{\uparrow}-\hat{H}^{\downarrow}$ and take only site-diagonal (local) part of this splitting. However, in LSDA, the off-diagonal elements of $\hat{H}$ can also depend on spin, giving rise to non-local contributions to the xc field [30]. Then, the use of the site-diagonal part alone will violate the sum rule because $\hat{\mathcal{R}}_{0}^{\uparrow \downarrow} \vec{b}$ with such $\vec{b}$ will no longer reproduce the ground-state magnetization $\vec{m}$ (see section 2.6). Another possibility is to reenforce the sum rule by defining the local xc field $\vec{b} = \hat{\mathcal{Q}}^{\uparrow \downarrow}_{0}\vec{m}$, which would reproduce the ground state magnetization $\vec{m}$. In the three-orbital (3o) $t_\mathrm{2g}$ model, this results only in minor change of the exchange interactions. However, starting from the five-orbital (5o) model, the difference becomes more significant (see table 2).

Similar analysis can be performed for the correlated model in the Hartree–Fock approximation (see appendix A). The non-interacting electron part of the model Hamiltonian in this case is evaluated in LDA. The corresponding electronic structure is shown in figures 2(c) and (d): since in CrCl₃ and CrI₃ the Cr 3d bands are well separated from the ligand ones, such model can be easily constructed for both compounds. The screened Coulomb interactions are evaluated within constrained RPA (cRPA) [93], as explained in [89]. The obtained (averaged) parameters of on-site Coulomb repulsion U, intra-atomic exchange interaction J, and nonsphericity 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 CrCl₃ (CrI₃). The screened U is not particularly large. Moreover, the screening is more efficient in CrI₃ due to proximity of the I 5d and Cr $t_\mathrm{2g}$ bands [89]. The corresponding densities of states (DOSs) are shown in figure 3. As expected [94, 95], the Coulomb U further increases the band gap in comparison with LSDA. The band gap is smaller in CrI₃ because U is smaller.

Zoom In Zoom Out Reset image size

Corresponding parameters of exchange interactions are summarized in tables 3 and 4 for CrCl₃ and CrI₃, respectively. The approximate method $\hat{b}$ systematically underestimates the FM interactions. For instance, all interactions in CrCl₃ are AFM, which clearly contradicts to the experimental situation. In CrI₃, only J₁ is weakly FM, which is not enough to stabilize the FM ground state [96]. The method M systematically improves the situation: the FM interactions clearly prevail and the FM ground state is stabilized both in CrCl₃ and CrI₃. ⁵ The corresponding Curie temperature, evaluated in RPA [97] (see appendix B for details), is also in reasonable agreement with the experimental data. The method $\hat{m}$ has a tendency to overestimate the interactions J₁, J₃, and J₆, making the FM ground state in CrI₃ unstable.

Table 3. Isotropic exchange interactions in CrCl₃ (in meV): results of the correlated Cr 3d model in the Hartree–Fock approximation. $\hat{b}$, $\hat{m}$, and M stand for the methods based on the infinitesimal rotations of, respectively, the xc field, the magnetization matrix, and the local spin moments. $T_{\mathrm{S}}$ is the corresponding spin transition temperature in RPA (in K). The methods $\hat{m}$ and M yield the ferromagnetic ground state, while the method $\hat{b}$ results in an incommensurate spin-spiral state. The notations of J_k are explained in figure 2(b).

Method	J1	J2	J3	J4	J5	J6	$T_{\mathrm{S}}$
$\hat{b}$	−1.69	−0.13	−0.28	−0.02	−0.03	−0.18	6
$\hat{m}$	$\phantom{-}6.66$	−0.06	−0.58	$\phantom{-}0.02$	$\phantom{-}0.02$	−0.35	11
M	$\phantom{-}3.76$	−0.02	−0.13	$\phantom{-}0.05$	$\phantom{-}0.05$	−0.14	27

Table 4. The same as table 3 but for CrI₃. The method M yields the ferromagnetic ground state, while the methods $\hat{b}$ and $\hat{m}$ result in an incommensurate spin-spiral state. (A) denotes the same parameters calculated in the antiferromagnetic state.

Method	J1	J2	J3	J4	J5	J6	$T_{\mathrm{S}}$
$\hat{b}$	0.76	−0.30	−0.39	0.04	0	−0.46	10
$\hat{m}$	5.87	−0.10	−0.21	0.33	0.16	−0.71	14
M	4.58	−0.08	−0.06	0.37	0.25	−0.37	51
M (A)	4.52	−0.07	−0.06	0.36	0.25	−0.41	49

Thus, the minimal model, which can capture the FM character of the exchange interactions in the chromium trihalides is the 5o model. Formally, the ferromagnetism can be obtained without invoking the ligand states. Nevertheless, it is crucial to consider the unoccupied $e_\mathrm{g}$ bands. Furthermore, it is crucially important to use the exact technique, formulated in terms of the inverse response function. The approximate method $\hat{b}$, which is linear in $\hat{\mathcal{R}}$, can lead to an incorrect magnetic ground state (see tables 3 and 4).

Generally, the interatomic exchange interactions, defined via infinitesimal rotations of spins, can depend on the magnetic state. In a number of cases, this dependence can be very strong, reflecting the dependence of the electronic structure on the magnetic arrangement [98, 99]. Considering how J_k in CrX₃ change depending on the method used for their calculations (for instance, M versus $\hat{b}$), it is reasonable to expect that the system is no longer in the strong-coupling limit and, therefore, these parameters can also depend on the magnetic state. However, this appears to be not the case. For instance, the exchange parameters calculated in the AFM state of CrI₃ are practically identical to the ones in the FM state (table 4). This observation is very important and means, for example, that the parameters derived in the ordered FM state can be also used to evaluate the spin transition temperature, $T_{\mathrm{S}}$. The main reason why the $\hat{b}$ and M methods yield different J_k is related to the fact that the crystal-field splitting, 10Dq, is comparable with the intra-atomic exchange splitting. However, this 10Dq does not depend on the magnetic state, and therefore does not contribute to the magnetic-state dependence of J_k .

Now we turn to the analysis of the most general model, which explicitly includes the contributions of the Cr 3d as well we the ligand Cl 3p or I 5p states.

3.3.1. Cr 3d and X p bands in LSDA.

We start with the analysis of general TB model, including both Cr 3d and X p bands, in LSDA. First, we note that the Cr 3d and ligand p states are antipolarized. In this case $M_{\mathrm{Cr}}$ exceeds the nominal value of $3\mu_{\mathrm{B}}$ (where $\mu_{\mathrm{B}}$ denotes the Bohr magneton). Nevertheless, it is compensated by the negative M_X on the ligand atoms, so that the total moment, $M_{\mathrm{tot}} = M_{\mathrm{Cr}} + 3 M_{X}$, is integer and equal to $3\mu_{\mathrm{B}}$. This antipolarization is a joint effect of the spin splitting on the Cr atoms and the hybridization between the occupied ligand states and unoccupied Cr $e_\mathrm{g}$ states. Since the $\uparrow$-spin Cr $e_\mathrm{g}$ states are closer to the ligand band, the hybridization is stronger, resulting in the stronger admixture of the $\uparrow$-spin Cr $e_\mathrm{g}$ states into the occupied ligand band (and transfer of the $\uparrow$-spin ligand states into the unoccupied Cr $e_\mathrm{g}$ band). Therefore, when the ligand band is added to the model, we have $M_{\mathrm{Cr}}\gt3\mu_{\mathrm{B}}$ but $M_{X}\lt0$, as schematically illustrated in figure 4. For instance, the LSDA yields $M_{\mathrm{Cr}} = 3.14$ (3.37) $\mu_{\mathrm{B}}$ and $M_{X} = -0.05$ (−0.12) $\mu_{\mathrm{B}}$ for CrCl₃ (CrI₃). As expected, the effect is stronger in CrI₃ where the Cr $e_\mathrm{g}$ states are closer to the I 5p band (see figure 2(d)). Moreover, the I 5p states are more extended in comparison with the Cl 3p ones, resulting in stronger hybridization.

Zoom In Zoom Out Reset image size

The new aspect of calculations of the exchange interactions is that now they can also depend on the strength of the Stoner coupling $\mathcal{I}_{X}$, which contributes to the second term of equation (37). Therefore, the first problem we have to solve is how to properly define $\mathcal{I}$. In fact, there are several possible definitions:

• (I)  
  $\mathcal{I}_{\mu} = -{\mathrm{Tr}}_{L} \{\hat{b}_{\mu} \hat{m}_{\mu} \}/(M_{\mu})^{2}$, where $\hat{b}_{\mu}$ is associated with the intra-atomic spin splitting of $\hat{H}^{\sigma}$ and $\hat{m}_{\mu}$ is the ground-state magnetization;
• (II)  
  $\mathcal{I}_{\mu} = -B_{\mu}/M_{\mu}$, where $B_{\mu} = \frac{1}{n} {\mathrm{Tr}}_{L} \{\hat{b}_{\mu} \}$ and $M_{\mu} = {\mathrm{Tr}}_{L} \{\hat{m}_{\mu} \}$, which is nothing but the spherically averaged version of (I);
• (III)  
  the same as (I), but taking $\hat{b}_{\mu}$ from the sum rule (25);
• (IV)  
  the same as (II), but taking B_µ from the sum rule $\vec{M} = \hat{\mathbb{R}}^{\uparrow \downarrow} \vec{B}$;
• (V)  
  the same as (I), but taking $\hat{m}_{\mu}$ and M_µ from the sum rule (and defining $\hat{b}_{\mu}$ as the intra-atomic spin splitting of $\hat{H}^{\sigma}$);
• (VI)  
  the same as (II), but taking M_µ from the sum rule.

The results are summarized in table 5. One can see that $\mathcal{I}_{\mathrm{Cr}}$ only weakly depends on the definition (thought it is different in CrCl₃ and CrI₃). However, we do not need this parameter in our calculations. On the other hand, $\mathcal{I}_{X}$ is very sensitive to the definition [76]. Apparently, one can discard the definitions V and VI as they yield incorrect $M_{\mathrm{tot}}$ violating the fundamental property $M_{\mathrm{tot}} = 3\mu_{\mathrm{B}}$. ⁶ Furthermore, it makes sense to enforce the sum rule by using the definitions III and IV instead of I and II. Finally, the definitions II and IV are more appropriate for the spherically averaged method M, while the definitions I and III should be used in the combination with the matrix form of the methods $\hat{m}$ and $\hat{b}$.

Table 5. The effective Stoner parameters in CrCl₃ and CrI₃ (in eV), obtained using six possible definitions (as explained in the text).

Definition	CrCl3		CrI3
	$\mathcal{I}_{\mathrm{Cr}}$	$\mathcal{I}_{\mathrm{Cl}}$	$\mathcal{I}_{\mathrm{Cr}}$	$\mathcal{I}_{\mathrm{I}}$
I	0.82	−3.08	0.76	−0.55
II	0.81	−0.75	0.75	−0.13
III	0.83	−0.28	0.77	$0.31$
IV	0.83	−3.32	0.77	−0.46
V	0.83	−2.50	0.77	−0.51
VI	0.83	−0.64	0.77	−0.13

While $\mathcal{I}_{\mathrm{Cr}}$ in CrX₃ is close to atomic values (${\sim}0.7$ eV) and can be interpreted as the intra-atomic exchange integral responsible for Hund's first rule [100–102], $\mathcal{I}_{X}$ is definitely not. In most of the cases $\mathcal{I}_{X}$ is negative (except scheme III for CrI₃). This is another manifestation of the fact that the moments M_X are merely induced by the hybridization with the Cr 3d states, which acts against B_X . According to equation (39), the negative $\mathcal{I}_{X}$ will tend to decrease the FM coupling, contrary to rather common believe that the ferromagnetism in CrX₃ is driven by the Hund's rule coupling on the ligand sites, as suggested by the phenomenological GKA rules [36–39].

The parameters of exchange interactions are summarized in tables 6 and 7, for CrCl₃ and CrI₃, respectively. We note the following. In comparison with the 5o model (section 3.2), the parameters J_k becomes mostly FM, except AFM J₆. The exchange interactions are sensitive to $\mathcal{I}_{X}$. This dependence is illustrated only for the method M, but very similar behavior is observed also for the methods $\hat{b}$ and $\hat{m}$. The use of the M method in the combination with $\mathcal{I}_{X}$ obtained in the spherically averaged scheme IV substantially improves the agreement with the experimental data for $T_{\mathrm{C}}$ (but not for the spin-wave dispersion, which will be considered in section 3.6). The exchange parameters obtained in the method $\hat{m}$ are unrealistically large [61] and not shown here. Furthermore, in the earlier work [61], the approximate method $\hat{b}$ was concluded to provide systematically worse description, especially when considers contributions of the ligand states. However, the situation crucially depends to two factors: (i) the proper choice for $\mathcal{I}_{X}$; (ii) the proper definition of the xc field $\hat{b}$. For instance, abilities of this method can be substantially improved by enforcing the sum rule in the definition of $\hat{b}$ and $\mathcal{I}_{X}$, and taking into account the aspherical contributions to $\mathcal{I}_{X}$, following the definition III.

Table 6. Isotropic exchange interactions in CrCl₃ (in meV): results of the ${\mathrm{Cr}} \, 3d+{\mathrm{Cl}} \, 3p$ model in LSDA. $\hat{b}$ and M stand for the methods based on the infinitesimal rotations of, respectively, the xc field and the local spin moments. The xc field was either associated with the on-site spin splitting of the TB Hamiltonian (H) or obtained from the sum rule (sr). $\mathcal{I}_{\mathrm{Cl}}$ is the value of the effective Stoner parameter used in the calculations (the results of the methods III and IV from table 5, in eV). In the ${\mathrm{L}}0$ method, the xc field on the ligand sites is set to be zero, while that on the Cr sites is set to reproduce the total magnetization $M_{\mathrm{tot}} = 3\mu_{\mathrm{B}}$ (further details are given in the text). $T_{\mathrm{C}}$ is corresponding Curie temperature in RPA (in K). The notations of J_k are explained in figure 2(b).

Method	$\mathcal{I}_{\mathrm{Cl}}$	J1	J2	J3	J4	J5	J6	$T_{\mathrm{C}}$
$\hat{b}$ (H)	−0.28	−0.51	0.26	0.34	0.16	0.15	−0.31	21
$\hat{b}$ (sr)	−0.28	$\phantom{-}4.21$	0.28	0.26	0.13	0.13	−0.34	49
M	−0.28	$\phantom{-}4.36$	0.28	0.25	0.16	0.16	−0.32	52
M	−3.32	$\phantom{-}0.64$	0.25	0.53	0.08	0.10	−0.42	20
${\mathrm{L}}0$	0	$\phantom{-}4.73$	0.29	0.25	0.17	0.17	−0.32	56

Table 7. Isotropic exchange interactions in CrI₃ (in meV): results of the ${\mathrm{Cr}} \, 3d+{\mathrm{I}} \, 5p$ model in LSDA (see table 6 for the notations).

Method	$\mathcal{I}_{\mathrm{I}}$	J1	J2	J3	J4	J5	J6	$T_{\mathrm{C}}$
$\hat{b}$ (H)	$\phantom{-}0.31$	3.93	0.77	0.75	0.87	0.58	−0.44	$\phantom{1}97$
$\hat{b}$ (sr)	$\phantom{-}0.31$	3.40	0.76	0.71	0.85	0.56	−0.46	$\phantom{1}88$
M	$\phantom{-}0.31$	4.29	0.87	0.89	1.03	0.67	−0.41	113
M	−0.46	1.32	0.75	0.77	0.82	0.54	−0.51	$\phantom{1}62$
${\mathrm{L}}0$	0	3.21	0.84	0.86	0.96	0.63	−0.46	$\phantom{1}96$

Recently, Ke and Katsnelson [103] reported the exchange parameters for CrI₃, which were also derived using the inverse response function (i.e. similar to our method M). In SDFT, they have found (using our notations, in meV): $J_{1} = 6.58$, $J_{2} = 0.38$, $J_{3} = 1.14$, $J_{4} = 0.62$, $J_{5} = 0.94$, and $J_{6} = -0.14$. These parameters are in reasonable agreement with our data, perhaps except J₁, which is systematically smaller in our case. Nevertheless, our experience shows that in order to obtain reasonable parameters J_k in CrI₃ it is absolutely crucial to consider the contributions of the ligand states (using the downfolding technique described in section 2.10). The bare exchange interactions between the Cr 3d states are strongly AFM and fails to reproduce the correct FM ground state [61]. Unfortunately, it is not clear how this problem was tackled by Ke and Katsnelson [103]. Furthermore, our parameters are sensitive to the choice of $\mathcal{I}_{\mathrm{I}}$, as seen in table 7.

The sensitivity of the exchange interactions to the parameter $\mathcal{I}_{X}$ can be regarded as the weak point of the downfolding technique. Nevertheless, one can propose an alternative option, which can be viewed as an extension of the method M for the FM insulators and half-metals. In the future, we will call it the ${\mathrm{L}}0$ scheme. The crucial observation is that M_X is induced by the hybridization between the Cr 3d and ligand p states (which is totally in line with the general idea of the downfolding technique, considered in section 2.10). Then, it is reasonable to enforce $B_{X} = 0$ and find the remaining parameter of the xc field, $B_{\mathrm{Cr}}$, from the equation $\vec{M} = \hat{\mathbb{R}}^{\uparrow \downarrow} \vec{B}$. This $B_{\mathrm{Cr}}$ induces the magnetic moments on the Cr sites as well as the ligand sites. Then, the parameter $B_{\mathrm{Cr}}$ can be chosen to reproduce $M_{\mathrm{tot}} = 3 \mu_{\mathrm{B}}$. For CrCl₃ (CrI₃) in LSDA, such procedure yields $M_{\mathrm{Cr}} = 3.18$ (3.40) $\mu_{\mathrm{B}}$ and $M_{X} = -0.06$ (−0.13) $\mu_{\mathrm{B}}$, which are pretty close to the values of spin magnetic moments in the LSDA ground state. The good aspect of this approximation is that $\mathcal{I}_{X} = 0$ and the exchange interactions are solely determined by the 1st term of equation (37) (but with the redefined $M_{\mathrm{Cr}}$ for $B_{X} = 0$). The corresponding values of J_k and $T_{\mathrm{C}}$ are also listed in tables 6 and 7. In comparison with the M scheme, these parameters somewhat overestimate $T_{\mathrm{C}}$ but otherwise fall within the range of typical estimates for J_k . We will continue to use this ${\mathrm{L}}0$ scheme in the future. It is especially good for semiquantitative estimates, which would illustrate the basic trend in the behavior of interatomic exchange interactions. The 2nd term in equation (37) also vanishes if $\hat{\unicode{x211A}}^{\uparrow \downarrow}_{\mathrm{LL}} = 0$, meaning that there is no ligand–ligand interactions. In this sense, there is certain similarity with the downfolding method proposed by Logemann et al [78], but reformulated for the scheme M.

3.3.2. Ligand states in correlated Cr 3d band.

In this section, we try to extend the correlated 5o model, by expanding its Wannier basis W into the pseudo-atomic Wannier basis, A, of the more general model, containing Cr 3d as well as the ligand p bands. Namely, we still consider only ten Cr 3d bands (for each spin), but transform $| C \rangle$ in equation (12) from the basis W to the basis A: $| C^\mathrm{W} \rangle \to | C^\mathrm{A} \rangle = \hat{T}| C^\mathrm{W} \rangle$, where $\hat{T}$ is the transformation matrix. Our intension here is to consider explicitly the X p states and all the contributions of these states to the exchange interactions in the 5o model. The magnetic moments in the A basis are $M_{\mathrm{Cr}} = 2.63$ (2.56) $\mu_{\mathrm{B}}$ and $M_{X} = 0.12$ (0.15) $\mu_{\mathrm{B}}$ for CrCl₃ (CrI₃). In the 5o model, the only occupied $\uparrow$-spin $t_\mathrm{2g}$ band provides three electrons, which are distributed between one Cr and three ligand atoms. Therefore, the Cr and ligand moments in the $t_\mathrm{2g}$ band are ferromagnetically coupled (see also figure 4), while in order to obtain the opposite polarization of these states, it is essential to consider a more general model, which would explicitly include the ligand p band.

Without the ligand bands, the inversion of the matrix $\hat{\mathbb{R}}^{\uparrow \downarrow}$ in the A basis becomes unstable as it contains small matrix elements associated with the ligand sites. While for calculations of the exchange interactions for the given $\mathcal{I}_{X}$ such inversion can be largely avoided using equations (38) and (39), the inverse matrix is still needed to evaluate $\mathcal{I}_{X}$ using the sum rule. Nevertheless, one can still try to estimate J_k using the method ${\mathrm{L}}0$, where it is sufficient to invert only the subblock of $\hat{\mathbb{R}}^{\uparrow \downarrow}$ in the subspace of the Cr sites. It yields the magnetic moments: $M_{\mathrm{Cr}} = 2.64$ (2.59) $\mu_{\mathrm{B}}$ and $M_{X} = 0.12$ (0.14) $\mu_{\mathrm{B}}$ for CrCl₃ (CrI₃), which are close to the ones reported above.

The corresponding exchange interactions and $T_{\mathrm{C}}$ for the 5o model reformulated in the pseudo-atomic basis, are clearly overestimated (see table 8, the rows denoted 5o). However, this is quite in line with general understanding. The interactions mediated by the tails of the Wannier functions spreading to the ligand sites and transferring there the FM magnetization can be viewed as an analog of the direct Heisenberg exchange [1]. In order to evaluate these direct exchange contributions numerically, one typically performs the six-dimensional integration in the real space [40, 104, 105]. Nevertheless, the downfolding method suggests how these contributions can be naturally evaluated within SDFT. The bare direct exchange integrals are typically large and need to be additionally scaled, in the spirit of cRPA, to account for the screening caused by other bands [40, 105]. The same situation is here: the attempt to include the ligand states into 5o model only worsen the description of exchange interactions. The discrepancy can be resolved by considering more general model, which would explicitly include the contributions of the ligand p band.

Table 8. Exchange interactions in CrCl₃ and CrI₃ (in meV) obtained in the scheme ${\mathrm{L}}0$ for the 5o model with the ligand states (5o) and after merging this model with the ligand p bands (5o + p). $T_{\mathrm{C}}$ is the corresponding Curie temperature in RPA (in K). The notations of J_k are explained in figure 2(b).

System	J1	J2	J3	J4	J5	J6	$T_{\mathrm{C}}$
CrCl3 (5o)	14.27	0.85	1.65	0.42	0.45	$\phantom{-}0.21$	242
CrCl3 (5o + p)	10.81	0.28	0.35	0.18	0.15	−0.27	102
CrI3 (5o)	14.03	2.37	3.63	1.89	1.59	$\phantom{-}0.76$	484
CrI3 (5o + p)	$6.87$	0.72	0.91	1.14	0.63	−0.31	119

The next important step is to merge the correlated Cr 3d bands with the ligand bands in the framework of the Cr $3d + X$ p model. The correlated 5o model is formulated in the Wannier basis for the Cr 3d bands in LDA. After solving this model in the Hartree–Fock approximation, we want to replace the original Cr 3d bands in the all-electron LDA band structure by these correlated bands in the pseudo-atomic A basis. Since such basis states of the 5o model and the ligand p bands are orthogonal to other states, such merging is not unique and an arbitrary energy shift of the Cr 3d and X p bands relative to each other will not change the magnetization (but will change the exchange interactions). Similar problem arises in the LDA $+\,U$ method [94], where the relative position of the transition-metal 3d and ligand p states is controlled by the empirical double-counting term [57, 95]. Therefore, we have to make some empirical conjecture about this merging and we choose it from the 'charge neutrality' condition, by requesting the Fermi energy of the correlated model to coincide with the one in LDA. If the system is insulating, Fermi energy is chosen in the middle of the gap. The electronic structure after such merging is shown in figure 5.

Zoom In Zoom Out Reset image size

We would like to emphasize that this electronic structure is rather artificial and would require different sets of the xc fields for the Cr 3d and X p bands (zero for the X p band, but finite for the Cr 3d one), which is hardly useful for our purposes. Therefore, in order to evaluate the exchange interactions, we again employ the ${\mathrm{L}}0$ technique. It yields the following magnetic moments: $M_{\mathrm{Cr}} = 3.35$ (3.49) $\mu_{\mathrm{B}}$ and $M_{X} = -0.12$ (−0.16) $\mu_{\mathrm{B}}$ for CrCl₃ (CrI₃). Thus, after adding the X p band, the ${\mathrm{L}}0$ method nicely capture the antipolarization of the Cr 3d and X p states, which is again stronger in CrI₃. However, it should be understood that this effect is 'mimicked' by the specific choice of the xc fields, while from the viewpoint of electronic structure itself, the X p band remains unpolarized and the Cr 3d band is polarized as in section 3.3.2. In the other words, in the ${\mathrm{L}}0$ method, the imperfection of the electronic structure is corrected by the physical choice of the xc fields (though the electronic structure itself was obtained not for these fields).

The corresponding exchange parameters are listed in table 8 (and denoted as 5o $+\,p$). We note a systematic improvement in comparison with the 5o model with the ligand states (denoted as 5o): all exchange interactions become smaller so as $T_{\mathrm{C}}$. However, this improvement is only partial as these parameters are still too large and clearly overestimate the tendency toward the ferromagnetism. This demonstrates the complexity of the problem: the ${\mathrm{L}}0$ method is only as the starting point of the self-consistent solution, which should take into account the feedback of correlated Cr 3d bands onto the ligand bands via the xc field and results in the magnetic polarization of the ligand bands. However, such field will inevitably mix up the Cr 3d and X p bands. After that, one has to redefine the Wannier basis for the correlated model and recalculate the model parameters. An alternative solution is LDA $+\,U$ [94, 95], which considers both transition-metal and ligand states, and takes into account the polarization of the ligand bands. Today there are many applications of this method for calculating the interatomic exchange interactions in various transition-metal oxides and other strongly correlated systems [26–28]. However, it should be understood that the LDA $+\,U$ approach is supplemented with additional approximations of purely empirical character, such as the form of the double counting as well as the choice of the Coulomb interaction parameters and the correlated subspace for the solution of the Hubbard-type model [57].

In this section, we briefly discuss the behavior of DM interactions in CrI₃, which are driven by the strong SO coupling of the heavy I atoms (so as other magnetic interactions of the relativistic origin [86]). The DM interactions in CrCl₃ are considerably weaker [49, 53] and not significant [83]. In addition to d^z , two other components of the DM vector, d^x and d^y , can be computed by rotating the coordinate frame and applying the same equation (28) (or similar equations for the $\hat{b}$ or $\hat{m}$ rotations). The symmetry of CrI₃ is such that two Cr 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 vanish [6, 7]. On the other hand, the interactions within the sublattices can exist and for each vector R connecting two Cr sites, the interactions in the sublattices I and II are related as: $\boldsymbol{d}^\mathrm{II}(\boldsymbol{R}) = - \boldsymbol{d}^\mathrm{I}(\boldsymbol{R})$. The strongest interaction is d ₃, which occurs in the 2nd coordination sphere of the honeycomb plane (together with the isotropic interaction J₃). In total, there are six bonds n connecting the central atom with the atoms in the 2nd coordination sphere. The corresponding DM vectors are $\boldsymbol{d}_{3,n} = ( d^{xy}\text{cos} ( \phi_{n}+\psi ), d^{xy}\text{sin} ( \phi_{n}+ \psi ), (-1)^{n}d^{z})$, where $\phi_{n} = (2n+1)\frac{\pi}{6}$ is the azimuthal angle specifying the direction of the bond n in the xy plane (relative to the bond along a in figure 2(a), corresponding to n = 1) and ψ specifies the direction of the DM vector in the plane relative to this bond. For the $R\overline{3}$ symmetry, the parameters d^xy , d^z , and ψ do not depend on n. They are listed in table 9.

Table 9. Parameters of Dzyaloshinskii–Moriya interactions in CrI₃ (d^xy and d^z are in meV and ψ is in degrees) in correlated five-orbital model (5o) and LSDA for the all-electron ${\mathrm{Cr}} \, 3d+{\mathrm{I}} \, 5p$ model.

Method	Model	$|d^{xy}|$	ψ	$|d^{z}|$
$\hat{b}$	5o	0.01	39	0.22
M	5o	0.02	19	0.20
$\hat{b}$	LSDA	0.08	−2	0.25
M	LSDA	0.08	0	0.28

We note that the schemes $\hat{b}$ and M provide very consistent description for the DM interactions. The small discrepancy for ψ in the correlated 5o model is probably due to the fact that the angle ψ is ill-defined when d^xy is small. There is also surprisingly good agreement between results of the correlated 5o model and the Cr $3d + \mathrm{I}\, 5p$ model in LSDA. However, such agreement is probably fortuitous because the models are very different ⁷ . The DM interactions in CrCl₃ are considerably weaker ($|d^{z}| \sim 0.02$ meV [53]). However, this is to be expected and the order of magnitude difference of the DM interactions in CrCl₃ and CrI₃ well correlates with the strength of the SO coupling on the ligand sites, which differs by the same order of magnitude: $\xi_{\mathrm{Cl}} = 95$ meV versus $\xi_{\mathrm{I}} = 881$ meV.

The experimental parameters of exchange interactions, derived from the inelastic neutron scattering data [81, 83, 84], are summarized in table 10. ⁸ One of the interesting features of the experimental spin-wave dispersion in CrI₃ is a ${\sim}4$ meV gap at the Dirac (K) point, indicating at the existence of large DM interaction d^z between 2nd neighbors in the honeycomb plane [81]. Quite expectably, no such feature was observed in CrCl₃ [83, 84], where this DM interaction is small. For the isotropic interactions, the agreement between theoretical and experimental data depends on the model and approximations employed for treating the Coulomb correlations and the ligand states, which we will discuss below. The theoretical DM interaction in CrI₃ appears to be underestimated by factor two, both in the correlated 5o model and LSDA for the ${\mathrm{Cr}} \, 3d+{\mathrm{I}} \, 5p$ model.

Table 10. Experimental parameters of exchange interactions in CrCl₃ and CrI₃.

Compound	Reference	J1	J2	J3	J6	dz
CrCl3	[83]	2.14	−0.02	0.05	−0.11	0.03
CrCl3	[84]	2.12		0.08	−0.16
CrI3	[81]	4.52	$\phantom{-}1.33$	0.36	−0.18	0.50

The theoretical spin-wave dispersion in CrCl₃ and CrI₃, calculated for some representative sets of parameters, is plotted in figures 6 and 7, respectively, in comparison with the experimental data. First, let us consider results of correlated 5o model (shown in the panels a and b of these figures), where there are no additional complications caused by the ligand states. Within this 5o model, the approximate scheme $\hat{b}$ severely underestimates the FM interactions, making the FM ground state unstable both in CrCl₃ and CrI₃. The situation is dramatically improved when we use the 'exact' scheme M: it strengthens the FM interactions and stabilizes the FM ground state in both the compounds. In CrCl₃, the dispersion is overestimated by factor two, while in CrI₃ we obtain an overall fair agreement with the experimental data, except for the band splitting in the point K, which is underestimated by factor two. Nevertheless, the situation changes significantly when we explicitly consider the contributions of the ligand states in the framework of the all-electron ${\mathrm{Cr}} \, 3d+X \, p$ model in LSDA (panels c and d). Now, the FM state becomes stable already in the scheme $\hat{b}$. The spin-wave dispersion in CrCl₃ remains overestimated by factor two. The agreement with the experimental data is better for CrI₃, again except for the band splitting in the point K. When we turn to the M scheme, which is expected to be more accurate, the theoretical spin-wave dispersion in CrCl₃ is substantially reduced, so that the lowest experimental band is reproduced pretty well. However, the behavior of theoretical upper bands is not satisfactory. This is related to the fact that the theoretical $J_{1} = 0.64$ meV (table 6) is underestimated by factor three. A similar situation occurs in CrI₃, where the theoretical $J_{1} = 1.32$ meV (table 7) is underestimated by the same amount, which worsens the agreement with the experimental data. The problem is partially related to our choice of $\mathcal{I}_{X}$, which appears to be more negative for the scheme M due to the additional spherical approximation, as discussed in section 3.3.1. For instance, if one uses the same $\mathcal{I}_{X}$ as in the scheme $\hat{b}$ (or the ${\mathrm{L}}0$ scheme, where we do not need any $\mathcal{I}_{X}$), one could get a better agreement with the experimental data for the upper bands (but not for the dispersion in the entire energy region).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Regarding the value of d^z and the band gap, first we note that there is an alternative mechanism of opening the band gap, related to long-range isotropic exchange interactions beyond the 6th coordination sphere [103]. These interactions were taken into account also in our calculations. ⁹ However, their effect is not particularly strong. Therefore, the band gap is mainly controlled by d^z and the discrepancy with the experimental data is caused by the underestimation of this interaction in the theoretical calculations, as was also reported by Kvashnin et al [49] and Olsen [106]. In this respect, we note that by explicitly considering the contributions of the ligand states in the correlated 5o model (section 3.3.2) and using the ${\mathrm{L}}0$ scheme for the evaluation of the exchange parameters, we can easily get $|d^{z}|$ as large as 1.10 meV, which exceeds the value obtained in the regular correlated 5o model by factor five and the experimental value by factor two [81]. Such enhancement is caused by the additional contribution of the basis functions, which explicitly takes into account the effect of the heavy I atoms. Nevertheless, when we try to combine this effect with another contribution stemming from the I 5p band itself (employing for these purposes the merging of the Cr 3d and I 5p bands, as discussed in section 3.4), again in the framework of the ${\mathrm{L}}0$ scheme, the DM interaction is reduced till $|d^{z}| = 0.16$ meV. Such strong reduction is apparently caused by 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 antipolarized, the fact of the cancellation itself is not surprising. However, such analysis clearly demonstrates the fragility of the situation, where the value of d^z appears to depend on the delicate balance of two large contributions arising from the Cr 3d and I 5p bands. The merging of these two bands considered in section 3.4 was probably too crude to explain the experimental situation.

To summarize this section, we would like to stress again the main points:

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

The half-metallic ferromagnetism is not very common in stoichiometric transition-metal oxides. Nevertheless, there are exceptions and CrO₂ is one of them, which is widely considered in various applications related to spintronics and magnetic recording [121]. One of the limitations of CrO₂ from this practical point of view is the relatively low $T_{\mathrm{C}} \sim 390$ K [121].

CrO₂ crystallizes in the rutile structure (the space group $P4_{2}/mnm$, figure 8(a)) [122]. The exchange interactions remain sizable up to at least the 8th coordination sphere. Moreover, since the rutile structure is nonsymmorphic, there are two types of interactions, J₇ and J₈, which are denoted by superscripts $\gt$ and $\lt$, as explained in figure 8(b). The calculations are performed for the experimental parameters of the crystal structure [122] on the mesh of the $10\times10\times16$ points both for k and q .

Zoom In Zoom Out Reset image size

The LDA band structure is featured by three separated bands: ${\mathrm{O}} \, 2p$ in the occupied part, ${\mathrm{Cr}} \, t_\mathrm{2g}$ near the Fermi level, and ${\mathrm{Cr}} \, e_\mathrm{g}$ in the unoccupied part. Therefore, one can consider two types of correlated models: the 3o model for the Cr $t_\mathrm{2g}$ bands and 5o model both for Cr $t_\mathrm{2g}$ and Cr $e_\mathrm{g}$ bands. In the former case, the on-site Coulomb and exchange interactions can be described in terms of two Kanamori parameters [123]: the intraorbital Coulomb interaction $\mathcal{U}$ and the exchange interaction $\mathcal{J}$, which can be evaluated within cRPA as 2.84 eV and 0.70 eV, respectively [120]. The third parameter of interorbital Coulomb interaction can be obtained from these two as $\mathcal{U}^{\prime} = \mathcal{U} - 2\mathcal{J}$ [123]. The parameters of 5o model are specified by U = 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 ${\mathrm{Cr}} \, 3d+{\mathrm{O}} \, 2p$ model in LSDA. The half-metallic character of the electronic structure is evident from the LSDA DOSs (figure 8(c)). The Coulomb interactions additionally split the $\uparrow$-spin $t_\mathrm{2g}$ band, placing the Fermi energy inside a pseudogap, ¹⁰ and shift the $\downarrow$-spin states to the higher energy region, away from the Fermi level.

Zoom In Zoom Out Reset image size

The parameters of exchange interactions are summarized in table 11. First, there is a large difference in the parameters obtained in the schemes $\hat{b}$ and M, especially in the Hartree–Fock approximations for correlated 3o and 5o models, where the more accurate method M strengthens the FM interactions practically in all the bonds. The situation is more modest in LSDA for the all-electron ${\mathrm{Cr}} \, 3d+{\mathrm{O}} \, 2p$ model: when going from the method $\hat{b}$ to the method M, the nearest FM interactions J₁ and J₂ increase only partially, while other longer-range interactions tend to decrease and become more AFM.

Table 11. Parameters of exchange interactions in CrO₂ (in meV), as obtained in correlated three- and five-orbital models (respectively, 3o and 5o), and in LSDA for the all-electron ${\mathrm{Cr}} \, 3d+{\mathrm{O}} \, 2p$ model using the schemes $\hat{b}$ and M (corresponding to the infinitesimal rotations of the xc fields and local spin moments, respectively). In the L0 scheme, xc field was redefined to enforce $\mathcal{I}_{\mathrm{O}} = 0$. The notations of exchange parameters are explained in figure 8(b). $D^{xx} = D^{yy}$ and D^zz are nonvanishing elements of the spin-stiffness tensor (in meV·Å²). $T_{\mathrm{C}}$ is the Curie temperature in RPA (in K).

	3o		5o		LSDA
	$\hat{b}$	M	$\hat{b}$	M	$\hat{b}$	M	L0
J1	9.18	15.76	5.24	17.47	30.87	36.86	37.63
J2	10.94	23.15	11.47	31.17	21.39	24.80	25.32
J3	1.12	1.55	0.33	0.44	3.04	2.54	2.57
J4	0.84	1.51	0.71	1.54	1.49	1.22	1.20
J5	−0.34	−0.38	−0.12	0.04	−0.79	−1.72	−1.83
J6	−1.92	−1.78	−1.75	−1.73	−3.58	−4.96	−5.14
$J_{7}^{\gt}$	−3.41	−3.36	−2.50	−2.79	−6.25	−8.39	−8.61
$J_{7}^{\lt}$	−1.09	−1.77	−0.38	−0.67	−2.09	−3.00	−3.11
$J_{8}^{\gt}$	−0.02	0.19	0.10	0.62	−1.50	−1.98	−2.08
$J_{8}^{\lt}$	−0.52	−0.64	−0.39	−0.30	−0.53	−0.66	−0.73
Dxx	37	258	113	532	123	113	109
Dzz	11	168	39	362	161	159	152
$T_{\mathrm{C}}$	310	1026	430	1563	868	880	875

In LSDA, the Cr and O atoms are antipolarized due to the joint effect of the hybridization and the intraatomic spin splitting, similar to CrX₃. The corresponding spin magnetic moments are $M_{\mathrm{Cr}} = 2.15\mu_{\mathrm{B}}$ and $M_{\mathrm{O}} = -0.08\mu_{\mathrm{B}}$. The effective Stoner parameters were calculated using the definitions III and IV, which should be used in the combination with the methods $\hat{b}$ and M, respectively (see section 3.3.1). As expected, $\mathcal{I}_{\mathrm{Cr}} = 0.87$ eV practically does not depend on the definition. On the other hand, the definitions III and IV yield different values of $\mathcal{I}_{\mathrm{O}}$: −0.78 eV and −0.48 eV, respectively. Nevertheless, the parameters are pretty small and the exchange interactions in CrO₂ appears to be less sensitive to the ligand states (at least, in comparison with CrX₃) [61]. For instance, the scheme L0, where we enforce $\mathcal{I}_{\mathrm{O}} = 0$, provides basically the same set of parameters as the M scheme with $\mathcal{I}_{\mathrm{O}} = -0.48$ eV (see table 11).

The spin-wave dispersions calculated with different sets of parameters for the 3o and 5o models is plotted in figure 10 (see appendix B for details). The nonvanishing xx and zz elements of the spin-wave stiffness tensor $\hat{D} = [D^{\alpha \beta}]$,

Zoom In Zoom Out Reset image size

$$\begin{align*} \omega_\mathrm{A}\left(\boldsymbol{q}\right) \approx \sum_{\alpha, \beta} D^{\alpha \beta} q_{\alpha} q_{\beta}, \end{align*}$$

describing the dispersion of the acoustic (A) spin-wave branch near the Γ point, are summarized in table 11. The experimental estimates for the averaged spin-wave stiffness in CrO₂ vary from 70 to 150 meV·Å² [125, 126]. The Curie temperature is about 390 K [121]. Then, the LSDA values of $T_{\mathrm{C}}$ seem to be overestimated, though the situation is rather subtle. On the one hand, the theoretical values of the spin-wave stiffness are within the experimental scatter. On the other hand, the exchange parameters derived for ground state are not supposed to reproduce $T_{\mathrm{C}}$ for this metallic system and more sophisticated methods, considering the temperature dependence of the exchange interactions, may be needed [22, 127, 128]. Anyways, in the light of well-known limitations of LSDA for the transition-metal oxides [94], such moderate disagreement would not be surprising.

Even more interesting situation is realized in correlated 3o and 5o models. At first glance, the approximate scheme $\hat{b}$ provides a much better description for the spin-wave stiffness and $T_{\mathrm{C}}$ in comparison with the scheme M, which is supposed to be more accurate. However, this is true only in the Hartree–Fock approximation, which has serious limitations for the metallic systems. If one goes beyond the Hartree–Fock approximation, the situation can change dramatically. Particularly, as expected for the half-metallic systems [119], the dynamic correlations lead to a strong redistribution of the electronic states. These effects can be treated in the framework of dynamical mean-field theory (DMFT) [129]. The latter can be regarded as an extension of SDFT in which the KS potential is replaced by the frequency-dependent self-energy $\Sigma^{\uparrow,\downarrow}(\omega)$ [130]. This method has been applied to CrO₂ in [120], using the same correlated 3o model. The exchange interactions were evaluated basically in the $\hat{b}$ scheme using equation (3), by considering the infinitesimal rotations of the frequency-dependent xc field $\Delta \Sigma(\omega) = \Sigma^{\uparrow}(\omega) - \Sigma^{\downarrow}(\omega)$ [131]. The dynamic correlations substantially reduce $\Delta \Sigma(\omega)$, resulting in much stronger AFM interactions J₇ and J₈, which overcome the effect of FM interactions J₁ and J₂, making the FM state unstable. The corresponding spin-wave dispersion is also plotted in figure 10, where the negative frequencies $\omega(\boldsymbol{q})$ along the directions Γ–X, Γ–${\mathrm{M}}$, and Γ–${\mathrm{M}}$ mean that the theoretical ground state should be in one of these points. Therefore, it was concluded in [120] that in order to reproduce the FM ground state of CrO₂, it is essential to extend the 3o model, by adding new ingredients such as the Cr $e_\mathrm{g}$ and O 2p bands. For instance, the spin-wave stiffness is systematically larger in the 5o model (table 11), which makes the FM state more stable. Nevertheless, the present analysis also suggests that the problem may be not in the 3o model itself, but in additional approximations underlying the scheme $\hat{b}$ for the exchange interactions. The method M systematically improves the stability of the FM ground state in CrO₂. This tendency may be overestimated in the Hartree–Fock approximation. However, a more rigorous treatment of correlation effects in the framework of DMFT, in the combination with the method M, could possibly bring the situation to a better agreement with the experimental data.

The DM interactions can take place between atoms in different Cr sublattices. The strongest ones are realized in the 2nd coordination sphere (in the combination with J₂, as shown in figure 8(b). If $\boldsymbol{\epsilon} = \frac{1}{\sqrt{2a^2+c^2}}(\pm a, \pm a, \pm c)$ are the directions of such bonds, the corresponding to them DM vectors have the following form $\boldsymbol{d} = d \epsilon^{z} [\boldsymbol{\epsilon} \times \boldsymbol{n}^{z}]$ (i.e. the vectors lie in the xy plane, are perpendicular to the bonds, and have opposite signs for the bonds in the directions ±z). The parameter d can be estimated as 0.85 and 1.38 meV, thus yielding $\boldsymbol{d} = (\pm 0.23, \pm 0.23, 0)$ and $(\pm 0.37, \pm 0.37, 0)$ meV in the Hartree–Fock approximation for correlated 5o model and all-electron LSDA, respectively (employing the scheme M in both cases). The interesting point is that the DM interactions are expected to be strong (and comparable with the ones in CrI₃) even though CrO₂ does not contain heavy elements. This is basically the consequence of two factors: (i) the partial filling of the $t_\mathrm{2g}$ states, opening the room for unquenched orbital magnetization; (ii) half-metallic character of the electronic structure, giving rise to the spin-current mechanism of the DM interactions in the $\uparrow$-spin $t_\mathrm{2g}$ band [58, 70]. Nevertheless, d is much smaller than the isotropic interaction J₂ operating in the same bonds.

Co₃Sn₂S₂ is a ferromagnet in which small spontaneous magnetization (about $0.3\mu_{\mathrm{B}}$ per Co atom) coexists with relatively high $T_{\mathrm{C}} = 177$ K [132]. It crystallizes in the rhombohedral $R\overline{3}m$ structure, hosting the kagome lattice of Co ions (figures 11(a)–(c)) [132]. The S atoms sit on the top of the Co₃ triangles, forming with them the network of alternating tetrahedra. There are two inequivalent Sn sites located in (Sn₂) and between (Sn₁) the kagome planes. Recently, Co₃Sn₂S₂ has attracted a lot of attention as a magnetic Weyl semimetal whose nontrivial topology of the electronic states gives rise to a large anomalous Hall effect [110–112]. Furthermore, Co₃Sn₂S₂ appears to be half-metallic, as was demonstrated in experimental [113] and theoretical [114, 115] studies.

Zoom In Zoom Out Reset image size

We have chosen Co₃Sn₂S₂ to illustrate the complexities one can face in the analysis of interatomic exchange interactions for this type of itinerant electron systems. At the present stage, we have in our disposition only all-electron ${\mathrm{Co}} \, 3d+{\mathrm{Sn}} \, 5p+{\mathrm{S}} \, 3p$ model, which was constructed in [115] within GGA [133], by employing the Vienna ab initio simulation package [134] for the electronic structure calculation and the maximal localization technique for the Wannier functions [135]. The corresponding DOS is shown in figure 11(h): the Fermi level crosses the $\uparrow$-spin band and falls inside the gap for the $\downarrow$-spin channel. The magnetic moments are $M_{\mathrm{Co}} = 0.33$, $M_{\mathrm{Sn}_{1}} = -0.03$, $M_{\mathrm{Sn}_{2}} = -0.04$, and $M_{\mathrm{S}} = 0.04\mu_{\mathrm{B}}$, so that the total moment per one unit cell is $1\mu_{\mathrm{B}}$. It is believed to be still an open question how good is GGA for Co₃Sn₂S₂. However, the attempts to construct a more compact model, which would also include the Coulomb correlations, were so far unsuccessful. Formally, the electronic structure in figure 11(h) can be divided into fully occupied 19 bands and remaining 8 bands, which are separated by a band gap. However, the construction of the Wannier basis for the upper eight-band manifold is not straightforward, as these Wannier functions will probably reside not on single atomic sites [136]. A compact model for Co₃Sn₂S₂ was constructed in [137], but using empirical arguments.

Co₃Sn₂S₂ provides an interesting example of itinerant magnetism. According to the GGA calculations [115], finite rotations of spins away from the FM ground state (for instance, by forcing three Co spins to form the umbrella texture) eventually leads to the collapse of magnetization, so that the system falls in the nonmagnetic state. ¹¹ A similar behavior was found in fcc Ni [138, 139] and several Ru-based oxides [140, 141]. Moreover, the size of the magnetic moments is expected to shrink due the temperature disorder, affecting both electronic structure and interatomic exchange interactions [115]. Therefore, the bilinear Heisenberg model cannot be defined in the global sense (for instance, for describing simultaneously low-temperature spin-wave dispersion and $T_{\mathrm{C}}$). Nevertheless, it can still be defined locally, for the analysis of local stability of the FM state with respect to the infinitesimal rotations of spins. The exchange interactions spread at least up to 9th coordination sphere around each Co site, as explained in figures 11(d)–(g). Moreover, some of these interactions, formally belonging to the same coordination sphere, can be inequivalent (as J₄ and $J_{4}^{\prime}$ in the plane, J₅ and $J_{5}^{^{\prime}}$, and J₉ and $J_{9}^{\prime}$ between the planes). The calculations have been performed on the mesh of the $20\times20\times20$ points both for k and q .

The effective Stoner parameters, calculated by using different definitions are summarized in table 12. As was discussed in section 3.3, the definitions II and IV should be used in combination with the spherically averaged scheme M, while the definitions I and III are more suitable for the matrix schemes $\hat{b}$ and $\hat{m}$. One can clearly see that all the parameters in Co₃Sn₂S₂ strongly depend on the definition. This holds even for $\mathcal{I}_{\mathrm{Co}}$, which can easily change by about 50%. This is clearly in contrast with other considered compounds, where the value of $\mathcal{I}$ on the magnetic transition-metal site practically did not depend on the definition. Nevertheless, the exchange interactions (37) between the Co sites do not explicitly depend on the choice of $\mathcal{I}_{\mathrm{Co}}$. The uncertainty with the choice of the parameters $\mathcal{I}_{\mathrm{Sn}}$ and $\mathcal{I}_{\mathrm{S}}$, which strongly depend on the way how they are defined, posses a more serious problem. The parameters $\mathcal{I}_{\mathrm{Sn}}$ are mainly negative (except $\mathcal{I}_{\mathrm{Sn}_{1}}$ in the case I) and tend to destabilize the FM interactions. On the other hand, positive $\mathcal{I}_{\mathrm{S}}$ strengthens the ferromagnetism. However, depending on the definition, $\mathcal{I}_{\mathrm{S}}$ can change by factor four, and $\mathcal{I}_{\mathrm{Sn}}$ changes by an order of magnitude. This dependence has a strong impact on the exchange interactions, which are summarized in table 13.

Table 12. The effective Stoner parameters (in eV) obtained using definitions I–IV (as explained in section 3.3.1). Sn₂ and Sn₁ are located, respectively, in and between kagome planes.

Definition	$\mathcal{I}_{\mathrm{Co}}$	$\mathcal{I}_{\mathrm{Sn}_{1}}$	$\mathcal{I}_{\mathrm{Sn}_{2}}$	$\mathcal{I}_{\mathrm{S}}$
I	1.30	$\phantom{-}0.17$	−0.10	1.22
II	0.98	−0.06	−0.20	0.88
III	1.55	−0.65	−0.65	2.84
IV	1.45	−3.44	−2.53	3.85

Table 13. Parameters of exchange interactions in Co₃Sn₂S₂ (in meV) calculated using the schemes $\hat{b}$, $\hat{m}$, and M (the infinitesimal rotations of, respectively, the xc field, magnetization matrix, and local spin moments). In the scheme $\hat{b}$, the xc field was either associated with the site-diagonal part of the TB Hamiltonian (H) or derived from the sum rule (sr). The roman number in the parentheses stands for the set of the parameters $\mathcal{I}_{\mathrm{Sn}}$ and $\mathcal{I}_{\mathrm{S}}$ from table 12, which was used in the calculations of J_k . In the L0 scheme, xc field was redefined to ensure $\mathcal{I} = 0$ on the sites Sn and S. The notations of J_k are explained in figures 11(d)–(g). $D^{xx} = D^{yy}$ and D^zz are nonvanishing elements of the spin-stiffness tensor (in meV·Å²).

	$\hat{b}$ (H,I)	$\hat{b}$ (sr,III)	$\hat{m}$ (III)	M (II)	M (IV)	L0
J1	$\phantom{-}1.60$	$\phantom{-}2.54$	$\phantom{-}7.71$	−0.08	−1.51	−0.09
J2	$\phantom{-}0.05$	$\phantom{-}0.08$	−0.14	−0.06	−0.25	−0.06
J3	$\phantom{-}0.09$	$\phantom{-}0.14$	$\phantom{-}0.22$	−0.06	−1.11	−0.07
J4	$\phantom{-}0.18$	$\phantom{-}0.22$	−0.01	$\phantom{-}0.34$	$\phantom{-}0.35$	$\phantom{-}0.41$
$J_{4}^{\prime}$	$\phantom{-}0.54$	$\phantom{-}0.77$	$\phantom{-}0.77$	$\phantom{-}0.61$	$\phantom{-}0.52$	$\phantom{-}0.73$
J5	$\phantom{-}0.42$	$\phantom{-}0.58$	$\phantom{-}0.88$	$\phantom{-}0.83$	−0.23	$\phantom{-}0.98$
$J_{5}^{\prime}$	$\phantom{-}0.65$	$\phantom{-}0.90$	$\phantom{-}1.56$	$\phantom{-}0.95$	$\phantom{-}0.92$	$\phantom{-}1.13$
J6	$\phantom{-}0.29$	$\phantom{-}0.41$	$\phantom{-}0.64$	$\phantom{-}0.44$	$\phantom{-}0.51$	$\phantom{-}0.53$
J7	$\phantom{-}0.10$	$\phantom{-}0.14$	$\phantom{-}0.08$	$\phantom{-}0.12$	$\phantom{-}0.18$	$\phantom{-}0.15$
J8	$\phantom{-}0.04$	$\phantom{-}0.06$	−0.13	$\phantom{-}0.10$	$\phantom{-}0.07$	$\phantom{-}0.11$
J9	−0.13	−0.19	−0.66	−0.27	−0.25	−0.31
$J_{9}^{\prime}$	$\phantom{-}0.07$	$\phantom{-}0.10$	$\phantom{-}0.03$	$\phantom{-}0.02$	$\phantom{-}0.16$	$\phantom{-}0.03$
Dxx	$\phantom{-}521$	$\phantom{-}727$	$\phantom{-}383$	$\phantom{-}490$	$\phantom{-}36$	$\phantom{-}582$
Dzz	$\phantom{-}246$	$\phantom{-}317$	$\phantom{-}73$	$\phantom{-}310$	$\phantom{-}48$	$\phantom{-}374$

As was pointed out before, the attempts to estimate the Curie temperature entirely from J_k , defined via the infinitesimal rotations of spins near the FM ground state, are pretty much meaningless in the case of Co₃Sn₂S₂, where proper theories for $T_{\mathrm{C}}$ should consider also the longitudinal change of the magnetization [115]. Nevertheless, these J_k can be used to evaluate the spin-wave dispersion and the spin stiffness. The latter has been measured experimentally and the nonvanishing elements of the spin-stiffness tensor are $D^{xx} = D^{yy} = 803 \pm 46$ and $D^{zz} = 237 \pm 13$ meV·Å² [142], which can be used for comparison with theoretical data.

First, we note that the results of the scheme $\hat{b}$ strongly depend on the definition of the xc field: the enforcement of the sum rule for $\hat{b}$ strengthens the FM interactions, bringing the spin stiffness to a good agreement with the experimental data. In this case, the exchange parameters are given mainly by the bare interactions, corresponding to the first term in equation (36). The corrections caused by the ligand states are small and do not play a decisive role. An interesting situation is realized in the scheme $\hat{m}$, based on the rigid rotations of the magnetization matrix: as expected, the individual parameters 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 exchange interactions in the scheme M are very sensitive to $\mathcal{I}_{\mathrm{Sn}}$ and $\mathcal{I}_{\mathrm{S}}$. If one uses the definition II, where the xc field in the expression for $\mathcal{I}$ is taken from the site-diagonal part of the TB Hamiltonian, $\mathcal{I}_{\mathrm{Sn}}$ and $\mathcal{I}_{\mathrm{S}}$ are relatively small so that the exchange interactions are given basically by the first term of equation (37). In this case, J₁–J₃ are weakly AFM, while other interactions are FM (except J₉, which is AFM in all considered methods). As the result, the FM state is stable, though D^xx is somewhat underestimated in comparison with the experiment.

The situation changes dramatically if one uses the definition IV, where the xc field is taken from the sum rule, which substantially strengthens both $\mathcal{I}_{\mathrm{Sn}}$ and $\mathcal{I}_{\mathrm{S}}$. Apparently, negative $\mathcal{I}_{\mathrm{Sn}}$ plays a dominant role and is responsible for the AFM character of interactions in some of the bonds: particularly, J₁–J₃ become prominently AFM and J₅ changes from strongly FM to AFM. Then, the FM state become unstable and the spins tend to form the 120° in the xy plane (at least, on the mean-field level). The corresponding spin stiffness is strongly underestimated compared to the experimental data. This is clearly an artifact of the model analysis, which is probably caused by unrealistic estimate for $\mathcal{I}_{\mathrm{Sn}}$. At least, the conclusion does not seem to be consistent with the brute-force GGA calculations, where the 120° alignment of spins in the xy plane leads to the collapse of the magnetic state [115].

The L0 approach allows us to eliminate the dependence of J_k on $\mathcal{I}_{\mathrm{Sn}}$ and $\mathcal{I}_{\mathrm{S}}$, assuming that all magnetic moments are induces by $B_{\mathrm{Co}}$, while $B_{\mathrm{Sn}} = B_{\mathrm{S}} = 0$ (so as $\mathcal{I}_{\mathrm{Sn}}$ and $\mathcal{I}_{\mathrm{S}}$). This results only in a small change of magnetic moments in comparison to the plain GGA values reported above: $M_{\mathrm{Co}} = 0.36$, $M_{\mathrm{Sn}_{1}} = -0.05$, $M_{\mathrm{Sn}_{2}} = -0.07$, and $M_{\mathrm{S}} = 0.02\mu_{\mathrm{B}}$ (thus, $M_{\mathrm{tot}}$ remains equal to $1\mu_{\mathrm{B}}$). The exchange interactions in these case reminiscent the ones obtained in the scheme M using the definition II for $\mathcal{I}_{\mathrm{Sn}}$ and $\mathcal{I}_{\mathrm{S}}$ with somewhat stronger tendency toward the ferromagnetism, which is reflected in larger values of D^xx and D^zz .

In principle, the schemes $\hat{b}$ and L0 both provide a reasonably good agreement with the experimental data for the spin-stiffness constants. Nevertheless, the behavior of exchange interactions is quite different. In the scheme $\hat{b}$, the strongest FM interaction is J₁. The longer-range interactions operating practically at the same distance in the 4th and 5th coordination spheres, in and between the kagome planes (see figures 11(e) and (g)), are also sizable but smaller than J₁. In the scheme L0, J₁ is weakly AFM, while the strongest FM interactions take place in the 4th and 5th coordination spheres. The corresponding spin-wave dispersions are plotted in figure 12. In the long wavelength limit $\boldsymbol{q} \to 0$, two methods provide very similar description, while the main difference is in the position of optical branches in the higher-energy region. The results of recent inelastic neutron scattering data (which were nevertheless limited by the behavior of the acoustic branch in the vicinity of $\boldsymbol{q} = 0$) where interpreted in terms of four interactions (in our notations): $J_{2} = -0.32$, $J_{3} = 0.64$, $\bar{J}_{4} = 1.66$, and $\bar{J}_{5} = 0.09$ meV [143], where $\bar{J}_{k} = \frac{1}{3}(J_{k} + 2J_{k}^{\prime})$ stands for the averaged exchange interactions in the 4th and 5th coordination spheres. ¹² There is certain similarity with the picture provided by the L0 scheme: at least J₁ is negligibly small, J₂ is AFM, and the main FM interactions occurs in the next coordination spheres. Nevertheless, there are also the differences. Especially, experimental $\bar{J}_{5}$ appears to be smaller than $\bar{J}_{4}$, while theoretical parameters are at least comparable (see table 13). According to the theoretical analysis, the values of inequivalent parameters J_k and $J_{k}^{\prime}$ can be very different and this difference was not considered in the fitting of the experimental spin-wave spectrum. On the other hand, it is not clear whether it can resolve the discrepancy between theoretical and experimental data. It is also unclear whether the experimental data available only in the vicinity of $\boldsymbol{q} = 0$ are enough to make a decisive conclusion about the complexity of the exchange interactions in Co₃Sn₂S₂.

Zoom In Zoom Out Reset image size

The Co₃Sn₂S₂ structure has three Co sublattices (which are transformed to each other by the threefold rotations). The spacial inversion transforms each Co atom to the same sublattice. Thus, there can be no DM interactions within the sublattices. Nevertheless, the DM interactions between the sublattices are permitted by the $R\overline{3}m$ space group. The strongest ones take place between nearest neighbors in the kagome plane. They are explained in figure 13. Using the convention, where the directions of the bonds ( ε ) in the triangles can be transformed to each other by the threefold rotations, the DM vectors can be presented in the form: $\boldsymbol{d} = d^{xy}[\boldsymbol{\epsilon} \times \boldsymbol{n}^{z}] + d^{z}\boldsymbol{n}^{z}$. If the FM moment is parallel to z (the experimental situation), d^z does not contribute to the energy change, while d^xy is responsible for the canting of spins and tends to form the umbrella texture. If $\boldsymbol{e}_{k} = (\theta \text{sin} \frac{\pi k}{3}, \theta \text{cos} \frac{\pi k}{3}, 1-\frac{\theta^2}{2})$ are the directions of spins in the triangle (k = 1, 2, or 3, as explained in figure 13), the energy gain due to the DM interaction is $\delta E_{\mathrm{DM}} = -\sqrt{3}d^{xy} \theta$ (per one Co site). The corresponding energy loss due to the isotropic exchange is $\delta E_{\mathrm{H}} = \frac{3}{4} J_{1,23} \theta^2$, where $J_{1,23}$ is the effective interaction of the atom in the 1st sublattice with the atoms of the 2nd and 3rd sublattices in all the coordination spheres. The values of the parameters obtained in the scheme L0 are $J_{1,23} = 6.04$ meV, $d^{xy} = 0.20$ meV, and $d^{z} = -0.29$ meV. ¹³ Thus, the angle $\theta = \frac{2d^{xy}}{\sqrt{3}J_{1,23}}$ can be estimates as 2°, which perfectly agrees with $\theta \sim 2^\circ$ obtained in the brute-force GGA calculations with the SO coupling [115].

Zoom In Zoom Out Reset image size

First, we investigate abilities of correlated 5o model: whether it can reproduce the main tendencies in the behavior of interatomic exchange interactions in LaMnO₃ and HoMnO₃, stabilizing the A-type AFM state in the former case and making it unstable with respect to the formation of the spin superstructures along the orthorhombic axis b in the latter one. Another important question is how good is the approximate scheme $\hat{b}$, which is frequently used for the analysis of interatomic exchange interactions in these orthorhombic manganites [42, 145], in comparison with the more accurate scheme M. The averaged on-site parameters of the Coulomb repulsion, intraatomic exchange interaction, and nonsphericity (see appendix A for definitions) calculated in cRPA are, respectively, 2.15 (2.16), 0.85 (0.85), and 0.09 (0.09) eV for LaMnO₃ (HoMnO₃) [42]. Thus, the Coulomb $U\sim2$ eV is not particularly strong. ¹⁴ This finding is very important as it naturally explains the existence of the long-range exchange interaction, $J_{2}^{b}$ and $J_{3}^{1}$, arising in higher orders of $\hat{H}/U$ beyond the conventional superexchange approximation [2], and responsible for the formation of the spin superstructures along b. On the other hand, the relatively small U also means that the strong-coupling limit is hardly satisfied. Therefore, it is reasonable to expect that the approximate scheme $\hat{b}$ may experience serious limitations for these orthorhombic manganites.

The DOSs obtained in the Hartree–Fock approximation for the A-type AFM order are displayed in figure 15. LaMnO₃ and HoMnO₃ have similar electronic structure. The only difference is slightly smaller $e_\mathrm{g}$ bandwidth in the case of HoMnO₃. Nevertheless, the behavior exchange interactions appears to be very different. These interactions are summarized in tables 15 and 16 (lines 5o) for LaMnO₃ and HoMnO₃, respectively.

Zoom In Zoom Out Reset image size

Table 15. Isotropic exchange interactions in LaMnO₃ (in meV) as obtained in the correlated five-orbital (5o) model in comparison with LSDA for the all-electron ${\mathrm{Mn}} \, 3d+{\mathrm{O}} \, 2p+{\mathrm{La}} \, 5d$ model. The notations of parameters are explained in figures 14(b) and (c). The AFM interactions $J_{2}^{a}$ and $J_{3}^{2}$ are considerably weaker and not shown here. $T_{\mathrm{N}}$ is the Néel temperature evaluated in RPA (in K) and $\boldsymbol{q} = (0,q_{b},0)$ is the theoretical ground state propagation vector (in units of reciprocal lattice translations, where $\boldsymbol{q} = 0$ corresponds to the A-type AFM state).

Method	$J_{1}^{\parallel}$	$J_{1}^{\perp}$	$J_{2}^{1}$	$J_{2}^{2}$	$J_{2}^{b}$	$J_{3}^{1}$	$T_{\mathrm{N}}$	qb
$\hat{b}$ in 5o	$\phantom{1}2.36$	−6.64	−0.97	−1.10	−1.00	−3.27	80	0.32
$M$ in 5o	17.23	$\phantom{-}3.51$	−1.42	−1.51	−1.13	−3.56	186	0
$\hat{b}$ in LSDA	15.10	$\phantom{-}7.37$	−3.14	−3.25	−1.78	−9.46	131	0
$M$ in LSDA	17.09	$\phantom{-}8.20$	−3.24	−3.27	−1.83	−9.21	185	0

Table 16. The same as table 15 but for HoMnO₃. The LSDA results are obtained in the ${\mathrm{Mn}} \, 3d+{\mathrm{O}} \, 2p+{\mathrm{Ho}} \, 5d$ model.

Method	$J_{1}^{\parallel}$	$J_{1}^{\perp}$	$J_{2}^{1}$	$J_{2}^{2}$	$J_{2}^{b}$	$J_{3}^{1}$	$T_{\mathrm{N}}$	qb
$\hat{b}$ in 5o	−6.53	−7.58	−0.54	−0.72	−1.44	−2.22	65	0.22
$M$ in 5o	$\phantom{-}4.02$	−1.44	−0.63	−0.82	−1.53	−2.36	53	0.26
$\hat{b}$ in LSDA	$\phantom{-}2.23$	−1.64	−1.72	−1.85	−2.57	−6.90	79	0.39
$M$ in LSDA	$\phantom{-}4.15$	−1.03	−1.71	−1.81	−2.46	−6.50	79	0.34

The scheme $\hat{b}$ yields strong interlayer coupling $J_{\phantom{1}}^{\perp} = J_{1}^{\perp} + 2J_{2}^{1} + 2J_{2}^{2} = -10.78$ (−10.1) meV for LaMnO₃ (HoMnO₃), which is even stronger than the intralayer one $J_{1}^{\parallel} = 2.36$ (−6.53) meV. This behavior is inconsistent with the experimental neutron-scattering data for LaMnO₃, indicating that $J_{\phantom{1}}^{\perp} \sim -4.7$ meV is weaker than $J_{1}^{\parallel} = 6.6$ meV [158, 159]. ¹⁵ Furthermore, the longer-range AFM interactions $J_{2}^{b}$ and $J_{3}^{1}$ in LaMnO₃ are comparable or even stronger than $J_{1}^{\parallel}$, making the experimental A-type AFM structure unstable.

On the contrary, the scheme M systematically improves the description of the interatomic exchange interactions. First, the interlayer coupling becomes considerably weaker: $J_{\phantom{1}}^{\perp} = -2.35$ (−4.34) meV for LaMnO₃ (HoMnO₃). Then, the itralayer interaction $J_{1}^{\parallel}$ in LaMnO₃ becomes strongly FM, as expected for the 'antiferro' orbital ordering [160], ¹⁶ and overcomes the AFM long-range interactions $J_{2}^{b}$ and $J_{3}^{1}$. As the result, the experimental A-type AFM phase becomes stable. The theoretical $T_{\mathrm{N}} = 186$ K, evaluated in RPA, is in fair agreement with the experimental value of 140 K [158, 159]. In HoMnO₃, $J_{1}^{\parallel}$ is significantly reduced due to the additional buckling of the Mn–O–Mn bonds to become comparable with $J_{2}^{b}$ and $J_{3}^{1}$. This makes the A-type AFM phase unstable. Regarding the direction of this instability, it is very important that $J_{2}^{a} = -0.35$ meV is much weaker than $J_{2}^{b} = -1.53$ meV. Therefore, the propagation vector for the expected theoretical ground state is parallel to b, in agreement with the experimental observation. This instability is further enhanced by the AFM interactions $J_{3}^{1}$. Nevertheless, in order to reproduce the experimental E phase with commensurate $\boldsymbol{q} = (0,\frac{1}{2},0)$, it is essential to consider other ingredients such as the exchange striction and single-ion anisotropy, which would lock the spin superstructure to the lattice [161–163]. In HoMnO₃, such lock-in transition occurs at 29 K, while the Néel temperature is about 41 K [149, 150], which is close to theoretical value of $T_{\mathrm{N}} = 53$ K. Other aspects of stability of the E-type AFM phase will be considered in section 5.4.

The all-electron model in LSDA provides an alternative description for the exchange interactions. We start with the analysis of effective Stoner parameters. As was discussed in section 3.3.1, among several possible definitions of the parameters $\mathcal{I}$, the most relevant seem to be III and IV, which should be considered in combination with the methods $\hat{b}$ and M, respectively. In the A-type AFM phase, the A and apical O atoms, located between the antiferromagnetically coupled layers, are nonmagnetic. Therefore, we set $\mathcal{I} = 0$ for them. The remaining parameters for Mn and planar O atoms are listed in table 17. $\mathcal{I}_{\mathrm{Mn}}$ is practically identical for LaMnO₃ and HoMnO₃, and does not depend on the definition. $\mathcal{I}_{\mathrm{O}}$ appears to be more sensitive to the environment and the definition. Nevertheless, the change of $\mathcal{I}_{\mathrm{O}}$ is rather modest. For both definitions, $\mathcal{I}_{\mathrm{O}}$ is large and positive, that will additionally strengthens the FM interactions.

Table 17. The effective Stoner parameters (in eV) for Mn and planar O atoms in LaMnO₃ and HoMnO₃, obtained using definitions III and IV (as explained in section 3.3.1). The A and apical O atoms are nonmagnetic in the A-type AFM phase and not considered here.

Definition	LaMnO3		HoMnO3
	$\mathcal{I}_{\mathrm{M}n}$	$\mathcal{I}_{\mathrm{O}}$	$\mathcal{I}_{\mathrm{M}n}$	$\mathcal{I}_{\mathrm{O}}$
III	0.98	3.71	0.97	4.20
IV	0.98	4.35	0.97	4.23

The exchange interactions are summarized in tables 15 and 16 (lines LSDA). Unlike for the correlated 5o model, the methods $\hat{b}$ and M in the all-electron LSDA provide a consistent description. On the one hand, there is an effect of the O 2p band, which is explicitly treated in this model. On the other hand, the underestimation of the FM interactions in the scheme $\hat{b}$ is partly compensated by $\mathcal{I}_{\mathrm{O}}$. Moreover, the correct definition of the xc field $\hat{b}$ using the sum rule is important. For instance, the use of $\hat{b}$ defined via site-diagonal elements of $\hat{H}^{\uparrow, \downarrow}$ underestimates the FM interactions and makes the experimental A-type AFM state unstable in LaMnO₃.

The exchange interactions are generally stronger than in correlated 5o model (except $J_{1}^{\parallel}$), partly due to the fact that in these insulating materials the exchange interactions are expected to increase when the effective Coulomb repulsion decreases [2]. Furthermore, the oxygen band can also contribute to the exchange interactions [22, 164]. Anyway, the total interlayer coupling in LaMnO₃, $J_{\phantom{1}}^{\perp} = -4.82$ (−5.41) in the scheme M ($\hat{b}$), is consistent with the experimental data [158, 159]. $J_{1}^{\parallel}$ is strongly FM, which appears to be sufficient to overcome the strong AFM interactions $J_{2}^{b}$ and $J_{3}^{1}$, and stabilize the experimental A-type AFM phase. In HoMnO₃, this $J_{1}^{\parallel}$ is substantially weaker, making the A-type AFM phase unstable with respect to an incommensurate spin structure propagating along b. The theoretical $T_{\mathrm{N}}$ is in fair agreement with experimental data.

All nn DM interactions in the orthorhombic planes, $\boldsymbol{d}^{\, \parallel}_{ij}$, are transformed to each other by the symmetry operations of the space group Pbnm [145]. The same holds for the interplane interactions $\boldsymbol{d}^{\perp}_{ij}$. Furthermore, since neighboring planes are connected by the mirror reflection, the z (c) component of $\boldsymbol{d}^{\perp}_{ij}$ is equal to zero. Therefore, there are five parameters describing all nn DM interactions: $d^{\, \parallel}_{a}$, $d^{\, \parallel}_{b}$, and $d^{\, \parallel}_{c}$ for the in-plane interactions, and $d^{\perp}_{a}$ and $d^{\perp}_{b}$ for the out-of-plane interactions. The corresponding DM vectors, attached to neighboring Mn–O–Mn bonds, are shown in figure 16. The values of the parameters, obtained in the scheme M, are summarized in table 18.

Zoom In Zoom Out Reset image size

Table 18. Parameters of nearest-neighbor DM interactions (in meV) as obtained in the correlated five-orbital (5o) model and all-electron LSDA using the scheme M. The results of [145], employing mixed perturbation theory, are shown for comparison.

Compound	Model	$d^{\, \parallel}_{a}$	$d^{\, \parallel}_{b}$	$d^{\, \parallel}_{c}$	$d^{\perp}_{a}$	$d^{\perp}_{b}$
LaMnO3	[145]	0.44	0.33	0.53	0.45	0.71
LaMnO3	LSDA	0.46	0.35	0.63	0.53	0.60
HoMnO3	LSDA	0.54	0.38	0.73	0.70	0.47
HoMnO3	5o	0.23	0.09	0.28	0.20	0.11

For LaMnO₃ in LSDA we note a good agreement with the results of [145], employing the mixed perturbation theory, where the rotation of the xc field on one Mn sites was combined with the SO coupling on another such site. This means that the 5d states of the heavy A atoms, which are located far in the unoccupied part of the spectrum (see figure 14), do not strongly contribute to the DM interactions. Partly, this may be due to the fact that the A sites remain nonmagnetic in the A-type AFM phase.

The parameters obtained in the correlated 5o model are generally smaller than in LSDA—similar to what was found for the isotropic interactions. However, this is not surprising and can be again explained by the Coulomb U in the denominator of the superexchange interactions [2].

The DM interactions in AMnO₃ mainly contribute to the spin canting. The magnetocrystalline anisotropy tends to align the spins parallel to the b axis [165–167]. In LaMnO₃, they order according to the A-type: $\boldsymbol{e}_{\mu} = (0,-1,0)$ for µ = 1 and 2 and $\boldsymbol{e}_{\mu} = (0,1,0)$ for µ = 3 and 4 in figure 16. This perfect AFM alignment will be further deformed by the DM interactions, which additionally rotate the spins along a and c by, respectively, $e_{a} = -d^{\, \parallel}_{c}/(J_{1}^{\parallel}-J_{2}^{1}-J_{2}^{2})$ and $e_{c} = d^{\, \perp}_{a}/(J_{1}^{\parallel}+2J_{2}^{1}+2J_{2}^{2})$ [62]. The a components will order according to the G-type, ¹⁷ while c components gives rise to the weak ferromagnetism [165, 166]. Using the obtained parameters of exchange interactions, $|e_{a}|$ and $|e_{c}|$ can be estimated as 0.019 and 0.110, respectively. Taking into account that in all-electron LSDA $M = 3.59\mu_{\mathrm{B}}$, the weak FM moment can be estimated as $0.4\mu_{\mathrm{B}}$ per Mn atom (being somewhat larger that the experimental $0.18\mu_{\mathrm{B}}$, derived from magnetization measurements [168]).

An interesting question is whether the multiferroicity associated with the E-type AFM phase in HoMnO₃ can coexist with the weak ferromagnetism. In the E-type AFM structure, each of the $J_{2}^{1}$ and $J_{2}^{2}$ will contribute to two FM bonds and two AFM ones. Therefore, these contributions cancel each other and the spin canting along c will be given by $e_{c} = \pm d^{\, \perp}_{a}/J_{1}^{\parallel}$ (for the sublattices with different spins in the ab plane). Using the parameters for the Pbnm structure of HoMnO₃, $|e_{c}|$ can be estimated as 0.169 and 0.050 in the all-electron LSDA and correlated 5o model, respectively, where the difference mainly comes from $d^{\, \perp}_{a}$. Nevertheless, this component will repeat the pattern of the E-type AFM phase in each ab plane and, therefore, there will be no weak ferromagnetism.

An interesting aspect of interatomic exchange interactions defined via infinitesimal rotations of spins near some equilibrium magnetic states is that these exchange interactions depend on the magnetic state in which they are calculated, reflecting the dependence of the electronic structure on the magnetic state. If it happens, the bilinear spin model (1) is ill-defined in the global sense, meaning, for instance, that the same set of the exchange parameters cannot be used for the analysis of the low-temperature spin-wave dispersion and the magnetic transition temperature, where the electronic structure is strongly modified by the spin disorder [22, 128]. Nevertheless, the model (1) can be defined locally, near each magnetic equilibrium. Then, in principle, one can expect some kind of self-organization phenomenon, when the change of the electronic structure in certain magnetic state itself can stabilize this magnetic state, at least locally. This is what happens at least with some of the exchange interactions in the 50% doped manganites, where the zigzag (CE-type) AFM alignment opens the band gap, induces the orbital ordering, and additionally stabilizes the FM coupling in the zigzag chains [98]. Can we expect a similar behavior in the E phase?

Intuitively, the reason for the formation of this noncentrosymmetric E-type AFM structure can be understood as follows: the AFM interactions $J_{2}^{b}$ and $J_{3}^{1}$ tend to align the corner spins, separated by the orthorhombic translations a and b, antiferromagnetically, as shown in figure 17(a). The central Mn atom is located in the inversion center. Then, due to the AFM alignment, the corner atoms should transform to each other by the symmetry operation $\hat{I}\hat{T}$ (where the spacial inversion $\hat{I}$ is combined with the time reversal $\hat{T}$), which is formally compatible with the Pbnm space group. However, the same symmetry operation would make the central Mn atom nonmagnetic. This is an eligible scenario from the symmetry point of view, but would lead to gigantic energy loss, $\frac{1}{4}\mathcal{I}_{\mathrm{Mn}}M_{\mathrm{Mn}}^2\sim 4$ eV, caused by the violation of the first Hund's rule for the atoms, which are potentially expected to be in the high-spin state. Therefore, the more favorable scenario is to keep the Mn atom magnetic, but break the inversion symmetry, which is quite common in multiferroic materials, especially those with high-spin ions.

Zoom In Zoom Out Reset image size

Then, the next question is whether the E-type AFM order is compatible with the change of the exchange interactions induced by this order.

Below we consider results of correlated 5o model (in the scheme M) for HoMnO₃. A qualitatively similar behavior has been found in LSDA for the all-electron model. The E-type AFM order results in the additional narrowing of the $t_\mathrm{2g}$ and $e_\mathrm{g}$ bands (figure 17), which is consistent with the decrease of the number of FM bonds around each Mn site (2 instead of 4 in the A-type AFM phase). The inversion symmetry breaking makes the nn interactions inequivalent, where $J^{\uparrow \uparrow}_{1}$ in the FM bond is generally different from $J^{\uparrow \downarrow}_{1}$ in the AFM one. In order to stabilize the E state, one would need at least $J^{\uparrow \uparrow}_{1} \gt J^{\uparrow \downarrow}_{1}$ (the FM coupling is stronger in the FM bond). However, we have found the opposite tendency: $J^{\uparrow \uparrow}_{1} = 2.99$ meV and $J^{\uparrow \downarrow}_{1} = 4.00$ meV, meaning that the E state corresponds to the energy maximum and any small rotations of spins will slide the system away from this equilibrium towards a new magnetic state. Furthermore, there is an intrinsic mechanism inducing these rotations of spins. Indeed, the E-type AFM order, producing some changes in the electronic structure, can induce the electric polarization even in the centrosymmetric Pbnm structure [161]. Then, it is reasonable to expect that the same changes in the electronic structure may lead to the appearance of new DM interactions, which would otherwise be forbidden in the Pbnm structure. Particularly, we have found finite interactions $\boldsymbol{d}_{2}^{b} = (\pm 0.01,0,\pm 0.02)$ meV and $\boldsymbol{d}_{3}^{1} = (\pm 0.01,\pm 0.01,\pm 0.02)$ meV (being analogs of isotropic $J_{2}^{b}$ and $J_{3}^{1}$, respectively, where the ± signs depend on the origin and direction of the bond). Although these interactions are not particularly strong, they are sufficient to shift the system of spins 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 HoMnO₃ is intrinsically unstable and cannot be stabilized by purely electronic mechanisms. For these purposes, it is essential to consider the exchange striction or the single-ion anisotropy (or both of them) [161–163].

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