Spin Hamiltonians in Magnets: Theories and Computations

Abstract

The effective spin Hamiltonian method has drawn considerable attention for its power to explain and predict magnetic properties in various intriguing materials. In this review, we summarize different types of interactions between spins (hereafter, spin interactions, for short) that may be used in effective spin Hamiltonians as well as the various methods of computing the interaction parameters. A detailed discussion about the merits and possible pitfalls of each technique of computing interaction parameters is provided.

1. Introduction

The utilization of magnetism can date back to ancient China when the compass was invented to guide directions. Since the relationship between magnetism and electricity was revealed by Oersted, Lorentz, Ampere, Faraday, Maxwell, and others, more applications of magnetism have been invented, which include dynamos (electric generators), electric motors, cyclotrons, mass spectrometers, voltage transformers, electromagnetic relays, picture tubes, and sensing elements. During the information revolution, magnetic materials were extensively employed for information storage. The storage density, efficiency, and stability were substantially improved by the discovery and applications of the giant magnetoresistance effect [1,2], tunnel magnetoresistance [3,4,5,6,7,8,9], spin-transfer torques [10,11,12,13], etc. Recently, more and more novel magnetic states such as spin glasses [14,15], spin ice [16,17], spin liquid [18,19,20,21,22], and skyrmions [23,24,25,26,27,28] were found, revealing both theoretical and practical significance. For example, hedgehogs and anti-hedgehogs can be seen as the sources (monopoles) and the sinks (antimonopoles) of the emergent magnetic fields of topological spin textures [29], while magnetic skyrmions have shown promise as ultradense information carriers and logic devices [24].

To explain or predict the properties of magnetic materials, many models and methods have been invented. In this review, we will mainly focus on the effective spin Hamiltonian method based on first-principles calculations, and its applications in solid-state systems. In Section 2, we will introduce the effective spin Hamiltonian method. Firstly, in Section 2.1, the origin and the computing methods of the atomic magnetic moments are presented. Then, from Section 2.2, Section 2.3, Section 2.4, Section 2.5 and Section 2.6, different types of spin interactions that may be included in the spin Hamiltonians are discussed. Section 3 will discuss and compare various methods of computing the interaction parameters used in the effective spin Hamiltonians. In Section 4, we will give a brief conclusion of this review.

2. Effective Spin Hamiltonian Models

Though accurate, first-principles calculations are somewhat like black boxes (that is to say, they provide the final total results, such as magnetic moments and the total energy, but do not give a clear understanding of the physical results without further analysis), and have difficulties in dealing with large-scale systems or finite temperature properties. In order to provide an explicit explanation for some physical properties and improve the efficiency of thermodynamic and kinetic simulations, the effective Hamiltonian method is often adopted. In the context of magnetic materials where only the spin degree of freedom is considered, it can also be called the effective spin Hamiltonian method. Typically, the effective spin Hamiltonian models need to be carefully constructed and include all the possibly important terms; then the parameters of the models need to be calculated based on either first-principles calculations (see Section 3) or experimental data (see Section 3.3). Given the effective spin Hamiltonian and the spin configurations, the total energy of a magnetic system can be easily computed. Therefore, it is often adopted in Monte Carlo simulations [30] (or quantum Monte Carlo simulations) for assessing the total energy of many different configurations so that the finite temperature properties of magnetic materials can be studied. If the effects of atom displacements are taken into account, the effective Hamiltonians can also be applied to the spin molecular dynamics simulations [31,32,33], which is beyond the scope of this review.

In this review, we mainly focus on the classical effective spin Hamiltonian method where atomic magnetic moments (or spin vectors) are treated as classical vectors. In many cases, these classical vectors are assumed to be rigid so that their magnitudes keep constant during rotations. This treatment significantly simplifies the effective Hamiltonian models, and it is usually a good approximation, especially when atomic magnetic moments are large enough. In this part, we shall first introduce the origin of the atomic magnetic moments as well as the methods of computing atomic magnetic moments. Then different types of spin interactions will be discussed. 2.1. Atomic Magnetic Moments

The origin of atomic magnetic moments is explained by quantum mechanics. Suppose the quantized direction is the z-axis, an electron with quantum numbers (n, l ,_ml, _ms) leads to an orbital magnetic moment_μl= −_μBland a spin magnetic moment_μs= −_{ge μB}s, with their z components_μlz= −_{ml μB}and_μsz= −_{ge ms μB}, where_μB= |e|ℏ2mis the Bohr magneton and_ge≈2is the g-factor for a free electron. The energy of a magnetic momentμin a magnetic fieldB(magnetic induction) along z-direction is−μ⋅B = −_μzB.

Considering the Russell-Saunders coupling (also referred to as L-S coupling), which applies to most multi-electronic atoms, the total orbital magnetic moment and the total spin magnetic moment are_μL= −_μBLand_μS= −_{ge μB}S, respectively, whereL = ∑i_liandS = ∑i_siare the summation over electrons. The quantum numbers of each electron can usually be predicted by Hund’s rules. Owing to the spin–orbit interaction,LandSboth precess around the constant vectorJ = L+S. The time-averaged effective total magnetic moment isμ = −_{gJ μB}J, where

_gJ= 1+J(J+1)+S(S+1)−L(L+1)2J(J+1)

The atomic magnetic moments discussed above are based on the assumption that the atoms are isolated. Taking the influence of other atoms and external fields into account, the orbital interaction theory, the crystal field theory [34], or the ligand field theory [35,36] may be a better choice for theoretically predicting atomic magnetic moments. Notice that half-filled shells lead to a totalL = 0, and that in solids and molecules, orbital moments of electrons are usually quenched, resulting in an effectiveL = 0 [37] (counterexamples may be found for 4f elements or for 3d⁷ configurations as in Co(II)). Therefore typicallyμ = _μS, with its z component_μz= −_{ge Sz μB}(_Szis restricted to discrete values:S, S−1,…,−S). Usually, a nonzero_μSresults from singly filled (localized) d or f orbitals, while s and p orbitals are typically either doubly filled or vacant so that they have no direct contribution to the atomic magnetic moments. Therefore, when referring to atomic magnetic moments, usually only the atoms in the d- and f-transition series need to be considered.

The atomic magnetic moments can also be predicted numerically employing first-principles calculations. Nevertheless, we should notice that traditional Kohn–Sham density functional theory (DFT) calculations [38,39] (based on single-electron approximation) are not reliable for predicting atomic magnetic moments, and hence require the consideration of strong correlation effect among electrons, especially when dealing with localized d or f orbitals. Based on the Hubbard model [40], such problem can often be remedied by introducing an intra-atomic interaction with effective on-site Coulomb and exchange parameters, U and J [41] (or only one parameter^Ueff= U−J in Dudarev’s approach [42]). This approach is the DFT+U method [41,42,43,44], including LDA+U (LDA: Local density approximation), LSDA+U (LSDA: Local spin density approximation), GGA+U (GGA: Generalized gradient approximation), and so forth, where “+U” indicates the Hubbard “+U” correction. The parameters U and J can be estimated according to experience or semi-empirically by seeking agreement with experimental results of some specific properties, which is convenient but not very reliable. Considering how the values of the parameters U and J affect the prediction of atomic magnetic moments and other physical properties, we may need to compute these parameters more rigorously. A typical approach is constrained DFT calculations [43,45,46,47], where the local d or f charges are constrained to different values in several calculations, so that the parameters U and J can be obtained. Another approach based on constrained random-phase approximation (cRPA) [48,49,50] allows for considering the frequency (or energy) dependence of the parameters. More methods of computing U and J are summarized in Ref. [43]. There are more accurate approaches for dealing with strong correlated systems like DFT + Dynamical Mean Field Theory (DFT+DMFT) [51,52,53,54,55] and Reduced Density Matrix Functional Theory (RDMFT) [56,57], but they are much more sophisticated and computationally demanding so that they may be impractical for large-scale calculations. Wave function (WF) methods, such as Complete Active Space Self-Consistent Field (CASSCF) [58,59,60,61], Complete Active Space second-order Perturbation Theory (CASPT2) [62,63,64], Complete Active Space third-order Perturbation Theory (CASPT3) [65], and Difference Dedicated Configuration Interaction (DDCI) [66,67,68], are also widely adopted by theoretical chemists for studying magnetic properties of materials (especially molecules), including atomic magnetic moments and magnetic interactions. These WF methods are also more accurate but more computationally demanding than the DFT+U method, more detailed discussions of which can be found in Ref. [69].

2.2. Heisenberg Model

The simplest effective spin Hamiltonian model is the classical Heisenberg model, which can be reduced to Ising model or XY model. The classical Heisenberg model can be written as

_Hspin = ∑i,j>i_{Jij Si}⋅_Sj,

where_Siand_Sjindicate the total spin vectors on atomsiandj, and the summation is over all relevant pairs(ij). Its form was suggested by Heisenberg, Dirac, and Van Vleck. Such an interaction comes from the energy splitting between quantized parallel (ferromagnetic, FM; triplet state) and antiparallel (antiferromagnetic, AFM; singlet state) spin configurations._Jij>0and_Jij<0prefer AFM and FM configurations, respectively. There may be a difference in the factor such as −1 and12between different definitions of_Hspin, which is also the case in other models as will be discussed.

