The DMI matrix D*ij* is antisymmetric and therefore has only three independent components, which can be expressed by a vector *Dij* with *Dij*,*x* = D*ij*,*yz*, *Dij*,*y* = D*ij*,*zx*, and *Dij*,*z* = D*ij*,*xy*. Thus, the DMI can be expressed by cross product: *S T i* D*ijSj* = *Dij* · *Si* × *Sj* . Such an interaction prefers the vectors *Si* and *Sj* to be orthogonal to each other, with a rotation (of *Sj* relative to *Si* ) around the direction of −*Dij*. Together with Heisenberg term *JijSi* · *Sj* , the preferred rotation angle between *Si* and *Sj* relative to the collinear state preferred by the Heisenberg term would be arctan|*Dij*| |*Jij*| . In Ref. [106], the DMI is shown to determine the chirality of the magnetic ground state of Cr trimers on Au(111). The DMI is also important in explaining the skyrmion states in many materials such as MnSi and FeGe [24,26,27,29,107–110]. Materials with skyrmion states induced by DMI usually have a large ratio of |*D*1| |*J*1| (typically 0.1~0.2) where the subscript "1" means nearest pairs [24,100]. In Ref. [100], strong enough DMI for the existence of helical cycloid phases and skyrmionic states are predicted in Cr(I,X)3 (X = Br or Cl) Janus monolayers (e.g., for Cr(I,Br)3, supposing |*Si* | = 3 2 for any *i*, the corresponding interaction parameters are computed as *J*1 = −1.800 meV and |*D*1| = 0.270 meV, thus |*D*1| |*J*1| = 0.150), though monolayers such as CrI3 only exhibit an FM state for lack of DMI. In Ref. [110], the nonreciprocal magnon spectrum (and the associated spectral weights) of MnSi, as well as its evolution as a function of magnetic field, is explained by a model including symmetric exchange, DMI, dipolar interactions, and Zeeman energy (related to the magnetic field).

The Kitaev matrix K*ij* has five independent components as a symmetric matrix with zero trace. For the specific cases when *Si* and *Sj* are parallel to each other (pointing in the same direction), *S T i* K*ijSj* would perform like *S T i* A*iSi* and show preference to the direction with the lowest eigenvalue of K*ij*; while when *Si* and *Sj* are antiparallel to each other, *S T i* K*ijSj* would prefer the direction with the highest eigenvalue of K*ij*. The difference between the highest and lowest eigenvalue of K*ij* can be defined as *Kij* (a scalar), which characterizes the anisotropic contribution of K*ij*. Generally, the favorite direction of the spins is decided by both SIA and Kitaev interactions. The long-range ferromagnetic order in monolayer CrI3 was explained by the anisotropic superexchange interaction since the Cr-I-Cr bond angle is close to 90◦ [111]. In Ref. [112], the interplay between the prominent Kitaev interaction and SIA was studied to explain the different magnetic behaviors of CrI3 and CrGeTe3 naturally. For CrI3, supposing |*Si* | = 3 2 for any *i*, the *Jij* and *Kij* parameters between nearest pairs are computed as −2.44 and 0.85 meV, respectively; while the only independent component A*i*,*zz* of A*i* is −0.26 meV. For CrGeTe3, these three parameters are calculated to be −6.64, 0.36, and 0.25 meV, respectively. These two kinds of interactions are induced by SOC of the heavy ligands (I or Te) in these two materials (rather than the commonly believed Cr ions). Among different types of quantum spin liquids (QSLs), the exactly solvable Kitaev model with a ground state being QSL (with Majorana excitations) [113] has attracted much attention. Materials that achieve the realization of such Kitaev QSLs as *α*-RuCl3 [114–116] and (Na1-*x*Li*x*)2IrO3 [117,118] (with an effective *S* = 1/2 spin value) with honeycomb lattices are discovered. A possible Kitaev QSL state is also predicted in epitaxially strained Cr-based monolayers with *S* = 3/2, e.g., CrSiTe3 and CrGeTe3 [119].

#### *2.4. Fourth-Order Interactions without SOC*

Sometimes, higher-order interactions are also crucial for explaining the magnetic properties of some materials, especially if the magnetic atoms have large magnetic moments or if the system is itinerant. As mentioned in Section 2.3, when SOC can be ignored, the

effective spin Hamiltonian should only include inner product terms of spins. Besides second-order Heisenberg terms like *JijSi* · *Sj* , the terms with the next lowest order (which is fourth-order) are biquadratic (exchange) terms like *Kij Si* · *Sj* 2 , three-body fourth-order terms like *Kijk Si* · *Sj* (*Si* · *Sk*), and four-spin ring coupling terms like *Kijkl Si* · *Sj* (*Sk* · *Sl*). That is to say, when SOC and external magnetic field are ignored, keeping the terms with orders no higher than fourth, the effective spin Hamiltonian can be expressed as

$$H\_{\rm spin} = \sum\_{i,j>i} I\_{ij} \mathbf{S}\_i \cdot \mathbf{S}\_j + \sum\_{i,j>i} K\_{ij} \left(\mathbf{S}\_i \cdot \mathbf{S}\_j\right)^2 + \sum\_{\substack{i,j>i\\k>j}} K\_{ijk} \left(\mathbf{S}\_i \cdot \mathbf{S}\_j\right) \left(\mathbf{S}\_i \cdot \mathbf{S}\_k\right) + \sum\_{\substack{i,j>i,\\k>j>i,\\l>k}} K\_{ijkl} \left(\mathbf{S}\_i \cdot \mathbf{S}\_j\right) \left(\mathbf{S}\_k \cdot \mathbf{S}\_l\right) \tag{10}$$

The biquadratic terms have been found to be important in many systems, such as MnO [120,121], YMnO3 [74], TbMnO3 [122], iron-based superconductor KFe1.5Se2 [123], and 2D magnets [124]. In the case of TbMnO3, besides the biquadratic terms, the four-spin couplings are also found to be important in explaining the non-Heisenberg behaviors [122]; the three-body fourth-order terms are also found to be important in simulating the total energies of different spin configurations [125] (a list of the fitted values of each important interaction parameter in TbMnO3 is provided in the supplementary material of Ref. [125]). According to Ref. [126], in a Heisenberg chain system constructed from alternating *S* > 1 2 and *S* = 1 2 site spins, the additional isotropic three-body fourth-order terms are found to stabilize a variety of partially polarized states and two specific non-magnetic states including a critical spin-liquid phase and a critical nematic-like phase. In Ref. [127], the four-spin couplings were found to have a large effect on the energy barrier preventing skyrmions (or antiskyrmions) collapse into the ferromagnetic state in several transitionmetal interfaces.

#### *2.5. Chiral Magnetic Interactions Beyond DMI*

Some high-order terms containing cross products of spins may also be necessary for fitting the models to the total energy or explaining some certain magnetic properties. Due to the chiral properties of these interactions like DMI, it is possible that they can also lead to, or explain, intriguing noncollinear spin textures such as skyrmions. In Ref. [128], topological–chiral interactions are found to be very prominent in MnGe, which includes chiral–chiral interactions (CCI) with the form

$$\kappa\_{ijkl'j'k'}^{\mathbf{CC}} \left[ \mathbf{S}\_i \cdot (\mathbf{S}\_j \times \mathbf{S}\_k) \right] \left[ \mathbf{S}\_{i'} \cdot \left( \mathbf{S}\_{j'} \times \mathbf{S}\_{k'} \right) \right] \tag{11}$$

whose local part has the form

$$\kappa\_{ijk}^{\rm CC} \left[ \mathbf{S}\_i \cdot \left( \mathbf{S}\_j \times \mathbf{S}\_k \right) \right]^2 \tag{12}$$

and spin–chiral interactions (SCI) with the form

$$\kappa\_{ijk}^{\mathbf{SC}} \left( \mathbf{r}\_{ijk} \cdot \mathbf{S}\_i \right) \left[ \mathbf{S}\_i \cdot \left( \mathbf{S}\_j \times \mathbf{S}\_k \right) \right] \tag{13}$$

where the unit vector *τijk* ∝ *Rj* − *Ri* × (*Rk* − *Ri*) is the surface normal of the oriented triangle spanned by the lattice sites *Ri* , *Rj* and *Rk* . The local scalar spin chirality *χijk* = *Si* · *Sj* × *Sk* among triplets of spins can be interpreted as a fictitious effective magnetic field *B e*ff ∝ *χijkτijk*, which leads to topological orbital moments *L TO* (TOM) [129–134] arising from the orbital current of electrons hopping around the triangles. The TOM is defined as

$$\mathbf{L}\_{i}^{TO} = \sum\_{(jk)} \mathbf{L}\_{ijk}^{TO} = \sum\_{(jk)} \kappa\_{ijk}^{TO} \chi\_{ijk} \mathbf{\sigma}\_{ijk} \tag{14}$$

where *κ TO ijk* is the local topological orbital susceptibility. CCI corresponds to the interaction between pairs of topological orbital currents (or TOMs), whose local part can be interpreted as the orbital Zeeman interaction *L TO i* · *B* eff. SCI arises from the SOC, which couples the where

TOM to single spin magnetic moments. An illustration of CCI and SCI is provided in Figure 2. Considerations of CCI and SCI improved the fitting of the total energy in MnGe (see details in Ref. [128]). Moreover, the authors showed the possibility that the CCI may lead to three-dimensional topological spin states and therefore may be vital in deciding the ground state of the spin configurations of MnGe, which was found to be a threedimensional topological lattice (possibly built up with hedgehogs and anti-hedgehogs) experimentally [135]. terpreted as the orbital Zeeman interaction ்ை ⋅ ୣ. SCI arises from the SOC, which couples the TOM to single spin magnetic moments. An illustration of CCI and SCI is provided in Figure 2. Considerations of CCI and SCI improved the fitting of the total energy in MnGe (see details in Ref. [128]). Moreover, the authors showed the possibility that the CCI may lead to three-dimensional topological spin states and therefore may be vital in deciding the ground state of the spin configurations of MnGe, which was found to be a threedimensional topological lattice (possibly built up with hedgehogs and anti-hedgehogs) experimentally [135].

tion between pairs of topological orbital currents (or TOMs), whose local part can be in-

*Molecules* **2021**, *26*, x FOR PEER REVIEW 10 of 26

்ை =

()

்ை

=

()

்ை is the local topological orbital susceptibility. CCI corresponds to the interac-

்ை

(14)

**Figure 2.** Schematic diagrams of chiral–chiral interactions (CCI) and spin–chiral interactions (SCI), as provided by S. Grytsiuk et al. in Ref. [128]. Spins and topological orbital moments (TOMs) are denoted as black arrows and blue arrows, respectively. CCI can be regarded as interactions between TOMs, while SCI can be interpreted as interactions between TOM and local spins. **Figure 2.** Schematic diagrams of chiral–chiral interactions (CCI) and spin–chiral interactions (SCI), as provided by S. Grytsiuk et al. in Ref. [128]. Spins and topological orbital moments (TOMs) are denoted as black arrows and blue arrows, respectively. CCI can be regarded as interactions between TOMs, while SCI can be interpreted as interactions between TOM and local spins.

A new type of chiral pair interaction ⋅ ( × ) ( ⋅ ), named as chiral biquadratic interaction (CBI), which is the biquadratic equivalent of the DMI, was derived from a microscopic model and demonstrated to be comparable in magnitude to the DMI in magnetic dimers made of 3*d* elements on Pt(111), Pt(001), Ir(111) and Re(0001) surface with strong SOC [136]. Similar but generalized chiral interactions such as ⋅ ( × ) ( ⋅ ) and ⋅ ( × ) ( ⋅ ), named four-spin chiral interactions, were discussed in Ref. [137], and they are found to be important in predicting a correct chirality for a spin spiral state of Fe chains deposited on the Re(0001) surface. A new type of chiral pair interaction *Cij* · *Si* × *Sj Si* · *Sj* , named as chiral biquadratic interaction (CBI), which is the biquadratic equivalent of the DMI, was derived from a microscopic model and demonstrated to be comparable in magnitude to the DMI in magnetic dimers made of 3*d* elements on Pt(111), Pt(001), Ir(111) and Re(0001) surface with strong SOC [136]. Similar but generalized chiral interactions such as *Dijjk* · *Sj* × *Sk Si* · *Sj* and *Dijkl* · (*Sk* × *Sl*) *Si* · *Sj* , named four-spin chiral interactions, were discussed in Ref. [137], and they are found to be important in predicting a correct chirality for a spin spiral state of Fe chains deposited on the Re(0001) surface.

When there is a magnetic field, for a nonbipartite lattice, the magnetic field can couple with the spin and produce a new term of the form ⋅ ൫ × ൯= [138,139], which can be termed the three-spin chiral interaction (TCI) [140]. Such a chiral term can induce a gapless line in frustrated spin-gapped phases, and a critical chiral strength can change the ground state from spiral to Néel quasi-long-range-order phase [138]. This chiral term is also found to produce a chiral spin liquid state [141], where the time-reversal symmetry is broken spontaneously by the emergence of long-range order of scalar chirality [142]. When there is a magnetic field, for a nonbipartite lattice, the magnetic field can couple with the spin and produce a new term of the form *JijkSi* · *Sj* × *Sk* = *Jijkχijk* [138,139], which can be termed the three-spin chiral interaction (TCI) [140]. Such a chiral term can induce a gapless line in frustrated spin-gapped phases, and a critical chiral strength can change the ground state from spiral to Néel quasi-long-range-order phase [138]. This chiral term is also found to produce a chiral spin liquid state [141], where the time-reversal symmetry is broken spontaneously by the emergence of long-range order of scalar chirality [142].

#### *2.6. Expansions of Magnetic Interactions*

In general, a complete basis can be used to expand the spin interactions. One example is the spin-cluster expansion (SCE) [143–145], where unit vectors denoting the directions of spins are used as independent variables and spherical harmonic functions are used in the expressions of the basis functions. When SOC and the external magnetic field are not important, as mentioned in Section 2.3, only inner products of spins need to be considered. Consequently, the expansion can be

$$H\_{\rm spin} = E\_0 + \sum\_{n=1}^{\infty} J\_n \left( \sum\_{\substack{\text{th nearest pairs } \langle k,l \rangle}}^{\prime} \mathbf{e}\_k \cdot \mathbf{e}\_l \right) + \sum\_{n'=1}^{\infty} J\_{n'}' \sum\_{\substack{\text{tr} \mathbf{l}' \text{ with type}}} (\mathbf{e}\_k \cdot \mathbf{e}\_l)(\mathbf{e}\_m \cdot \mathbf{e}\_{\mathbf{l}}) + \cdots \tag{15}$$

as used in Ref. [125]. When using expansions of spin vectors (or directions of the spins), suitable truncations on interaction distances and interaction orders are needed. Otherwise, the number of terms would be infinite, and as a result, the problem would be unsolvable.

#### **3. Methods of Computing the Parameters of Effective Spin Hamiltonian Models**

In this part, we mainly discuss the methods of computing the spin interaction parameters based on first-principles calculations of crystals, where periodic boundary conditions are tacitly supposed. These methods include different kinds of energy-mapping analysis (see Section 3.1) and Green's function method based on magnetic-force linear response theory (see Section 3.2). Discussions on the rigid spin rotation approximation and other assumptions are provided in Section 3.3. The cases of clusters, where periodic boundary conditions do not exist, will be briefly discussed in Section 3.3. Methods of obtaining spin interaction parameters from experiments will also be briefly mentioned in Section 3.3.

