Methods of Modeling of Strongly Correlated Electron Systems

The discovery of high-Tc superconductivity in cuprates in 1986 moved strongly correlated systems from exotic worlds interesting only for pure theorists to the focus of solid-state research. In recent decades, the majority of hot topics in condensed matter physics (high-Tc superconductivity, colossal magnetoresistance, multiferroicity, ferromagnetism in diluted magnetic semiconductors, etc.) have been related to strongly correlated transition metal compounds. The highly successful electronic structure calculations based on density functional theory lose their predictive power when applied to such compounds. It is necessary to go beyond the mean field approximation and use the many-body theory. The methods and models that were developed for the description of strongly correlated systems are reviewed together with the examples of response function calculations that are needed for the interpretation of experimental information (inelastic neutron scattering, optical conductivity, resonant inelastic X-ray scattering, electron energy loss spectroscopy, angle-resolved photoemission, electron spin resonance, and magnetic and magnetoelectric properties). The peculiarities of (quasi-) 0-, 1-, 2-, and 3- dimensional systems are discussed.


Introduction
The strongly correlated transition metal compounds remain a focus of attention in the condensed matter scientific community since they may show high-T c superconductivity in quasi-two-dimensional cuprates [1,2], frustrated magnetism in low-dimensional cuprates and other materials [3][4][5][6], colossal magnetoresistance and observation of Griffiths phase in manganites [7][8][9][10], multiferroism [11,12], or spin liquid behavior [13][14][15], as well as many other interesting properties. The interpretation of the vast experimental information depends on the calculation of various response functions. The task becomes highly nontrivial for the transition metal compounds because of the strong correlations in the dor f -electron shell of the transition metal ions. The difficulty is due to the necessity to take into account the many-body effects. These effects cannot be described by finite orders of the perturbation theory because of high (macroscopically large) degeneracy of the unperturbed state. Thus, the summation of terms up to the infinite order of perturbation theory or non-perturbative theoretical methods should be used [16].
The commonly used approach has serious limitations: for an interpretation of an experiment, the response function is calculated for a model, the parameters of the model being taken from a fit to experimental data. The model thus describes the dynamics of the studied system in an energy range that is relevant for the calculated response. However, the interrelation between responses on different energy scales is missing in this approach. In particular, the models describing the optical response (the energy range of several eV) has usually nothing in common with the spin-Hamiltonians that describe the magnetic properties of the same compound (the energy range of several meV).
The aim of this review is to outline a systematic theoretical approach for the calculation of physical properties and response functions for different kinds of transition metal As we mentioned in the Introduction, the description of transition metal compounds demands the account of correlation in the motion of d and f -electrons [16]. The simplest system that demonstrates the role of correlation effects is the electronic structure of the hydrogen molecule. We adopt the Born-Oppenheimer approximation, where the nuclei positions are fixed in the points R 1 , R 2 . So, we consider the motion of two electrons that experience an attraction by protons and a repulsion from each other.
In the mean field Hückel theory [17][18][19][20], it is assumed that every electron moves in a self-consistent external field V SCF (r) that is a sum of the attraction to the nuclei and the repulsion from the second electron averaged over the positions of this electron. On the basis of orthogonalized functions of the hydrogen atom ground state φ i (r) = φ(r − R i ), the mean-field theory Hamiltonian has the form (in the second quantization) where the diagonal and non-diagonal matrix elements are (t > 0) operator a † i,s creates an electron with orbital function φ i (r) and spin index s =↑, ↓. Without loss of generality, we may assume the orthogonality of basis functions φ i (r)|φ j (r) = δ ij ; then creation and annihilation operators satisfy usual fermionic commutation relations a i,s a † j,s + a † j,s a i,s = δ ij (for the generalization for the non-orthgonal basis case see chapter 2 of the P. Fulde book [16]).
Within the approximations made, the HamiltonianĤ SCF (1) has (for each value of the spin projection s) two eigenfunctions (bonding and antibonding orbitals, see Figure 1) with energies e ± = ε SCF ± t In the ground state, both electrons occupy the lowest energy bonding orbital |ψ − . The electrons should have opposite spins. So, two electron ground state function has the form |g SCF = a † −,↑ a † −,↓ |vac = 1 √ 2 (|s 1 + |s 2 ) (4) [φ 1 (r 1 )φ 2 (r 2 ) + φ 2 (r 1 )φ 1 (r 2 ) + φ 1 (r 1 )φ 1 (r 2 ) + φ 2 (r 1 )φ 2 (r 2 )](α 1 β 2 − β 1 α 2 ), The spinors α, β refer to spin up and spin down, respectively. We see that electrons move independently in the mean field approximation. The probability of finding them on the same site (it is given by the square of the coefficient of the so-called ionic configurations that are given by the terms of the form a † i,↑ a † i,↓ |vac ) equals the probability to finding the electrons on different sites (square of the coefficient of the terms a † 1,↑ a † 2,↓ |vac ). This means the absence of the correlations in the electron motion in the mean field approximation.
In order to consider the correlations, we have to go beyond the mean-field approximation and take into account so-called residual interaction, which is the difference between the bare electron-electron Coulomb repulsion ijkl,s,s w ijkl a † i,s a † j,s a l,s a k,s , w ijkl = d 3 rd 3 r ϕ * i (r)ϕ * j (r ) 1 |r − r | ϕ k (r)ϕ l (r ) (6) and its part that was accounted for in the mean field [16]. The residual interaction is much more localized compared to the bare Coulomb interaction (6). That is why Hubbard proposed to use the basis of localized atomic-like Wannier functions and to take into account only the largest terms having i = j = k = l [21]. Then, the HamiltonianĤ SCF (1) is supplemented by the term where n ≡ g SCF | ∑ s a † i,s a i,s |g SCF . The one-particle HamiltonianĤ c subtracts the part of interactionÛ, which was accounted for in the self-consistent field. As a result, the average g SCF |Ĥ H,res |g SCF = 0. In our problem, n = 1. The general form ofĤ H,res for the full Coulomb interaction (6) is given in Eq. (2.3.35) of Ref. [16]. In our notations [the Coulomb matrix elements of Ref. [16] are connected with ours (6) via relation V ikjl = w ijkl ], it iŝ where the bond order P ij ≡ g SCF | ∑ s a † i,s a j,s |g SCF is introduced. Now, the Hamiltonian for the hydrogen molecule acquires the form of the Hubbard HamiltonanĤ where we have dropped the constant term Un 2 /4. The diagonal matrix element is now ε = ε SCF − Un/2. We take it as the zero of energy ε ≡ 0. We shall solve the problem by the Heitler-London approach [22], which uses the many-body function basis. The Hamiltonian conserves the total spin S of the system. We can find the solution separately for singlet (S = 0) and triplet (S = 1) sectors of the Hilbert space. In addition, the Hamiltonian conserves the parity of the wave functions. Triplet sector has three degenerate wave functions The energy of the triplet is E t = 2ε = 0.
In the singlet sector, we have three basis functions, one odd with the eigenenergy E 3 = U, and two even: |s 1 and |s 2 (5). We obtain a 2 × 2 problem for even singletsĤ H |s 1 = −2t|s 2 ,Ĥ H |s 2 = U|s 2 − 2t|s 1 . Its eigenvalues are where ν = ±1. The last approximate equality (12) is valid in so-called strongly correlated limit U t. In this limit, The eigenvectors are Thus, for the ground state |g S = |g − , we have Comparing it with |g SCF (4), we see that the weight of ionic states |s 2 where electrons are found on the same site is strongly suppressed in the correlated wave function.
In the mean field approximation, the lowest excited state corresponds to one-electron excitation from bonding to antibonding orbital. It is separated from the ground state by the energy ∆E SCF = 2t. In contrast to the mean field theory, the many-body approach obtains the first excited state of the system as magnetic excitation that flips the spin of one electron and transfers the singlet ground state |g S to one of the triplet states |t, m , m = 0, ±1. The excitation energy is The set of four lowest states |g S , |t, m is separated from other states by an energy ∆E U ∼ U. It is easy to show that this set may be described by a low-energy effective Hamiltonian which is just the antiferromagnetic Heisenberg Hamiltonian.