The spatial wave function of two electrons should possess the form of_ψ± = 12[_ψa(_r1)_ψb(_r2)±_ψb(_r1)_ψa(_r2)], where_ψaand_ψbare any single-electron spatial wave functions. A parallel triplet spin state and an antiparallel singlet spin state should correspond to an antisymmetric (_ψ−) and a symmetric (_ψ+) spatial wave function, respectively. For given_ψaand_ψb, the expectation value for the total energy of_ψ−can be different from that of_ψ+, which gives a preference to the AFM or FM spin configuration. Whether the AFM or FM spin configuration is preferred depends on the circumstances, and their energy difference can be described as a Heisenberg term_{J12 S1}⋅_S2.

In the simple case of an H₂ molecule, an AFM singlet spin state is preferred, whose symmetric spatial wave function_ψ+ corresponds to a bonding state [70,71]. However, this leads to the total magnetic moment of zero because the two antiparallel electrons share the same spatial state. In another simple case, where_ψaand_ψbstands for two degenerate and orthogonal orbitals of the same atom, an FM triplet spin state is preferred, which is in agreement with Hund’s rules. Consider a set of orthogonal Wannier functions with_ϕnλ(r −_rα)resembling the nth atomic orbital with spinλcentered at theαth lattice site, and suppose there are Nh electrons each localized on one of the N lattice sites, each ion possessing h unpaired electrons. If these h electrons have the same exchange integrals with all the other electrons, the interaction resulting from the antisymmetrization of the wave functions can be expressed as

_Hex = ∑α^α′n^n′_Jn^n′(_rα,_r^α′ )[14+S(_rα)⋅S(_r^α′ )]

which is called the Heisenberg exchange interaction [72]. After removing the constant terms, we can see such an interaction has the form of_Hspin = ∑i,j>i_{Jij Si}⋅_Sj.

Based on molecular orbital analysis using_ϕaand_ϕbto denote the singly filled d orbitals of the two spin-12 magnetic ions (i.e., d⁹ ions), Hay et al. [73] showed that the exchange interaction between the two ions can be approximately expressed as

J = −2_Kab+^∆2Ueff = _JF+_JAF

where

_Kab∝⟨_ϕa(1)_ϕb(2)|1_r12|_ϕb(1)_ϕa(2)⟩ = ∫_ϕa*(_r1)_ϕb*(_r2)1_r12ϕb(_r1)_ϕa(_r2)>0

^Ueff =_Jaa−_Jab∝⟨_ϕa(1)_ϕa(2)|1_r12|_ϕa(1)_ϕa(2)⟩−⟨_ϕa(1)_ϕb(2)|1_r12|_ϕa(1)_ϕb(2)⟩>0

and∆indicates the energy gap between the bonding state and antibonding state constructed by_ϕaand_ϕb. The two components of J have opposite signs, i.e.,_JF = −2_Kab<0and_JAF = ^∆2Ueff>0, which give preference to FM and AFM spin configurations, respectively. For general dⁿ cases, more orbitals should be considered. Therefore the expression ofJ will be more complicated, but the exchange interaction can still be similarly decomposed into FM and AFM contributions [73]. An application of this analysis is that when calculating exchange parameterJusing Dudarev’s approach of DFT+U with a parameter^Ueff, the calculated value ofJshould vary with^Ueffapproximately asJ = _JF+^∆2Ueffwith_JFand^∆2 to be fitted [74]. However, this is no longer correct when^Ueff→0.

Another mechanism that leads to FM spin configurations is the double exchange, in which the interaction between two magnetic ions is induced by spin coupling to mobile electrons that travels from one ion to another. A mobile electron has lower energy if the localized spins are aligned. Such a mechanism is essential in metallic systems containing ions with variable charge states [75,76].

The superexchange is another important indirect exchange mechanism, where the interaction between two transition-metal (TM) ions is induced by spin coupling to two electrons on a non-magnetic ligand (L) ion that connects them, forming an exchange path of TM-L-TM type. Different mechanisms were proposed to explain the superexchange interaction. In Anderson’s mechanism [77,78], the superexchange results from virtual processes in which an electron is transferred from the ligand to one of the neighboring magnetic ions, and then another electron on the ligand couples with the spin of the other magnetic ion through exchange interaction. In Goodenough’s mechanism [79,80], the concept of semicovalent bonds was invented, where only one electron given by the ligand predominates in a semicovalent bond. Because of the exchange forces between the electrons on the magnetic ion and the electron given by the ligand, the ligand electron with its spin parallel to the net spin of the magnetic ion will spend more time on the magnetic ion than that with an antiparallel spin if the d orbital of the magnetic ion is less than half-filled, and vice versa. The magnetic atom and the ligand are supposed to be connected by a semicovalent bond or a covalent bond when they are near, or by an ionic bond (or possibly a metallic-like bond) otherwise. The superexchange interaction with semicovalent bonds existing is also called semicovalent exchange interaction. Kanamori summarized the dependence of the sign of the superexchange parameter (whether FM or AFM) on bond angle, bond type and number of d electrons (in different mechanisms), which is often referred as Goodenough–Kanamori (GK) rules [80,81,82]. For the 180° (bond angle) case, generally, AFM interaction is expected between cations of the same kind (counterexamples may exist for d⁴ cases such as Mn³⁺-Mn³⁺, where the sign depends on the direction of the line of superexchange), and FM interaction is expected between two cations with more-than-half-filled and less-than-half-filled d-shells, respectively [81]. For the 90° case, the results are usually the opposite [81]. A schematic diagram of superexchange interactions (between cations both with more-than-half-filled d-shell) is given in Figure 1. More details of the discussions can be found in Ref. [81] and Ref. [82].

A counterexample of the GK rules can be found in the layered magnetic topological insulator MnBi₂Te₄, which possesses intrinsic ferromagnetism [83]. In contrast, the prediction of the GK rules leads to a weak AFM exchange interaction between Mn ions. In Ref. [84], the presence of Bi³⁺ was found to be essential for explaining this anomaly: d⁵ ions in TM-L-TM spin-exchange paths would prefer FM coupling if the empty p orbitals of a nonmagnetic cation M (which is Bi³⁺ ion in the case of MnBi₂Te₄) hybridize strongly with those of the ligand L (but AFM coupling otherwise). Oleś et al. [85] pointed out that the GK rules may not be obeyed in transition metal compounds with orbital degrees of freedom (e.g., d¹ and d² electronic configurations) due to spin-orbital entanglement.

Exchange interactions between two TM ions also take place through the exchange paths of TM-L…L-TM type [86], referred to as super-superexchanges, where TM ions do not share a common ligand. Each TM ion of a solid forms a TML_n polyhedron (typically, n = 3–6) with the surrounding ligands L, and the unpaired spins of the TM ion are accommodated in the singly filled d-states of TML_n. Since each d-state has a d-orbital of TM combined out-of-phase with the p-orbitals of L, the unpaired spin of TM does not reside solely on the d-orbital of TM, as assumed by Goodenough and Kanamori, but is delocalized into the p-orbitals of the surrounding ligands L. Thus, TM-L…L-TM type exchanges occur and can be strongly AFM when their L…L contact distances are in the vicinity of the van der Waals distance so that the ligand p-orbitals overlap well across the L…L contact.

Another mechanism is the indirect coupling of magnetic moments by conduction electrons, referred to as Ruderman–Kittel–Kasuya–Yosida (RKKY) interaction [87,88,89,90]. This kind of interaction between two spins_S1and_S2is also proportional to_S1⋅_S2with an expression

HRKKY∝∑qχ(q)eiq⋅r21 S1⋅S2.

The magnetic dipole–dipole interaction (between magnetic moments_μ1and_μ2located at different atoms) with energy

V = 1^R3[_μ1⋅_μ2−3(R^⋅_μ1)(R^⋅_μ2)]

also has contributions to the bilinear term, but it is typically much weaker than the exchange interactions in most solid-state materials such as iron and cobalt. The characteristic temperature of dipole–dipole interaction (or termed as “dipolar interaction”) in magnetic materials is typically of the order of 1 K, above which no long-range order can be stabilized by such an interaction [37]. However, in some cases, such as in several single-molecule magnets (SMMs), the exchange interactions can be so weak that they are comparable to or weaker than dipolar interactions, thus the dipolar interactions must not be neglected [91].

For most magnetic materials, the Heisenberg interaction is the most predominant spin interaction. As a result, the simple classical Heisenberg model is able to explain the magnetic properties such as the ground states of spin configurations (FM or AFM) and the transition temperatures (Curie temperature for FM states or Néel temperature for AFM states) for many magnetic materials.