#### *3.1. Energy-Mapping Analysis*

In an energy-mapping analysis, we do several first-principles calculations to assess the total energies of different spin configurations. Then, we use the effective spin Hamiltonian to provide the expressions of the total energies of these spin configurations (with the expressions containing several undetermined parameters). By mapping the total energy expressions given by the effective spin Hamiltonian model to the results of first-principles calculations, the values of the undetermined parameters can be estimated. There are several types of energy-mapping analysis. For the first type, a minimal number of configurations are used, and a concrete expression for calculating the parameters can be given in advance. An example is by mapping between the eigenvalues and eigenfunctions of exact Hamiltonians and the effective spin Hamiltonian models (typically the Heisenberg model) to estimate the exchange parameters for relatively simple systems [146,147]. Several broken symmetry (BS) approaches are also of this type, where broken-symmetry states (instead of eigenstates of exact Hamiltonians) are adopted for energy mapping between the models and results of first-principles calculations [147,148]. A typical example of BS approach is the four-state method [148,149] where four special states are chosen for calculating each component of the parameters, which will be introduced in Section 3.1.1. For the second type, more configurations are used, and the parameters in the supposed effective spin Hamiltonian model are determined by employing least-squares fitting, which will be introduced in Section 3.1.2. The third type is similar to the second one, but the concrete form of the effective spin Hamiltonian model is not determined in advance. In the beginning, one includes many terms in the mode Hamiltonians. The relevance of each individual term depends on the fitting performance with respect to first-principles calculations. While selecting the relevant terms for a model Hamiltonian, it is important to search for the minimal Hamiltonian for a given magnetic system, namely, the one with the minimal number of parameters that capture its essential physics. This type of energy-mapping analysis will be introduced in Section 3.1.3.

In this section, we will mainly focus on the applications in solid state systems with periodic boundary conditions. The total energy of a designated configuration (which is usually a broken-symmetry state) is typically provided by first-principles calculations (e.g., DFT+U calculations) with constrained directions of magnetic moments.

#### 3.1.1. Four-State Method

The energy-mapping analysis based on four ordered spin states [148,149], also referred to as the four-state method, assumes the effective spin Hamiltonian include only second-

order terms (i.e., isotropic Heisenberg terms, DMI terms, Kitaev terms, and SIA terms). Each component of the parameters like *Jij*, *Dij*,*x* , A*i*,*xy*, (A*i*,*yy* − A*i*,*xx*) and (A*i*,*zz* − A*i*,*xx*) can be obtained by first-principles calculations for four specified spin states [148]. Taking the isotropic Heisenberg parameter *Jij* for example, with the spin-orbit coupling (SOC) switched off during the first-principles calculations, we use *Eij*,*αβ* (*α*, *β* =↑, ↓ ) to denote the energy of the configuration where spin *i* is parallel or antiparallel to the *z* direction (if *α* =↑ or ↓ , respectively), spin *j* is parallel or antiparallel to the *z* direction (if *β* = ↑ or ↓ , respectively), and all the spins except *i* and *j* are kept unchanged in the four states (which will be referred to as the "reference configuration", usually chosen to be a low-energy collinear state). Then *Jij* can be expressed as ferred to as the four-state method, assumes the effective spin Hamiltonian include only second-order terms (i.e., isotropic Heisenberg terms, DMI terms, Kitaev terms, and SIA terms). Each component of the parameters like , ,௫, ,௫௬, (,௬௬ − ,௫௫) and (,௭௭ − ,௫௫) can be obtained by first-principles calculations for four specified spin states [148]. Taking the isotropic Heisenberg parameter for example, with the spin-orbit coupling (SOC) switched off during the first-principles calculations, we use ,ఈఉ (, =↑, ↓) to denote the energy of the configuration where spin *i* is parallel or antiparallel to the *z* direction (if =↑ or ↓, respectively), spin *j* is parallel or antiparallel to the *z* direction (if = ↑ or ↓, respectively), and all the spins except *i* and *j* are kept unchanged in the four states (which will be referred to as the "reference configuration", usually chosen to be a low-energy collinear state). Then can be expressed as

In this section, we will mainly focus on the applications in solid state systems with periodic boundary conditions. The total energy of a designated configuration (which is usually a broken-symmetry state) is typically provided by first-principles calculations

The energy-mapping analysis based on four ordered spin states [148,149], also re-

(e.g., DFT+U calculations) with constrained directions of magnetic moments.

$$J\_{i\dot{j}} = \frac{E\_{i\dot{j},\uparrow\uparrow} + E\_{i\dot{j},\downarrow\downarrow} - E\_{i\dot{j},\uparrow\downarrow} - E\_{i\dot{j},\downarrow\uparrow}}{4S^2} \tag{16}$$

The schematic diagrams of these four states are shown in Figure 3. The schematic diagrams of these four states are shown in Figure 3.

*Molecules* **2021**, *26*, x FOR PEER REVIEW 12 of 26

3.1.1. Four-State Method

**Figure 3.** Schematic diagrams of the four states adopted in the four-state method for computing . Except atoms *i* and *j*, all the other atoms are omitted in this diagram (which are kept the same in the four states). Spins pointing up and down are indicated with orange and blue balls, respectively. The total energies given by first-principles calculations corresponding to these four states are denoted as (**a**) ,↑↑, (**b**) ,↑↓, (**c**) ,↓↑, and (**d**) ,↓↓, respectively. **Figure 3.** Schematic diagrams of the four states adopted in the four-state method for computing *Jij*. Except atoms *i* and *j*, all the other atoms are omitted in this diagram (which are kept the same in the four states). Spins pointing up and down are indicated with orange and blue balls, respectively. The total energies given by first-principles calculations corresponding to these four states are denoted as (**a**) *Eij*,↑↑, (**b**) *Eij*,↑↓, (**c**) *Eij*,↓↑, and (**d**) *Eij*,↓↓, respectively.

In general, the second-order effective spin Hamiltonian takes the form of In general, the second-order effective spin Hamiltonian takes the form of

$$H\_{\rm spin} = \sum\_{i,j>i} \mathbf{S}\_i^T \mathbb{J}\_{ij} \mathbf{S}\_j + \sum\_i \mathbf{S}\_i^T \mathbb{A}\_i \mathbf{S}\_i \tag{17}$$

and each component of the matrix and can also be obtained by first-principles calculations for four specified spin states. To compute , (, = , , ), we use ,,ఈఉ and each component of the matrix J*ij* and A*i* can also be obtained by first-principles calculations for four specified spin states. To compute J*ij*,*ab* (*a*, *b* = *x*, *y*, *z*), we use *Eij*,*ab*,*αβ* ( *α*, *β* = ↑, ↓ ) to denote the energy of the configuration where spin *i* is parallel or antiparallel to the *a* direction (if *α* = ↑ or ↓ , respectively), spin *j* is parallel or antiparallel to the *b* direction (if *β* = ↑ or ↓ , respectively), and all the spins except for *i* and *j* are kept unchanged and parallel to the *c*-axis (*c* = *x*, *y* or *z*, *c* 6= *a*, *c* 6= *b*) with an appropriate reference configuration and kept unchanged. Then, J*ij*,*ab* can be expressed as

$$\mathbb{L}\_{ij,ab} = \frac{E\_{ij,ab,\uparrow\uparrow} + E\_{ij,ab,\downarrow\downarrow} - E\_{ij,ab,\uparrow\downarrow} - E\_{ij,ab,\downarrow\uparrow}}{\mathbf{4S}^2} \tag{18}$$

To compute A*i*,*ab* (*a*, *b* = *x*, *y*, *z* with *a* 6= *b*), we use *Ei*,*ab*,*αβ* (*α*, *β* = ↑, ↓ ) to denote the energy of the configuration where spin *i* is parallel to the direction whose *a* component is ± √ 2 2 (for *α* = ↑ or ↓ , respectively), *b* component is ± √ 2 2 (for *β* = ↑ or ↓ , respectively), and the other component is 0. Here all the spins except for *i* and *j* are parallel to the *c*-axis (*c* = *x*, *y* or *z*, *c* 6= *a*, *c* 6= *b*) with an appropriate reference configuration. Then, A*i*,*ab* can be expressed as

$$\mathbb{A}\_{i,ab} = \frac{E\_{i,ab,\uparrow\uparrow} + E\_{i,ab,\downarrow\downarrow} - E\_{i,ab,\uparrow\downarrow} - E\_{i,ab,\downarrow\uparrow}}{4S^2} \tag{19}$$

To compute (A*i*,*aa* − A*i*,*bb*) (*a*, *b* = *x*, *y*, *z* with *a* 6= *b*), we use *Ei*,*ab*,*αβ* ( *α*, *β* = ↑, ↓ or 0) to denote the energy of the configuration where spin *i* is parallel to the direction whose *a* component is ±1 or 0 (if *α* = ↑, ↓ or 0, respectively), *b* component is ±1 or 0 (if *β* = ↑, ↓ or 0, respectively) and the other component is 0, while all the spins except *i* are parallel to the *c*-axis (*c* = *x*, *y* or *z*, *c* 6= *a*, *c* 6= *b*) with an appropriate reference configuration. Then (A*i*,*aa* − A*i*,*bb*) can be expressed as

$$\mathbb{A}\_{i,aa} - \mathbb{A}\_{i,bb} = \frac{E\_{i,ab,\uparrow0} + E\_{i,ab,\downarrow0} - E\_{i,ab,0\uparrow} - E\_{i,ab,0\downarrow}}{4S^2} \tag{20}$$

It is easy to verify that, by employing this four-state method, each component of the J*ij* and A*i* can be obtained with the effects of other second-order terms entirely cancelled. Now we take the effects of fourth-order terms (without SOC) into account and check if the algorithms for computing J*ij* and A*i* are still rigorous. For A*i* , we can find out that the effects of all these terms are correctly cancelled. For J*ij*, the effects of biquadratic terms, three-body fourth-order terms, and most of the four-spin ring coupling terms are perfectly cancelled, while only the terms like *Si* · *Sj* (*Sk* · *Sl*) (*k*, *l* 6= *i*, *j*) will interfere with the calculation of *Jij* because (*Sk* · *Sl*) is constant during the calculation and therefore mixed with the contribution of *Jij Si* · *Sj* . The error of the calculated *Jij* originated from four-spin ring coupling terms

$$\sum\_{\substack{i,j>i,\\k>i,\\|k>i,\\l>k}} \mathbb{K}\_{ijkl} (\mathbf{S}\_i \cdot \mathbf{S}\_j)(\mathbf{S}\_k \cdot \mathbf{S}\_l) \tag{21}$$

is given by

$$\sum\_{\substack{k \neq i \text{ or } j \\ l \neq i \text{ or } j \\ (l > k)}} \mathbb{K}\_{ijkl}(\mathbf{S}\_k \cdot \mathbf{S}\_l) \tag{22}$$

while there is no easy way to get rid of this problem perfectly. Other parts of the J*ij* (including D*ij* and K*ij*) are not affected by these fourth-order terms (without SOC), but by other types of fourth-order terms (like four-spin chiral interactions) if SOC is taken into account (because of the similar reason).

In Ref. [122], the four-spin ring coupling interaction is found to be important in TbMnO3, and therefore leads to instability in calculating the Heisenberg parameters using the four-state method when changing the reference configurations. This problem is remedied by calculating the Heisenberg parameter *Jij* twice with the four-state method using the FM and A-type AFM (A-AFM, see Figure 4c) as the reference configurations, and use their average value as the final estimation of *Jij*. The parameter of the vital ring coupling interaction is obtained by calculating the difference between the two calculated *Jij* values (with FM and A-AFM reference configurations, respectively). This effective remedy is based on the assumption that only one kind of ring coupling interaction is essential. However, if there are more kinds of significant ring coupling or if we do not know which ring coupling is essential in advance, such a method of calculating *Jij* is still not very trustworthy. Nevertheless, it is found that by taking the average of the calculated *Jij* with four-state

method with FM and G-type AFM (G-AFM, see Figure 4a) reference configurations, the influences of *Kijkl Si* · *Sj* (*Sk* · *Sl*) with *k* and *l* being nearest pairs are eliminated. The terms *Kijkl Si* · *Sj* (*Sk* · *Sl*) with non-nearest pairs of *k* and *l* still interfere with the calculation of *Jij*, but they are generally very weak. Therefore, such a remedy to calculate *Jij* using the four-state method should work well in most cases. In cases when a G-AFM state (in which all the nearest pairs of spins are antiparallelly aligned) cannot be defined (e.g., a triangular or a Kagomé lattice), there may be more than two reference configurations to use in the four-state method so as to eliminate the effects of *Kijkl Si* · *Sj* (*Sk* · *Sl*) terms with non-nearest pairs of *k* and *l*. These reference configurations need to be designed carefully according to the specific circumstances to get rid of the effects of *Kijkl Si* · *Sj* (*Sk* · *Sl*) terms as much as possible. state method with FM and G-type AFM (G-AFM, see Figure 4a) reference configurations, the influences of ( ⋅ )( ⋅ ) with *k* and *l* being nearest pairs are eliminated. The terms ( ⋅ )( ⋅ ) with non-nearest pairs of *k* and *l* still interfere with the calculation of , but they are generally very weak. Therefore, such a remedy to calculate using the four-state method should work well in most cases. In cases when a G-AFM state (in which all the nearest pairs of spins are antiparallelly aligned) cannot be defined (e.g., a triangular or a Kagomé lattice), there may be more than two reference configurations to use in the four-state method so as to eliminate the effects of ( ⋅ )( ⋅ ) terms with non-nearest pairs of *k* and *l*. These reference configurations need to be designed carefully according to the specific circumstances to get rid of the effects of ( ⋅ )( ⋅ ) terms as much as possible.

the FM and A-type AFM (A-AFM, see Figure 4c) as the reference configurations, and use their average value as the final estimation of . The parameter of the vital ring coupling interaction is obtained by calculating the difference between the two calculated values (with FM and A-AFM reference configurations, respectively). This effective remedy is based on the assumption that only one kind of ring coupling interaction is essential. However, if there are more kinds of significant ring coupling or if we do not know which ring coupling is essential in advance, such a method of calculating is still not very trustworthy. Nevertheless, it is found that by taking the average of the calculated with four-

*Molecules* **2021**, *26*, x FOR PEER REVIEW 14 of 26

**Figure 4.** Schematic diagrams of (**a**) G-type AFM, (**b**) C-type AFM, and (**c**) C-type AFM states. Spins pointing to two opposite directions (e.g., up and down) are indicated with orange and blue balls, respectively. **Figure 4.** Schematic diagrams of (**a**) G-type AFM, (**b**) C-type AFM, and (**c**) C-type AFM states. Spins pointing to two opposite directions (e.g., up and down) are indicated with orange and blue balls, respectively.