Temperature Dependence of Optical Conductivity
It is clear that the triplet states will be populated at temperatures k B T ∼ J. The temperature dependence of some response functions of a strongly correlated system becomes observable at these temperatures. The values of effective exchange integral J/k B may vary from one to hundreds of Kelvins. Figure 2 illustrates the dependence of the charge response of the system on the magnetic initial state: the transition over the gap is allowed in the singlet state (left panel) and is prohibited in the triplet state (right panel) due to the Pauli exclusion principle. The population of the triplet states by temperature will evidently affect the response. Let us illustrate this statement by the calculation of the optical conductivity for an ensemble of molecules having a temperature T. Optical conductivity for an ensemble of finite systems is given by (polarization of light is along the x-axis parallel to the molecule) is the thermodynamic weight of the initial state |µ , ω νµ = E ν − E µ is the energy of the transition, V is the volume per one system,P x and x = −i P x ,Ĥ H are the polarization and the current operators. The absorptive real part of the optical conductivity is where σ µ (ω) is the contribution to the optical conductivity of transitions from the state |µ . In our restricted basis for a two site system, only optical transitions with charge transfer between sites is possible. The current operator iŝ d is the distance between sites, and e is the electron charge. Non-zero matrix elements exist only between singlets |s 1 and |s 3 . We thus have s 1 ||s 3 = −2ited, and g ν ||s 3 = −2itedα ν and we are able to calculate the optical conductivity analytically for any temperature and parameters of Hamiltonian (9). Figure 3 shows the optical conductivity for strongly correlated regime U/t = 10 (a finite imaginary part η = 0.1t was added to ω in order to visualize the δ-function). The typical value of t is ∼ 1 eV=11604 K·k B . Thus, for any realistic temperature, only the transition from the ground state (15) will be observed both in strongly (U t) and in weakly (U t) correlated limits. We may introduce the weight of the transition (the coefficient before δ ω − ω νµ ) The temperature dependence of W 3g is given in Figure 4. We see that the peculiarity of a strongly correlated system is a strong variation of W(T) with temperature. As the transition energy ω 3g ∼ U T , we may neglect w 3 (T). In addition, we may neglect all terms ∼ exp(−U/k B T) in the denominator of Equation (19), which is the partition function. Then, This means that the characteristic temperature of optical response variation is the magnetic energy J. This is the general property of strongly correlated systems [23][24][25].
In the next subsection, we will show that the same conclusion holds for the temperature dependence of resonant inelastic X-ray scattering (RIXS) spectra for another simple finite system.    We consider the three-site "Cu-O-Cu" cluster, which is described by a many-band Hubbard model (see Figure 5) [26] where the operator p † σ creates a particle with the spin projection index σ =↑, ↓ on the uncorrelated "O" site, is the Hubbard projection operator that creates a particle with the spin projection index σ on the site i = 1, 2, where strong correlations prohibit the double occupancy of the site (see Appendix C). Thus, we consider a limiting case of the Emery (multi-band Hubbard) model [27] (later in Section 4.2 we will discuss the Emery model for cuprates more in detail) when 1/U d = U p = 0, U d (U p ) being on-site repulsion on the "Cu"("O") site. The remaining energetic parameters are the charge-transfer energy ∆ (positive for hole representation), and the hopping t. The typical values for cuprates are ∆ ∼ 4, t ∼ 1 eV, and below we give the analytic formulas for exact values together with the expansion over t/∆ 1. We consider the insulating case when we have two holes in the cluster. The Classification of States. We may characterize the wave function by total spin value, its projection on the z-axis, and parity with respect to exchange of "Cu" sites 1 ↔ 2. In the two-particle sector, we have ( Figure 6) (i) S = 0 (1) even: where is the Zhang-Rice singlet state [28,29] formed by holes on neighboring Cu-O sites. (2) odd: (ii)S = 1, S z = 1 (1) even: (2) odd ( Figure 7): The wave functions with other values of S z may be obtained by the action of operator S − on the above states. Summary of the spectrum. The Hamiltonian (25) has 13 two-particle eigenstates (See Appendix A.1 for the details): 4 singlets (31), (A6), and 3 triplets (32) ,(A12), and (A13). In the low-energy states | f = 0 (A6) and |Gt (A12), the dominant contribution comes from the states |s d (27) and |t d (33), both having two particles occupying different "Cu" sites. In this low-energy part of the spectrum, the system has only spin degrees of freedom which are described by an effective Heisenberg Hamiltonian with the superexchange parameter J = E t,0 − E 0 ≈ 4t 4 /∆ 3 t, ∆. The constant shift 2e d = −2t 2 /∆ comes from the hybridization contribution to the "crystal field" on the "Cu" site.
The excited states corresponding to a transfer of charge between the "Cu" and the "O" sites lie higher in energy by the value about ∆.