If some pairs of spins favor FM spin configurations while other pairs favor AFM configurations, frustration may occur, leading to more complicated and more interesting noncollinear spin configurations. For example, the FM effects of double exchange resulting from mobile electrons in some antiferromagnetic lattices give rise to a distortion of the ground-state spin arrangement and lead to a canted spin configuration [92]. A magnetic solid with moderate spin frustration lowers its energy by adopting a noncollinear superstructure (e.g., a cycloid or a helix) in which the moments of the ions are identical in magnitude but differ in orientation or a collinear magnetic superstructure (e.g., a spin density wave, SDW) in which the moments of the ions differ in magnitude but identical in orientation [93,94]. For a cycloid formed in a chain of magnetic ions, each successive spin rotates in one direction by a certain angle, so there are two opposite ways of rotating the successive spins hence producing two cycloids opposite in chirality but identical in energy. When these two cycloids occur with equal probability below a certain temperature, their superposition leads to a SDW [93,94]. On lowering the temperature further, the electronic structure of the spin-lattice relaxes to energetically favor one of the two chiral cycloids so that one can observe a cycloid state. The latter, being chiral, has no inversion symmetry and gives rise to ferroelectricity [95]. The spin frustration is also a potential driving force for topological states like skyrmions and hedgehogs [29].

2.3. The J Matrices and Single-Ion Anisotropy

The classical Heisenberg model can be generalized to a matrix form to include all the possible second-order interactions between two spins (or one spin itself):

_Hspin = ∑i,j>i_{SiT Jij Sj}+∑i_{SiT Ai Si}

where_Jijand_Aiare 3 × 3 matrices called the J matrix and single-ion anisotropy (SIA) matrix. The_Jijmatrix can be decomposed into three parts: The isotropic Heisenberg exchange parameter_Jij = (_Jij,xx+_Jij,yy+_Jij,zz)/3as in the classical Heisenberg model, the antisymmetric Dzyaloshinskii–Moriya interaction (DMI) matrix_Dij = (_Jij−_JijT)/2 [96,97,98], and the symmetric (anisotropic) Kitaev-type exchange coupling matrix_Kij = (_Jij+_JijT)/2−_JijI(whereIdenotes a 3 × 3 identity matrix). Thus_Jij =_JijI+_Dij+_Kij [99,100].

Now we analyze the possible origin of these terms by means of symmetry analysis. When considering interaction potential between (or among) spins, we should notice that the total interaction energy should be invariant under time inversion ({_Sk}→{−_Sk}). Therefore, any odd order term in the spin Hamiltonian should be zero unless an external magnetic field is present when a term−∑i_μi⋅B= ∑i_{ge μB}B⋅_Sishould be added to the effective spin Hamiltonian. Ignoring the external magnetic field, the spin Hamiltonian should only contain even order terms, with the lowest order of significance being the second-order (the zeroth-order term is a constant and therefore not necessary). If the spin-orbit coupling (SOC) is negligible, the total effective spin Hamiltonian_Hspinshould be invariant under any global spin rotations, therefore_Hspinshould be expressed by only inner product terms of spins like terms proportional to_Si⋅_Sj,(_Si⋅_Sj)(_Sk⋅_Sl)and so on. That is to say, when SOC is negligible, the second-order terms in the_Hspinshould only include the classical Heisenberg term∑i,j>i_{Jij Si}⋅_Sj, which implies that those interactions described by_Ai,_Dij, and_Kijmatrices all originate from SOC (_HSO=λS^⋅L^). What is more, if the spatial inversion symmetry is satisfied by the lattice,_Jijshould be equal to_JijTso that there will be no DMI (_Dij=0). That is to say, the DMI can only exist where the spatial inversion symmetry is broken.

The SIA matrix_Aihas only six independent components and is usually assumed to be symmetric. If we suppose the magnitude of the classical spin vector_Sito be independent of its direction, the isotropic part13(_Ai,xx+_Ai,yy+_Ai,zz)Iwould be of no significance, and therefore_Aiwould have only five independent components after subtracting the isotropic part from itself. It is evident that the_{SiT Ai Si}prefers the direction of_Sialong the eigenvector of_Aiwith the lowest eigenvalue. If this lowest eigenvalue is two-fold degenerate, the directions of_Sifavored by SIA will be those belonging to the plane spanned by the two eigenvectors that share the lowest eigenvalue, in which case we say the ion i has easy-plane anisotropy. On the contrary, if the lowest eigenvalue is not degenerate while the higher eigenvalue is two-fold degenerate, we say the ion i has easy-axis anisotropy. In these two cases (easy-plane or easy-axis anisotropy), by defining the direction of z-axis parallel to the nondegenerate eigenvector, the_{SiT Ai Si}part would be simplified to_Ai,zz ^(_Siz)2 with only one independent component. The easy-axis anisotropy has been found to be helpful in stabilizing the long-range magnetic order and enhancing the Curie temperature in two-dimensional or quasi-two-dimensional systems [101]. The easy-plane anisotropy in three-dimensional ferromagnets can lead to the effect called “quantum spin reduction”, where the mean spin at zero temperature has a value lower than the maximal one due to the quantum fluctuations [101,102]. Recently, several materials with unusually large easy-plane or easy-axis anisotropy were found [103,104,105], which, as single-ion magnets (SIMs), are promising for applications such as high-density information storage, spintronics, and quantum computing.

The DMI matrix_Dijis antisymmetric and therefore has only three independent components, which can be expressed by a vector_Dijwith_Dij,x=_Dij,yz,_Dij,y=_Dij,zx, and_Dij,z=_Dij,xy. Thus, the DMI can be expressed by cross product:_{SiT Dij Sj}=_Dij⋅(_Si×_Sj). Such an interaction prefers the vectors_Siand_Sjto be orthogonal to each other, with a rotation (of_Sjrelative to_Si) around the direction of−_Dij. Together with Heisenberg term_{Jij Si}⋅_Sj, the preferred rotation angle between_Siand_Sjrelative to the collinear state preferred by the Heisenberg term would bearctan|_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,108,109,110]. Materials with skyrmion states induced by DMI usually have a large ratio of|_D1||_J1| (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)₃ (X = Br or Cl) Janus monolayers (e.g., for Cr(I,Br)₃, supposing|_Si|=32for any i, the corresponding interaction parameters are computed as_J1=−1.800meV and|_D1|=0.270meV, thus|_D1||_J1|=0.150 ), though monolayers such as CrI₃ 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_Kijhas five independent components as a symmetric matrix with zero trace. For the specific cases when_Siand_Sjare parallel to each other (pointing in the same direction),_{SiT Kij Sj}would perform like_{SiT Ai Si}and show preference to the direction with the lowest eigenvalue of_Kij; while when_Siand_Sjare antiparallel to each other,_{SiT Kij Sj}would prefer the direction with the highest eigenvalue of_Kij. The difference between the highest and lowest eigenvalue of_Kijcan be defined as_Kij(a scalar), which characterizes the anisotropic contribution of_Kij. Generally, the favorite direction of the spins is decided by both SIA and Kitaev interactions. The long-range ferromagnetic order in monolayer CrI₃ was explained by the anisotropic superexchange interaction since the Cr-I-Cr bond angle is close to90° [111]. In Ref. [112], the interplay between the prominent Kitaev interaction and SIA was studied to explain the different magnetic behaviors of CrI₃ and CrGeTe₃ naturally. For CrI₃, supposing|_Si|=32for any i, the_Jijand_Kijparameters between nearest pairs are computed as −2.44 and 0.85 meV, respectively; while the only independent component_Ai,zzof_Ai is −0.26 meV. For CrGeTe₃, 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α -RuCl₃ [114,115,116] and (Na_1-xLi_x)₂IrO₃ [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., CrSiTe₃ and CrGeTe₃ [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_{Jij Si}⋅_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

_Hspin=∑i,j>i_{Jij Si}⋅_Sj+∑i,j>i_Kij ^(_Si⋅_Sj)2+∑i,j,k>j_Kijk(_Si⋅_Sj)(_Si⋅_Sk)+∑i,j>i,k>i,l>k_Kijkl(_Si⋅_Sj)(_Sk⋅_Sl)

The biquadratic terms have been found to be important in many systems, such as MnO [120,121], YMnO₃ [74], TbMnO₃ [122], iron-based superconductor KFe_1.5Se₂ [123], and 2D magnets [124]. In the case of TbMnO₃, 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 TbMnO₃ is provided in the supplementary material of Ref. [125]). According to Ref. [126], in a Heisenberg chain system constructed from alternatingS>12andS=12 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 transition-metal 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

_{κijk^{i′ j′ k′}CC}[_Si⋅(_Sj×_Sk)][_S^i′ ⋅(_S^j′ ×_S^k′ )]

whose local part has the form

_κijkCC ^{[_Si⋅(_Sj×_Sk)]2}

and spin–chiral interactions (SCI) with the form

_κijkSC(_τijk⋅_Si)[_Si⋅(_Sj×_Sk)]

where the unit vector_τijk∝(_Rj−_Ri)×(_Rk−_Ri)is the surface normal of the oriented triangle spanned by the lattice sites_Ri,_Rjand_Rk. The local scalar spin chirality_χijk=_Si⋅(_Sj×_Sk)among triplets of spins can be interpreted as a fictitious effective magnetic field^Beff∝_{χijk τijk}, which leads to topological orbital moments^LTO (TOM) [129,130,131,132,133,134] arising from the orbital current of electrons hopping around the triangles. The TOM is defined as

_LiTO=∑(jk)_LijkTO=∑(jk)_{κijkTO χijk τijk}

where_κijkTOis 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_LiTO⋅^Beff . 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 three-dimensional topological lattice (possibly built up with hedgehogs and anti-hedgehogs) experimentally [135].

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 3d 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_{Jijk Si}⋅(_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,144,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

_Hspin=_E0+∑n=1∞_Jn(∑nth nearest pairs ⟨k,l⟩_ek⋅_el)+∑^n′=1∞_J^n′′∑^n′th type (_ek⋅_el)(_em⋅_en)+⋯

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,_Ai,xy, (_Ai,yy−_Ai,xx) and (_Ai,zz−_Ai,xx ) can be obtained by first-principles calculations for four specified spin states [148]. Taking the isotropic Heisenberg parameter_Jijfor 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_Jijcan be expressed as

_Jij=_Eij,↑↑+_Eij,↓↓−_Eij,↑↓−_Eij,↓↑4^S2

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

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

_Hspin=∑i,j>i_{SiT Jij Sj}+∑i_{SiT Ai Si}

and each component of the matrix_Jijand_Aican also be obtained by first-principles calculations for four specified spin states. To compute_Jij,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≠a,c≠b) with an appropriate reference configuration and kept unchanged. Then,_Jij,abcan be expressed as

_Jij,ab=_{Eij,ab,↑↑}+_{Eij,ab,↓↓}−_{Eij,ab,↑↓}−_{Eij,ab,↓↑}4^S2

To compute_Ai,ab(a,b=x,y,zwitha≠b), we use_Ei,ab,αβ(α,β= ↑,↓) to denote the energy of the configuration where spin i is parallel to the direction whose a component is±22(forα= ↑or↓, respectively), b component is±22(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≠a,c≠b) with an appropriate reference configuration. Then,_Ai,abcan be expressed as

_Ai,ab=_Ei,ab,↑↑+_Ei,ab,↓↓−_Ei,ab,↑↓−_Ei,ab,↓↑4^S2

To compute(_Ai,aa−_Ai,bb)(a,b=x,y,zwitha≠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≠a,c≠b) with an appropriate reference configuration. Then(_Ai,aa−_Ai,bb)can be expressed as

_Ai,aa−_Ai,bb=_Ei,ab,↑0+_Ei,ab,↓0−_Ei,ab,0↑−_Ei,ab,0↓4^S2

It is easy to verify that, by employing this four-state method, each component of the_Jijand_Aican 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_Jijand_Aiare still rigorous. For_Ai, we can find out that the effects of all these terms are correctly cancelled. For_Jij, 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≠i,j) will interfere with the calculation of_Jijbecause(_Sk⋅_Sl)is constant during the calculation and therefore mixed with the contribution of_Jij(_Si⋅_Sj). The error of the calculated_Jijoriginated from four-spin ring coupling terms