The main advantages of the four-state method are its relatively small amount of firstprinciples calculations and its relatively good cancellations of other relevant terms. A weakness is that it cannot analyze the uncertainties of the parameters, so that we do not know how precise those estimated values are. Another weakness, which is also shared with other energy-mapping analysis methods, is that the computed between and is actually the sum over ᇲ with any lattice vector =ᇲ − . Therefore, to get rid The main advantages of the four-state method are its relatively small amount of first-principles calculations and its relatively good cancellations of other relevant terms. A weakness is that it cannot analyze the uncertainties of the parameters, so that we do not know how precise those estimated values are. Another weakness, which is also shared with other energy-mapping analysis methods, is that the computed J*ij* between *Si* and *Sj* is actually the sum over J*ij*0 with any lattice vector *R* = *rj* 0 − *rj* . Therefore, to get rid of the effects of other spin pairs, a relatively large supercell is needed.

of the effects of other spin pairs, a relatively large supercell is needed. The four-state method can also be generalized to compute biquadratic parameters, where the SOC needs to be switched off during the first-principles calculations. For calculating in the term ( ⋅ )ଶ, we can let pointing to the (1,0,0) direction, pointing to the (1,0,0), (−1,0,0), ቀ ଵ √ଶ , ଵ √ଶ , 0ቁ and ቀ− ଵ √ଶ , − ଵ √ଶ , 0ቁ directions, with other spins parallel to the *z*-axis. The corresponding total energies are denoted as ଵ, ଶ, ଷ, and ସ, respectively. Then the can be expressed as The four-state method can also be generalized to compute biquadratic parameters, where the SOC needs to be switched off during the first-principles calculations. For calculating *Kij* in the term *Kij Si* · *Sj* 2 , we can let *Si* pointing to the (1, 0, 0) direction, *Sj* pointing to the (1, 0, 0), (−1, 0, 0), √ 1 2 , √ 1 2 , 0 and − √ 1 2 , − √ 1 2 , 0 directions, with other spins parallel to the *z*-axis. The corresponding total energies are denoted as *E*1, *E*2, *E*3, and *E*4, respectively. Then the *Kij* can be expressed as

$$K\_{ij} = E\_1 + E\_2 - E\_3 - E\_4 \tag{23}$$

It can be easily checked that the effects of other terms not higher than fourth order are totally eliminated. Therefore, this approach of calculating should be relatively rigorous theoretically. It can be easily checked that the effects of other terms not higher than fourth order are totally eliminated. Therefore, this approach of calculating *Kij* should be relatively rigorous theoretically.

Note that the four-state methods [148,149] could also give the derivatives of exchange interactions with respect to the atomic displacements without doing additional first-principles calculations due to the Hellmann-Feynman theorem. These derivatives are useful for the study of spin-lattice coupling related phenomena.

#### 3.1.2. Direct Least Squares Fitting

Another type of energy-mapping analysis, instead of the four-state method, uses more first-principles calculations with different spin configurations and fits the results to the effective spin Hamiltonian using the least-squares method to estimate the parameters [74,122–124]. Ways of choosing spin configurations can depend on which parameters to estimate.

In Ref. [74], for the four Mn sites in a unit cell of YMnO3, the polar and azimuthal angles (*θ*, ϕ) of their spins are given by (0, 0), (0, 0), *θ*, 3*π* 2 , and *θ*, *π* 2 , respectively. By changing the *θ* from 0◦ to 180◦ , different configurations are produced. If the effective spin Hamiltonian only contains Heisenberg terms, there will be a systematic deviation between the predicted value given by the effective Hamiltonian and the calculated value given by first-principles calculations. Such a deviation is well remedied by adding biquadratic exchange interactions into the effective Hamiltonian model. Thus, the biquadratic parameters can be fitted. Similar approaches were adopted by others to calculate and show the importance of biquadratic parameters and topological chiral–chiral contributions [122–124,128]. Apart from using an angular variable for generating spin configurations, using two or more variables is also practicable, or randomly chosen directions [106] can also be considered. Thus, more diverse configurations will be produced. The least-squares fitting will also work, but whether a systematic deviation exists will not be as apparent as the case when only one angular variable is used for generating configurations, and the fitting task may be more laborious.

The main virtues of this method are that the reliability of the model can be checked by the fitting performance and that the uncertainties of the parameters can be estimated if needed. This method is especially suitable for calculations of biquadratic parameters and topological CCI. However, when talking about calculations of Heisenberg parameters, this method needs more first-principles calculations and is thus less efficient. Furthermore, the fitted result of the Heisenberg parameter *Jij* is vulnerable to the effects of other fourthorder interactions such as terms as *Si* · *Sj* (*Si* · *Sk*) and *Si* · *Sj* (*Sk* · *Sl*). Therefore, the estimations of the Heisenberg parameters may not be very reliable if any of such fourthorder interactions are essential. This problem can be remedied by adding the related terms into the effective Hamiltonian model, while it is not easy to decide which terms to include in the model beforehand.

A possible way to get rid of the effects of other high-order terms is to perform artificial calculations where most of the magnetic ions are replaced by similar but nonmagnetic ions (e.g., substituting Fe3+ ions with nonmagnetic Al3+ ions) except for one or more ions to be studied [150]. For example, when calculating SIA, only one magnetic ion is not substituted, and by rotating this magnetic ion and calculating the total energy, the SIA can be studied. When studying two-body interactions between *Si* and *Sj* , only two magnetic ions (at site *i* and *j*) are not substituted, and by rotating the spins (or magnetic moments) of these two magnetic ions, the interactions between them can be studied. Such a technique of substituting atoms can be applied to energy mapping analysis based on either the four-state method or least-squares fitting. In this way, effects from interactions involving other sites are effectively avoided. Nevertheless, this substitution method can make the chemical environments of the remaining magnetic ions different from those in the system with no substitution. This may make the calculations of the interaction parameters untrustworthy.

#### 3.1.3. Methods Based on Expansions and Selecting Important Terms

The traditional energy-mapping analysis needs to construct an effective spin Hamiltonian first and then fit the undetermined parameters. However, it is not easy to give a perfect guess, especially when high-order interactions are essential. Such problems can be solved by considering almost all the possible terms utilizing some particular expansion with appropriate truncations. Usually, there are too many possible terms to be considered, so a direct fitting is impracticable; it requires at least as many first-principles calculations as the number of terms to determine, but leads to over-fitting problems due to too many parameters to determine. So, it is necessary to decide whether or not to include each term into the effective spin Hamiltonian on the basis of their contributions to the fitting performance.

In Ref. [145], SCE is adopted for the expansion of spin interactions of bcc and fcc Fe. After truncations based on the interacting distance and interaction orders, they considered 154 (179) possible different interaction terms in bcc (fcc) Fe. They randomly generated 3954 (835) different spin configurations in a 2 × 2 × 2 supercell for fitting. Their method of choosing terms is as follows: starting from the effective Hamiltonian model with only a constant term, try adding each possible term into the temporary model and accept the one providing the best fitting performance, in which way the terms are added to the model one by one. This method is the forward selection in variable selection problems, which is simple and straightforward, and it works well in most cases. However, this method may include unnecessary interactions to the effect Hamiltonian.

In Ref. [125], a machine learning method for constructing Hamiltonian (MLMCH) is proposed, which is more efficient and more reliable than the traditional forward selection method. Firstly, a testing set is used to avoid over-fitting problems. Secondly, not only adding terms but also deleting and substituting terms are considered during the search for the appropriate model. Thus, if an added term is judged to be unnecessary, it can still be removed from the model afterward. A penalty factor *p ˘* (*λ* ≥ 1), where *p* is the number of parameters in the temporary model and *λ* is a given parameter, is used together with the loss function *σ* 2 (the fitting variance) to select models with fewer parameters. Several techniques are used to reduce the search space and enhance the search efficiency to select important terms out of tens of thousands of possible terms. The flow charts of this method of variable selection as well as the forward selection method are shown in Figure 5. *Molecules* **2021**, *26*, x FOR PEER REVIEW 17 of 26

**Figure 5.** Flow charts of several methods of variable selection. In these flow charts, the fitting variance ଶ estimated by the training set and the testing set are denoted as ଶ and ᇱଶ, respectively. (**a**) The traditional forward selection method as adopted in Ref. [145]. (**b**) The forward selection method using a testing set to check if over-fitting problems occur. (**c**) A simplified flow chart of the algorithm used in MLMCH [125], where some details are omitted. A testing set is used to check if over-fitting problems occur. The criterion for a better model is a smaller ൫ଶ ⋅ ൯ with ≥1. There are ఒ values of , which are set in advance, saved in the array "lambdas(1: ఒ)", whose components are usually arranged in descending order of magnitude. For each value of , the best model is selected by using the criterion ൫ଶ ⋅ ൯. This method is advantageous in two ways: (a) Constructing the effective spin Ham-**Figure 5.** Flow charts of several methods of variable selection. In these flow charts, the fitting variance *σ* 2 estimated by the training set and the testing set are denoted as ˆ*σ* 2 and ˆ *σ* 02 , respectively. (**a**) The traditional forward selection method as adopted in Ref. [145]. (**b**) The forward selection method using a testing set to check if over-fitting problems occur. (**c**) A simplified flow chart of the algorithm used in MLMCH [125], where some details are omitted. A testing set is used to check if over-fitting problems occur. The criterion for a better model is a smaller ˆ*σ* 2 · *λ p* with *λ* ≥ 1. There are *Nλ* values of *λ*, which are set in advance, saved in the array "lambdas(1: *Nλ*)", whose components are usually arranged in descending order of magnitude. For each value of *λ*, the best model is selected by using the criterion ˆ*σ* 2 · *λ p* .

have some room for improvement.

calized orbitals. By defining

interaction terms; (b) this method is general, so it can be applied to most magnetic materials. The least-squares fitting needed in this approach can also provide the estimations for the uncertainties of the parameters. The flaw is that it needs lots of (typically hundreds of) first-principles calculations, which could be impracticable when a very large supercell is needed (especially when the material is metallic so that long-range interactions are essential). The way to generate spin configurations (typically randomly distributed among all possible directions, sometimes deviating only moderately from the ground state) may

In Green's function method based on magnetic-force linear response theory [151– 159], we need localized basis functions ఙ() (, , indicating the site, orbital, and spin indices, respectively) based on the tight-binding model. The localized basis functions can be provided by DFT codes together with Wannier90 [160,161] or codes based on lo-

*3.2. Green's Function Method Based on Magnetic-Force Linear Response Theory*

This method is advantageous in two ways: (a) Constructing the effective spin Hamiltonians is carried out comprehensively, which makes it less likely to miss some critical interaction terms; (b) this method is general, so it can be applied to most magnetic materials. The least-squares fitting needed in this approach can also provide the estimations for the uncertainties of the parameters. The flaw is that it needs lots of (typically hundreds of) first-principles calculations, which could be impracticable when a very large supercell is needed (especially when the material is metallic so that long-range interactions are essential). The way to generate spin configurations (typically randomly distributed among all possible directions, sometimes deviating only moderately from the ground state) may have some room for improvement.

#### *3.2. Green's Function Method Based on Magnetic-Force Linear Response Theory*

In Green's function method based on magnetic-force linear response theory [151–159], we need localized basis functions *ψimσ*(*r*) (*i*, *m*, *σ* indicating the site, orbital, and spin indices, respectively) based on the tight-binding model. The localized basis functions can be provided by DFT codes together with Wannier90 [160,161] or codes based on localized orbitals. By defining

$$\mathbb{H}\_{imjm'\sigma\sigma'}(\mathbf{R}) = \langle \psi\_{im\sigma}(\mathbf{r}) | H | \psi\_{im\sigma}(\mathbf{r} + \mathbf{R}) \rangle \tag{24}$$

$$\mathbb{S}\_{imjm'\sigma\sigma'}(\mathbf{R}) = \langle \psi\_{im\sigma}(\mathbf{r}) | \psi\_{im\sigma}(\mathbf{r} + \mathbf{R}) \tag{25}$$

$$\mathbb{H}(\mathbf{k}) = \sum\_{\mathbf{R}} \mathbb{H}(\mathbf{R}) e^{i\mathbf{k} \cdot \mathbf{R}} \tag{26}$$

$$\mathbb{E}\mathbb{S}(\mathbf{k}) = \sum\_{\mathbf{R}} \mathbb{S}(\mathbf{R}) e^{i\mathbf{k}\cdot\mathbf{R}} \tag{27}$$

the Green's function in reciprocal space and real space are defined as

$$\mathbb{G}(\mathbf{k},\varepsilon) = \left(\varepsilon \mathbb{G}(\mathbf{k}) - \mathbb{H}(\mathbf{k})\right)^{-1} \tag{28}$$

and

$$\mathbb{G}(\mathbf{R}, \varepsilon) = \int\_{\text{BZ}} \mathbb{G}(\mathbf{k}, \varepsilon) e^{-i\mathbf{k} \cdot \mathbf{R}} \mathbf{d}k \tag{29}$$

Based on the magnetic force theorem [162], the total energy variation due to a perturbation (which is the rotation of spins in this case) from the ground state equals the change of single-particle energies at the fixed ground-state potential:

$$
\delta E = \int\_{-\infty}^{E\_F} \varepsilon \delta n(\varepsilon) \mathbf{d}\varepsilon = -\int\_{-\infty}^{E\_F} \delta \mathbf{N}(\varepsilon) \mathbf{d}\varepsilon \tag{30}
$$

where

$$m(\varepsilon) = -\frac{1}{\pi}Im\text{Tr}(\mathbb{G}(\varepsilon))\tag{31}$$

and

$$N(\varepsilon) = -\frac{1}{\pi}Im\text{Tr}(\varepsilon - \mathbb{H})\tag{32}$$

where traces are taken over orbitals. By defining

$$\mathbb{P}\_{\text{i}} = \mathbb{H}\_{\text{ii}}(\mathbf{R} = \mathbf{0}) = \mathbb{I}\_{\text{i}}^{0} \Leftarrow + \stackrel{\rightarrow}{\mid}\_{\text{i}} \cdot \stackrel{\rightarrow}{\sigma} \tag{33}$$

with its component

$$\mathbb{P}\_{imm'} = p\_{imm'}^0 \prec + p\_{imm'} \overset{\rightarrow}{e}\_{imm'} \cdot \overset{\rightarrow}{\sigma} \tag{34}$$

where → *σ* is the vector composed of Pauli matrices. By defining