O K RIXS Spectrum for Finite Temperature
General Expression The O K RIXS process has the following stages: (i) in a system being in an initial state |g , an X-ray quantum excites an electron from an oxygen core 1s-state to a 2p-state on the same site R; (ii) the valence electron system propagates in the presence of the immovable core-hole at site R; (iii) the 2p-electron recombines on the same site R with the core 1s-hole, another X-ray quantum is emitted, and the system is left in a finite state | f .
The RIXS spectrum intensity at finite temperature is given by (see e.g., [30][31][32][33]) where · · · T denote the statistical average over initial states |g for a temperature T, k ≈ 1/11604 eV/K is the Boltzmann constant, ( ) is the polarization vectors of incident(emitted) photons, and Ω and ω are the incident and emitted photon energies. For a given |g , the intensity of the O K RIXS signal is (see e.g., [31][32][33]) whereT µ,R = ∑ σ s † Rσ p Rµσ + h.c., s † Rσ is the creation operator of the O 1s-hole with spin projection σ at site R, p † Rµσ creates a 2p-hole at the same site, µ and ν are Cartesian indices of 2p-orbitals, and |m, R is the eigenstate of the Hamiltonian of the system in the presence of the core hole at site R (in the stage ii) In Equation (38), the first termĤ pd is the generalized many-band Hubbard Hamiltonian that describes the valence electron system of cuprates;Ĥ C,R describes the O 1s hole and its interaction with valence p-holes which is assumed to be reduced to local Coulomb repulsion. The sum in Equation (37) runs over sites R where the core hole is created at stage (i), rests at stage (ii), and annihilates at stage (iii), cf. Equation (3) of Ref. [33] where apparently a triple sum over R is present. In fact, the expression (38) implies that the core hole does not move and is annihilated at the same site where it was created. This reduces the triple sum to a single one.
With the assumptions inherent in Equation (38), the role of the core hole is reduced to the change of the on-site energy of valence p-states at the stage (ii) of the RIXS process. This may be shown in the following way: let us recall that the core hole is absent in the initial and final states, i.e., we can write |f = | f ⊗ |0 C , |0 C being the vacuum for core states. Then, we have with z ≡ E g + Ω − ıΓ. In the derivation of (39), we have used the relation The substitution of (39) into (37) gives Equation (41). We see that the O K RIXS spectral function (within the approximation made) is defined by the dynamics of valence electrons only. Thus, we obtain where p † Rµσ creates a 2p-hole with spin projection σ at oxygen site R, µ = x, y, z.Ĥ pd is the generalized many-band Hubbard Hamiltonian that describes the valence electron system, ε s is the energy of 1s hole level, and U C is the interaction strength between the 1s-and valence 2p-holes. The interaction is assumed to be reduced to local Coulomb repulsion.
The expressions (36)- (42) involves only valence states. The stages i)-iii) may be reformulated as: (i ) in a system being in the N-hole ground state |g , a hole in a 2p-state on the site R is annihilated; (ii ) the N − 1 hole system is perturbed by the increase of site energy on the site R by the value U C ; (iii ) the 2p−hole on the same site R is created, and the system is left in an excited state | f . Application to the Three-Site Model As we have already mentioned, the charge-transfer excitations have energies of the order of several eV. They will never be populated at temperatures reachable in an RIXS experiment (T < 0.1 eV). So, we should make the statistical average only over the low-energy states where w s = 1/Q(T) and w t = 3 exp(−J/kT)/Q(T) are statistical weights of the lowest singlet and triplet states, and Q(T) = 1 + 3 exp(−J/kT) is the partition function.
In the intermediate state, our system has only one hole. The eigenenergies of even states are where here sin γ and cos γ are given by expressions similar to Equations (A14) with the change R t → r, and ∆ → ∆ + U C . An odd state |σ, a = 1 2 |vac has the energy E a = 0. The calculated RIXS and XAS spectra are shown in Figures 8-11. We clearly see the resonant character of the spectra shown on Figure 8. The resonsnce occurs at different incident energies for singlet and triplet initial states. The occupation of different initial states depends on temperature. This leads to temperature dependence of the RIXS spectra. The changes of temperature on the scale of magnetic interaction value J leads to drastic changes of spectra on a much larger scale t J as shown in Figure 10.
The strong temperature dependence of the RIXS spectrum was first observed for Li 2 CuO 2 and CuGeO 3 edge-shared cuprate compounds in Ref. [34] The XAS spectrum also depends on temperature, but this dependence is weak. . T-dependence of XAS spectrum in a "Cu-O-Cu" molecule.

Resolvent Method (Löwdin Downfolding)
This method provides the simplest way to obtain a low-energy effective Hamiltonian from a full Hamiltonian of a system. It was proposed in a series of works of P.-O. Löwdin (e.g. Refs. [35,36]), where he called this method the "Partitioning technique". Here, we briefly review the technique.
We assume that a full Hilbert space of states of a system described by a Hamiltonian H may be divided into two parts A and B, A being an "interesting low-energy part" in some sense. Then, we may write the Hamiltonian matrix and an eigenvector symbolically Then, the secular equation Hc = Ec becomes We substitute the expression for b found from the second equation into the first one and obtain Equation (49) Note that H e f f (E) depends on the eigenenergy. So, Equation (49) is a non-linear equation. We should write E = E (0) + E (2) + · · · , and solve the equation iteratively.
Let us pay attention to an outstanding feature of the expansion (50). It never diverges if the states in the subspace A are separated by an energy gap from the states in B. This is the case when we derive an effective magnetic Hamiltonian for a strongly correlated insulating system. Then, the subspace B contains the states with charge excitations which are separated by a charge-transfer or a Hubbard gap for the charge-transfer or Mott-Hubbard insulators, respectively (for the classification of correlated systems see Section 4.4).