∑i,j>i,k>i,l>k_Kijkl(_Si⋅_Sj)(_Sk⋅_Sl)

is given by

∑k≠i or jl≠i or j(l>k)_Kijkl(_Sk⋅_Sl)

while there is no easy way to get rid of this problem perfectly. Other parts of the_Jij(including_Dijand_Kij) 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 TbMnO₃, 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_Jijvalues (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_Jijis 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_Jijusing 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.

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_Jijbetween_Siand_Sjis actually the sum over_Ji^j′with any lattice vectorR=_r^j′ −_rj. Therefore, to get rid 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_Kijin the term_Kij ^(_Si⋅_Sj)2, we can let_Sipointing to the(1,0,0)direction,_Sjpointing to the(1,0,0),(−1,0,0),(12,12,0)and(−12,−12,0)directions, with other spins parallel to the z-axis. The corresponding total energies are denoted as_E1,_E2,_E3, and_E4, respectively. Then the_Kijcan be expressed as

_Kij=_E1+_E2−_E3−_E4

It can be easily checked that the effects of other terms not higher than fourth order are totally eliminated. Therefore, this approach of calculating_Kijshould 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,123,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 YMnO₃, 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,123,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_Jijis vulnerable to the effects of other fourth-order 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 fourth-order 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 Fe³⁺ ions with nonmagnetic Al³⁺ 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_Siand_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), wherepis 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.

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,152,153,154,155,156,157,158,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

_{ℍimj^m′σ^σ′}(R)=⟨_ψimσ(r)|H|_ψimσ(r+R)⟩

_{Simj^m′σ^σ′}(R)=⟨_ψimσ(r)|_ψimσ(r+R)⟩

ℍ(k)=∑Rℍ(R)^eik⋅R

S(k)=∑RS(R)^eik⋅R

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

G(k,ε)=^{(εS(k)−ℍ(k))−1}

and

G(R,ε)=_∫BZG(k,ε)^e−ik⋅Rdk

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:

δE=_∫−∞EFεδn(ε)dε=−_∫−∞EFδN(ε)dε

where

n(ε)=−1πImTr(G(ε))

and

N(ε)=−1πImTr(ε−ℍ)

where traces are taken over orbitals. By defining

_ℙi=_ℍii(R=0)=_??i0??+_??→i⋅σ→

with its component

_ℙim^m′=_pim^m′0??+_{pim^m′ e→im^m′}⋅σ→

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

_Gim,j^m′=_Gim,j^m′0??+_{G→im,j^m′}⋅σ→

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

δ_Eij=−2π_∫∞EFImTr(δ_ℍiGδ_ℍjG)dε

withδ_ℍi=δ_ϕi×_pi. After mathematical simplification, the expression ofδ_Eijcan be mapped to that given by the effective Hamiltonian model

_Hspin=∑i_{SiT Ai Si}+∑i,j>i[_{Jij Si}⋅_Sj+_Dij⋅(_Si×_Sj)+_{SiT Kij Sj}]+∑i,j>i_Kij ^(_Si⋅_Sj)2

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

_Jij=Im(_Aij00−_Aijxx−_Aijyy−_Aijzz−2_{Aijzz Siref}⋅_Sjref)

_Kijuv=Im(_Aijuv+_Aijvu)

_Diju=Re(_Aij0u−_Aiju0)

_Bij=Im(_Aijzz)

where

_Aijuv=1π_∫−∞EFTr{_{??iz Giju ??jz Gjiv}}dε

withu,v∈{0,x,y,z}(the trace is also taken over orbitals), and_Sirefindicates 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_Jijis 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δEto 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_{Jij Si}⋅_Sjcan be expressed as_{Jij Si}⋅_Sj=_{Jij Si ei}⋅_{Sj ej}; the magnitudes of the spins_Siand_Sjdepend on the angles between these two spins or neighboring spin directions (e.g.,_ek). Therefore,

_{Jij Si}⋅_Sj=_{J˜0, ij ei}⋅_ej+_C1 ^(_ei⋅_ej)2+_C2(_ei⋅_ek)(_ei⋅_ej)+_C3(_ej⋅_ek′)(_ei⋅_ej)+⋯.

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., single-molecule 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,164,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,169,170,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, four-spin 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 higher-order 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).

Figure 1. A schematic diagram of superexchange interactions between transition-metal (TM) ions both with more-than-half-filled d-shell. According to Goodenough-Kanamori (GK) rules, the 180° and the 90° cases favor antiferromagnetic (AFM) and ferromagnetic (FM) arrangements of TM ions, respectively. The main difference is whether the two electrons of L occupy the same p orbital, leading to different tendencies for the alignments of the two electrons of L that interact with two TM ions.

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 3. Schematic diagrams of the four states adopted in the four-state method for computingJij . 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.

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 5. Flow charts of several methods of variable selection. In these flow charts, the fitting varianceσ2estimated by the training set and the testing set are denoted asσ2^andσ'2^ , 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 areNλ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).

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.; writing-review 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

1. Baibich, M.N.; Broto, J.M.; Fert, A.; Vandau, F.N.; Petroff, F.; Eitenne, P.; Creuzet, G.; Friederich, A.; Chazelas, J. Giant magnetoresistance of (001)Fe/(001) Cr agnetic superlattices. Phys. Rev. Lett. 1988, 61, 2472-2475.

2. Binasch, G.; Grunberg, P.; Saurenbach, F.; Zinn, W. Enhanced magnetoresistance in layered magnetic-structures with antiferromagnetic interlayer exchange. Phys. Rev. B 1989, 39, 4828-4830.

3. Julliere, M. Tunneling between ferromagnetic-films. Phys. Lett. A 1975, 54, 225-226.

4. Butler, W.H.; Zhang, X.G.; Schulthess, T.C.; MacLaren, J.M. Spin-dependent tunneling conductance of Fe vertical bar MgO vertical bar Fe sandwiches. Phys. Rev. B 2001, 63, 054416.

5. Mathon, J.; Umerski, A. Theory of tunneling magnetoresistance of an epitaxial Fe/MgO/Fe(001) junction. Phys. Rev. B 2001, 63, 220403.

6. Bowen, M.; Cros, V.; Petroff, F.; Fert, A.; Boubeta, C.M.; Costa-Kramer, J.L.; Anguita, J.V.; Cebollada, A.; Briones, F.; de Teresa, J.M.; et al. Large magnetoresistance in Fe/MgO/FeCo(001) epitaxial tunnel junctions on GaAs(001). Appl. Phys. Lett. 2001, 79, 1655-1657.

7. Yuasa, S.; Nagahama, T.; Fukushima, A.; Suzuki, Y.; Ando, K. Giant room-temperature magnetoresistance in single-crystal Fe/MgO/Fe magnetic tunnel junctions. Nat. Mater. 2004, 3, 868-871.