$$\mathbf{G}\_{im,jm'} = \mathbf{G}\_{im,jm'}^{0} \prec + \stackrel{\rightarrow}{\mathbf{G}}\_{im,jm'} \cdot \stackrel{\rightarrow}{\sigma} \tag{35}$$

the energy variation due to the two-spin interaction between sites *i* and *j* is

$$\delta E\_{\rm ij} = -\frac{2}{\pi} \int\_{\infty}^{E\_F} Im \text{Tr} \{ \delta \mathbb{H}\_i \mathbb{G} \delta \mathbb{H}\_j \mathbb{G} \} \, \text{d}\varepsilon \tag{36}$$

with *δ*H*i* = *δφi* × *pi* . After mathematical simplification, the expression of *δEij* can be mapped to that given by the effective Hamiltonian model

$$H\_{\rm spin} = \sum\_{i} \mathbf{S}\_{i}^{T} \mathbf{A}\_{i} \mathbf{S}\_{i} + \sum\_{i,j>i} \left[ I\_{\rm ij} \mathbf{S}\_{i} \cdot \mathbf{S}\_{j} + \mathbf{D}\_{\rm ij} \cdot \left( \mathbf{S}\_{i} \times \mathbf{S}\_{j} \right) + \mathbf{S}\_{i}^{T} \mathbb{K}\_{i\uparrow} \mathbf{S}\_{j} \right] + \sum\_{i,j>i} \mathbf{K}\_{ij} \left( \mathbf{S}\_{i} \cdot \mathbf{S}\_{j} \right)^{2} \tag{37}$$

(including all the second-order terms and a biquadratic term) to obtain the expressions of the parameters:

$$J\_{\vec{ij}} = \operatorname{Im} \left( A\_{\vec{ij}}^{00} - A\_{\vec{ij}}^{\text{xx}} - A\_{\vec{ij}}^{yy} - A\_{\vec{ij}}^{zz} - \mathfrak{A} A\_{\vec{ij}}^{zz} \mathbf{S}\_{\vec{i}}^{\text{ref}} \cdot \mathbf{S}\_{\vec{j}}^{\text{ref}} \right) \tag{38}$$

$$\mathbb{K}\_{\rm ij}^{\mu\upsilon} = \operatorname{Im} \left( A\_{\rm ij}^{\mu\upsilon} + A\_{\rm ij}^{\mu\iota} \right) \tag{39}$$

$$D\_{ij}^{\mu} = \text{Re}\left(A\_{ij}^{0\mu} - A\_{ij}^{\mu 0}\right) \tag{40}$$

$$B\_{i\bar{j}} = \operatorname{Im} \left( A\_{i\bar{j}}^{z\bar{z}} \right) \tag{41}$$

where

$$A\_{ij}^{\mu\upsilon} = \frac{1}{\pi} \int\_{-\infty}^{\text{E}\_{\text{F}}} \text{Tr} \left\{ \mathbf{i}\_{i}^{\text{Z}} \mathbf{G}\_{ij}^{\text{u}} \mathbf{j}\_{j}^{\text{Z}} \mathbf{G}\_{ji}^{\upsilon} \right\} \text{d}\varepsilon \tag{42}$$

with *u*, *v* ∈ {0, *x*, *y*, *z*} (the trace is also taken over orbitals), and *S* ref *i* indicates the unperturbed vector *Si* . An *xyz* average strategy can be adopted so that some components inaccessible from one first-principles calculation can be obtained [157].

The main advantages of this method are that it only requires one or three DFT calculations to obtain all the parameters of second-order terms and biquadratic terms between different atoms, using only a small supercell (with a dense enough *k*-point sampling) to obtain interaction parameters between spins far away from each other. Therefore, it saves the computational cost compared to the energy-mapping analysis, especially when long-range interactions are essential. It also avoids the difficulties of reaching self-consistent-field convergence in DFT calculations for high-energy configurations, which may occur in the energy-mapping analysis. This method is good at describing states near the ground state but may not be so good at describing high-energy states. A limitation is that this method cannot obtain SIA parameters, and its calculations for biquadratic parameters are not very trustworthy [157]. The calculated Heisenberg parameter *Jij* is mixed with the contributions of other fourth-order interactions such as terms like *Si* · *Sj* (*Si* · *Sk*) and *Si* · *Sj* (*Sk* · *Sl*). Therefore, the results may be unreliable if any of such fourth-order interactions are essential. Another little flaw is the noise of a typical order of magnitude of a few *µ*eV introduced by the process of obtaining the Wannier orbitals [157]. In addition, the uncertainties of the parameters cannot be obtained by this method.

A recent study [140] considered the rotations of more than two spins as the perturbation and mapped the *δE* to the corresponding quantity given by an effective Hamiltonian with more types of interactions, including terms proportional to *Si* · *Sj* × *Sk* , *Si* · *Sj* (*Sk* · *Sl*) and *Si* × *Sj* (*Sk* · *Sl*). This enables one to obtain the expressions needed for calculating the associated parameters. The derivations and forms of the expressions are much more complicated compared with those from the second-order interactions discussed above. This generalization of the traditional approach for calculating second-order interaction parameters remedied the problem of ignoring the effects of other high-order interactions to some extent. However, it is still a challenging task to get rid of the effects of high-order interactions on calculating the Heisenberg parameters (and other second-order parameters). The noise introduced by Wannier orbitals, the inability to determine SIA and

the uncertainties of the resulting parameters are still the drawbacks of this approach. In addition, this method cannot give the derivatives of exchange interactions with respect to the atomic displacements, in contrast to the four-state method discussed above.

#### *3.3. More Discussions on Calculating Spin Interaction Parameters*

We should notice that, for all the methods discussed above, a rigid spin rotation approximation is used. The latter is equivalent to the supposition that the magnitudes of the spins should be constant in different configurations. However, this is not always true. For example, the magnitudes of the spins may be a little different in FM and AFM states. In Ref. [128], the energy-mapping analysis based on direct least-squares fitting (with configurations generated with different *θ* values) is adopted to study the spin interactions of MnGe and FeGe, to find that the agreement between the calculations and the model is enhanced by allowing the magnitudes of the spins to depend on the parameter *θ* (which decides the configurations) instead of using the fixed magnitudes (see Supplementary Materials of Ref. [128]). It is possible to obtain the relationship between the magnitudes of the spins and the parameter *θ* with an appropriate fitting or interpolation so that for a configuration with a new value of *θ*, the magnitudes of spins and the total energy can be predicted. Nevertheless, for a general spin configuration that cannot be described by a single *θ*, the prediction for the magnitudes of the spins can be very difficult. This is why one commonly employs the effective spin Hamiltonian by assuming that the magnitudes of the spins are constant.

Another perspective for the rigid spin rotation approximation, as implied in Ref. [143], is that even if the magnitudes of the spins are highly relevant to the configurations, the total energy can be fitted by using the directions of spins, instead of the spin vectors themselves, as the independent variables (which is mathematically equivalent to supposing the magnitudes of the spins to be constant) and considering an appropriate expansion of these variables (spin directions). For example, supposing SOC can be ignored (supposing SOC is switched off during first-principles calculations), the Heisenberg term *JijSi* · *Sj* can be expressed as *JijSi* · *Sj* = *JijSiei* · *Sjej* ; the magnitudes of the spins *Si* and *Sj* depend on the angles between these two spins or neighboring spin directions (e.g., *ek* ). Therefore,

$$J\_{\mathbf{i}\mathbf{j}}\mathbf{S}\_{\mathbf{i}}\cdot\mathbf{S}\_{\mathbf{j}} = \tilde{J}\_{0,\ i\mathbf{j}}\mathbf{e}\_{\mathbf{i}}\cdot\mathbf{e}\_{\mathbf{j}} + \mathbb{C}\_{1}\left(\mathbf{e}\_{\mathbf{i}}\cdot\mathbf{e}\_{\mathbf{j}}\right)^{2} + \mathbb{C}\_{2}\left(\mathbf{e}\_{\mathbf{i}}\cdot\mathbf{e}\_{\mathbf{k}}\right)\left(\mathbf{e}\_{\mathbf{i}}\cdot\mathbf{e}\_{\mathbf{j}}\right) + \mathbb{C}\_{3}\left(\mathbf{e}\_{\mathbf{j}}\cdot\mathbf{e}\_{\mathbf{k}\mathbf{t}}\right)\left(\mathbf{e}\_{\mathbf{i}}\cdot\mathbf{e}\_{\mathbf{j}}\right) + \cdots \quad\text{(43)}$$

That is to say, the relevance of the magnitudes of spins to the configurations can be transferred to higher-order interactions when supposing the magnitudes of spins to be constant. These "artificial" higher-order terms, only emerging for compensating for such configuration dependency, are not very physical but can somewhat improve the fitting performances (when such dependency is prominent).

All the methods discussed above assumes, besides the rigid spin rotation approximation, that magnetic moments are localized on the atoms. In addition, we note that the DFT calculation results depend on the chosen exchange-correlation functional and the value of DFT+U parameters [157].

In the above discussions, we have supposed the periodic boundary conditions, for we mainly focus on the studies of crystals. When dealing with clusters (e.g., singlemolecule magnets), we can still arrange a cluster in a crystal (using periodic boundary conditions) [159] with enough vacuum space to prevent the interactions between two clusters belonging to different periodic cells. If no periodic boundary conditions exist, the energy-mapping analysis can still work, while Green's function method based on magnetic-force linear response theory will fail because the reciprocal space is not defined. For the cases without periodic boundary conditions, theoretical chemists have developed several other approaches (such as wave-function based quantum chemical approaches) for studying magnetic interactions [69,146,147,163–165], detailed discussions of which are beyond the scope of this review.

The spin interaction parameters can also be obtained from comparing the experimental results of observable quantities such as transition temperatures, magnetization [166], specific heat [166], magnetic susceptibility [166,167], and magnon spectrum (given by inelastic neutron scattering measurements) [110,124,168–171] with the corresponding predictions given by the effective Hamiltonian model, which is similar to the idea of energy-mapping analysis based on least squares fitting. While the transition temperatures can only be used to roughly estimate the major interaction (typically the Heisenberg interaction between nearest pairs), the magnon spectrum can provide more detailed information and thus widely adopted for obtaining interaction parameters. These experimental results can also be used for checking the reliability of effective spin Hamiltonian models and the corresponding parameters obtained from first-principles calculations [172].

#### **4. Conclusions**

In this review, we summarized different types of spin interactions that an effective spin Hamiltonian may include. Recent studies have shown the importance of several kinds of high-order terms in some magnetic systems, especially biquadratic terms, fourspin ring interactions, topological chiral interactions, and chiral biquadratic interactions. In addition, we discussed in some detail the advantages and disadvantages of various methods of computing interaction parameters of the effective spin Hamiltonians. The energy-mapping analysis is easier to use, and it is less vulnerable to the effects of higherorder interactions (if carefully treated). Compared with the energy-mapping analysis, Green's function method requires less first-principle calculations and a relatively small supercell. The energy-mapping analysis usually gives a relatively good description of many kinds of states with diverse energies, while Green's function method provides a more accurate description of states close to the ground state (or the reference state). Both methods usually provide similar results and are both widely adopted in the studies of magnetic materials. We expect that first-principles based effective spin Hamiltonian will continue to play a key role in the investigation of novel magnetic states (e.g., quantum spin liquid and magnetic skyrmions).

**Author Contributions:** X.L. and H.Y. contributed equally to this work. Topic selection, X.L., F.L., J.F., and H.X.; literature search, X.L., H.Y., J.F., F.L., and H.X.; analysis and discussions, X.L., H.Y., J.F., H.X., and M.-H.W.; graphs, X.L., H.Y., and F.L.; writing-original draft preparation, X.L.; writingreview and editing, H.Y., F.L., H.X., M.-H.W., and J.F.; supervision, H.X. The manuscript was written through contributions of all authors. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is supported by NSFC (11825403, 11991061), the Program for Professor of Special Appointment (Eastern Scholar), the Qing Nian Ba Jian Program. J. S. Feng acknowledges the support from Anhui Provincial Natural Science Foundation (1908085MA10) and the Opening Foundation of State Key Laboratory of Surface Physics Fudan University (KF2019\_07).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

**Myung-Hwan Whangbo 1,2,\* , Hyun-Joo Koo 1 and Reinhard K. Kremer 3**

**Abstract:** In this review on spin exchanges, written to provide guidelines useful for finding the spin lattice relevant for any given magnetic solid, we discuss how the values of spin exchanges in transition metal magnetic compounds are quantitatively determined from electronic structure calculations, which electronic factors control whether a spin exchange is antiferromagnetic or ferromagnetic, and how these factors are related to the geometrical parameters of the spin exchange path. In an extended solid containing transition metal magnetic ions, each metal ion M is surrounded with main-group ligands L to form an MLn polyhedron (typically, n = 3–6), and the unpaired spins of M are represented by the singly-occupied d-states (i.e., the magnetic orbitals) of MLn. Each magnetic orbital has the metal d-orbital combined out-of-phase with the ligand p-orbitals; therefore, the spin exchanges between adjacent metal ions M lead not only to the M–L–M-type exchanges, but also to the M–L . . . L–M-type exchanges in which the two metal ions do not share a common ligand. The latter can be further modified by d0 cations A such as V5+ and W6+ to bridge the L . . . L contact generating M–L . . . A . . . L–M-type exchanges. We describe several qualitative rules for predicting whether the M–L . . . L–M and M–L . . . A . . . L–M-type exchanges are antiferromagnetic or ferromagnetic by analyzing how the ligand p-orbitals in their magnetic orbitals (the ligand p-orbital tails, for short) are arranged in the exchange paths. Finally, we illustrate how these rules work by analyzing the crystal structures and magnetic properties of four cuprates of current interest: α-CuV2O6 , LiCuVO4 , (CuCl)LaNb2O7 , and Cu3 (CO3 )2 (OH)2 .

**Keywords: Keywords**: spin exchange; magnetic orbitals; ligand p-orbital tails; M–L–M exchange; M–L . . . L–M exchange; α-CuV2O6 ; LiCuVO4 ; (CuCl)LaNb2O7 ; Cu3 (CO3 )2 (OH)2

#### **1. Introduction**

An extended solid consisting of transition metal magnetic ions has closely packed energy states (Figure 1a,b) so that, at a given non-zero temperature, the ground state as well as a vast number of the excited states can be thermally occupied. The thermodynamic properties such as the magnetic susceptibility and the specific heat of a magnetic system represents the weighted average of the properties associated with all thermally occupied states, with their Boltzmann factors as the weights. Such a quantity is difficult to calculate if all states were to be determined by first principle electronic structure calculations.

To generate the states of a given magnetic system and subsequently calculate the thermally-averaged physical property, a model Hamiltonian (also called a toy Hamiltonian) is invariably employed [1–3]. A typical model Hamiltonian used for this purpose is the Heisenberg-type spin Hamiltonian, Hspin, expressed as:

$$\mathbf{H}\_{\rm spin} = \sum\_{\mathbf{i} > \mathbf{j}} \mathbf{J}\_{\mathbf{i}\mathbf{j}} \mathbf{\hat{S}}\_{\mathbf{i}} \cdot \mathbf{\hat{S}}\_{\mathbf{j}} \tag{1}$$

**Citation:** Whangbo, M.-H.; Koo, H.-J.; Kremer, R.K. Spin Exchanges between Transition Metal Ions Governed by the Ligand p-Orbitals in Their Magnetic Orbitals. *Molecules* **2021**, *26*, 531. https://doi.org/ 10.3390/molecules26030531

Academic Editors:

Alessandro Stroppa and Igor Djerdj Received: 19 December 2020 Accepted: 15 January 2021 Published: 20 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

where the energy spectrum of a magnetic system as the sum of the pairwise spin exchange interactions JijSˆ i · Sˆ j is approximated. The spin operators Sˆ i and Sˆ j (at the spin sites i and j, respectively) can be treated as the spin vectors → S i and → S j , respectively, unless they operate on spin states. If all magnetic ions of a given system are identical with spin S, each term → S i · → S j can be written as → S i · → S j = S 2 cos θij where θij is the angle between the two spin vectors. In such a case, Equation (1) is rewritten as:

$$\mathbf{H}\_{\rm spin} = \sum\_{\mathbf{i} > \mathbf{j}} \mathbf{J}\_{\mathbf{i}\mathbf{j}} \mathbf{S}^2 \cos \theta\_{\mathbf{i}\mathbf{j}} \tag{2}$$

**Figure 1.** (**a**,**b**) Allowed energy states of a magnetic solid. Between the lowest-lying excited state and the ground state, there is no energy gap in (**a**), but a non-zero energy gap in (**b**). (**c**–**e**) Examples of simple spin lattices: a uniform chain in (**c**); isolated spin dimers in (**d**); and an alternating chain in (**e**). Here, all nearest-neighbor spins are antiferromagnetically coupled.

In a collinearly ordered spin arrangement of a magnetic solid, every two spin arrangements are either ferromagnetic (FM, i.e., parallel (θij = 0 ◦ )) or antiferromagnetic (AFM, i.e., antiparallel (θij = 180◦ )). With the definition of a spin Hamiltonian as in Equation (1), AFM and FM spin exchanges are represented by positive and negative Jij values, respectively. For any collinearly ordered spin arrangement, the total energy is readily written as a function of the various spin exchanges Jij.

The energy spectrum allowed for a magnetic system, and hence its magnetic properties, depend on its spin lattice. The latter refers to the repeat pattern of predominant spin exchange paths, i.e., those with large |Jij| values. For example, between the ground and the excited states, a uniform half-integer spin AFM chain (Figure 1c) has no energy gap (Figure 1a), whereas an isolated AFM dimer (Figure 1d) and an alternating AFM chain (Figure 1e) have a non-zero energy gap (Figure 1b). The spin lattice of a given magnetic solid is determined by its electronic structure, which makes it interesting how to identify the spin lattice of a magnetic solid on the basis of its atomic and electronic structures. Strong AFM exchanges between magnetic ions are often termed magnetic bonds, in contrast to chemical bonds determined by strong chemical bonding. Thus, Figure 1c represents a uniform AFM chain, and Figure 1d represents isolated AFM dimers. The magnetic bonds do not necessarily follow the geometrical pattern of the magnetic ion arrangement dictated by chemical bonding.

In a magnetic solid, transition metal ions M often share a common ligand L to form M–L–M bridges (Figure 2a). The spin exchange between the magnetic ions in an M–L–M bridge has been termed "superexchange" [4–7]. Whether the spin arrangement between two metal ions becomes FM or AFM, as described by the Goodenough–Kanamori rules formulated in the late 1950s, depends on the geometry of the M–L–M bridge [5–9]. Ever since, the Goodenough–Kanamori rules have greatly influenced the thinking of inorganic and solid-state chemists dealing with magnetic systems. Since the early 2000s, it has become increasingly clear that the magnetic properties of certain compounds cannot be conclusively explained unless also one takes into consideration the spin exchanges of the M–L . . . L–M (Figure 2b) or the M–L . . . A . . . L–M (Figure 2c) types, termed super-superexchanges [1–3],

in which the metal ions do not share a common ligand. In the late 1950s, it was impossible to imagine that spin exchanges could take place in such paths in which the M . . . M distances were very long, because the prevailing concept of chemical bonding at that time, mainly based on the valence bond picture [10], suggested that an unpaired electron of a magnetic ion is accommodated in a pure d-orbital of M. In an extended solid, each transition metal ion M is typically surrounded with main-group ligand atoms L to form an MLn polyhedron (n = 3–6), and each unpaired electron of MLn resides in a singly-occupied d-state (referred to as a magnetic orbital) of MLn, in which the d-orbital of M is combined out-of-phase with the p-orbitals of L. In this molecular orbital picture, the unpaired spin density is already delocalized from the d-orbital of M (the "magnetic orbital head") to the p-orbitals of L (the "magnetic orbital tails") [1–3]. Thus, it is quite natural to think that the spin exchange between the metal ions in an M–L . . . L–M path takes place through the overlap between the ligand p-orbital tails present in the L . . . L contact. This interaction between the p-orbital tails can be modified by the empty dπ orbitals of the d0 cation A (e.g., V5+ and W6+) in an M–L . . . A . . . L–M exchange.

**Figure 2.** Three types of spin exchange paths associated with two magnetic ions: (**a**) M–L–M, (**b**) M–L . . . L–M, and (**c**) M–L . . . A . . . L–M, where A represents a d0 cation such as V5+ or W6+. The M, L, and A are represented by red, blue, and green circles, respectively.

Given that spin exchanges between magnetic ions are determined by the p-orbital tails of their magnetic orbitals, it is not surprising that spin exchanges of the M–L . . . L–M or M–L . . . A . . . L–M type can be stronger than those of the M–L–M type, and that spin exchanges are generally strong only along a certain direction of the crystal structure. Although this aspect has been repeatedly pointed out in review articles over the years [1–3], it is still not infrequent to observe that experimental results are incorrectly interpreted simply because spin lattices have been deduced considering only the M–L–M type spin exchanges. Such an unfortunate mishap is akin to providing solutions in search of a problem. In this review, written as a tribute to John B. Goodenough for his long and illustrious scientific career culminating with the Nobel Prize in 2019, we review what electronic factors govern the nature of the M–L–M, M–L . . . L–M and M–L . . . A . . . L–M type spin exchanges, with an ultimate goal to provide several qualitative rules useful for finding the spin lattice relevant for any given magnetic system.

Our work is organized as follows: Section 2 examines how to quantitatively determine the values of spin exchanges by carrying out the energy-mapping analysis based on electronic structure calculations. Section 3 explores the electronic factors controlling whether a spin exchange is AFM or FM. In Section 4, we derive several qualitative rules that enable one to predict whether the M–L . . . L–M and M–L . . . A . . . L–M type spin exchanges are AFM or FM by analyzing how their ligand p-orbitals are arranged in the exchange paths. To illustrate how the resulting structure–property relationships operate, in Section 5 we examine the crystal structures and magnetic properties of α-CuV2O6, LiCuVO4, (CuCl)LaNb2O7 and Cu3(CO3)2(OH)2. Our concluding remarks are presented in Section 6.

#### **2. Energy Mapping Analysis for Quantitative Evaluation of Spin Exchanges**

This section probes how to define the relevant spin Hamiltonian for a given magnetic system on a quantitative level. This requires the determination of the spin exchanges to include in the spin Hamiltonian. For various collinearly ordered spin states of a given

magnetic system, one finds the expressions for their relative energies in terms of the spin exchange parameters Jij to evaluate, performs DFT+U [11] or DFT+hybrid [12] electronic structure calculations for the ordered spin states to determine the numerical values for their relative energies, and finally maps the two sets of relative energies to find the numerical values of the exchanges.

#### *2.1. Using Eigenstates*

To gain insight into the meaning and the nature of a spin exchange, we examine a spin dimer made up of two S = 1/2 ions (Figure 3a). The spin Hamiltonian describing the energies of this dimer is given by:

$$\mathbf{H}\_{\rm spin} = \mathbf{J} \mathbf{\hat{S}}\_{\rm i} \cdot \mathbf{\hat{S}}\_{\rm j} \tag{3}$$

**Figure 3.** (**a**) Spin dimer made up of two S = 1/2 ions at sites i and j. Each site has one unpaired spin. The magnetic orbitals at the sites i and j are represented by φi and φj , respectively. (**b**) The bonding and antibonding states, Ψ1 and Ψ2 , resulting from the interactions between φi and φj , are split in energy by ∆e. (**c**) The singlet and triplet electron configurations resulting from Ψ1 and Ψ2 .

The spin states allowed for this dimer are the singlet and triplet states, ΨS and ΨT, respectively.

$$\begin{array}{l} \Psi\_{\mathsf{T}} = |\uparrow\uparrow\rangle, \ |\downarrow\downarrow\rangle, \ (|\uparrow\downarrow\rangle + |\downarrow\uparrow\rangle)/\sqrt{2} \\ \Psi\_{\mathsf{S}} = (|\uparrow\downarrow\rangle - |\downarrow\uparrow\rangle)/\sqrt{2} \end{array} \tag{4}$$

It can be readily shown that these states are the eigenstates of the spin Hamiltonian by rewriting Equation (3) as:

$$\mathbf{H}\_{\rm spin} = \mathbf{J}\mathbf{\hat{S}}\_{\rm iz}\mathbf{\hat{S}}\_{\rm jz} + (\mathbf{J}/2)[\mathbf{\hat{S}}\_{\rm i}\mathbf{\hat{S}}\_{\rm j-} + \mathbf{\hat{S}}\_{\rm i-}\mathbf{\hat{S}}\_{\rm j+}] \tag{5}$$

where Sˆmz (m = i, j) is the z-component of Sˆm, while Sˆm+ and Sˆm− are the raising and lowering operators associated with Sˆm, respectively. Then, it is found that:

$$\begin{aligned} \mathsf{H}\_{\mathrm{spin}}\Psi\_{\mathrm{T}} &= (\mathsf{J}/4)\Psi\_{\mathrm{T}}\\ \mathsf{H}\_{\mathrm{spin}}\Psi\_{\mathrm{S}} &= (-\mathsf{3J}/4)\Psi\_{\mathrm{S}} \end{aligned} \tag{6}$$

Thus:

$$
\Delta \mathbf{E}\_{\text{eigen}} = \mathbf{E}\_{\text{T}} - \mathbf{E}\_{\text{S}} = \langle \mathbf{\varPsi}\_{\text{T}} | \mathbf{H}\_{\text{spin}} | \mathbf{\varPsi}\_{\text{T}} \rangle - \langle \mathbf{\varPsi}\_{\text{S}} | \mathbf{H}\_{\text{spin}} | \mathbf{\varPsi}\_{\text{S}} \rangle = \mathbf{J} \tag{7}
$$

Therefore, the spin exchange J is related to the energy difference between the singlet and triplet states of the spin dimer, as illustrated in Figure 4a.

triplet states of the spin dimer, as illustrated in Figure 4a.

**Figure 4.** Relationships of the spin exchange J to the energy difference between two spin states of a spin dimer made up of two S = 1/2 ions in terms of (**a**) the eigenstates and (**b**) the brokensymmetry states. The legends HS and BS in (**b**) refer to the high-symmetry and broken-symmetry **Figure 4.** Relationships of the spin exchange J to the energy difference between two spin states of a spin dimer made up of two S = 1/2 ions in terms of (**a**) the eigenstates and (**b**) the broken-symmetry states. The legends HS and BS in (**b**) refer to the high-symmetry and broken-symmetry states, respectively.

E E E H H J Δ eigen = T − S = ψT spin ψT − ψS spin ψS = (7)

Therefore, the spin exchange J is related to the energy difference between the singlet and

#### states, respectively. *2.2. Using Broken-Symmetry States*

*2.2. Using Broken-Symmetry States*  In general, it is not an easy task to find the eigenstates of a general spin Hamiltonian (e.g., Equation (1)), which makes it difficult to relate the spin exchanges to the energy differences between the eigenstates of a magnetic system. However, the energies for the broken-symmetry (BS) states (i.e., the non-eigenstates) of a spin Hamiltonian are easy to evaluate. For the magnetic dimer in Section 2.1, the high-spin (HS) state, or ψHS = ↑↑ ↓↓ In general, it is not an easy task to find the eigenstates of a general spin Hamiltonian (e.g., Equation (1)), which makes it difficult to relate the spin exchanges to the energy differences between the eigenstates of a magnetic system. However, the energies for the brokensymmetry (BS) states (i.e., the non-eigenstates) of a spin Hamiltonian are easy to evaluate. For the magnetic dimer in Section 2.1, the high-spin (HS) state, ΨHS = |↑↑i or |↓↓i, which has an FM spin arrangement, is an eigenstate of the spin Hamiltonian (Equation (3)). However, the low-spin (LS) state can be expressed as:

$$\Psi\_{\rm LS} = |\uparrow\downarrow\rangle \text{ or } |\downarrow\uparrow\rangle$$

ψLS = ↑↓ or ↓↑ This state has an AFM spin arrangement and is the BS state of the spin Hamiltonian. Using Equation (5), it is found that:

H ( J/ 4) (J/ 2)

↑↓ = − ↑↓ + ↓↑

$$\begin{array}{c} \mathbf{H}\_{\text{spin}}|\uparrow\downarrow\rangle = (-\mathbf{J}/\mathbf{4})\left|\uparrow\downarrow\right\rangle + (\mathbf{J}/\mathbf{2})\left|\downarrow\uparrow\right\rangle\\ \mathbf{H}\_{\text{spin}}|\downarrow\uparrow\rangle = (-\mathbf{J}/\mathbf{4})\left|\downarrow\uparrow\right\rangle + (\mathbf{J}/\mathbf{2})\left|\uparrow\downarrow\right\rangle \end{array} \tag{8}$$

Therefore:

Therefore:

$$
\langle \uparrow | \mathbf{H}\_{\text{spin}} | \uparrow \downarrow \rangle = \langle \downarrow \uparrow | \mathbf{H}\_{\text{spin}} | \downarrow \uparrow \rangle = -\mathbf{J}/4 \tag{9}
$$

From Equations (6) and (9), we obtain:

spin

spin

$$
\Delta \mathbf{E\_{BS}} = \langle \Psi\_{\rm HS} | \mathbf{H\_{spin}} | \Psi\_{\rm HS} \rangle - \langle \Psi\_{\rm BS} | \mathbf{H\_{spin}} | \Psi\_{\rm BS} \rangle = \mathbf{J}/2 \tag{10}
$$

E H H J/ 2 Δ BS = ψHS spin ψHS − ψBS spin ψBS = (10)

From Equations (6) and (9), we obtain: The spin exchange J is related to the energy difference between the states including the BS states (Figure 4b).

#### *2.3. Energy Mapping*

The spin exchange J is related to the energy difference between the states including the BS states (Figure 4b). *2.3. Energy Mapping*  For a general spin lattice, the energy difference between any two states (involving BS and HS states) can be readily determined by using the spin Hamiltonian of Equation (3), which is expressed as a function of unknown constants Jij. To evaluate the spin exchanges Jij of a general spin Hamiltonian, it is necessary to numerically determine the relative energies of various ordered spin states. This is done by performing DFT+U [11] or For a general spin lattice, the energy difference between any two states (involving BS and HS states) can be readily determined by using the spin Hamiltonian of Equation (3), which is expressed as a function of unknown constants Jij. To evaluate the spin exchanges Jij of a general spin Hamiltonian, it is necessary to numerically determine the relative energies of various ordered spin states. This is done by performing DFT+U [11] or DFT+hybrid [12] electronic structure calculations for the collinearly ordered spin states of a magnetic system on the basis of the electronic Hamiltonian, Helec. These types of calculations ensure that the electronic structures calculated for various ordered spin states have a bandgap as expected for a magnetic insulator. Suppose that N different spin exchange paths Jij are considered to describe a given magnetic solid. If one considers N + 1 ordered spin states, for example, one can determine N different relative energies ∆Espin(i) (i = 1, 2,···, N) expressed in terms of N different spin exchanges Jij (in principle, one could use more spin states, but at least N+1 are necessary; using more would produce error bars and increase the precision of the analysis). By carrying out electronic structure calculations for the N+1 ordered spin states of the magnetic system, one obtains the numerical values for the N different relative energies ∆Eelec(i) (i = 1, 2,···, N). Then, by equating the ∆Espin(i) (i = 1, 2, ···, N) values

to the corresponding ∆Eelec(i) (i = 1, 2, ···, N) values, the N different spin exchanges Jij are obtained.

$$
\Delta \mathbf{E}\_{\text{spin}}(\mathbf{i}) \text{ (i = 1, 2, \dots, N)} \leftrightarrow \Delta \mathbf{E}\_{\text{elec}}(\mathbf{i}) \text{ (i = 1, 2, \dots, N)}\tag{11}
$$

The spin exchanges discussed so far are known as Heisenberg exchanges. There are other variants of interactions between spins which, though weaker than Heisenberg exchanges in strength, are needed to explain certain magnetic properties not covered by the symmetrical Heisenberg exchanges. They include Dzyaloshinskii–Moriya exchanges (or antisymmetric exchanges) and asymmetric exchanges [2,3]. The energy-mapping analysis based on collinearly ordered spin states allows one to determine the Heisenberg spin exchanges only. To evaluate the Dzyaloshinskii–Moriya and asymmetric spin exchanges, the energy-mapping analysis employs the four-state method [2,13], in which non-collinearly ordered broken-symmetry states are used. A further generalization of this energy-mapping method was developed to enable the evaluation of other energy terms that one might include in a model spin Hamiltonian [14].

#### **3. Qualitative Features of Spin Exchange**

A spin Hamiltonian appropriate for a given magnetic system is one that consists of the predominant spin exchange paths. Such a spin Hamiltonian can be determined by evaluating the values of various possible spin exchanges for the magnetic system by performing the energy-mapping analysis as described in Section 2. If the spin lattice is chosen without quantitatively evaluating its spin exchanges, one might inadvertently choose a spin lattice irrelevant for the interpretation of the experimental data. When simulating the thermodynamic properties using a chosen set of spin exchanges, the values of the spin exchanges are optimized until they provide the best possible simulation even if the chosen spin lattice is incorrect from the viewpoint of electronic structure. Thus, in principle, more than one spin lattice might provide an equally good simulation. In interpreting the experimental results of a magnetic system with a correct spin lattice, it is crucial to know what electronic and structural factors control the signs and the magnitudes of spin exchanges.

#### *3.1. Parameters Affecting Spin Exchanges*

To examine what energy parameters govern the sign and magnitude of a spin exchange, we revisit the spin exchange of the spin dimer (Figure 3a) by explicitly considering the electronic structures of its singlet and triplet states. For simplicity, we represent each spin site with one magnetic orbital. As will be discussed in Section 5, the nature of the magnetic orbital plays a crucial role in determining the sign and magnitude of a spin exchange. We label the magnetic orbitals located at the spin sites i and j as φi and φj , respectively (Figure 3a). These orbitals overlap weakly, and hence interact weakly, to form the in-phase and out-of-phase states, Ψ1 and Ψ2, respectively, with energy split ∆e between the two (Figure 3b). The overlap integral hφi |φj E between φi and φj is small for magnetic systems, so Ψ1 and Ψ2 are well approximated by:

$$\Psi\_1 \approx \frac{\left(\Phi\_{\mathbf{i}} + \Phi\_{\mathbf{j}}\right)}{\sqrt{2}}, \Psi\_2 \approx \frac{\left(\Phi\_{\mathbf{i}} - \Phi\_{\mathbf{j}}\right)}{\sqrt{2}} \tag{12}$$

For simplicity, it is assumed here that the in-phase combination (i.e., the bonding combination) is described by the plus combination, which amounts to the assumption that the overlap integral hφi |φj E is positive. The energy split ∆e is approximately proportional to the overlap integral:

$$
\Delta \mathbf{e} \propto \langle \phi\_{\mathbf{i}} | \phi\_{\mathbf{j}} \rangle \tag{13}
$$

In understanding the qualitative features describing the electronic energy difference between the singlet and triplet states, ∆Eelec, and hence the spin exchange J, it is necessary to consider two other quantities. One is the on-site repulsion Uii:

$$\mathbf{U\_{ii}} = \langle \phi\_{\mathbf{i}}(\mathbf{1})\phi\_{\mathbf{i}}(\mathbf{2})|\mathbf{1}/\mathbf{r}\_{12}|\,\phi\_{\mathbf{i}}(\mathbf{1})\phi\_{\mathbf{i}}(\mathbf{2})\rangle\ (\mathbf{i} = \mathbf{1}, \mathbf{2})\tag{14}$$

where the product φi (m)φi (m) is the electron density ρii(m) associated with the orbital φi (m) occupied by electron m (=1, 2). Thus, Uii is the self-repulsion when the orbital φi is occupied by two electrons. If the spin sites i and j are identical in nature, the on-site repulsion Ujj at the site j is the same as Uii, therefore it is convenient to use the symbol U to represent both Uii and Ujj (i.e., U = Uii = Ujj). The other quantity of interest is the exchange repulsion Kij between φi and φj :

$$\mathbf{K}\_{\mathbf{ij}} = \left\langle \phi\_{\mathbf{i}}(\mathbf{1})\phi\_{\mathbf{j}}(\mathbf{2}) | \mathbf{1}/\mathbf{r}\_{12}| \,\phi\_{\mathbf{j}}(\mathbf{1})\phi\_{\mathbf{i}}(\mathbf{2}) \right\rangle \tag{15}$$

This is the self-repulsion arising from the overlap electron density:

$$
\rho\_{\vec{\text{ij}}}(\mathbf{m}) \, \, = \phi\_{\vec{\text{i}}}(\mathbf{m}) \phi\_{\vec{\text{j}}}(\mathbf{m}) \, \left(\mathbf{m} \, \, \, = \, \, 1 \, \, \mathcal{D}\right) \tag{16}
$$

To illustrate the difference between the overlap integral and overlap electron density, we consider the px and py atomic orbitals located at a same atomic site (Figure 5a,b). The product px py represents the overlap electron density ρxy, which consists of four overlapping regions (Figure 5c); two regions of positive electron density (colored in pink) and two regions of negative electron density (colored in cyan). The overlap integral hpx |py E is the sum of these four overlap electron densities, which adds up to zero.

**Figure 5.** The overlap density resulting from two p-orbital at a given atomic site: (**a**) a px orbital; (**b**) a py orbital; and (**c**) the overlap density between the two orbitals, ρxy = px py . The pink and cyan regions have positive and negative values, respectively.

The exchange repulsion between px and py is written as:

$$\mathbf{K\_{xy}} = \left< \mathbf{p\_x}(1)\mathbf{p\_x}(2) \middle| 1/\mathbf{r\_{12}} \middle| \mathbf{p\_y}(1)\mathbf{p\_y}(2) \right> \tag{17}$$

Equation (17) is given by the sum of the self-repulsion resulting from each overlapping region of the overlap electron density ρxy (Figure 5c). Each overlapping region, be it positive or negative, leads to a positive repulsion; therefore, the exchange repulsion Kxy is positive.

#### *3.2. Two Competing Components of Spin Exchange*

In Figure 3c, the configuration ΨT represents the triplet state of the dimer while the configurations ΨS1 and ΨS2 each represent a singlet state. For a magnetic system, for which ∆e is very small, the singlet state is not well described by ΨS1 alone and is represented by a linear combination of ΨS1 and ΨS2, i.e., ΨS = c1ΨS1 + c2ΨS2, where the mixing coefficients c1 and c2 are determined by the interaction between ΨS1 and

ΨS2, namely, hΨS1|Helec|ΨS2i = Kij, as well as by the energies of ΨS1 and ΨS2, that is, hΨS1|He*l*ec|ΨS1i = E1 and hΨS2|Helec|ΨS2i = E2. After some lengthy manipulations under the condition, Kij >> |E1−E2|, which is satisfied for magnetic systems, the energy difference between the triplet and singlet states, and hence the spin exchange J, is expressed as [15]:

$$\mathbf{J} = \mathbf{E}\_{\rm T} - \mathbf{E}\_{\rm S} \approx -2\mathbf{K}\_{\rm ij} + \frac{\left(\Delta \mathbf{e}\right)^{2}}{\mathbf{U}} \tag{18}$$

Note that the spin exchange J consists of two components, J = JF + JAF, where:

$$\begin{array}{l} \mathrm{J}\_{\mathrm{F}} = -2\mathrm{K}\_{\mathrm{ij}} < 0\\ \mathrm{J}\_{\mathrm{AF}} = \frac{\left(\Delta \mathrm{e}\right)^{2}}{\mathrm{U}} > 0 \end{array} \tag{19}$$

The magnitude of the FM component | JF| increases as the exchange repulsion Kij increases, namely, as the overlap density ρij = φiφj  increases. The strength of the AFM component JAF increases with increasing the energy split ∆e, i.e., with increasing the overlap integral hφi |φj E , while JAF decreases as the on-site repulsion U increases.

As already discussed in Section 2, the quantitative values of spin exchanges can be accurately determined by the energy-mapping analysis based on first-principles DFT+U or DFT+hybrid calculations. The purpose of Equations (18) and (19) is not to determine the numerical value of any spin exchange, but to show that each spin exchange J consists of two competing components, JF and JAF, that the overall sign of J is determined by which component dominates, and which electronic parameters govern the strength of each component.

#### **4. Spin Exchanges Determined by the Ligand p-Orbitals in the Magnetic Orbitals**

To illustrate how spin exchanges of transition metal magnetic ions are controlled by the ligand p-orbitals in their magnetic orbitals, we consider various spin exchanges involving Cu2+ (d9 , S = 1/2) ions as an example, which typically form axially elongated CuL6 octahedra (Figure 6a). The energies of their d-states are split as (xz, yz) < xy < 3z2–r2 < x2–y2 (Figure 6b) so that the magnetic orbital of each Cu2+ ion is represented by the x 2–y2 state, in which the Cu x2–y2 orbital induces σ-antibonding with the p-orbitals of four equatorial ligands L (Figure 6c). In this section, we examine how the various types of spin exchanges associated with Cu2+ ions are controlled by the ligand p-orbitals of their x2–y2 states. The major component of the magnetic orbital of a Cu2+ ion (Figure 6d) is the Cu d-orbital (i.e., the magnetic orbital "head"), and the minor component the ligand p-orbitals (i.e., the magnetic orbital "tail"). In this section, we probe how the nature and strengths of M–L–M, M–L . . . L–M and M–L . . . A . . . L–M-type exchanges are determined by how the ligand p-orbital tails of their magnetic orbitals are arranged in their exchange paths.

**Figure 6.** (**a**) An axially elongated CuL6 octahedron. (**b**) The electron configuration of a Cu2+(d9 ) ion at an axially elongated octahedral site. (**c**) The magnetic orbital of the Cu2+(d9 ) ion at an axially elongated octahedral site, which is contained in the CuL4 equatorial plane. (**d**) The head and tails of the magnetic orbital.

cation.

#### *4.1. M–L* . . . *L–M and M–L* . . . *A* . . . *L–M Spin Exchanges* overlap electron density resulting from the p-orbital tails will be practically zero. Then, the spin exchange J would be zero. If the L…L contact of such an M**–**L…L**–**M exchange

vanishes (i.e., Δe = 0) (Figure 9c) so that JAF = 0. In addition, JF should vanish because the

through-space interaction in the M**–**L…L**–**M exchange leads to a weak overall M**–** L…A…L**–**M spin exchange due to the effect of the through-bond interaction (Figure 8c). An example of very weak Cu**–**L…L**–**Cu exchange is shown in Figure 9a, in which the

*Molecules* **2019**, *24*, x FOR PEER REVIEW 10 of 29

The next-nearest-neighbor (nnn) spin exchange Jnnn (Figure 7a) that occurs in a CuL2 (L = O, Cl, Br) ribbon chain, is obtained by sharing the opposite edges of CuL4 square planes. This spin exchange is an example of a strong Cu–L . . . L–Cu exchange (Figure 7b), when the L . . . L contact distance is in the vicinity of the van der Waals distance. The latter will be assumed to be the case in what follows. The two magnetic orbitals interact across the L . . . L contacts through the overlap of their p-orbital tails. This through-space interaction leads to the in-phase and out-of-phase combinations (Ψ+ and Ψ−, respectively) of the magnetic orbitals (Figure 7c), with energy split ∆e between the two (Figure 7d). This makes the AFM component JAF non-zero. The overlap electron density associated with the interacting p-orbital tails is nearly zero, so the FM component JF is practically zero. As a result, Jnnn becomes AFM. path is bridged by a d0 transition metal cation A to form an M**–**L…A…L**–**M spin exchange (Figure 9d), the ψ− state of the M**–**L…L**–**M path interacts with the empty dπ orbital of A, thereby lowering its energy (Figure 9e) while that of the ψ+ state is unchanged. Thus, when the M**–**L…L**–**M exchange has a very weak through-space interaction, the throughbond interaction induces the large energy split Δe, so that the overall M**–**L…A…L**–**M spin exchange becomes strong (Figure 9f). In short, a strong M**–**L…L**–**M exchange becomes a weak M**–**L…A…L**–**M exchange when the L…L linkage is bridged by d0 cations, while a weak M**–**L…L**–**M exchange becomes a strong M**–**L…A…L**–**M exchange when the L…L linkage is bridged by a d0

**Figure 7.** (**a**) The next-nearest-neighbor spin exchange Jnnn in a CuL2 ribbon chain made up of edgesharing CuL4 square planes. (**b**) A case of strong Cu**–**L…L**–**Cu exchange, which represents the next-nearest-neighbor exchange Jnnn. (**c**) The in-phase and out-of-phase combinations of the two magnetic orbitals, ψ+ and ψ− , respectively, associated with the Cu**–**L…L**–**Cu exchange. (**d**) The large energy split Δe resulting from the through-space (TS) interaction in the Cu**–**L…L**–**Cu exchange becomes small in the Cu**–**L…A…L**–**Cu exchange as a result of the through-bond (TB) interaction that occurs with the ψ− state. (**e**) A Cu**–**L…A…L**–**Cu exchange generated when each **Figure 7.** (**a**) The next-nearest-neighbor spin exchange Jnnn in a CuL2 ribbon chain made up of edge-sharing CuL4 square planes. (**b**) A case of strong Cu–L . . . L–Cu exchange, which represents the next-nearest-neighbor exchange Jnnn. (**c**) The in-phase and out-of-phase combinations of the two magnetic orbitals, Ψ+ and Ψ−, respectively, associated with the Cu–L . . . L–Cu exchange. (**d**) The large energy split ∆e resulting from the through-space (TS) interaction in the Cu–L . . . L–Cu exchange becomes small in the Cu–L . . . A . . . L–Cu exchange as a result of the through-bond (TB) interaction that occurs with the Ψ− state. (**e**) A Cu–L . . . A . . . L–Cu exchange generated when each L . . . L contact is bridged by a d0 metal cation A. (**f**) The bonding interaction of the Ψ− state with the dπ orbital of A.

L…L contact is bridged by a d0 metal cation A. (**f**) The bonding interaction of the ψ− state with the dπ orbital of A. In the M–L . . . L–M spin exchange Jnnn discussed above, there are two equivalent exchange paths due to the ribbon structure. Suppose that each L . . . L contact of the M–L . . . L–M exchange path is bridged by a d0 transition metal cation A (e.g., V5+ and W6+) to form an M–L . . . A . . . L–M spin exchange (Figure 7e). We now analyze the relative strengths of the M–L . . . L–M and M–L . . . A . . . L–M spin exchanges. Across the L . . . L contact of the M–L . . . L–M path, the p-orbital tails of L are combined in-phase in Ψ+, but out-of-phase in Ψ−. Thus, the empty dπ orbital of the cation A interacts in-phase with Ψ− (Figure 7f) to lower the energy of Ψ−, but it does not interact with Ψ+, therefore the energy of Ψ+ is unaffected. This selective interaction of the bridging d0 cation A with the L . . . L contact of the M–L . . . L–M spin exchange has a dramatic consequence on the strength of an M–L . . . A . . . L–M spin exchange. When the M–L . . . L–M exchange has a strong through-space interaction, the through-bond interaction reduces the large energy split ∆e to a small value, thereby weakening the overall M–L . . . A . . . L–M spin exchange (Figure 7d).

Another example of a strong M–L . . . L–M exchange occurs when the p-orbital tails of L are pointed to each other along the L . . . L contact (Figure 8a). The through-space interaction between the magnetic orbitals leads to the in-phase and out-of-phase combinations, Ψ+ and Ψ−, respectively (Figure 8b), with a large energy split ∆e between the two (Figure 8c) and a negligible overlap electron density between the interacting p-orbital tails. As a result, the M–L . . . L–M spin exchange becomes AFM. If the L . . . L contact of an M–L . . . L–M exchange path is bridged by a d0 transition metal cation A to form an M–L . . . A . . . L–M spin exchange (Figure 8d), only the Ψ− state of the M–L . . . L–M path interacts effectively with one of the empty dπ orbitals of A (Figure 8e). Thus, a strong through-space interaction in the M–L . . . L–M exchange leads to a weak overall M–L . . . A . . . L–M spin exchange due to the effect of the through-bond interaction (Figure 8c).

**Figure 8.** (**a**) A case of strong Cu–L . . . L–Cu spin exchange where the two Cu–L bonds leading to the L . . . L contacts are linear. (**b**) The in-phase and out-of-phase combinations (Ψ+ and Ψ−, respectively) of the two magnetic orbitals, resulting from the through-space (TS) interactions. (**c**) The energy split ∆e between Ψ+ and Ψ− is large when the overlap between the p-orbital tails is large. (**d**) A Cu–L . . . A . . . L–Cu exchange generated when the L . . . L contact is bridged by a d0 metal cation A. (**e**) The bonding interaction of the Ψ− state with the dπ orbital of A. (**f**) The large energy split ∆e resulting from the through-space (TS) interaction in the Cu–L . . . L–Cu exchange becomes small in the Cu–L . . . A . . . L–Cu exchange as a result of the through-bond (TB) interaction that occurs primarily with the Ψ− state.

An example of very weak Cu–L . . . L–Cu exchange is shown in Figure 9a, in which the two magnetic orbitals are arranged such that the p-orbital tails are orthogonal to each other (Figure 9b), and their overlap vanishes so that the energy split between Ψ+ and Ψ− vanishes (i.e., ∆e = 0) (Figure 9c) so that JAF = 0. In addition, JF should vanish because the overlap electron density resulting from the p-orbital tails will be practically zero. Then, the spin exchange J would be zero. If the L . . . L contact of such an M–L . . . L–M exchange path is bridged by a d0 transition metal cation A to form an M–L . . . A . . . L–M spin exchange (Figure 9d), the Ψ− state of the M–L . . . L–M path interacts with the empty dπ orbital of A, thereby lowering its energy (Figure 9e) while that of the Ψ+ state is unchanged. Thus, when the M–L . . . L–M exchange has a very weak through-space interaction, the through-bond interaction induces the large energy split ∆e, so that the overall M–L . . . A . . . L–M spin exchange becomes strong (Figure 9f).

**Figure 9.** (**a**) A case of weak Cu–L . . . L–Cu spin exchange, where the two Cu–L bonds leading to the L . . . L contacts are orthogonal. (**b**) The in-phase and out-of-phase combinations (Ψ+ and Ψ−, respectively) of the two magnetic orbitals. (**c**) The vanishing energy split ∆e resulting from the through-space (TS) interaction in the Cu–L . . . L–Cu exchange is enhanced in the Cu–L . . . A . . . L–Cu exchange as a result of the through-bond (TB) interaction that occurs primarily with the Ψ− state. (**d**) A Cu–L . . . A . . . L–Cu exchange generated when the L . . . L contact is bridged by a d0 metal cation A. (**e**) The in-phase interaction of the Ψ− state with the dπ orbital of A.

In short, a strong M–L . . . L–M exchange becomes a weak M–L . . . A . . . L–M exchange when the L . . . L linkage is bridged by d0 cations, while a weak M–L . . . L–M exchange becomes a strong M–L . . . A . . . L–M exchange when the L . . . L linkage is bridged by a d 0 cation.

#### *4.2. M–L–M Spin Exchanges*

The Goodenough–Kanamori rules cover these types of spin exchanges [5–9]. For the sake of completeness, we discuss these types of spin exchanges from the viewpoint of the ligand p-orbital tails on the basis of Equations (18) and (19). Let us consider a Cu2L6 dimer resulting from two CuL4 square planes obtained by sharing an edge (Figure 10a), where the ligand L can be O, Cl, or Br. The two magnetic orbitals associated with the nearest-neighbor (nn) spin exchange Jnn, presented in Figure 10b, interact at the bridging ligands L of the M–L–M paths. If the CuL4 units have an ideal square planar shape, the ∠M–L–M angle becomes 90◦ so that the two p-orbital tails at the bridging ligands L are orthogonal to each other. Thus, as discussed in Section 3, the overlap integral between them is zero. Therefore, the in-phase and out-of-phase combinations of the two magnetic orbitals (Figure 10b) is not split in energy, so ∆e = 0 (Figure 10c) and JAF = 0. However, the overlap electron density between the two p-orbital tails at the bridging ligand L is not zero, i.e., JF > 0. Thus, the M–L–M spin exchange becomes FM. When the ∠M–L–M angle deviates from 90◦ (Figure 10d), the two p-orbital tails at the bridging ligands L are no longer orthogonal to each other so the overlap integral between them is non-zero. Therefore, the in-phase and out-of-phase combinations of the two magnetic orbitals (Figure 10e) differ in energy, so ∆e > 0 (Figure 10f) and JAF is non-zero. The overlap electron density between the two p-orbital tails is non-zero, so JF is non-zero. Thus, whether the spin exchange is FM or AFM depends on which component, JF or JAF, dominates, which in turn depends on the ∠M–L–M angle φ. Typically, the angle φ where FM changes to AFM is slightly greater than 90◦ due to the involvement of the ligand s-orbital [15], which is commonly neglected for simplicity.

**Figure 10.** (**a**–**c**) A Cu–L–Cu spin exchange with ∠M–L–M angle φ = 90◦ : This occurs in a Cu2L6 dimer shown in (**a**), which is made up of two coplanar CuL4 square planes by sharing an edge. The in-phase and out-of-phase combinations of the two magnetic orbitals in (**b**), and the energy split ∆e between the two in (**c**). (**d**–**f**) A Cu–L–Cu spin exchange with ∠M–L–M angle φ > 90◦ : This occurs when two non-coplanar CuL4 square planes are corner-shared as in (**d**). The in-phase and out-of-phase combinations of the two magnetic orbitals in (**e**), and the energy split ∆e between the two in (**f**). (**g**–**k**) Effect of the molecular anions OH− on the spin exchange of the edge-sharing Cu2O6 dimer: the shared edge consists of two OH− ions in (**g**); the p-orbital tail of one OH− ligand in a CuO4 square plane expected if the OH− is treated as O2− in (**h**); the three oxygen lone pairs associated with an isolated OH− in (**i**); the tilting of the O 2p-orbital tail by mixing two oxygen lone pairs in (**j**); and the arrangement of the tilted O 2p-tails at the bridging O atoms of the Cu2O6 dimer in (**k**).

So far in our discussion, it has been implicitly assumed that each main-group ligand L exists as a spherical anion (e.g., each O as an O2− anion, and each Cl atom as a Cl− anion). However, this picture is not quite accurate when the ligand atom makes a strong covalent bonding with another main-group element to form a molecular anion such as OH−. Suppose that each O atom on the shared edge of the CuO6 (L = O) dimer is not a O2− but an OH− anion (Figure 10g). Then, the ligand p-orbital tail on that O cannot be the p-orbital pointed along one lobe of the Cu x2–y2 orbital (Figure 10h) because it is incompatible with the O–H bonding, which has three directional O lone pairs depicted in Figure 10i. To satisfy both the strong covalent-bonding with H and the weak covalentbonding with Cu, the O lone pair of OH− tilts slightly toward one lobe of the x2–y2 orbital (Figure 10j). As a result, the ligand p-orbital tails, arising from the two magnetic orbitals at the bridging O atoms, are not orthogonal as in Figure 10b but become more parallel to each other (Figure 10k). As a result, the spin exchange between the two Cu2+ ions in Figure 10g becomes AFM (see below). Another molecular anion of interest is the carbonate ion CO3 2−, in which each O atom makes a strong covalent bond with C; therefore, the O atoms of CO3 2− should not be treated as isolated O2− anions in their coordination with transition metal cations. In general, the presence of molecular anions such as OH− and CO3 2− in a magnetic solid makes it difficult to deduce, on a qualitative reasoning, what its spin lattice would be. This is where the quantitative energy-mapping analysis is indispensable, because it does not require any qualitative reasoning.

#### *4.3. Qualitative Rules for Spin Exchanges Based on the p-Orbital Tails of Magnetic Orbitals*

In a magnetic orbital of an MLn polyhedron, the d-orbital of M dictates by its symmetry which p-orbitals of the ligands L become the p-orbital tails. From the viewpoint of orbital interaction, a spin exchange between magnetic ion is none other than the interaction

between their magnetic orbitals. The latter is caused by the interaction between their p-orbital tails, not by that between their d-orbital heads. In other words, a spin exchange is not a "head-to-head" interaction but a "tail-to-tail" interaction. By considering these tail-to-tail interactions described above, we arrive at the following four qualitative rules governing the nature and strengths of the M–L . . . L–M and M–L . . . A . . . L–M type exchanges under the assumption that the L . . . L contact distance is in the vicinity of the van der Waals distance:

These rules on the M–L . . . L–M and M–L . . . A . . . L–M exchanges should be used together with the Goodenough–Kanamori rules in choosing a proper set of spin exchanges to evaluate using the energy-mapping analysis based on DFT+U or DFT+hybrid calculations. In principle, this analysis can provide quantitative values for any possible exchanges of a given magnetic system. However, even this quantitative tool cannot determine the value of any spin exchange unless it is included in the set of spin exchanges for the energymapping analysis. It is paramount to consider in detail the structural features governing the strengths of spin exchanges in order not to miss exchange paths crucial for defining the correct spin lattice of a given solid.

#### **5. Representative Examples**

In Section 4, we analyzed the structural features governing the nature of the three types of spin exchanges, i.e., M–L–M, M–L . . . L–M and M–L . . . A . . . L–M, which occur in various magnetic solids. This section will discuss the occurrence of these exchanges in actual magnetic solids by analyzing the crystal structures and magnetic properties of four representative magnetic solids, α-CuV2O6, LiCuVO4, (CuCl)LaNb2O7 and Cu3(CO3)3(OH)3. α-CuV2O6, LiCuVO4, and (CuCl)LaNb2O7 were chosen to show that correct spin lattices can be readily predicted by the qualitative rules of Section 4.3, although they have to be confirmed by performing the energy-mapping analyses. Azurite Cu3(CO3)2(OH)2 was chosen to demonstrate that the spin lattice of a certain magnetic system cannot be convincingly deduced solely on the basis of the qualitative rules. The magnetic ions of such a system are coordinated with molecular anions in which the first-coordinate main-group ligands L make strong covalent bonds with other main-group elements (e.g., H in the OH− ion, and C in the CO3 2− ion). In such a case, use of the energy-mapping analysis is the only recourse with which to find the spin lattice correct for a given system.

#### *5.1. Two-Dimensional Behavior of α-CuV2O6*

The magnetic properties of α-CuV2O6 were initially analyzed in terms of a onedimensional (1D) spin S = 1/2 Heisenberg chain model with uniform nearest-neighbor AFM spin exchange [16–18]. However, a rather high Néel temperature of ~22.4 K indicated the occurrence of a substantial interchain spin exchange of the order of 50% of the intrachain exchange, casting serious doubts on the applicability of a simple chain description. α-CuV2O6 consists of CuO4 chains, made up of edge-sharing CoO6 octahedra, which run along the a-direction (Figure 11a). If the two axially elongated Cu–O bonds are removed from each CuO6 octahedron to identify its CuO4 equatorial plane containing the magnetic orbital, one finds that each CuO4 chain of edge-sharing CoO6 octahedra becomes a chain of stacked CuO4 square planes (Figure 11b). In each stack-chain along the a-direction, adjacent CuO4 square planes are parallel to each other such that the adjacent magnetic

orbitals generate neither a non-zero overlap nor a non-zero overlap electron density. The same is true between adjacent CuO4 square planes along the b-direction (Figures 11c and 12a). However, between adjacent CuO4 square planes along the c-direction, an almost linear Cu–O . . . O–Cu contact occurs with O . . . O distance of 2.757 Å (Figures 11c and 12b), slightly shorter than the van der Waals distance of 2.80 Å. Thus, this Cu–O . . . O–Cu spin exchange along the c-direction (Jc) should be substantial.

**Figure 11.** Crystal structure of α-CuV2O6 , where the blue spheres represent the Cu atoms, and the large and small red spheres the V and O atoms, respectively: (**a**) A CuO4 chain along the a-direction, which is made up of edge-sharing, axially-elongated, CuO6 octahedra. (**b**) A stack of CuO4 square planes along the a-direction, which results from the chain of edge-sharing CuO6 octahedra by removing the axial Cu-O bonds. (**c**) Stacks of CuO4 square planes running along the c-direction. (**d**) One sheet of CuO4 stack chains parallel to the a–b plane condensed with chains of corner-sharing VO4 tetrahedra on one side of the sheet. (**e**) Stacking of CuV2O6 layers forming α-CuV2O6 .

**Figure 12.** Spin exchange paths and spin lattice of α-CuV2O6 : (**a**) Spin exchange along the b-direction, Jb ; (**b**) Spin exchange along the c-direction, Jc; (**c**) Spin exchange along the a-direction, Ja; (**d**) Spin exchange along the (a + b)-direction, Ja+b; (**e**) Rectangular spin lattice made up of Jc and Ja+b.

Each sheet of CuO4 stack-chains parallel to the a–b plane is corner-shared with V2O6 chains, above and below the sheet (Figure 11d), to form layers of composition CuV2O6 (Figure 11e). Each V2O6 chain is made up of corner-sharing VO4 tetrahedra containing V5+ (d0 , S = 0) ions. As a consequence, adjacent CuO4 square planes are corner-shared with V2O6 chains (Figure 12c). The adjacent CuO4 square planes are not bridged by a single VO4 tetrahedron; therefore, the spin exchange along the a-direction, Ja, is expected to be weak. Along the (a + b)-direction, two adjacent CuO4 square planes are bridged by VO4 to form two Cu–O . . . V 5+ . . . O–Cu paths, in which the two Cu–O bonds in each path have a near-orthogonal arrangement.

As already discussed (Figure 10), this Cu–O . . . V 5+ . . . O–Cu spin exchange should be substantial due to the through-bond effect of the V5+ cation. Consequently, the spin lattice of α-CuV2O6 must be described by a two-dimensional (2D) rectangular lattice defined by Ja+b and Jc. In support of this analysis, the energy-mapping analysis based on DFT+U calculations with Ueff = 4 eV show that Ja+b and Jc are the two dominant spin exchanges, and are nearly equal in magnitude, namely, Ja+b = 86.8 K and Jc/Ja+b = 0.88 [19]. In agreement with this finding, one re-investigation of the magnetic properties of α-CuV2O6 clearly attested a 2D S = 1/2 rectangular spin lattice model with an anisotropy ratio of 0.7 [19]. The magnetic structure determined from neutron powder diffraction data was in best agreement with these DFT+U calculations. In terms of chemical bonding, α-CuV2O6 consists of CuV2O6 layers stacked along the c-direction. There is no chemical bonding between adjacent layers; only van der Waals interactions. In terms of magnetic bonding, however, α-CuV2O6 consists of 2D spin lattices parallel to the (a + b)-c plane. There is negligible magnetic bonding perpendicular to this plane.

#### *5.2. One-Dimensional Chain Behavior of LiCuVO4*

LiCuVO4 consists of axially elongated CuO6 octahedra, which form edge-sharing CuO4 chains along the b-direction, which are corner-shared with VO4 tetrahedra containing V 5+ (d0 , S = 0) ions (Figure 13a). When the axial Cu–O bonds are deleted, one finds CuVO4 layers in which the CuO2 ribbon chains are corner-shared by VO4 tetrahedra (Figure 13b). A perspective view of a single CuVO4 layer approximately along the c-direction (Figure 13c) shows that each VO4 tetrahedron bridges two neighboring CuO2 ribbon chains. Thus, the spin exchanges of interest are the nearest-neighbor exchanges Jnn of the Cu–O–Cu type and the next-nearest-neighbor spin exchange Jnnn of the Cu–O . . . O–Cu type within each CuO2 ribbon chain as well as the interchain spin exchange Ja along the a-direction of the Cu-O . . . V5+ . . . O-Cu type (Figure 13d).

In the absence of the V5+ ion, the exchange Ja would be similar in strength to Jnnn. However, the O . . . V 5+ . . . O bridges will weaken the strength of Ja, as discussed in Figure 7. Furthermore, there is no spin exchange path between adjacent CuVO4 layers. Consequently, the spin lattice of LiCuVO4 is a 1D chain running along the b-direction, as is the 1D ribbon chain. The major cause for the occurrence of the 1D chain character is the Cu–O . . . V 5+ . . . O–Cu spin exchange, which nearly vanishes because the effect of the through-space interaction is canceled by that of the through-bond interaction. In agreement with this reasoning, DFT+U calculations with Ueff = 4 eV show that Jnnn is strongly AFM (i.e., 208.7 K), while Jnn and Ja are weakly FM (i.e., Jnn/Jnnn = −0.12, and Ja/Jnnn = −0.08) [20].

As in the case of α-CuV2O6, the magnetic properties of LiCuVO4 were initially analyzed in terms of a 1D chain with uniform nearest-neighbor Heisenberg spin exchange [21–23]. The magnetic susceptibility of LiCuVO4 showed a typical broad maximum at about 28 K, characteristic for 1D behavior with strong intrachain spin exchanges, whereas long-range AFM ordering was detected only below ~2.2 K. The need to modify this simple 1D description was brought about by Gibson et al., who determined the magnetic structure of LiCuVO4 from single crystal neutron diffraction [24]. They found that the spins of each CuO2 ribbon chain had a cycloid structure (Figure 14a) with adjacent Cu2+ moments making an angle of slightly less than 90◦ . The incommensurate cycloid structure is explained by spin frustration due to competing exchange Jnn, which is FM, and Jnnn, which is AFM [20,24–26]. If all the spins of a CuO2 ribbon chain were to be collinear, the ribbon chain cannot satisfy all nearest-neighbor exchanges FM and all next-nearest-neighbors AFM simultaneously. Thus, the spin arrangement of the CuO2 ribbon chain is spin-frustrated. To reduce the extent of this spin frustration, the spins of the ribbon chain adopt a noncollinear spin arrangement. The noncollinear spin arrangement observed for LiCuVO4 has a cycloid structure in each ribbon chain, in which the nearest-neighbor spins are nearly orthogonal to each other, while the next-nearest-neighbor spins generate a near-AFM arrangement (Figure 14a) [24,25]. In a cycloid, each successive spin in the CuO2 ribbon chain rotates in one direction by a certain angle. Thus, a cycloid structure is chiral in nature, which means that the alternative cycloid structure opposite in chirality but identical in energy is equally probable (Figure 14b) [27].

**Figure 13.** (**a**) The crystal structure of LiCuVO4 viewed approximately along the b-direction, where the blue and green spheres represent the Cu and Li atoms, respectively, and the large and small red spheres represent the V and O atoms, respectively. (**b**) The CuVO4 lattice resulting from LiCuVO4 by removing the axial Cu–O bonds and Li atoms. (**c**) A perspective view of a single CuVO4 layer approximately along the c-direction. (**d**) The spin exchange paths present in a single CuVO4 layer, where the labels nn, nnn and a represent the spin exchanges Jnn, Jnnn and Ja, respectively.

**Figure 14.** (**a**,**b**) Two cycloids, which are opposite in chirality but are identical in energy.

In general, when the temperature is lowered below a certain temperature, TSDW, a moderately spin-frustrated magnetic system gives rise to two cycloids of opposite chirality with equal probability. The resulting superposition of the two (Figure 15a–c) leads to a state known as a spin density wave (SDW) [27,28]. The latter becomes transverse if the preferred spin orientation at each magnetic ion is perpendicular to the SDW propagation direction, but becomes longitudinal if the spin orientation prefers the SDW propagation direction

(Figure 15d–f). When the temperature is lowered further below TSDW, the electronic structure of the spin-lattice may relax to energetically favor one of the two chiral cycloids so that one can observe a cycloid state at a temperature slightly below TSDW. The repeat unit of a cycloid is determined by the spin frustration present in the magnetic system, therefore a cycloid phase is typically incommensurate.

**Figure 15.** The superposition of the two chiral cycloids in (**a**); (**b**) leads to the transverse SDW in (**c**). The superposition of the two chiral cycloids in (**d**); (**e**) leads to the longitudinal SDW in (**f**). For ease of illustration, the angle of the spin rotation was taken to be 45◦ .

In the cycloid state, LiCuVO4 exhibits ferroelectricity [29–31], because a cycloid structure lacks inversion symmetry. The polarization of LiCuVO4 can be switched with magnetic and electric fields [32–35]. LiCuVO4 has attracted special attention for the possibility of inducing new phases by applying external magnetic fields. Saturation of the Cu moments occurs above ~45 T, depending on the orientation of the crystal [36]. The competition between Jnn and Jnnn in LiCuVO4 is considered a promising setting for unusual bond nematicity and spin nematic phases close to the full magnetic saturation [37–45].

#### *5.3. Spin Gap Behavior from the CuO2Cl2 Perovskite Layer of (CuCl)LaNb2O7*

(CuCl)LaNb2O7 consists of CuClO2 and LaNb2O7 layers, which alternate along the c-direction by sharing their O corners [46]. Each LaNb2O7 layer represents two consecutive layers (Figure 16a) of LaNbO3 perovskite. The building blocks of the CuClO2 layer are the CuCl4O2 octahedra with their O atoms at the apical positions forming a linear O–Cu–O bond aligned along the c-direction. Suppose that the CuCl4 equatorial planes are square in shape with four-fold rotational symmetry around the O–Cu–O axis (Figure 16b,c), and such CuCl4O2 octahedra corner-share their Cl atoms to form a perovskite layer CuClO2. Then, the highest-occupied d-states of each CuCl4O2 octahedron are degenerate and have three electrons to accommodate, causing a Jahn–Teller instability of each CuCl4O2 octahedron. In the CuClO2 layer, the CuCl4O2 octahedra must undergo a cooperative Jahn–Teller distortion. In the tetragonal structure (SG, P4/mmm) of (CuCl)LaNb2O7 [46], each Cl site is split into four positions (Figure 16d). A Jahn–Teller distortion available to such a CuCl4O2 octahedron is an axial-elongation, in which one linear Cl–Cu–Cl bond is shortened while lengthening the other linear Cl–Cu–Cl bond, as shown in Figure 16e, which can be simplified as in Figure 16f. By choosing one of the four split positions from each Cl site, it is possible to construct the CuClO2 layer with a cooperative Jahn–Teller distortion, as presented in Figure 17a, which shows the Cu–Cl–Cu–Cl zigzag chains running along the b-direction with the plane of each CuCl2O2 rhombus perpendicular to the a–b plane. Each Cl site has four split positions when all these four possibilities occur equally.

With the local coordinate axes of each CuCl2O2 rhombus taken as in Figure 17a, then the magnetic orbital of the Cu2+ ion can be described as the x2–z2 state in which Cu x2–z2 orbital makes σ-antibonding interactions with the p-orbitals of O and Cl. Extended Hückel tight binding calculations [47] for the CuCl2O2 rhombus show that in the magnetic orbital

of the CuCl2O2 rhombus, the Cu2+ ion is described by 0.646 (3z2–r2 ) − 0.272 (x2–y2 ). The latter is rewritten as:

$$\begin{array}{rcl} 0.646(3\mathbf{z}^2 - \mathbf{r}^2) - 0.272(\mathbf{x}^2 - \mathbf{y}^2) & \approx & [2(3\mathbf{z}^2 - \mathbf{r}^2) - (\mathbf{x}^2 - \mathbf{y}^2)]/3 \\ & \approx & (\mathbf{z}^2 - \mathbf{x}^2) + 0.25(\mathbf{x}^2 - \mathbf{y}^2) \approx & (\mathbf{z}^2 - \mathbf{x}^2) \end{array}$$

**Figure 16.** (**a**) One LaNb2O7 layer of tetragonal (CuCl)LaNb2O7 , where the La atoms at the 12 coordinate sites formed by eight corner-sharing NbO6 octahedra are not shown for simplicity. (**b**) Perspective and (**c**) projection views of a CuCl4O2 octahedron with 4-fold rotational symmetry around the O–Cu–O axis aligned along the c-direction. (**d**) A projection view of a CuClO2 layer of tetragonal (CuCl)LaNb2O7 in which every Cl site is split into four positions. (**e**,**f**) Projection views of the Jahn–Teller distorted CuCl4O2 octahedron, where the axially elongated Cu–Cl bonds are shown in (**e**), but are not shown in (**f**) for simplicity.

**Figure 17.** The structures of the CuClO2 layer with cooperative Jahn–Teller distortion in (**a**) the tetragonal and (**b**) orthorhombic phases of (CuCl)LaNb2O7 . The Cl–Cu–Cl unit is linear in the tetragonal structure, but is slightly bent in the orthorhombic structure. The latter has an important consequence on the spin exchanges.

Namely, it is dominated by the (z2–x2 ) character. This is consistent with the NMR/NQR study of (CuCl)LaNb2O7, which showed that the d-state of the Cu2+ ion is mostly characterized by 3z2–r2 with some contribution of x2–y2 [48].

We note that the CuCl2O2 rhombuses of each Cu–Cl–Cu–Cl zigzag chain make Cu–Cl . . . Cl–Cu contacts of Cl . . . Cl = 4.03 Å with its adjacent zigzag chains (Figure 17a). The Cl p-orbitals of the (z2–x2 ) magnetic orbital (Figure 18a) are pointed approximately

along the Cl . . . Cl contact, so their overlap is substantial. In all other spin exchange paths, the Cl p-orbitals in their magnetic orbitals are not arranged to overlap well. Thus, only the exchange along the Cu–Cl . . . Cl–Cu direction is expected to be substantially AFM. Then the spin lattice of the CuClO2 layer would be 1D Heisenberg uniform AFM chains running along the Cu–Cl . . . Cl–Cu directions, e.g., the (a + 2b)- and (a − 2b)-directions (Figure 17a). Uniform AFM chains do not have a spin gap, but the magnetic properties of (CuCl)LaNb2O7 reveal a spin gap behavior [49]. Thus, the tetragonal structure of (CuCl)LaNb2O7 is inconsistent with experiment. If one discards the possibility of cooperative Jahn-Teller distortions, one can generate several nonuniform clusters made up of distorted CuCl4O2 octahedra [50]. However, the latter would lead to several different spin gaps rather than that observed by Kageyama et al. [49], who showed that the magnetic susceptibility of (CuCl)LaNb2O7 can be approximated by an isolated spin dimer model with the intradimer distance of approximately 8.8 Å, which corresponds to the fourth-nearest-neighbor Cu . . . Cu distance. The spin-gap behavior of (CuCl)LaNb2O7 was surprising, given the belief of the square lattice arrangement of Cu2+ ions is spin-frustrated [49,51,52]. This led to several DFT studies designed to find the precise crystal structure of (CuCl)LaNb2O7 [53–55], leading to the conclusion that an orthorhombic structure of space group Pbam is correct for (CuCl)LaNb2O7 [54,55].

**Figure 18.** (**a**) The magnetic orbital lying in the CuCl2O2 rhombus plane, which is best described as the x2–z2 state. (**b**–**d**) Three different types of the magnetic orbital arrangements found in the CuClO2 layers of orthorhombic CuCl(LaV2O7 ). In this idealized representation of the CuCl2O2 rhombus, the slight bending of the Cl–Cu–Cl linkage is neglected for clarity.

The cause for the spin gap behavior of (CuCl)LaNb2O7 is found from its orthorhombic structure (SG, *Pbam*) (Figure 17b) [56,57], in which the arrangement of the Cu2+ ions are no longer tetragonal so that adjacent Cu–Cl–Cu–Cl zigzag chains have two kinds of Cl . . . Cl contacts (i.e., 3.83 and 4.23 Å), and the Cu–Cl . . . Cl–Cu chains become alternating with shorter and longer Cl . . . Cl contacts. Furthermore, although the Cl–Cu–Cl unit of each CuCl2O2 rhombus is slightly bent in the orthorhombic structure, the latter has an important consequence on the spin exchanges (see below).

The six spin exchange paths J1–J6 of the CuClO2 layer are depicted in Figure 19a. In the spin exchanges J1 and J2, the Cl p-orbital tails are approximately pointed toward each other (Figure 18b). J3 is a Cu–Cl–Cu exchange with ∠Cu–Cl–Cu angle somewhat greater than 90◦ (namely, 109.0◦ ). In J4 and J5, the Cl p-orbital tails are approximately orthogonal to each other (Figure 18c), but they differ due to the bending in the Cl–Cu–Cl units. In the J6 path, the Cl p-orbital tails are not pointed toward to each other but are approximately