The Effective Hamiltonian after Fourth-Order Canonical Transform
Let us consider another way to obtain an effective Hamiltonian. We denotê where the eigenvalues E m and eigenvectors |m of an unperturbed HamiltonianĤ 0 are assumed to be known. Without loss of generality, we may assume that the perturbationV contains only non-diagonal terms. The operator has the property Then, up to the fourth order, the canonical transformation giveŝ The explicit calculation giveŝ The advantage of this approach is that it gives an energy-independent effective Hamiltonian. However, we can see a substantial difference of the Hamiltonian (54) from (50). The denominators of the third and fourth orders (58) and (59) contain the energy differences between intermediate states. These terms may diverge if the states are (quasi-)degenerate even if these states are well separated by an energy gap from the states in the subspace A. In fact, the generatorŴ has excluded the non-diagonal terms due to the property (53) only up to the second order. For the derivation of the fourth order effective Hamiltonian, we need to also exclude non-diagonal terms inĤ 2 . So, we perform a second transform with the generatorŴ Then, the fourth-order term becomeŝ After some algebra (see Appendix B for the details), we obtain Now, compare Equation (61) with G njklm given by (64) and the Löwdin result (50). We see that the first term of (64) produces the form similar to (50), if we assume that both |m , |n belong to the subspace A. The second term of (64) looks different. However, we can note the following: i) it vanishes when E m = E n , i.e., when the subspace A is degenerate; ii) we have written the transformation that tries to remove all non-diagonal terms in (51), whereas only the part of them, namely H AB , are removed in the Löwdin approach. If we divideV in (51) and make the transformation that removes onlyV AB , we will have the form similar to (50), but in practice, it is difficult to make the decomposition (65). As we mentioned above, one of the advantage of the canonical transform is that it gives an energy-independent effective Hamiltonian. The second important advantage is that it easily gives the transformed form of any operator to the same order: the transformed form of the wave function may be found as well.