8. Parkin, S.S.P.; Kaiser, C.; Panchula, A.; Rice, P.M.; Hughes, B.; Samant, M.; Yang, S.H. Giant tunnelling magnetoresistance at room temperature with MgO (100) tunnel barriers. Nat. Mater. 2004, 3, 862-867.

9. Ikeda, S.; Hayakawa, J.; Ashizawa, Y.; Lee, Y.M.; Miura, K.; Hasegawa, H.; Tsunoda, M.; Matsukura, F.; Ohno, H. Tunnel magnetoresistance of 604% at 300 K by suppression of Ta diffusion in CoFeB/MgO/CoFeB pseudo-spin-valves annealed at high temperature. Appl. Phys. Lett. 2008, 93, 082508.

10. Ralph, D.C.; Stiles, M.D. Spin transfer torques. J. Magn. Magn. Mater. 2008, 320, 1190-1216.

11. Jonietz, F.; Muehlbauer, S.; Pfleiderer, C.; Neubauer, A.; Muenzer, W.; Bauer, A.; Adams, T.; Georgii, R.; Boeni, P.; Duine, R.A.; et al. Spin Transfer Torques in MnSi at Ultralow Current Densities. Science 2010, 330, 1648-1651.

12. Pai, C.-F.; Liu, L.; Li, Y.; Tseng, H.W.; Ralph, D.C.; Buhrman, R.A. Spin transfer torque devices utilizing the giant spin Hall effect of tungsten. Appl. Phys. Lett. 2012, 101, 122404.

13. Mellnik, A.R.; Lee, J.S.; Richardella, A.; Grab, J.L.; Mintun, P.J.; Fischer, M.H.; Vaezi, A.; Manchon, A.; Kim, E.A.; Samarth, N.; et al. Spin-transfer torque generated by a topological insulator. Nature 2014, 511, 449-451.

14. Edwards, S.F.; Anderson, P.W. Theory of Spin Glasses. J. Phys. F-Met. Phys. 1975, 5, 965-974.

15. Binder, K.; Young, A.P. Spin-glasses-experimental facts, theoretical concepts, and open questions. Rev. Mod. Phys. 1986, 58, 801-976.

16. Bramwell, S.T.; Gingras, M.J.P. Spin ice state in frustrated magnetic pyrochlore materials. Science 2001, 294, 1495-1501.

17. Castelnovo, C.; Moessner, R.; Sondhi, S.L. Magnetic monopoles in spin ice. Nature 2008, 451, 42-45.

18. Anderson, P.W. Resonating valence bonds: A new kind of insulator? Mater. Res. Bull. 1973, 8, 153-160.

19. Shimizu, Y.; Miyagawa, K.; Kanoda, K.; Maesato, M.; Saito, G. Spin liquid state in an organic Mott insulator with a triangular lattice. Phys. Rev. Lett. 2003, 91, 107001.

20. Balents, L. Spin liquids in frustrated magnets. Nature 2010, 464, 199-208.

21. Yan, S.; Huse, D.A.; White, S.R. Spin-Liquid Ground State of the S = 1/2 Kagome Heisenberg Antiferromagnet. Science 2011, 332, 1173-1176.

22. Han, T.-H.; Helton, J.S.; Chu, S.; Nocera, D.G.; Rodriguez-Rivera, J.A.; Broholm, C.; Lee, Y.S. Fractionalized excitations in the spin-liquid state of a kagome-lattice antiferromagnet. Nature 2012, 492, 406-410.

23. Skyrme, T.H.R. A unifield field theory of mesons and baryons. Nucl. Phys. 1962, 31, 556-569.

24. Fert, A.; Cros, V.; Sampaio, J. Skyrmions on the track. Nat. Nanotechnol. 2013, 8, 152-156.

25. Muehlbauer, S.; Binz, B.; Jonietz, F.; Pfleiderer, C.; Rosch, A.; Neubauer, A.; Georgii, R.; Boeni, P. Skyrmion Lattice in a Chiral Magnet. Science 2009, 323, 915-919.

26. Neubauer, A.; Pfleiderer, C.; Binz, B.; Rosch, A.; Ritz, R.; Niklowitz, P.G.; Boeni, P. Topological Hall Effect in the A Phase of MnSi. Phys. Rev. Lett. 2009, 102, 186602.

27. Pappas, C.; Lelievre-Berna, E.; Falus, P.; Bentley, P.M.; Moskvin, E.; Grigoriev, S.; Fouquet, P.; Farago, B. Chiral Paramagnetic Skyrmion-like Phase in MnSi. Phys. Rev. Lett. 2009, 102, 197202.

28. Roessler, U.K.; Bogdanov, A.N.; Pfleiderer, C. Spontaneous skyrmion ground states in magnetic metals. Nature 2006, 442, 797-801.

29. Fujishiro, Y.; Kanazawa, N.; Tokura, Y. Engineering skyrmions and emergent monopoles in topological spin crystals. Appl. Phys. Lett. 2020, 116, 090501.

30. Hastings, W.K. Monte-carlo sampling methods using markov chains and their applications. Biometrika 1970, 57, 97-109.

31. Ma, P.-W.; Woo, C.H.; Dudarev, S.L. Large-scale simulation of the spin-lattice dynamics in ferromagnetic iron. Phys. Rev. B 2008, 78, 024434.

32. Ma, P.-W.; Dudarev, S.L.; Woo, C.H. Spin-lattice-electron dynamics simulations of magnetic materials. Phys. Rev. B 2012, 85, 184301.

33. Omelyan, I.P.; Mryglod, I.M.; Folk, R. Algorithm for molecular dynamics simulations of spin liquids. Phys. Rev. Lett. 2001, 86, 898-901.

34. Van Vleck, J.H. Theory of the variations in paramagnetic anisotropy among different salts of the iron group. Phys. Rev. 1932, 41, 208-215.

35. Ballhausen, C.J. Introduction to Ligand Field Theory; McGraw-Hill: New York, NY, USA, 1962.

36. Griffith, J.S. The Theory of Transition-Metal Ions; Cambridge University Press: Cambridge, UK, 1964.

37. Mattis, D.C. The Theory of Magnetism Made Simple: An. Introduction to Physical Concepts and to Some Useful Mathematical Methods; World Scientific Publishing Company: Singapore, 2006.

38. Hohenberg, P.; Kohn, W. Inhomogeneous electron gas. Phys. Rev. B 1964, 136, B864.

39. Kohn, W.; Sham, L.J. Self-consistent equations including exchange and correlation effects. Phys. Rev. 1965, 140, A1133.

40. Hubbard, J. Electron correlations in narrow energy bands. Proc. R. Soc. Lond. Ser. A-Math. Phys. Sci. 1963, 276, 238-257.

41. Liechtenstein, A.I.; Anisimov, V.I.; Zaanen, J. Density-functional theory and strong-interactions-orbital ordering in mott-hubbard insulators. Phys. Rev. B 1995, 52, R5467-R5470.

42. Dudarev, S.L.; Botton, G.A.; Savrasov, S.Y.; Humphreys, C.J.; Sutton, A.P. Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study. Phys. Rev. B 1998, 57, 1505-1509.

43. Himmetoglu, B.; Floris, A.; de Gironcoli, S.; Cococcioni, M. Hubbard-Corrected DFT Energy Functionals: The LDA + U Description of Correlated Systems. Int. J. Quantum Chem. 2014, 114, 14-49.

44. The VASP Manual-Vaspwiki. Available online: https://www.vasp.at/wiki/index.php/The_VASP_Manual (accessed on 3 February 2021).

45. Dederichs, P.H.; Blugel, S.; Zeller, R.; Akai, H. Ground-states of constrained systems-application to cerium impurities. Phys. Rev. Lett. 1984, 53, 2512-2515.

46. Hybertsen, M.S.; Schluter, M.; Christensen, N.E. Calculation of coulomb-interaction parameters for la2cuo4 using a constrained-density-functional approach. Phys. Rev. B 1989, 39, 9028-9041.

47. Gunnarsson, O.; Andersen, O.K.; Jepsen, O.; Zaanen, J. Density-functional calculation of the parameters in the anderson model-application to Mn in CdTe. Phys. Rev. B 1989, 39, 1708-1722.

48. Springer, M.; Aryasetiawan, F. Frequency-dependent screened interaction in Ni within the random-phase approximation. Phys. Rev. B 1998, 57, 4364-4368.

49. Kotani, T. Ab initio random-phase-approximation calculation of the frequency-dependent effective interaction between 3d electrons: Ni, Fe, and MnO. J. Phys. Condens. Matter 2000, 12, 2413-2422.

50. Aryasetiawan, F.; Imada, M.; Georges, A.; Kotliar, G.; Biermann, S.; Lichtenstein, A.I. Frequency-dependent local interactions and low-energy effective models from electronic structure calculations. Phys. Rev. B 2004, 70, 195104.

51. Metzner, W.; Vollhardt, D. Correlated lattice fermions in dinfinity dimensions. Phys. Rev. Lett. 1989, 62, 324-327.

52. Georges, A.; Kotliar, G. Hubbard-model in infinite dimensions. Phys. Rev. B 1992, 45, 6479-6483.

53. Georges, A.; Kotliar, G.; Krauth, W.; Rozenberg, M.J. Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions. Rev. Mod. Phys. 1996, 68, 13-125.

54. Pavarini, E.; Koch, E.; Vollhardt, D.; Lichtenstein, A. The LDA+ DMFT approach to strongly correlated materials. Reihe Modeling Simul. 2011, 1, 2-13.

55. Embedded DMFT Functional Tutorials. Available online: http://hauleweb.rutgers.edu/tutorials/ (accessed on 3 February 2021).

56. Yang, W.T.; Zhang, Y.K.; Ayers, P.W. Degenerate ground states and a fractional number of electrons in density and reduced density matrix functional theory. Phys. Rev. Lett. 2000, 84, 5172-5175.

57. Sharma, S.; Dewhurst, J.K.; Lathiotakis, N.N.; Gross, E.K.U. Reduced density matrix functional for many-electron systems. Phys. Rev. B 2008, 78, 201103.

58. Roos, B.O.; Taylor, P.R.; Siegbahn, P.E.M. A complete active space scf method (casscf) using a density-matrix formulated super-ci approach. Chem. Phys. 1980, 48, 157-173.

59. Siegbahn, P.E.M.; Almlof, J.; Heiberg, A.; Roos, B.O. The complete active space scf (casscf) method in a newton-raphson formulation with application to the HNO molecule. J. Chem. Phys. 1981, 74, 2384-2396.

60. Knowles, P.J.; Sexton, G.J.; Handy, N.C. Studies using the casscf wavefunction. Chem. Phys. 1982, 72, 337-347.

61. Snyder, J.W., Jr.; Parrish, R.M.; Martinez, T.J. alpha-CASSCF: An Efficient, Empirical Correction for SA-CASSCF To Closely Approximate MS-CASPT2 Potential Energy Surfaces. J. Phys. Chem. Lett. 2017, 8, 2432-2437.

62. Roos, B.O.; Linse, P.; Siegbahn, P.E.M.; Blomberg, M.R.A. A simple method for the evaluation of the 2nd-order perturbation energy from external double-excitations with a casscf reference wavefunction. Chem. Phys. 1982, 66, 1-2, 197-207.

63. Andersson, K.; Malmqvist, P.A.; Roos, B.O.; Sadlej, A.J.; Wolinski, K. Second-order perturbation theory with a CASSCF reference function. J. Phys. Chem. 1990, 94, 5483-5488.

64. Andersson, K.; Malmqvist, P.A.; Roos, B.O. 2nd-order perturbation-theory with a complete active space self-consistent field reference function. J. Chem. Phys. 1992, 96, 1218-1226.

65. Werner, H.J. Third-order multireference perturbation theory-The CASPT3 method. Mol. Phys. 1996, 89, 645-661.

66. Miralles, J.; Daudey, J.P.; Caballol, R. Variational calculation of small energy differences-The singlet-triplet gap in Cu2Cl6 2. Chem. Phys. Lett. 1992, 198, 555-562.

67. Miralles, J.; Castell, O.; Caballol, R.; Malrieu, J.P. Specific ci calculation of energy differences-transition energies and bond-energies. Chem. Phys. 1993, 172, 33-43.

68. Malrieu, J.P.; Daudey, J.P.; Caballol, R. Multireference self-consistent size-consistent singles and doubles configuration-interaction for ground and excited-states. J. Chem. Phys. 1994, 101, 8908-8921.

69. Malrieu, J.P.; Caballol, R.; Calzado, C.J.; de Graaf, C.; Guihery, N. Magnetic Interactions in Molecules and Highly Correlated Materials: Physical Content, Analytical Derivation, and Rigorous Extraction of Magnetic Hamiltonians. Chem. Rev. 2014, 114, 429-492.

70. Heitler, W.; London, F. Wechselwirkung neutraler Atome und homöopolare Bindung nach der Quantenmechanik. Z. Für Phys. 1927, 44, 6-7, 455-472.

71. Levine, D.S.; Head-Gordon, M. Clarifying the quantum mechanical origin of the covalent chemical bond. Nat. Commun. 2020, 11, 4893.

72. White, R.M.; White, R.M.; Bayne, B. Quantum Theory of Magnetism; Springer: Berlin/Heidelberg, Germany, 1983.

73. Hay, P.J.; Thibeault, J.C.; Hoffmann, R. Orbital interactions in metal dimer complexes. J. Am. Chem. Soc. 1975, 97, 4884-4899.

74. Novak, P.; Chaplygin, I.; Seifert, G.; Gemming, S.; Laskowski, R. Ab-initio calculation of exchange interactions in YMnO(3). Comput. Mater. Sci. 2008, 44, 79-81.

75. Zener, C. Interaction between the d-shells in the transition metals. 2. ferromagnetic compounds of manganese with perovskite structure. Phys. Rev. 1951, 82, 403-405.

76. Anderson, P.W.; Hasegawa, H. Considerations on double exchange. Phys. Rev. 1955, 100, 675-681.

77. Anderson, P.W. Antiferromagnetism-Theory of superexchange interaction. Phys. Rev. 1950, 79, 350-356.

78. Kramers, H. L'interaction entre les atomes magnétogènes dans un cristal paramagnétique. Physica 1934, 1, 182-192.

79. Goodenough, J.B.; Loeb, A.L. Theory of ionic ordering, crystal distortion, and magnetic exchange due to covalent forces in spinels. Phys. Rev. 1955, 98, 391-408.

80. Goodenough, J.B. Theory of the role of covalence in the perovskite-type manganites la,m(ii) mno3. Phys. Rev. 1955, 100, 564-573.

81. Kanamori, J. Superexchange interaction and symmetry properties of electron orbitals. J. Phys. Chem. Solids 1959, 10, 87-98.

82. Goodenough, J.B. Magnetism and the Chemical Bond; Interscience Publishers: New York, NY, USA, 1963; Volume 1.

83. Gong, C.; Li, L.; Li, Z.; Ji, H.; Stern, A.; Xia, Y.; Cao, T.; Bao, W.; Wang, C.; Wang, Y.; et al. Discovery of intrinsic ferromagnetism in two-dimensional van der Waals crystals. Nature 2017, 546, 265-269.

84. Li, J.; Ni, J.Y.; Li, X.Y.; Koo, H.J.; Whangbo, M.H.; Feng, J.S.; Xiang, H.J. Intralayer ferromagnetism between S = 5/2 ions in MnBi2Te4: Role of empty Bi p states. Phys. Rev. B 2020, 101, 201408.

85. Oles, A.M.; Horsch, P.; Feiner, L.F.; Khaliullin, G. Spin-orbital entanglement and violation of the Goodenough-Kanamori rules. Phys. Rev. Lett. 2006, 96, 147205.

86. Whangbo, M.H.; Koo, H.J.; Dai, D. Spin exchange interactions and magnetic structures of extended magnetic solids with localized spins: Theoretical descriptions on formal, quantitative and qualitative levels. J. Solid State Chem. 2003, 176, 417-481.

87. Ruderman, M.A.; Kittel, C. Indirect exchange coupling of nuclear magnetic moments by conduction electrons. Phys. Rev. 1954, 96, 99-102.

88. Zener, C. Interaction between the d-shells in the transition metals. Phys. Rev. 1951, 81, 440-444.

89. Kasuya, T. A Theory of metallic ferromagnetism and antiferromagnetism on zeners model. Prog. Theor. Phys. 1956, 16, 45-57.

90. Yosida, K. Magnetic properties of cu-mn alloys. Phys. Rev. 1957, 106, 893-898.

91. Zhang, L.; Zhang, Y.-Q.; Zhang, P.; Zhao, L.; Guo, M.; Tang, J. Single-Molecule Magnet Behavior Enhanced by Synergic Effect of Single-Ion Anisotropy and Magnetic Interactions. Inorg. Chem. 2017, 56, 7882-7889.

92. Degennes, P.G. Effects of double exchange in magnetic crystals. Phys. Rev. 1960, 118, 141-154.

93. Gordon, E.E.; Derakhshan, S.; Thompson, C.M.; Whangho, M.-H. Spin-Density Wave as a Superposition of Two Magnetic States of Opposite Chirality and Its Implications. Inorg. Chem. 2018, 57, 9782-9785.

94. Koo, H.-J.; PN, R.S.; Orlandi, F.; Sundaresan, A.; Whangbo, M.-H. On Ferro- and Antiferro-Spin-Density Waves Describing the Incommensurate Magnetic Structure of NaYNiWO6. Inorg. Chem. 2020, 59, 17856-17859.

95. Xiang, H.J.; Whangbo, M.H. Density-functional characterization of the multiferroicity in spin spiral chain cuprates. Phys. Rev. Lett. 2007, 99, 257203.

96. Dzyaloshinsky, I. A thermodynamic theory of "weak" ferromagnetism of antiferromagnetics. J. Phys. Chem. Solids 1958, 4, 241-255.

97. Moriya, T. Anisotropic superexchange interaction and weak ferromagnetism. Phys. Rev. 1960, 120, 91-98.

98. Moriya, T. New mechanism of anisotropic superexchange interaction. Phys. Rev. Lett. 1960, 4, 228-230.