Hubbard Model
One of the simplest models that shows the peculiarities of the physics of strongly correlated electron systems was introduced by J. Hubbard in Ref. [21] (implicitly this model was used by P.W. Anderson for consideration of the superexchange [37]). The Hubbard Hamiltonian readŝ where the summation goes over the sites R of a lattice (here we consider an infinite crystal lattice), a † R,σ creates an electron in a state with a wave function φ(r − R) localized at a site R with spin projection σ,n R,σ = a † R,σ a R,σ is the on-site operator of the number of electrons with spin projection σ, and vector g joins nearest neighbors.
In the mean field approximation, the Hamiltonian reduces to a single-band tightbinding HamiltonianĤ where the operator annihilating an electron in a band state is given by the Fourier transform where N is the number of sites in the lattice. Our system is translationally invariant; thus, the on-site average of the electron number ∑ σnR,σ = n does not depend on R. Here, we consider a non-spin-polarized case, when also n R,σ =n = n/2 does not depend on σ and R.
The simplest approximation that shows the strongly correlated behavior for is the so-called Hubbard-I approximation [21], which is a decoupling scheme for the two-time Green's function technique.
The aim is the calculation of the retarded Green's function where [A, B] η ≡ AB − ηBA, the time dependence of an operatorÂ(t) is given byÂ(t) = e itĤÂ e −itĤ , and the angular brackets denote the thermodynamic average creates an electron with spin projection σ at point r; We consider the system in a non-magnetic state and will drop the spin index of the Green's function. The equation of motion for the Green's function reads The Hubbard-I decoupling is introduced in the equation for higher-order function (76) ( In a non-magnetic state, the average does not depend on σ. Then, Equation (77) becomes and we find This equation may be solved using Fourier transform where From the last Equation (80), we obtain The Green's function (82) is always diagonal in k-space for translationally invariant Hamiltonians. The momentum-dependent spectral density is the main characteristic of the electronic structure of strongly correlated systems. It contains information both about the quasiparticle energy dispersion (given by poles of G k (z)) and about the incoherent bands (corresponding to the branch cuts of the Green's function). It is proportional to ARPES intensity in the so-called direct-transitions limit (see Section 5.4). Let us study the property of the Green's function (82). We rewrite it in the form It is clear that the self-energy Σ(ω) is due to the interaction. Only the static (ωindependent) part of it, Σ SCF = Un, is taken into account in the mean-field approximation. Then, the Green's function has a simple pole form the pole position being the one-particle energy of the mean-field Hamiltonian (68). The spectral density has a single delta-functional peak with a unit weight for each spin direction. The correlations are responsible for the dynamic part of the self-energy, which is local (k-independent) in the Hubbard-I approximation. J. Hubbard was the first to show that the correlations split the single mean-field band into two subbands, which are now called the low-and upper-Hubbard bands. Indeed, the Green's function (84) has two poles where ω 1,2 are the solutions of the equation [G k (ω)] −1 = 0, which gives This expression may be simplified in the strong correlation limit (70). Up to the terms of order tU −1 , we may write [16] Now, it is clear that the two bands are separated by a gap of the order of U.

Anderson and Emery Models
The single impurity Anderson model (SIAM) Hamiltonian readŝ where the first term ofĤ SI AM describes a band of uncorrelated electrons, the second term H f is a generalized single-site Hubbard term, which is the Hamiltonian of a transition metal impurity that has a localized degenerate level with the single-partical energy ε f and strong Coulomb repulsion U. Operator f † m creates an electron in a localized state, m is the set of quantum numbers that characterize the state (e.g. combination of orbital and spin projection numbers for d-electrons or total moment j projection for f -electrons),n f m ≡ f † m f m . The last termV represents the hybridization between the localized and delocalized states. The HamiltonianĤ SI AM (91) was first introduced in Ref. [42] for the explanation of the existence of localized magnetic moments in dilute magnetic alloys (see also Ref. [43]). The moments are localized on d-ion impurities in non-magnetic metals and are due to strong correlations within the d-shell of the ions. The model was intensively studied and allowed to explain magnetic and transport properties of the dilute magnetic alloys. In the limit of small s-d mixing, it was shown to be equivalent to the model introduced by Kondo [44]. The equivalence was proved by Schrieffer and Wolff [45] by means of canonical transform (see Section 3). In the beginning of the 1980s, exact solutions for the Anderson and Kondo models were found (see the review of Tsvelick and Wiegmann [46]). These solutions are valid only for the impurities in "good" metals, where the Fermi energy of uncorrelated electrons E F is the largest energy parameter of the model, and the electron spectrum may be linearized ε(k) ≈ v F (k − k F ). Then, the model is reduced to a one-dimensional problem and the Bethe Ansatz approach may be applied [46]. In Ref. [47], the Anderson model was applied for the description of d-ion impurities in semiconductors. In the following works, it was widely used for the description of diluted magnetic semiconductors.
The interest in SIAM increased considerably when the dynamical mean-field theory (DMFT) [48,49] was formulated. The authors of Ref. [50] showed that the Hubbard model in infinite dimensions may be exactly mapped onto a single-impurity Anderson model.
An important generalization of the Anderson model is the periodic Anderson model (PAM), where transition metal ions form a periodic sublattice. It is given by the Hamiltonian where i enumerates the sites of the transition metal sublattice, operatorsĤ f ,i ,V i have the form of Equation (92) with the substitution f † m → f † m,i . The discovery of high-T c cuprate superconductors (HTSC) [1] immediately made the low-dimensional strongly correlated systems the focus of scientific community attention. P.W. Anderson realized the importance of correlations for the physics of cuprates [51,52]. The model for the description of hole motion in the CuO 2 planes of HTSC was proposed by Emery [27] as a generalization of the Hubbard model where i labels a copper or an oxygen site, the operator a † i,σ creates a hole with spin index σ in the Cu(d x 2 −y 2 ) or O(p x,y ) orbitals, which are the ones most strongly hybridized. Only site diagonal terms ( p,d , U p,d ) and nearest neighbor hopping ( i,j = ±t) and interaction (U i,j = V) terms were taken into account. The parameter regime is relevant for HTSC (here B A should be understood as B/A 2.5). Later, it was found that the account of the next-nearest neighbor oxygen-oxygen hopping t pp t is necessary for a realistic description of HTSC [53,54]. The average number of holes in the unit cell of the Emery model n d + 2n p = 1 + x, |x| < 1. Positive values of x correspond to holedoped HTSC, whereas negative values describe electron-doped HTSC. The model with x = 0 describes parent compounds, which are antiferromagnetic insulators. In this case, the low-energy spectrum of the Emery model may be described by an effective isotropic Heisenberg Hamiltonian (see, e.g., Refs. [55,56]). If one neglects correlations on oxygen sites (U p ≈ 0), the Emery model becomes a special case of the periodic Anderson model (93).

Spin-Fermion and t − J Models
The downfolding of the Emery model in the regime given by Equation (95) allows obtaining low-energy models with a reduced number of degrees of freedom. The 'minimal' Emery model that exhibits the essential properties of layered cuprates (1/U d = U p = t pp = 0) reads (in hole notation) where the Fermi operatorp r,γ annihilates a hole at site r of the oxygen sublattice with spin projection index γ, and the Hubbard projection operatorZ 0γ R , Equation (26), (see also Appendix C) annihilates a hole with spin index γ on a singly occupied copper site. The double occupancy of copper sites is thus excluded from (96). The first term,Ĥ 0 , includes the on-site energies (∆ = p − d , d is taken as zero of energy),V is the p-d hybridization, α = x, −x, y, −y characterizes the direction of a nearest-neighbor vector a, and the phase factors inV are absorbed into the definition of the operatorsp r,σ ,Z 0γ R . In the limit t/∆ 1, further downfolding by means of a canonical transformation of operators of the form leads to the model (see also Ref. [54] for the notation). Here, p and Z mean transformed operators,Ĵ s is the AFM copper-copper superexchange interaction, and g points to neighboring copper sites. The parameters are τ = t 2 /∆, and the AFM exchange J ∝ t 4 /∆ 3 . The model (96) is called the spin-fermion model [29,57,58]. As we have mentioned in a previous subsection, in the absence of doping, the Emery model is equivalent to the nearest-neighbor AFM Heisenberg modelĴ s . An extra hole on the oxygen site forms a Zhang-Rice singlet and triplet states with a neighboring Cu site [28,29], the triplet state being ∼ 8τ higher in energy than the singlet. Exclusion of the triplet states leads to the t − J model [28]

Classification of Strongly Correlated Systems
In the seminal work of J. Zaanen, G. A. Sawatzky, and J. W. Allen [64], the transition metal compounds were classified according to the relations between the energetic parameters: The hopping integral t between a ligand and transition metal ion.
According to this work (see Figure 3 of Ref. [64]), the transition metal strongly correlated compounds may behave as:

Many-Band Generalization of the Models
For a realistic description of a specific compound, the generalizations of the above models are necessary. In most compounds, several correlated states per site and the dependence of the V k,m,σ matrix element in Equation (92) on the symmetry of m-th orbital should be taken into account (see e.g., [65]). An account of the geometry of bonds and the symmetry of anion's ligand orbitals is also necessary for the quantitative description of the transition metal compounds.
In Ref. [66], the five-band p − d model was introduced for the electronic structure of so-called edge-shared compounds (see also [33,34,67]). The model is used for the unified consideration of magnetic properties and the optic and RIXS spectra of the compounds (e.g., [24,25,34,[68][69][70][71]). The orbital basis (Figure 12) consists of a single 3d xy orbital on each Cu site and the 2p x and 2p y orbitals on each oxygen site. It is p † x,l,σ p x,l,σ p † y,l,σ p y,l,σ x,l,↑ p y,l,↑ p † x,l,↓ p y,l,↓ + p † y,l,↑ p x,l,↑ p † y,l,↓ p x,l,↓ ) where m, m are Cu site indices, l, l are oxygen site indices, α = x, y are orbital indices for the 2p x,y orbitals, . . . is a sum over nearest neighbors, and n d , n p α are the usual number operators for the Cu and O orbitals. Besides the one-particle on-site energies d,m , p,l,α , hoppings t ml d,p α , t ll p α ,p α , the Hamiltonian (100) accounts for the Hubbard terms on Cu and O sites with parameters U d and U p , the Hund coupling on O site (K p ), the direct Coulomb (U pd ) and exchange (K pd ) interactions between neighboring Cu and O sites, and Coulomb interaction (U dd ) between neighboring Cu sites. When electron-lattice coupling is strong due to the Jahn-Teller effect, the orbital degrees of freedom come into play. Then, the relevant low-energy model is the so-called Kugel-Khomskii Hamiltonian [72]. Its main feature is the appearance of pseudo-spin values that characterize the orbital degrees of freedom. There exists a close analogy between a hole moving in the antiferromagnetic background and an electron moving in the alternating orbital environment of double exchange ferromagnets. The term "orbital polaron" was originally introduced by R. Kilian and G. Khaliullin for a quasi-particle for which the charge degree of freedom is not only coupled to orbital fluctuations, but also to the lattice [73].
Later, the term was used in studies on effective, low-energy t − J such as Hamiltonians in the field of manganites for orbital quasi-particles [74][75][76][77][78].

Response Functions Calculations and Spectroscopies
A theoretical interpretation of experimental data demands calculations of response functions. Above, we have given some examples of the calculations for finite systems, Equations (18) and (37). Below, we give more examples of the response function calculations for extended strongly correlated systems.

Ab initio Ligand Field Theory to Determine Electronic Multiplet Properties
There is a fundamental problem in electronic structure theory of solids, namely the proper description of multiplet effects of local magnetic centers built up of d or f electrons, which are intrinsically many-body states, in translational invariant settings. The many-electron multiplet levels are characterized by strong Coulomb interactions, electron correlations, and spin orbit coupling. The multiplets have been well understood for many years in atomic physics. Such multiplets persist in solids, either as sharp levels in the gap of insulators or semiconductors or as resonances in metals and small gap semiconductors. The influence of the surrounding crystal on the d or f electron shell of an ion is described by a few new parameters that are traditionally called crystal field (CF) or ligand field (LF) [79,80] parameters. The knowledge of LF parameters allows describing the splittings and mixing of single-ion many-body states in a crystal and calculating the response functions, which are determined by local multiplet effects.
In the literature, one can find several approaches to calculate LF parameters. First, there are wave function quantum chemistry methods [81]. However, it is difficult for these methods to treat a periodic crystal, and they become numerically expansive for heavy ions and large systems. There exist numerous attempts in the scientific literature to calculate multiplets and LF parameters in an ab initio style and based on density functional theory (DFT) [82][83][84]. The authors of Refs. [12,85] follow a much simpler way by starting with the non-spin-polarized calculation using GGA functional [86]. Then, they obtain the LF parameters by a Wannier fit to the non-magnetic band structure. In Ref. [85], the LF parameters serve as input for an exact diagonalization computer program ELISA (electrons localized in single atom) to calculate the response functions sensible to local multiplet effects, i.e., electron paramagnetic resonance (EPR), optical spectroscopy, inelastic neutron scattering (INS), X-ray absorption, and X-ray magnetic circular dichroism (XAS and XMCD) as well as resonant inelastic X-ray scattering.
When RIXS experiment exploits a resonance on a core level of a transition metal ion, its intensity is given by the formula where ω = ω in − ω out is the energy transfer, the indices i and f denote initial and final states, respectively, and with the scattering amplitude where we sum over all intermediate states m, and E in , E out are polarization vectors of incoming and outgoing X-ray radiation.
The optical absorption spectra are calculated by using the approach of Sugano and Tanabe [87], where the d-d transitions between two states a and b become possible by combining a parity changing perturbation V odd with the dipole operator P = q · r to give the transition probability by where ∆E is the energy difference between the given configuration with an incomplete d shell and the first excited configuration with odd parity.
To calculate the optical and X-ray spectra, the dipole transition probabilities are calculated in the ELISA code as it was published for XAS and XMCD [88].

Spin-Hamiltonians and Magnetic Response.
As we have already mentioned, in strongly correlated systems, charge and spin degrees of freedom are well separated in energy. For insulators, the spin and phonon excitations have the lowest energy and thus are responsible for thermodynamic properties. Magnetic properties are described by effective spin-Hamiltonians. Theoretical determination of parameters of the spin-Hamiltonians is the application of Löwdin downfolding to the generalized Hubbard model. The downfolding for magnetic impurities in non-magnetic materials is performed in three steps [89]: First, the virtual hoppings of electrons between the impurity ion and surrounding ligands are eliminated, and one obtains an effective ligand-field Hamiltonian. In the second step, one takes into account the fact that the largest part of the LF Hamiltonian (usually, the cubic splitting) is smaller than the remaining Coulomb interactions. Thus, an effective Hamiltonian for the lowest multiplet is obtained. Finally, the couplings of the ground state manifold with higher levels due to the smallest low-symmetry LF terms and by the spin-orbit interaction are eliminated. Thus, it is possible to obtain an analytical closed expression which connects the parameters of the microscopic Hamiltonian of the generalized Hubbard model with the parameters of the effective spin Hamiltonian [85,89,90].
Analogous downfolding is possible for the extended systems where magnetic ions form a regular sublattice. We have already mentioned the equivalence of the Heisenberg model and the low-energy behavior of the half-filled Hubbard model [38][39][40][41] and of the undoped Emery model (e.g., [55,91,92] and references therein). The account of spin-orbit interaction allows obtaining anisotropic terms [93][94][95]. Now, it becomes a standard to obtain the parameters of a spin-Hamiltonian from the spin-density-functional calculations. Comparison of total energies of different magnetic configuration allows finding the exchange values. Here, we cite only a few examples of such determination of parameters [96][97][98] just for the illustration of the method. In most cases, the density-functional theory (DFT) calculations do not provide exact values of the exchange interactions, but they allow establishing the hierarchy of the interactions. This is very important for the compounds with competing frustrated interactions. The refinement of the parameters is then possible by comparison with an experiment. Fitting of the inelastic neutron scattering (INS) spectra within the linear spin-wave theory may be ambiguous. It was the case of the quasi-one-dimensional LiCu 2 O 2 edge-shared cuprate compound. First, the INS spectrum of LiCu 2 O 2 was interpreted in terms of an antiferromagnetic J 1 − J 2 model [99], but the DFT calculations and other experimental evidence has shown that the main interactions along the chain are ferromagnetic nearest neighbor J 1 and antiferromagnetic next nearest neighbor J 2 [100,101]. Both sets of parameters explain the measured INS spectrum.
Once the spin-Hamiltonian is established, all the magnetic properties of the compound may be described. For example, to model the temperature dependence of the magnetic susceptibility χ(T), the open-source program code HTE10 may be used [102][103][104]. It provides the tenth-order high-temperature expansion (HTE) of a general Heisenberg model with up to four different exchange parameters J 1 , J 2 , J 3 , and J 4 . The tenth-order HTE is indispensable for systems where the scale of the exchange interactions JS(S + 1) is comparable to or exceeds the scale of thermal energy in the entire range of measurement temperatures [71,105,106]. Since the maximal measurement temperature is several hundred Kelvins under normal laboratory conditions, it has to be used in the cases where the exchange energy values (i.e., J/k B ∼ 100 K) are of comparable order and where the susceptibility does not obey the Curie-Weiss law that follows from the second-order HTE. The Curie-Weiss law fitting of experimental data, which is a common practice, may lead to false estimates of magnetic interactions in such a material [105]. The program HTE10 calculates the exact coefficients c n of the normalized susceptibility per one spin calculated in the tenth-order of HTE, and Padé approximants (ratios of two polynomials of m-th and n-th order, P m (T) and P n (T), respectively) χ HTE (T) ≈ [m, n] = P m (T)/P n (T), m + n ≤ 10. The Padé approximants allow extending the region of validity of the HTE [103].

Electron Energy Loss Spectroscopy
What is actually measured in transmission electron energy loss spectroscopy (EELS) experiments is the partial cross section [107,108] that may be decomposed into an amplitude factor and a dynamic structure factor The dynamic structure factor characterizes the linear response of the whole electronic system on longitudinal electric fields with the momentum q and frequency ω (the ionic contribution may be neglected for the considered frequency range of the order of several eV). Pronounced peaks in S(q, ω) are related with charge excitations: plasmons and excitons (sometimes all of them are called "excitons" [109]). The dynamic structure factor is related to the density-density correlation function where β is the inverse temperature, is the electronic density operator in the localized basis, the summation runs over all lattice sites r and orbital sorts s, and . . . means the thermodynamic average. For βω 1, we have is the retarded Green's function that defines the inverse dielectric function with v c being the volume of the unit cell, and e is the electronic charge. The function N(q, ω) describes the response to the unscreened external potential. The response to the total, screened potential is given by the function [110] N s (q, ω) = ε(q, ω)N(q, ω), In the diagrammatic language, the linear response to the total field may be expressed by the polarization operator where only irreducible graphs (which do not contain the contribution of the macroscopic electric field) should be taken into account [111,112]. Combining Equations (108) and (109), we express the dielectric permittivity via N s (q, ω) Substituting ε(q, ω) from Equation (110) to Equation (109), we obtain the relation which is exact for q → 0, as it was shown in Ref. [111]. The density response function N H (q, ω) calculated within a generalized Hubbard model is an approximation to N s (q, ω) [113]. In other words, it describes the motion of transverse (or "mechanical" by terminology introduced in Section 2.2.2 of the Agranovich and Ginzburg book [109]) excitons. The transverse ("mechanical") excitons are excitations that correspond to poles of dielectric permittivity, Equation (110), zeros of the inverse dielectric function, Equation (108), are determined by short-range interactions.
Using the spectral representation, we may write Here, we bear in mind that the Hubbard model contributes to transitions in the low frequency region ω < ω 0 with ω 0 of the order of the bandwidth, and the electrons of the rest of the solid are excited only at higher energies. In zero approximation, we may assume that in the frequency region ω > ω 0 , the electronic polarization of the rest of the solid follows the external field immediately N ∞ (q, z) ≈ N ∞ (q, 0). In other words, the Hubbard model is embedded into the medium with dielectric permeability ε ∞ (q) = 1 − 4πe 2 v c q 2 N ∞ (q, 0). In fact, ε ∞ may have its own dispersion and may be quite anisotropic for a layered or quasione-dimensional compound. In principle, it should be taken from, e.g., LDA calculations (we have assumed that the rest of the solid is uncorrelated) or from the experiment. It is obvious that the peak positions of the loss function L(q, ω) ≡ −Im ε −1 (q, ω) (113) and their intensity strongly depend on the value of ε ∞ (q). Usually, one neglects the qdependence and the anisotropy of ε ∞ , but it is a crude approximation, as well as another one which assumes ε(q, 0) = const. For a quantitative description of EELS experiments, the detailed knowledge of ε ∞ (q) is necessary. Then, the total dielectric function and its inverse are In Ref. [56], the problem of dielectric response in the strong coupling regime of a charge-transfer insulator was considered. An approach that starts from the correlated paramagnetic ground state with strong antiferromagnetic fluctuations was proposed. A set of coupled equations of motion for the two-particle Green's function was obtained and approximately solved by means of the projection technique. The solution is expressed by a two-particle basis that includes the excitonic states with electron and hole separated at various distances. The method was applied to the multiband Hubbard (Emery) model that describes layered cuprates. It was shown that strongly dispersive branches exist in the excitonic spectrum of the 'minimal' Emery model (1/U d = U p = t pp = 0). For this purpose, the downfolding to the spin-fermion model, Equation (98), was performed using the canonical transform, Equation (97), for the Hamiltonian and for the density operator, Equation (106). Then, the motion of electrons and holes in the effective Hamiltonian was considered. The exciton spectrum dependencies were analyzed on finite oxygen hopping t pp and on the value of on-site repulsion on oxygen, U p .

Angle-Resolved Photoemission Spectroscopy
It is commonly believed that the intensity of the angle-resolved photoemission spectra (ARPES) is proportional to the one-particle spectral function [114], which is an imaginary part of the retarded Green's function [cf. Equation (84)] divided by −π where a k,s,α annihilates an electron in a bulk Bloch state, Â ,B ≡ÂB +BÂ, the timedependent operatorÂ(t) isÂ(t) = exp(itĤ)Â exp(−itĤ), and the angular brackets denote the ground state or thermodynamic average, Equation (72). In Ref. [115], it was shown that this is the case only for one-and two-dimensional systems with a negligible dispersion normal to the surface. However, the actual crystals are three-dimensional, and the ARPES intensity (i.e., the steady radial photocurrent of electrons emerging from the solid along the observation direction defined by the unit vector q with energies between E and E + dE) is proportional to the spectral function of a more complicated Green's function where the operatorĈ creates an electron in a state with the wave function Here, ϕ > is the low-energy electron diffraction (LEED) wave function,Ô(x) is the operator of electron-light coupling. The function χ(r) decays into the solid owing to the spatial decay of the LEED function, and, at the same time, it rapidly vanishes in the vacuum owing to the confinement of the initial states.
It was shown how the spectra depend on physical properties of the initial and final states of the photoemission process. Both kinds of states are solutions of the Schrödinger equation with the same Hamiltonian. For the initial states, it is necessary to find the Green's function of the semi-infinite crystal. In the description of final states, the inelastic scattering due to electron-electron interaction in the propagation of the outgoing electron may be taken into account phenomenologically by introducing an absorbing optical potential into the effective Schrödinger equation for the function ϕ > (x,q, E) [116][117][118].

Application of the Methods to Specific Material Families
In this section, we briefly outline applications of the above methods for studies of strongly correlated materials.

High-T c Cuprate Superconductors
In cuprate materials, the most important features of the electronic structure are a large hybridization of O2p and Cu3d states in the pdσ-band and a strong local Coulomb repulsion on Cu3d states in the CuO 2 plane. As we have already mentioned in Section 4.2, the generalized Hubbard model was proposed by Emery [27], Equation (94). The hopping parameters t, t pp were taken from DFT band-structure calculations that were regarded as equivalent to applying mean-field theory to the model. Hubbard repulsion parameters were derived from photoemission and optical experiments in Refs. [119][120][121]. Later, the calculations with the constrained-density-functional approach [122,123] confirmed the obtained values. The importance of the account of Cu-O direct exchange K pd was pointed out in Ref. [124].
The t − J model (see Section 4.3) obtained by downfolding of the Emery model allowed describing most of the low-energy physics of the high-T c superconductors. A comprehensive review is given in the book [2].

Edge-Shared Cuprates
The electronic structure close to Fermi energy of these compounds is defined by Cu 3d and oxygen 2p states in CuO 2 chains ( Figure 12). As we have noted in Section 4.5, the electronic structure of the chain is described by the five-band p − d model, Equation (100). The magnetic part of its spectrum is well described by the one-dimensional spin-1/2 J 1 -J 2 Heisenberg model.
The ESC compounds represent a particular class of quantum magnets in which the local geometry gives rise to competing nearest ferromagnetic (FM) or antiferromagnetic (AFM) exchange coupling J 1 and frustrating antiferromagnetic next-nearest neighbor J 2 superexchange couplings.
In this rich family, one of the most studied compounds is Li 2 CuO 2 . The INS studies reported in Ref. [68] supplemented by the DFT calculations and exact diagonalization studies of Cu n O 2n+2 -clusters (n = 5, 6) allowed establishing a set of consistent parameters of the five-band p − d model that describe optical, EELS, O 1s XAS-spectral data [24,129], RIXS spectra [34,70], value of the magnetic saturation field [130], and temperature dependence of magnetic susceptibility [105].

Concluding Remarks: Building of a Microscopic Model for a Description of a Specific Material
We have outlined a realistic strategy in the description of transition metal compounds. The common steps of the model building are outlined in Figure 13. First, state-ofthe-art density functional theory calculations should be performed for the given composition and the crystal structure. Numerous computer codes are available for this purpose. We mention here the most popular codes: the Vienna ab initio simulation package (VASP) [131,132], WIEN2k [133], and FPLO [134]. A comparison of accuracy and an extensive list of the DFT codes may be found in Refs. [135,136]. If it is difficult to make a DFT calculation (e.g., because of a large unit cell or in the case of modeling of an impurity), Harrison's model (see Appendix D) [137] may be used.
Photoemission spectra provide an additional input that allow estimating the value of the Hubbard repulsion U for the transition metal ions.
On the base of the DFT calculations and photoemission data, the hierarchy of interactions in the compound may be established. Then, the generalized many-band Hubbard model is formulated. This model allows describing the electronic structure of the system on the energy scale of several eV. The general features of charge response (approximate energies of charge-transfer transitions) may be estimated on this step.
For the calculation of magnetic response and of thermodynamic properties, one needs to pass to the low-energy scale (about onetenth or one one-hundredth of eV) using Löwdin downfolding, canonical transform or a mapping of models using energy spectra of small clusters. In previous Sections 2.1 and 2.2, we have given examples of mapping of lowenergy spectra of Hubbard models onto the Heisenberg model. These methods allow connecting the parameters of the low-energy model with the parameters of the initial Hubbard model.
We have also shown that the detailed calculations of charge response should take into account the magnetic state of the system.
where A 2 f = ∑ 2 j=0 P j (E f ) 2 , the basis states |j = 0, 1, 2 are given by Equations (27)-(29) correspondingly, and P j (E) are polinoms, given by the recursion EP n (E) = a n P n (E) + b n P n−1 (E) + b n+1 P n+1 (E), with the initial conditions P 0 = 1, P −1 = 0; a n = H s,n,n are diagonal, and b n = H s,n,n+1 are off-diagonal elements of the matrix (A1), thus Up to the second order, we have for the eigenvector matrix here y ≡ t/∆. The odd singlet (31) has the energy E ZRSa = ∆ its eigenvalues are where R t ≡ 1 + 8(t/∆) 2 . The eigenvectors are |Gt = cos φ t |t d + sin φ t |ZRT , (A12) where

Appendix B. Simplification of the Result of Fourth-Order Canonical Transform
Let us simplify the expression (63) We recall that E mj ≡ E m − E j , then Now, the last term we rewrite as following Combining terms with similar denominators, we have