99. Xu, C.; Xu, B.; Dupe, B.; Bellaiche, L. Magnetic interactions in BiFeO3: A first-principles study. Phys. Rev. B 2019, 99, 104420.

100. Xu, C.; Feng, J.; Prokhorenko, S.; Nahas, Y.; Xiang, H.; Bellaiche, L. Topological spin texture in Janus monolayers of the chromium trihalides Cr(I, X)(3). Phys. Rev. B 2020, 101, 060404.

101. Spirin, D.V. Magnetization of two-layer ferromagnets with easy-plane and/or easy-axis single-site anisotropy. J. Magn. Magn. Mater. 2003, 264, 121-129.

102. Onufrieva, F.P. Magnetic property anomalies of ferromagnets with mono-ionic anisotropy. Fiz. Tverd. Tela 1984, 26, 3435-3437.

103. Zadrozny, J.M.; Xiao, D.J.; Atanasov, M.; Long, G.J.; Grandjean, F.; Neese, F.; Long, J.R. Magnetic blocking in a linear iron(I) complex. Nat. Chem. 2013, 5, 577-581.

104. Liu, J.; Chen, Y.-C.; Liu, J.-L.; Vieru, V.; Ungur, L.; Jia, J.-H.; Chibotaru, L.F.; Lan, Y.; Wernsdorfer, W.; Gao, S.; et al. A Stable Pentagonal Bipyramidal Dy(III) Single-Ion Magnet with a Record Magnetization Reversal Barrier over 1000 K. J. Am. Chem. Soc. 2016, 138, 5441-5450.

105. Deng, Y.-F.; Wang, Z.; Ouyang, Z.-W.; Yin, B.; Zheng, Z.; Zheng, Y.-Z. Large Easy-Plane Magnetic Anisotropy in a Three-Coordinate Cobalt(II) Complex Li(THF)(4) Co(NPh2)(3). Chem. A Eur. J. 2016, 22, 14821-14825.

106. Antal, A.; Lazarovits, B.; Udvardi, L.; Szunyogh, L.; Ujfalussy, B.; Weinberger, P. First-principles calculations of spin interactions and the magnetic ground states of Cr trimers on Au(111). Phys. Rev. B 2008, 77, 174429.

107. Bak, P.; Jensen, M.H. Theory of helical magnetic-structures and phase-transitions in mnsi and fege. J. Phys. C Solid State Phys. 1980, 13, L881-L885.

108. Yu, X.Z.; Kanazawa, N.; Onose, Y.; Kimoto, K.; Zhang, W.Z.; Ishiwata, S.; Matsui, Y.; Tokura, Y. Near room-temperature formation of a skyrmion crystal in thin-films of the helimagnet FeGe. Nat. Mater. 2011, 10, 106-109.

109. Huang, S.X.; Chien, C.L. Extended Skyrmion Phase in Epitaxial FeGe(111) Thin Films. Phys. Rev. Lett. 2012, 108, 267201.

110. Weber, T.; Waizner, J.; Tucker, G.S.; Georgii, R.; Kugler, M.; Bauer, A.; Pfleiderer, C.; Garst, M.; Boeni, P. Field dependence of nonreciprocal magnons in chiral MnSi. Phys. Rev. B 2018, 97, 224403.

111. Lado, J.L.; Fernandez-Rossier, J. On the origin of magnetic anisotropy in two dimensional CrI3. 2d Materials 2017, 4, 035002.

112. Xu, C.; Feng, J.; Xiang, H.; Bellaiche, L. Interplay between Kitaev interaction and single ion anisotropy in ferromagnetic CrI3 and CrGeTe3 monolayers. Npj Comput. Mater. 2018, 4, 57.

113. Kitaev, A. Anyons in an exactly solved model and beyond. Ann. Phys. 2006, 321, 2-111.

114. Banerjee, A.; Bridges, C.A.; Yan, J.Q.; Aczel, A.A.; Li, L.; Stone, M.B.; Granroth, G.E.; Lumsden, M.D.; Yiu, Y.; Knolle, J.; et al. Proximate Kitaev quantum spin liquid behaviour in a honeycomb magnet. Nat. Mater. 2016, 15, 733-740.

115. Baek, S.H.; Do, S.H.; Choi, K.Y.; Kwon, Y.S.; Wolter, A.U.B.; Nishimoto, S.; van den Brink, J.; Buechner, B. Evidence for a Field-Induced Quantum Spin Liquid in alpha-RuCl3. Phys. Rev. Lett. 2017, 119, 037201.

116. Zheng, J.; Ran, K.; Li, T.; Wang, J.; Wang, P.; Liu, B.; Liu, Z.-X.; Normand, B.; Wen, J.; Yu, W. Gapless Spin Excitations in the Field-Induced Quantum Spin Liquid Phase of alpha-RuCl3. Phys. Rev. Lett. 2017, 119, 227208.

117. Cao, G.; Qi, T.F.; Li, L.; Terzic, J.; Cao, V.S.; Yuan, S.J.; Tovar, M.; Murthy, G.; Kaul, R.K. Evolution of magnetism in the single-crystal honeycomb iridates (Na1-xLix)(2)IrO3. Phys. Rev. B 2013, 88, 220414.

118. Manni, S.; Choi, S.; Mazin, I.I.; Coldea, R.; Altmeyer, M.; Jeschke, H.O.; Valenti, R.; Gegenwart, P. Effect of isoelectronic doping on the honeycomb-lattice iridate A(2)IrO(3). Phys. Rev. B 2014, 89, 245113.

119. Xu, C.; Feng, J.; Kawamura, M.; Yamaji, Y.; Nahas, Y.; Prokhorenko, S.; Qi, Y.; Xiang, H.; Bellaiche, L. Possible Kitaev Quantum Spin Liquid State in 2D Materials with S = 3/2. Phys. Rev. Lett. 2020, 124, 087205.

120. Harris, E.A.; Owen, J. Biquadratic exchange between Mn2+ ions in MgO. Phys. Rev. Lett. 1963, 11, 9.

121. Rodbell, D.S.; Harris, E.A.; Owen, J.; Jacobs, I.S. Biquadratic exchange and behavior of some antiferromagnetic substances. Phys. Rev. Lett. 1963, 11, 10.

122. Fedorova, N.S.; Ederer, C.; Spaldin, N.A.; Scaramucci, A. Biquadratic and ring exchange interactions in orthorhombic perovskite manganites. Phys. Rev. B 2015, 91, 165122.

123. Zhu, H.-F.; Cao, H.-Y.; Xie, Y.; Hou, Y.-S.; Chen, S.; Xiang, H.; Gong, X.-G. Giant biquadratic interaction-induced magnetic anisotropy in the iron-based superconductor A(x)Fe(2-y)Se(2). Phys. Rev. B 2016, 93, 024511.

124. Kartsev, A.; Augustin, M.; Evans, R.F.L.; Novoselov, K.S.; Santos, E.J.G. Biquadratic exchange interactions in two-dimensional magnets. Npj Comput. Mater. 2020, 6, 150.

125. Li, X.-Y.; Lou, F.; Gong, X.-G.; Xiang, H. Constructing realistic effective spin Hamiltonians with machine learning approaches. New J. Phys. 2020, 22, 053036.

126. Ivanov, N.B.; Ummethum, J.; Schnack, J. Phase diagram of the alternating-spin Heisenberg chain with extra isotropic three-body exchange interactions. Eur. Phys. J. B 2014, 87, 226.

127. Paul, S.; Haldar, S.; von Malottki, S.; Heinze, S. Role of higher-order exchange interactions for skyrmion stability. Nat. Commun. 2020, 11, 4756.

128. Grytsiuk, S.; Hanke, J.P.; Hoffmann, M.; Bouaziz, J.; Gomonay, O.; Bihlmayer, G.; Lounis, S.; Mokrousov, Y.; Bluegel, S. Topological-chiral magnetic interactions driven by emergent orbital magnetism. Nat. Commun. 2020, 11, 511.

129. Hoffmann, M.; Weischenberg, J.; Dupe, B.; Freimuth, F.; Ferriani, P.; Mokrousov, Y.; Heinze, S. Topological orbital magnetization and emergent Hall effect of an atomic-scale spin lattice at a surface. Phys. Rev. B 2015, 92, 020401.

130. Hanke, J.P.; Freimuth, F.; Nandy, A.K.; Zhang, H.; Bluegel, S.; Mokrousov, Y. Role of Berry phase theory for describing orbital magnetism: From magnetic heterostructures to topological orbital ferromagnets. Phys. Rev. B 2016, 94, 121114.

131. Dias, M.D.S.; Bouaziz, J.; Bouhassoune, M.; Bluegel, S.; Lounis, S. Chirality-driven orbital magnetic moments as a new probe for topological magnetic structures. Nat. Commun. 2016, 7, 13613.

132. Hanke, J.-P.; Freimuth, F.; Bluegel, S.; Mokrousov, Y. Prototypical topological orbital ferromagnet gamma-FeMn. Sci. Rep. 2017, 7, 41078.

133. Dias, M.D.S.; Lounis, S. Insights into the orbital magnetism of noncollinear magnetic systems. Proc. SPIE 2017, 10357, 103572A.

134. Lux, F.R.; Freimuth, F.; Bluegel, S.; Mokrousov, Y. Engineering chiral and topological orbital magnetism of domain walls and skyrmions. Commun. Phys. 2018, 1, 60.

135. Tanigaki, T.; Shibata, K.; Kanazawa, N.; Yu, X.; Onose, Y.; Park, H.S.; Shindo, D.; Tokura, Y. Real-Space Observation of Short-Period Cubic Lattice of Skyrmions in MnGe. Nano Lett. 2015, 15, 5438-5442.

136. Brinker, S.; Dias, M.D.S.; Lounis, S. The chiral biquadratic pair interaction. New J. Phys. 2019, 21, 083015.

137. Laszloffy, A.; Rozsa, L.; Palotas, K.; Udvardi, L.; Szunyogh, L. Magnetic structure of monatomic Fe chains on Re(0001): Emergence of chiral multispin interactions. Phys. Rev. B 2019, 99, 184430.

138. Parihari, D.; Pati, S.K. Effect of chiral interactions in frustrated magnetic chains. Phys. Rev. B 2004, 70, 180403.

139. Kostyrko, T.; Bulka, B.R. Canonical perturbation theory for inhomogeneous systems of interacting fermions. Phys. Rev. B 2011, 84, 035123.

140. Mankovsky, S.; Polesya, S.; Ebert, H. Extension of the standard Heisenberg Hamiltonian to multispin exchange interactions. Phys. Rev. B 2020, 101, 174401.

141. Wen, X.G.; Wilczek, F.; Zee, A. Chiral spin states and superconductivity. Phys. Rev. B 1989, 39, 11413-11423.

142. Bauer, B.; Cincio, L.; Keller, B.P.; Dolfi, M.; Vidal, G.; Trebst, S.; Ludwig, A.W.W. Chiral spin liquid and emergent anyons in a Kagome lattice Mott insulator. Nat. Commun. 2014, 5, 5137.

143. Drautz, R.; Fahnle, M. Spin-cluster expansion: Parametrization of the general adiabatic magnetic energy surface with ab initio accuracy. Phys. Rev. B 2004, 69, 104404.

144. Singer, R.; Faehnle, M. Construction of basis functions for the spin-cluster expansion of the magnetic energy on the atomic scale in rotationally invariant systems. J. Math. Phys. 2006, 47, 113503.

145. Singer, R.; Dietermann, F.; Faehnle, M. Spin Interactions in bcc and fcc Fe beyond the Heisenberg Model. Phys. Rev. Lett. 2011, 107, 017204.

146. Reinhardt, P.; Habas, M.P.; Dovesi, R.; Moreira, I.D.; Illas, F. Magnetic coupling in the weak ferromagnet CuF2. Phys. Rev. B 1999, 59, 1016-1023.

147. Illas, F.; Moreira, I.D.R.; de Graaf, C.; Barone, V. Magnetic coupling in biradicals, binuclear complexes and wide-gap insulators: A survey of ab initio wave function and density functional theory approaches. Theor. Chem. Acc. 2000, 104, 265-272.

148. Xiang, H.; Lee, C.; Koo, H.-J.; Gong, X.; Whangbo, M.-H. Magnetic properties and energy-mapping analysis. Dalton Trans. 2013, 42, 823-853.

149. Xiang, H.J.; Kan, E.J.; Wei, S.-H.; Whangbo, M.H.; Gong, X.G. Predicting the spin-lattice order of frustrated systems from first principles. Phys. Rev. B 2011, 84, 224429.

150. Weingart, C.; Spaldin, N.; Bousquet, E. Noncollinear magnetism and single-ion anisotropy in multiferroic perovskites. Phys. Rev. B 2012, 86, 094413.

151. Liechtenstein, A.I.; Katsnelson, M.I.; Antropov, V.P.; Gubanov, V.A. Local spin-density functional-approach to the theory of exchange interactions in ferromagnetic metals and alloys. J. Magn. Magn. Mater. 1987, 67, 65-74.

152. Katsnelson, M.I.; Lichtenstein, A.I. First-principles calculations of magnetic interactions in correlated systems. Phys. Rev. B 2000, 61, 8906-8912.

153. Katsnelson, M.I.; Kvashnin, Y.O.; Mazurenko, V.V.; Lichtenstein, A.I. Correlated band theory of spin and orbital contributions to Dzyaloshinskii-Moriya interactions. Phys. Rev. B 2010, 82, 100403.

154. Lounis, S.; Dederichs, P.H. Mapping the magnetic exchange interactions from first principles: Anisotropy anomaly and application to Fe, Ni, and Co. Phys. Rev. B 2010, 82, 180404.

155. Szilva, A.; Costa, M.; Bergman, A.; Szunyogh, L.; Nordstrom, L.; Eriksson, O. Interatomic Exchange Interactions for Finite-Temperature Magnetism and Nonequilibrium Spin Dynamics. Phys. Rev. Lett. 2013, 111, 127204.

156. Korotin, D.M.; Mazurenko, V.V.; Anisimov, V.I.; Streltsov, S.V. Calculation of exchange constants of the Heisenberg model in plane-wave-based methods using the Green's function approach. Phys. Rev. B 2015, 91, 224405.

157. He, X.; Helbig, N.; Verstraete, M.J.; Bousquet, E. TB2J: A Python Package for Computing Magnetic Interaction Parameters. Available online: https://arxiv.org/abs/2009.01910 (accessed on 3 February 2021).

158. Wan, X.; Yin, Q.; Savrasov, S.Y. Calculation of magnetic exchange interactions in Mott-Hubbard systems. Phys. Rev. Lett. 2006, 97, 266403.

159. Boukhvalov, D.W.; Dobrovitski, V.V.; Katsnelson, M.I.; Lichtenstein, A.I.; Harmon, B.N.; Kogerler, P. Electronic structure and exchange interactions in V-15 magnetic molecules: LDA+U results. Phys. Rev. B 2004, 70, 054417.

160. Mostofi, A.A.; Yates, J.R.; Lee, Y.-S.; Souza, I.; Vanderbilt, D.; Marzari, N. wannier90: A tool for obtaining maximally-localised Wannier functions. Comput. Phys. Commun. 2008, 178, 685-699.

161. Pizzi, G.; Vitale, V.; Arita, R.; Bluegel, S.; Freimuth, F.; Geranton, G.; Gibertini, M.; Gresch, D.; Johnson, C.; Koretsune, T.; et al. Wannier90 as a community code: New features and applications. J. Phys. Condens. Matter 2020, 32, 165902.

162. Wang, X.D.; Wang, D.S.; Wu, R.Q.; Freeman, A.J. Validity of the force theorem for magnetocrystalline anisotropy. J. Magn. Magn. Mater. 1996, 159, 337-341.

163. Flygare, W.H. Magnetic-interactions in molecules and an analysis of molecular electronic charge distribution from magnetic parameters. Chem. Rev. 1974, 74, 653-687.

164. Brown, J.M.; Sears, T.J. Reduced form of the spin-rotation hamiltonian for asymmetric-top molecules, with applications to ho2 and nh2. J. Mol. Spectrosc. 1979, 75, 111-133.

165. Lahti, P.M.; Ichimura, A.S. Semiempirical study of electron exchange interaction in organic high-spin pi-systems-classifying structural effects in organic magnetic molecules. J. Org. Chem. 1991, 56, 3030-3042.

166. Chattopadhyay, S.; Lenz, B.; Kanungo, S.; Panda, S.K.; Biermann, S.; Schnelle, W.; Manna, K.; Kataria, R.; Uhlarz, M.; Skourski, Y.; et al. Pronounced 2/3 magnetization plateau in a frustrated S = 1 isolated spin-triangle compound: Interplay between Heisenberg and biquadratic exchange interactions. Phys. Rev. B 2019, 100, 094427.

167. Takano, K.; Sano, K. Determination of exchange parameters from magnetic susceptibility. J. Phys. Soc. Jpn. 1997, 66, 1846-1847.

168. Turek, I.; Kudrnovsky, J.; Drchal, V.; Bruno, P. Exchange interactions, spin waves, and transition temperatures in itinerant magnets. Philos. Mag. 2006, 86, 1713-1752.

169. Lancon, D.; Ewings, R.A.; Guidi, T.; Formisano, F.; Wildes, A.R. Magnetic exchange parameters and anisotropy of the quasi-two-dimensional antiferromagnet NiPS3. Phys. Rev. B 2018, 98, 134414.

170. Chen, L.; Chung, J.-H.; Gao, B.; Chen, T.; Stone, M.B.; Kolesnikov, A.I.; Huang, Q.; Dai, P. Topological Spin Excitations in Honeycomb Ferromagnet CrI3. Phys. Rev. X 2018, 8, 041028.

171. Oba, M.; Nakamura, K.; Akiyama, T.; Ito, T.; Weinert, M.; Freeman, A.J. Electric-Field-Induced Modification of the Magnon Energy, Exchange Interaction, and Curie Temperature of Transition-Metal Thin Films. Phys. Rev. Lett. 2015, 114, 107202.

172. Nishino, M.; Yoshioka, Y.; Yamaguchi, K. Effective exchange interactions and magnetic phase transition temperatures in Prussian blue analogs: A study by density functional theory. Chem. Phys. Lett. 1998, 297, 51-59.