1 Introduction

Strong electron–electron correlations can give rise to surprising co-operative phenomena, typically very hard to explain and even harder to predict. Paradigmatic examples are high-temperature superconductivity in the cuprates, other types of unconventional superconductivity, Mott insulating or heavy fermion behavior, spin and orbital ordering, spin and orbital liquid phenomena. Understanding many-body effects in materials is, by all means, a grand challenge since the early days of quantum mechanics. And yet, what are, exactly, strong electronic correlations? In principle, one could say that condensed matter physics is all about electron–electron interactions. Indeed, the general electronic Hamiltonian that describes all possible systems, in the non-relativistic limit and the Born–Oppenheimer approximation, is

$$\begin{aligned} {\hat{H}}_e= \underbrace{\sum _i^{N_e} \bigg ({{\hat{T}}_i +\sum _\alpha ^{N_n} {\hat{V}}_{i\alpha }} \bigg )}_{ {\hat{T}}_e+{\hat{V}}_{en}} +\underbrace{\frac{1}{2}\sum _{\alpha \ne \alpha '}^{N_n} {\hat{V}}_{\alpha \alpha '}}_{{\hat{V}}_{nn}} + \underbrace{\frac{1}{2}\sum _{i\ne i'}^{N_e} {\hat{V}}_{ii'} } _{{\hat{V}}_{ee}}, \end{aligned}$$
(1)

where \(N_e\) is the number of electrons, \(N_n\) the number of nuclei, \({\hat{T}}_e\) the kinetic energy of the electrons, \({\hat{V}}_{en}\) the electron–nuclei attraction, \({\hat{V}}_{nn}\) the nuclei–nuclei repulsion, and \({\hat{V}}_{ee}\) the electron–electron repulsion. If the latter could be neglected (\({\hat{V}}_{ee}=0\)), the Hamiltonian would be separable, and the electrons independent. In this case the solution of the many-body problem would be straightforward. For a single electron, eigenvalues and eigenvectors could be found solving the equation

$$\begin{aligned} {\hat{H}}^0_i\phi _n(\mathbf{r}_i)=\varepsilon _n \phi _n(\mathbf{r}_i), \end{aligned}$$
(2)

where

$$\begin{aligned} {\hat{H}}^0_i= {\hat{T}}_i + \sum _\alpha ^{N_n} {\hat{V}}_{i\alpha }+ \frac{1}{2}\sum _{\alpha \ne \alpha '}^{N_n} {\hat{V}}_{\alpha \alpha '}. \end{aligned}$$
(3)

The many-electron states would then be single Slater determinants, trivially constructed from the complete basis \(\{\phi _n(\mathbf{r})\}\) just obtained. For example, in the case of a half-filled band described by the dispersion relation \(\varepsilon ({\mathbf{k}})\) and Bloch states \(\phi _{\mathbf{k}\uparrow } (\mathbf{r})\), such a Slater determinant has the form

$$\begin{aligned} \varPsi (\mathbf{r}_1,\mathbf{r}_2,\dots ,\mathbf{r}_{N_e})=\frac{1}{\sqrt{N_e!}} \; \begin{array}{|cccc|} \phi _{\mathbf{k}_1\uparrow } (\mathbf{r}_1) &{}\phi _{\mathbf{k}_1\uparrow } (\mathbf{r}_2) &{}\dots &{} \phi _{\mathbf{k}_1\uparrow } (\mathbf{r}_{N_e}) \\ \phi _{\mathbf{k}_1\downarrow }(\mathbf{r}_1) &{}\phi _{\mathbf{k}_1\downarrow }(\mathbf{r}_2) &{}\dots &{} \phi _{\mathbf{k}_1\downarrow }(\mathbf{r}_{N_e}) \\ \vdots &{} \vdots &{}\vdots &{} \vdots \\ \phi _{\mathbf{k}_{\frac{N_e}{2}}\uparrow } (\mathbf{r}_1) &{}\phi _{\mathbf{k}_{\frac{N_e}{2}}\uparrow } (\mathbf{r}_2) &{}\dots &{} \phi _{\mathbf{k}_{\frac{N_e}{2}}\uparrow } (\mathbf{r}_{N_e}) \\ \phi _{\mathbf{k}_{\frac{N_e}{2}}\downarrow } (\mathbf{r}_1) &{}\phi _{\mathbf{k}_{\frac{N_e}{2}}\downarrow } (\mathbf{r}_2) &{}\dots &{} \phi _{\mathbf{k}_{\frac{N_e}{2}}\downarrow } (\mathbf{r}_{N_e}) \\ \end{array}. \end{aligned}$$
(4)

In real materials, however, the electronic Coulomb repulsion is both large and long ranged, and there is no obvious reason to neglect it. This can be already seen from the average bare Coulomb energy (Hartree term)

$$\begin{aligned} {E}_{H}=\frac{1}{2}\int d \mathbf{r} \int d \mathbf{r}' \frac{n(\mathbf{r})n(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|}, \end{aligned}$$
(5)

where \(n(\mathbf{r})\) is the electron density.Footnote 1 Even if electronic charges are localized at different atomic sites (labeled with \({{\mathbf {R}}}_{\alpha }\)), i.e., if

$$\begin{aligned} n(\mathbf{r})\approx \sum _{\alpha }\delta (\mathbf{r}-{{\mathbf {R}}}_{\alpha }), \end{aligned}$$
(6)

the integrand still decays as the inverse of the nuclei–nuclei distance. Since the electron–electron repulsion cannot be neglected, a many-body eigenstate can be, in principle, a combination of infinite Slater determinants. It is thus remarkable that, de facto, for understanding several phenomena, Hamiltonian (1) can be replaced to a good approximation with an effective non-interacting model for quasi electrons

$$\begin{aligned} {\hat{H}}_e\rightarrow \tilde{{H}}_e^0= \sum _i^{N_e} \tilde{{H}}^0_i, \end{aligned}$$
(7)

whose eigenstates are single Slater determinants. The independent-electron approximation is, e.g., sufficient to explain the key differences between transition metals (Ni, Cu, Ag, ...), alkali metals (Li, Na, K,...), semiconductors (Si, Ge, GaAs, ...) and band insulators (diamond,...). It is so successful that modern solid-state physics text books still devote a substantial volume to results obtained starting from it. Remarkably, an independent–electron model is also used in a different context. This is the Kohn–Sham (KS) construction in density-functional theory (DFT), the most popular tool for electronic-structure calculations. In DFT, the KS Hamiltonian is, however, merely an auxiliary model, instrumental in obtaining the exact ground-state energy and density—it has, in principle, no physical meaning. Nevertheless, the practice of using it as an effective model has proven successful for entire categories of systems and problems. The KS band-structure is thus an established tool for studying, understanding and predicting the electronic properties of materials. In the light of such an unexpected success, stepping outside the rigorous principles of DFT, one could view also the KS Hamiltonian as a (very effective) mean-field mapping of type (7).

Mapping a complicated many-body Hamiltonian into a simpler effective auxiliary problem has by and large proven a powerful technique in physics. Mean-field approaches capture, e.g., essential features of magnetic broken symmetry phases, although the problem of quantum magnetism is vastly more complex. It is thus tempting to conclude that, if the mapping in Eq. (7) is an excellent approximation in many cases, it should work for all systems, once we found the correct non-interacting effective particles. Strong correlation phenomena are those for which such an Ansatz epically fails, and a simple effective non-interacting model is not expected to be found for reasons of principle.

The article is organized as follows. First we will briefly discuss the nature of the strong correlation problem. We will then recall the basics of density-functional theory, the standard method used to perform electronic structure calculations in materials. In this context we will discuss the Kohn–Sham construction and, leaving aside the principles of DFT, the strength and the advantages of the Kohn–Sham picture. We will stress the fundamental reasons for its failure for strongly-correlated compounds, illustrating the two main strategies used for overcoming the problem. The first consists in improving the approximation of the exchange–correlation functional—even if at the price of introducing ad hoc corrections or free parameters, and thus operating outside the strict principles of density-functional theory. This scheme is used, e.g., in the DFT+U approach. The limitations of such a tactic will be emphasized. The second strategy consist in using the KS orbital as a basis for building materials-specific many-body models. This is the most effective strategy so far, and the one used in the DFT+DMFT approach. Successes and open problems will be reviewed. The differences between static (DFT+U) and dynamical (DFT+DMFT) mean-field methods will be elucidated. Towards the end of the article we will take a glimpse to extensions of the DFT+DMFT approach and future developments.

2 The quantum many-body problem

In the reductionist viewpoint, knowing the fundamental interactions, here given in Eq. (1), is alone sufficient to reconstruct the behavior of any system. There is no doubt that reductionism led to crucial advances in science. It also harbors a key misunderstanding, however, as Anderson pointed out in the well-known article More is different. [1]. Taking it to the extreme, reductionism implies that solving the many-body problem described by Hamiltonian (1) is a mere exercise, perhaps a complicated one, an enterprise however from which nothing fundamentally new can be learned. If this was true, phenomena like superconductivity should have been predicted, rather that found experimentally, to be understood only several decades later. Even more, in the era of supercomputers, all phenomena in chemistry and condensed matter physics would be unraveled. Saying it with Anderson [1]

The constructionist hypothesis breaks down when confronted with the twin difficulties of scale and complexity. The behavior of large and complex aggregates of elementary particles, it turns out, is not to be understood in terms of a simple extrapolation of the properties of a few particles. Instead, at each level of complexity entirely new properties appear, and the understanding of the new behaviors requires research which I think is as fundamental in its nature as any other.

Fig. 1
figure 1

The dimension D of the Hilbert space for a \(N_n\)-sites chain with 4 states per site

In fact, only problems involving a very small number of particles can be solved exactly in an exhaustive manner. Exact diagonalization becomes quickly impracticable increasing the system size, since the Hilbert space grows combinatorially. This can be already understood by assuming that the Hilbert space of a single atom is merely made by 4 states, the vacuum, the two single-electron states, and one two-electron state, given below with their energy, \(E(N_e)\)

$$\begin{aligned} \begin{array}{lrccl} |0\rangle &{}&{}&{}\quad \quad &{} E(0)=0\\ |1,\sigma \rangle &{}=&{}c_{d\sigma }^\dagger |0\rangle &{}\quad &{} E(1)=\varepsilon _d\\ |2\rangle &{}=&{}c_{d\uparrow }^\dagger c_{d\downarrow }^\dagger |0\rangle &{}\quad &{} E(2)=2\varepsilon _d+U. \end{array} \end{aligned}$$

For a chain made by \(N_n\) of such atoms, the dimension is \(D=4^{N_n}\); the faster than exponential increase of D as a function of \(N_n\) is shown in Fig. 1. In the last decades, exploiting the advances in algorithms and supercomputers, quantum chemistry made humungous progress in finding efficient ways of studying the ground-state properties of systems for which a very large numbers of Slater determinants is needed, thanks, e.g., to coupled-cluster theory or full configuration interaction quantum Monte Carlo [2,3,4,5]. The density-matrix renormalization group [6,7,8] is also becoming a powerful tool in this respect [9,10,11,12]. Despite the impressive progress, the thermodynamic convergence remains slow, however. Staying with our toy model, in the thermodynamic limit a Hamiltonian matrix of dimension \(D=4^{N_n}\) cannot be stored on any classical supercomputer, not to speak about diagonalizing it and using the eigenvectors afterwards to calculate interesting properties. The problem becomes obviously even harder in real materials, where the Hilbert space of the atomic building blocks is not limited to four states. Truncations of the basis can lead to large finite size errors and/or, in the worse case, to entirely miss phenomena whose description requires high accuracy in specific energy windows. While quantum computers might change the game in the future, the second fundamental problem with the growing size of the Hilbert space is that it becomes rapidly impossible to make sense of the results. Indeed, the final aim of all simulations is not to merely reproduce experiments, but to give us answer to the truly interesting questions. Why does a systems behave the way it does? With an exponentially large number of degrees of freedom, even if we had the means of obtaining the exact solution, we might be left with an ocean of information impossible to disentangle. Having the solution would then not necessary lead to real scientific progress. An example of this paradox comes from the classical gravitational N-body problem. In this case, an exact solution in terms of power series has been obtained. First it was done for the three-body problem [13] and later on it was extended to the N-body case [14]. Nevertheless, limited progress in understanding the N-body problem comes from power-series solutions themselves. This is due to the fact that they converge slowly and simulations can become quickly impossibly long [15]. In a similar way, most progress in unravelling the quantum many-body problem in materials comes from methods which construct minimal effective models capturing the essence of a given phenomenon with sufficient degree of sophistication, and from solutions of such models, perhaps approximate, but allowing us to answer to essential “why” questions. This category includes a whole range of techniques, from those for solving paradigmatic models (describing only the generic features of a given phenomenon) all the way, I will argue, to the so-called first-principles methods, discussed in the next section.

3 The standard model: density-functional theory

Density-functional theory [16, 17] has revolutionized the way we deal with the many-electron problem. For this reason, DFT can be to some extent viewed as the standard model of condensed-matter physics. Given such a status, there are many reviews and/or books devoted to DFT, covering the principles, the applications, the history, or all of this together. A selection—definitely non-exhausting—can be found in Refs. [18,19,20,21,22,23,24]. Here I will merely recall the aspects most relevant for strongly-correlated systems. DFT shifts the focus from the electronic wave-function, which depends on \(N_e\) three dimensional coordinates, \(\varPsi (\mathbf{r}_1,\dots \mathbf{r}_{N_e})\), to the electron density \(n(\mathbf{r})\), a function of a single coordinate. This changes completely the perspective. For this reason, in 1998, in his Nobel lecture [18], Walter Kohn described the most important contributions of DFT starting as follows

[..] The first is in the area of fundamental understanding. Theoretical chemists and physicists, following the path of the Schrödinger equation, have become accustomed to think in a truncated Hilbert space of single particle orbitals. The spectacular advances achieved in this way attest to the fruitfulness of this perspective. However, when very high accuracy is required, so many Slater determinants are required (in some calculations up to \(\sim 10^9\)) that comprehension becomes difficult. DFT provides a complementary perspective. It focuses on quantities in the real, three-dimensional coordinate space, principally on the electron density \(n(\mathbf{r})\) of the ground state. [..] These quantities are physical, independent of representation and easily visualisable even for large systems.

Decades later this conclusion remains unchallenged. Still, the rise of DFT as the state-of-the art approach for describing actual materials is to a large extent due to a second step, the Kohn–Sham construction. In practical implementations, the electronic ground-state density \(n(\mathbf{r})\) is typically calculated by mapping the many-body Hamiltonian onto an auxiliary non-interacting KS Hamiltonian, \({\hat{H}}^{\mathrm{KS}}_e\), with the same electronic ground-state density,

$$\begin{aligned} n(\mathbf{r})=n_0(\mathbf{r})=\sum _{l}^{\mathrm{occ}} |\phi ^{\mathrm{KS}}_{l}(\mathbf{r})|^2. \end{aligned}$$

Here \(\{\phi _l(\mathbf{r})\}\) are Kohn–Sham orbitals; the associated Kohn–Sham eigenvalues are Lagrange multipliers, obtained solving the equation

$$\begin{aligned} {\hat{h}}^{\mathrm{0}}_{e} (\mathbf{r})\; \phi ^{\mathrm{KS}}_{l}(\mathbf{r})= \big (-\frac{1}{2} \nabla ^2 +v_{\mathrm{KS}}(\mathbf{r})\big )\phi ^\mathrm{KS}_{l}(\mathbf{r})=\varepsilon _{l}\phi ^{\mathrm{KS}}_{l}(\mathbf{r}) \end{aligned}$$
(8)

with effective potential given by

$$\begin{aligned} v_{\mathrm{KS}}(\mathbf{r})=\underbrace{- \sum _{\alpha } \frac{Z_\alpha }{|\mathbf{r}-\mathbf{R}_\alpha |}}_{v_{en}{(\mathbf{r})}}+ \underbrace{\int d \mathbf{r}^\prime \frac{n(\mathbf{r}^\prime )}{|\mathbf{r}-\mathbf{r}^\prime |}}_{v_{\mathrm{H}}(\mathbf{r})} +\underbrace{\frac{\delta E_{xc}[n]}{\delta n(\mathbf{r})}}_{v_{\mathrm{xc}}(\mathbf{r})}. \end{aligned}$$
(9)

In this expression the first term (\(v_{en}{(\mathbf{r})}\)) is the external potential, the second and third are the Hartree (\(v_\mathrm{H}{(\mathbf{r})}\)) and the exchange-correlation (\(v_{\mathrm{xc}}(\mathbf{r})\)) potential; \(E_{xc}[n]\) is the exchange-correlation functional. The exact ground-state energy of the system can be expressed as

$$\begin{aligned} E_G[n]=\sum _l^{\mathrm{occ}} \varepsilon _l -E_H[n]+E_{\mathrm{xc}}[n] -\int d \mathbf{r}\; \frac{\delta E_{\mathrm{xc}}[n]}{\delta n(\mathbf{r})} n(\mathbf{r}), \end{aligned}$$
(10)

where

$$\begin{aligned} E_H[n]&=\frac{1}{2}\sum _{ll'}\int d\mathbf{r}_1 \int d \mathbf{r}_2\; \overline{\phi ^{\mathrm{KS}}_{l}(\mathbf{r}_1)} \overline{\phi ^{\mathrm{KS}}_{l'}(\mathbf{r}_2)} \frac{1}{|\mathbf{r}_1-\mathbf{r}_2|}\phi ^{\mathrm{KS}}_{l}(\mathbf{r}_1) \phi ^{\mathrm{KS}}_{l'}(\mathbf{r}_2) \end{aligned}$$
(11)
$$\begin{aligned}&=\frac{1}{2}\int d \mathbf{r} \int d \mathbf{r}'\; n(\mathbf{r}) \frac{1}{|\mathbf{r}-\mathbf{r}'|} n(\mathbf{r}'). \end{aligned}$$
(12)

is the Hartree energy. The main difficulty is to find good approximations to the unknown functional \(E_{xc}[n]\). One of the best known is the local-density approximation (LDA), in which \(E_{xc}[n]\) is replaced by its expression for a homogeneous interacting electron gas (HEG) with electron density equal to \(n(\mathbf{r})\)

$$\begin{aligned} E_{xc}[n] \sim \int d\mathbf{r} \; \epsilon _{xc}^{\mathrm{HEG}}(n(\mathbf{r})) \; n(\mathbf{r}). \end{aligned}$$
(13)

The interacting electron gas problem has been solved numerically via quantum Monte Carlo (QMC) [27], and analytical expressions of \(\epsilon _{xc}^{\mathrm{HEG}}(n(\mathbf{r}))\), obtained fitting the numerical data, are known [28]. Beside the LDA, a large number of alternative functionals are used in practice, often devised to cure shortcomings of the LDA in specific areas. In fact, the LDA has been placed at the bottom of Jacob’s ladder of density functional approximations [29, 30]. For strongly-correlated system, however, none of the known functionals solve the key problems,Footnote 2 as we will discuss later; hence we will refer often to the LDA as representative functional.

Fig. 2
figure 2

The KS band structure (LDA) of the insulating perovskite LaMnO\(_3\). The Mn d bands are yellow. In the KS picture, the system not is only metallic, but also is in the electronic configuration \(t_{2g}^4\), instead than \(t_{2g}^3e_g^1\). Results discussed in Refs. [25, 26]

The Kohn–Sham construction is in principle only a scheme for calculating the ground-state energy and density. To the general surprise, in the early days of DFT, it became quickly clear that the Kohn–Sham eigenvalues are in many cases excellent approximations to the actual eigenenergies of a given material. Since then, Kohn–Sham Hamiltonians have been very successfully used for studying and explaining all kind of properties and systems. Together with the successes, problems emerged. There are entire classes of materials for which this practice fails qualitatively due to strong local electron–electron repulsion effects. An example is shown in Fig. 2. Many transition-metal oxides have phases in which they behave experimentally as a paramagnetic insulator, but they are instead metallic in the Kohn–Sham picture, if the Kohn–Sham bands are calculated with the LDA, the GGA or similar functionals. An insulator is only obtained in the magnetically ordered phase, i.e., if spin-polarized calculations are performed. The origin of this failure of the KS picture is fundamental in nature. It is the same reason why a partially filled band, if no interaction breaks Kramers degeneracy, always describes a metallic system, no matter what the shape of the band is.

Strictly speaking, it should not be a surprise that the KS picture fails in some cases, since the KS eigenvalues were never supposed to be interpreted as elementary excitations. The purist approach, in this situation, consists in returning to the foundations of DFT, viewing the Kohn–Sham eigenvalues as the Lagrange parameters they are, and focusing on the total energy and other quantities that can indeed be calculated exactly, if the exact exchange-correlation functional is known. While this is a fully valid approach, it also limits the type of problems that can be addressed to those that can be solved by calculating the total energy and electron density of the ground state. Thus this strategy will not be discussed further in this article.

Here we will focus instead on less rigorous but more flexible methods. The beauty and simplicity of the KS picture and its great successes create indeed hope that the KS eigenvalues and eigenvectors represent a good starting point, even when they are not sufficient alone. In this view, it should be possible to describe the missing effects as corrections of some kind. Which, however? There are two different strategies that can be adopted to make progress. The first consists in finding, even stepping outside the rigorous principles of DFT, corrections to the exchange-correlation functional, for example in order to obtain a KS gap in Mott insulators; this approach is adopted, e.g., to the DFT+U technique, which we will discuss in Sect. 6. This choice has the advantage that the mapping to an effective one-electron model, with all its simplicity, is preserved. If it works, it is very effective. Unfortunately only some aspects of strong correlations can be recast into an independent quasi-electron problem. The second strategy consists in interpreting the KS orbitals as an optimal basis for constructing the general many-body Hamiltonian. The crucial step is building from it minimal many-body models, approximate but as realistic as possible, that can be solved with state-of-the-art many-body methods.Footnote 3 This is the viewpoint taken in the DFT+DMFT method and its extensions, illustrated in Sects. 8 and 13.

4 Strongly-correlated compounds

In his Lectures on Physics [31], Feynman writes

If, in some cataclysm, all of scientific knowledge were to be destroyed, and only one sentence passed on to the next generations of creatures, what statement would contain the most information in the fewest words? I believe it is the atomic hypothesis (or the atomic fact, or whatever you wish to call it) that all things are made of atoms–little particles that move around in perpetual motion, attracting each other when they are a little distance apart, but repelling upon being squeezed into one another. In that one sentence, you will see, there is an enormous amount of information about the world, if just a little imagination and thinking are applied.

Atoms are not only one of the crucial discoveries in science, but also the simplest systems in which strong correlations play a fundamental role. The non-interacting hydrogen-like atomic problem can be solved analytically, and this is sufficient to explain the main trends in the periodic table, progressively filling the atomic shells. This is not the full story, however. The electronic Hamiltonian for an atom, in the non-relativistic limit, reads

$$\begin{aligned} {\hat{H}}_e=\sum _{L\sigma } \varepsilon _L c^\dagger _{L\sigma } c_{L\sigma } +\frac{1}{2}\sum _{\sigma \sigma '} \sum _{LL'L''L'''} U_{LL'L''L'''} c^{\dagger }_{L\sigma } c^{\dagger }_{L' \sigma '} c_{L'''\sigma '} c_{L''\sigma }. \end{aligned}$$
(14)

Here \(L=nlm\) is a shortcut for the quantum numbers n (principal) and lm (angular momentum), \(\varepsilon _L=-{Z^2}/{2n^2}\), where Z is the atomic number, and \(U_{LL'L''L'''}\) is the Coulomb interaction tensor

$$\begin{aligned} U_{LL'L''L'''} =\int d \mathbf{r}_1 \int d \mathbf{r}_2\; \overline{\phi _{L} (\mathbf{r}_1)} \overline{\phi _{L'} (\mathbf{r}_2)} \frac{1}{|\mathbf{r}_1-\mathbf{r}_2|} {\phi _{L''} (\mathbf{r}_1)} {\phi _{L'''} (\mathbf{r}_2)}. \end{aligned}$$
(15)

The latter gives rise to different effects. The average static mean-field Coulomb interaction merely modifies the Coulomb potential and hence the energy levels \(\varepsilon _L\). This has an interest of its own, because it can break the “accidental” O(4) symmetry of the hydrogen atom, so that s and p levels with the same principal quantum number n can have different energies, unlike in the case of hydrogen. It does not change the overall picture, however. The qualitative new aspect brought in by the Coulomb interaction is the formation of multiplet states. The latter, if we neglect relativistic effect, are characterized by a specific total angular momentum \(L_{\mathrm{T}}\) and a specific total spin \(S_{\mathrm{T}}\), and can be labeled as \(^{2S_{\mathrm{T}}+1}L_{\mathrm{T}}\). Consequences of the Coulomb interaction (15) are, e.g., the first and second Hund’s rule determining the ground multiplet. The existence of atomic/ionic ground multiplets is essential to many co-operative phenomena in strongly-correlated lattices, from magnetic and orbital order to the metal–insulator transition and the Kondo effect. The emergence of local spins can be seen already in the case of an idealized one-level atom

$$\begin{aligned} {\hat{H}}_e=\varepsilon _d \sum _\sigma c^\dagger _{d\sigma } c_{d\sigma } +U{\hat{n}}_{d\uparrow } {\hat{n}}_{d\downarrow }, \end{aligned}$$
(16)

where \(\varepsilon _d\) is the on-site energy, which we can take as the energy zero for simplicity, and U the Coulomb repulsion. This Hamiltonian can be rewritten as

$$\begin{aligned} {\hat{H}}_e=\varepsilon _d \sum _\sigma c^\dagger _{d\sigma } c_{d\sigma } +U\bigg (-{\hat{S}}^2_z+\frac{{\hat{N}}_e^2}{4}\bigg ), \end{aligned}$$
(17)

where \({\hat{S}}_{z} =({\hat{n}}_{d\uparrow }-{\hat{n}}_{d\downarrow })/2\) is the total spin and \({\hat{N}}_e ={\hat{n}}_{d\uparrow }+{\hat{n}}_{d\downarrow }\) is the total number of electrons operator. This shows clearly that at half filling (\(N_e=1\)) such an idealized atom behaves as an effective spin 1/2.

Strongly-correlated materials are typically systems which retain the complexity of atomic (or molecular) physics in the crystalline state. As a rule of thumb, this happens when the building blocks are ions with localized open shells (usually f and d, sometimes p), and the crystal structure is such that the bands stemming from these shells are narrow. There are different families of strongly-correlated materials. The first group includes 4f and 5f electron compounds exhibiting either heavy fermion or Kondo behavior [32,33,34,35,36,37,38]. These systems can host fermionic quasi-particles orders of magnitude heavier than in conventional metals [32]. They can present quantum critical points and non-Fermi liquid behavior [39], and can become unconventional superconductors at low temperature [40]. Examples are CeAl\(_3\), CeCu\(_2\)Si\(_2\), CeCu\(_{6-x}\)Au\(_x\), YbRh\(_2\)Si\(_2\), UBe\(_{13}\) and UPt\(_3\). The second group includes 3d, 4d, and sometimes 5d compounds with partially occupied d bands. The most well known systems in this class are high-temperature superconducting cuprates. Other representative cases are titanates, manganites, cuprates, and ruthenates, among which LaTiO\(_3\), LaMnO\(_3\), KCuF\(_3\) and Sr\(_2\)RuO\(_4\). Additional systems of this kind are iron-based superconductors, iridates and nichelates. All these materials exhibit fingerprints of Mott physics, ranging from enhanced quasi-particle masses to a fully blown Mott insulating phase [41,42,43,44]. They can show in addition Hund’s metal behavior, spin, charge and orbital order, spin and orbital frustration, spin- and orbital-liquid phenomena. Finally, a third group of correlated systems collects molecular crystals, such as C\(_{60}\) compounds [45]; here the narrow bands stem from the small inter-molecular hopping integrals. The three families of systems just described have in common the fact that their unexpected electronic properties arise from the competition between hopping integrals, favoring electron delocalization, and on-site Coulomb interaction, fostering atomic-like behavior. The anomalous properties of correlated materials are usually very sensitive to small changes in external parameters, such as temperature and pressure. In the rest of the article we will focus on the representative class of transition-metal compounds.

5 It is not all about the gap

A classical example of strong-correlation effects is the Mott metal–insulator transition, characteristic feature of many transition–metal compounds. From the point of view of the mechanism, the Mott transition is well captured by the single-band Hubbard model

$$\begin{aligned} {\hat{H}}_e=-\sum _\sigma \sum _{ ii'} t^{i,i'} c^\dagger _{i \sigma } c_{i^\prime \sigma } +U\sum _i {\hat{n}}_{i\uparrow }{\hat{n}}_{i\downarrow }, \end{aligned}$$
(18)

where \(t^{i,i'} \) are the hopping integrals and U is the on-site Coulomb term. At half filling, for \(U=0\), this model describes a conventional metal with band-width W. Increasing the ratio U/W, the system becomes a strongly-correlated electron liquid. When the on-site Coulomb repulsion U is large with respect to W the electrons tend to stay as far apart as possible and the expectation value of double occupations, \(\langle n_{i\uparrow } n_{i\downarrow }\rangle \), correspondingly decreases. Eventually, for U larger than a critical \(U_c\), the system becomes an insulator. When \(W{=}0\), the model describes an insulating collection of decoupled atoms.

Since the most remarkable effect described by (18) is the formation, for \(U>U_c\) and at half filling, of a paramagnetic insulating state in a system which, in the independent-electron picture, should be metallic, the focus of the discussion often goes to the gap in the spectral function. The emphasis on the charge gap can be, however, misleading. Indeed, the gap is only one of the very surprising electronic properties characterizing a system described by Hamiltonian (18). At half filling, the list includes low-energy quasi-particles with enhanced masses and short lifetimes, along with their co-existence with Hubbard bands. It extends to local-moment behavior and low-temperature super-exchange-driven magnetic order, and the presence of both atomic and delocalized electron signatures in response functions. Away from half-filling, a doped Mott insulator is a metallic system with unusual transport and magnetic properties, it can exhibit superconductivity, Nagaoka-like ferromagnetism, and much more. Although we cannot solve (18) in the general case, we know from approximate solutions in specific regimes and/or limit cases that the complex behavior emerging from it is quite different from the one expected for a simple band metal or from a conventional doped band insulator.

Fig. 3
figure 3

The Born–Mayer repulsion is key for understanding the orbitally-ordered structure of KCuF\(_3\). Symbols indicate structures at specific temperatures, pressures or for different cations. The figure shows the change in total energy as a function of the short Cu–F bond (s) for structures with different lattice constants a. The minimum \(s_{\mathrm{min}}\) is practically independent of the temperature so that the experimentally observed increase in \(\varDelta =(l-s)/2\), where l is the long Cu–F bond, is merely an effect of thermal expansion. Reprinted figure with permission from Sims et al. [46]. Copyright (2017) by the American Physical Society

The properties of real correlated systems are, of course, in general not captured by the simple version of the Hubbard model given in Eq. (18). First, in a given material there can be more than one strongly correlated band. Second, other interactions beside the on-site Coulomb repulsion can play an important role. The behavior actually observed in experiments is typically the result of the interplay of several of them: electron–lattice coupling, crystal-field splittings, spin–orbit interaction, just to mention a few. This is of course not a surprise—it is true for any material, independently on the degree of correlation. For weakly correlated systems, such a complexity can typically be taken into account and disentangled in full, however, exploiting the Kohn Sham picture. In a correlated system, where the Kohn Sham picture fails at the qualitative level, it is instead often very tricky to identify the true causes of a specific phenomenon and/or to explain a given observation. A metal–insulator transition could, for example, arise from a pure Mott mechanism, but could also be driven by lattice distortion reducing the symmetry, or could be the co-operative effect of all this together. An example is the case of single layered ruthenates [47,48,49]. Often chicken-and-egg problems arise, as in the case of orbital-ordering. For the latter, it was long debated if the ordering arises from the electron–lattice interaction (Jahn–Teller coupling) or from a purely electronic many-body super-exchange mechanism. This riddle was solved only recently, after decades of hot discussions [25, 46, 50,51,52]. It was shown that super-exchange yields a large transition temperature, but is not sufficient to explain alone the persistence of distortions till very high temperatures; in the case of ionic systems such as KCuF\(_3\), a new mechanism was identified [46], with the Born–Mayer repulsion playing a key role in determining the actual experimental structure (Fig. 3). Given this complexity, the classification of a system as strongly correlated is justified not by a single observation, e.g., the mere existence of a Mott-like gap, but via a series of experimental results which coherently fit within the phenomena described by the Hubbard model or its generalizations. It is this coherent picture, built, e.g., collecting strong-correlation fingerprints in families of similar compounds or in different energy and temperature regimes, which makes it clear that it is necessary to take the local Coulomb interaction explicitly into account—not the fact that one specific isolated experiment in a specific system does not agree with calculations based on the KS Hamiltonian. This makes also clear that descriptions trying to circumvent the introduction of the Hubbard U, in order to become valid alternatives, must, of course, satisfy the same requirements, and capture the global picture.

Fig. 4
figure 4

Quasiparticles and Hubbard bands in the doped 2D Hubbard model (square lattice) with dispersion \(\varepsilon (\mathbf{k})=-2t (\cos k_x+ \cos k_y)+4t'\cos k_x \cos k_y\), obtained with the dynamical mean-field theory approach [53, 54]. The special points are \(\varGamma =(0,0,0)\), \(X=(0,\pi /a,0)\), \(M=(\pi /a,\pi /a,0)\), and Z = \((0,0,\pi /a)\). Left panels: \(t'/t=0.2\). Right panels: \(t'/t=0.4\). Calculations were performed for \(U=7\) eV, t = 0.4 eV, 290 K. The quantum impurity solver adopted is the Hirsch–Fye quantum Monte Carlo [55] method in the implementation of Ref. [53]. The number of holes is indicated with x

The complexity of Hubbard-U effects is illustrated, e.g., in Fig. 4 for the single-band Hubbard model. The figure shows the evolution of the \(\mathbf{k}\)-resolved spectral function as a function of hole doping for the 2-dimensional square lattice. The non-interacting dispersion adopted is

$$\begin{aligned} \varepsilon _\mathbf{k}=-2t(\cos k_x+\cos k_y)+4t'\cos k_x \cos k_y. \end{aligned}$$
(19)

This is a good approximation of the generic \(x^2-y^2\) band stemming from CuO\(_2\) planes in high-temperature superconducting cuprates [56, 57]. The figure shows the effect of introducing x holes in the CuO\(_2\) planes. For \(x\ll 1\) the spectral function shows a quasi-particle band close to the Fermi level as well as Hubbard bands at higher energy; the latter lose intensity with increasing x. The finite lifetime of quasi-particles is reflected in the intensity of the lines.

Going back to rigorous DFT, the theory gives in principle the exact ground-state gap, the difference between ionization energy and electron affinity

$$\begin{aligned} E_{\mathrm{gap}}=E(N_e+1)+E(N_e-1)-2E(N_e), \end{aligned}$$
(20)

if the exact exchange-correlation functional is used. The Kohn–Sham picture does not, however. The KS gap is defined as

$$\begin{aligned} E^{\mathrm{KS}}_\mathrm{gap}=\varepsilon _{\gamma +1}(N_e)-\varepsilon _{\gamma } (N_e), \end{aligned}$$
(21)

where \(\gamma \) labels the highest occupied Kohn–Sham orbital and \(\gamma +1\) the lowest unoccupied Kohn–Sham orbital. It has been shown long ago [58, 59] that the KS gap can differ sizably from the exact value, due to the discontinuities in the reference potential. In fact

$$\begin{aligned} E_{\mathrm{gap}}=E^{\mathrm{KS}}_{\mathrm{gap}}+\underbrace{ \frac{\delta E_{xc}}{\delta n}\Bigg |_{N_e+\delta } -\frac{\delta E_{xc}}{\delta n}\bigg |_{N_e-\delta }}_{\varDelta V_{xc}}. \end{aligned}$$
(22)

The extra term \(\varDelta V_{xc}\) is large in semiconductors [60], and can be the dominant contribution in strongly correlated systems, as it was shown explicitly for, e.g., a chain of idealized hydrogen atoms [61] or the two-site Hubbard model [62]. A characteristic of Mott insulators is that \(E^\mathrm{KS}_{\mathrm{gap}}=0\) in the LDA, the GGA, or similar approximations. The gap is zero even for the exact exchange-correlation functional –if Kramers degeneracy is not broken in the ground state. In discussing the KS gap, it is in fact important to remember that Mott insulators have typically a magnetically ordered ground state. In broken symmetry phases, a finite KS gap, even if often too small, is typically already obtained in calculations based on the local-spin-density approximation (LSDA) or other spin-density functionals, and thus it is likely to be obtained also with the exact functional. Hence, if one focuses on the broken symmetry ground-state only, the KS Hamiltonian, calculated with a well chosen or corrected functional, can give an insulator; this is the basis of the success of the DFT+U approach, discussed in the next section. In order to understand the complex electronic properties of Mott insulators this is not sufficient, however. In the KS description the gap is present only in the broken symmetry phase; furthermore the effects of hole/electron doping are the same as in regular band insulators and the magnetic excitation spectrum is qualitatively wrong. The real problem is thus the fact that the Kohn–Sham description fails to provide the global picture; the absence of a gap in the paramagnetic phase is not the only issue, it is merely a paradigmatic example of this failure, the one to which is perhaps simpler to point to.

6 Static mean-field theory and the DFT+U approach

The DFT+U method [63,64,65,66,67] is one of the approaches based on the strategy of improving the effective Kohn–Sham potential, given in Eq. (9). The correction is based on the physical insight that the shortcomings of the KS Hamiltonian, calculated with LDA, its spin-polarized extension, LSDA, or similar functionals arise from the inadequate treatment of the local Coulomb interaction, the Hubbard U. In this view, the KS Hamiltonian \({\hat{H}}^{\mathrm{KS}}_e\) should be augmented by a Hubbard-like correction

$$\begin{aligned} {\hat{H}}^{\mathrm{KS}}_e\longrightarrow {\hat{H}}^\mathrm{KS}_e+\varDelta {\hat{H}}_e^U \end{aligned}$$
(23)

where \(\varDelta {\hat{H}}_e^U\) is a Hubbard-like augmentation operator. Assuming for simplicity that the latter acts only on the electrons with quantum numbers m in, e.g., shell d,

$$\begin{aligned} \varDelta {\hat{H}}_e^U=\underbrace{\frac{1}{2} U\sum _i\!\!\!\sum _{m\sigma \ne m^\prime \sigma ^\prime }{\hat{n}}_{i m \sigma } {\hat{n}}_{i m^\prime \sigma ^\prime } }_{{\hat{H}}_U} -{\sum _i{\hat{H}}^i_{\mathrm{DC}}}. \end{aligned}$$
(24)

The term \({\hat{H}}_{\mathrm{DC}}\) is a double-counting (DC) correction, subtracting the effects of \({{\hat{H}}_U}\) already included in \({\hat{H}}^{\mathrm{KS}}_e\) via the Hartree and exchange-correlation potential. Here, for simplicity, we approximate it with the fully localized limit (FLL) expression

$$\begin{aligned} {\hat{H}}^i_{\mathrm{DC}}\sim \frac{1}{2}U N_d {\hat{N}}_{id}-\mu _{\mathrm{AT}} {\hat{N}}_{id}, \end{aligned}$$
(25)

where \(\hat{N_{id}}=\sum _{m\sigma } {\hat{n}}_{i m\sigma }\) yields \(N_d\), the number of d electrons at site i, \(\mu _{\mathrm{AT}}=U/2\) is the atomic chemical potential at half filling and the expectation value \(\frac{1}{2}U N_d N_d\) is the Hartree energy.

In DFT+U, the many-body Hubbard term in Eq. (23) is approximated with an effective single-electron operator, i.e., a modification of the parameters of the KS Hamiltonian. This is done via static mean-field decoupling.Footnote 4 The static mean-field Hamiltonian obtained from Eqs. (23) and (24) in this way is

$$\begin{aligned} {\hat{H}}^{\mathrm{MF}}_e= {\hat{H}}_e^{\mathrm{KS}}+\sum _{im\sigma } {\varDelta \varepsilon }_{im\sigma }{\hat{n}}_{im\sigma }, \quad \mathrm{with} \quad {\varDelta \varepsilon }_{im\sigma }=U\left( \frac{1}{2}-\langle {\hat{n}}_{im\sigma }\rangle \right) . \end{aligned}$$
(26)

In \({\hat{H}}^{\mathrm{MF}}_e\) the original levels \(\varepsilon _{im\sigma }\) are shifted by \({\varDelta \varepsilon }_{im\sigma }\), and

$$\begin{aligned} \varDelta \varepsilon _{im\sigma }= \left\{ \begin{array}{cc} -\frac{U}{2} &{}\quad \mathrm{if } \;\langle {\hat{n}}_{im\sigma }\rangle =1 \\ +\frac{U}{2} &{}\quad \mathrm{if }\; \langle {\hat{n}}_{im\sigma }\rangle =0 . \end{array}\right. \end{aligned}$$
(27)

Assuming that one spin state is empty and the other full, the associated shifts correspond to the poles of the atomic Green function at half filling

$$\begin{aligned} E(N){-}E(N{-}1)-\mu =&-U/2 \end{aligned}$$
(28)
$$\begin{aligned} E(N{+}1){-}E(N)-\mu =&+U/2, \end{aligned}$$
(29)

and thus to the center of the Hubbard bands. The modified (U-dependent) total energy “functional” giving Eq. (26) is

$$\begin{aligned} E_{\mathrm{LDA+U}}[n]=E_{\mathrm{LSDA}}[n]+ \sum _i \left[ \frac{1}{2} U \sum _{m\sigma \ne m^\prime \sigma ^\prime }\langle {\hat{n}}_{im\sigma }\rangle \langle {\hat{n}}_{im^\prime \sigma '}\rangle -E_{\mathrm{DC}}\right] , \end{aligned}$$
(30)

where

$$\begin{aligned} E_{\mathrm{DC}}=\frac{1}{2}U N_d(N_d-1). \end{aligned}$$
(31)

One can verify with Janak’s theorem [68] that indeed

$$\begin{aligned} \varepsilon _{im\sigma }^{\mathrm{LDA+U}}=\frac{\partial E_\mathrm{LDA+U}}{\partial \langle {\hat{n}}_{im\sigma }\rangle }= \varepsilon _{im\sigma }^{\mathrm{LSDA}}+U\left( \frac{1}{2}-\langle {\hat{n}}_{im\sigma }\rangle \right) = \varepsilon _{im\sigma }^\mathrm{LSDA}+\varDelta \varepsilon _{im\sigma }. \end{aligned}$$
(32)

Since U is a parameter, not determined univocally by the density, calculations based on the DFT+U approach are not fully ab initio and do not follow rigorously the principles of DFT. The correction yields however a better KS description of the magnetic ground state of Mott insulators, with respect to the one obtained via more rigorous DFT-based approaches. In DFT+U the gap calculated from the eigenvalues of the modified Kohn–Sham Hamiltonian (26) increases linearly with U as in the exact solution in the large U/W limit; thus the term \(\varDelta V_{\mathrm{xc}}\) in Eq. (22) is less relevant. The DFT+U approach has been generalized to account for additional elements of the Coulomb-interaction tensor, e.g., the Hund’s rule coupling J or the coupling between first neighbors, V. Furthermore, different types of DC corrections are used in different regimes/situations [69,70,71,72,73,74,75]. Hence, in practice, the name DFT+U collects a plethora of different Hartree–Fock-like corrections to the LSDA or similar functionals, all sharing however the principles just outlined.

The main strength of DFT+U, the ability of taking into account key effects of the local Coulomb repulsion via a mere modification of the parameters of the KS Hamiltonian, is also its main limit. The correction corresponds to a static mean field treatment of the Hubbard U. Thus, while it well describes the magnetic ground state of correlated materials, is not suitable for studying the metal-insulator transition itself. This failure can be understood in a simple way. In the absence of static local spin and orbital polarization, i.e., if

$$\begin{aligned} \langle {\hat{n}}_{im\sigma }\rangle =\frac{n_i}{2N_d}, \end{aligned}$$
(33)

the DFT+U correction is merely a level shift, identical for all d electrons. Hence, if the original KS Hamiltonian describes a metal, everything else remaining the same (in particular space group, magnetic group and primitive cell), so does the DFT+U Hamiltonian. Conversely, when an insulating state is obtained in DFT+U, it is of Slater- (and not of Mott) type. It requires a reduction of symmetry, e.g., via spin or charge disproportionation, and/or, depending on the system, orbital disproportionation, long-range spin and orbital order or, generalizing, spin-glass-like behavior.

As we have already pointed out, Mott insulators have typically a magnetically ordered ground state; it is perfectly justified to use DFT+U to describe it. Problems arise if the approach is used to analyze the behavior of a Mott insulator above the magnetic transition temperature, usually a paramagnetic insulating phase. An (ideal) local-moment paramagnet is a system in which local magnetic moments fluctuate in time. Hence \(\langle S_z^i\rangle =\langle S_x^i\rangle =\langle S_y^i\rangle =0\) at each site, but at the same time \(\langle \mathbf{S}_i\cdot \mathbf{S}_i\rangle \sim S(S+1)\); the latter is the local moment appearing, e.g., in the expression of the Curie-Weiss susceptibility. A static potential, even if orbital, site and spin dependent, fails to capture this behavior (and several associated phenomena). In fact, a paramagnetic insulator is very different from, e.g., a system characterized by spatial fluctuations of \(\langle S_z^i\rangle \) in a supercell (no matter how large this cell is or how quasi-random the fluctuations are).Footnote 5 This remains true even if \(\sum _i\langle S_z^i\rangle =0\) for the unit cell.Footnote 6 Furthermore, an ideal paramagnet and an ideal disordered spin system can be distinguished experimentally. While disordered systems do exist also among strongly-correlated materials, in order to capture the correct behavior of a paramagnetic insulator one needs a method which can explicitly account for dynamical fluctuations.

Let us now examine the interplay between the DC correction and charge self-consistency. In order to keep it simple, always with transition-metal oxides in mind, we consider a toy model with two sites, the first representing a transition metal ion and the second an oxygen ion. The model is

$$\begin{aligned} \hat{H_e}=&\hat{H_0}+\hat{H_U}\nonumber \\ \hat{H_0}=&\varepsilon _d \sum _{\sigma } c^\dagger _{d\sigma } c_{d\sigma } +\varepsilon _p \sum _{\sigma } c^\dagger _{p\sigma } c_{p\sigma } -t_{pd}\sum _\sigma (c^\dagger _{d\sigma }c_{p\sigma }+c^\dagger _{p\sigma }c_{d\sigma })\nonumber \\ \hat{H_U}=&Un_{d \uparrow } n_{d\downarrow } \end{aligned}$$
(34)

with \(\varepsilon _p=\varepsilon _d-\varDelta _{pd}\). The difference \(\varDelta _{pd}>0\) is the charge-transfer energy. For three electrons the exact ground doublet of this Hamiltonian is

$$\begin{aligned} |G,\sigma \rangle = a_1\, c^\dagger _{d\sigma } c^\dagger _{p\uparrow }c^\dagger _{p\downarrow }|0\rangle +a_2\, c^\dagger _{d\uparrow }c^\dagger _{d\downarrow } c^\dagger _{p\sigma }|0\rangle , \end{aligned}$$
(35)

where \(a_2^2+a_1^2=1\) and

$$\begin{aligned} a_1^2=\frac{t_{pd}^2}{ t_{pd}^2 +\bigg (\frac{\varDelta _{pd}+U}{2}-\sqrt{\big (\frac{\varDelta _{pd}+U}{2}\big )^2+t_{pd}^2}\bigg )^2 }. \end{aligned}$$
(36)

Its energy is

$$\begin{aligned} E_G(3)=3\varepsilon _d-2\varDelta _{pd} +\frac{\varDelta _{pd}+U}{2}-\sqrt{\bigg (\frac{\varDelta _{pd}+U}{2}\bigg )^2+t_{pd}^2}, \end{aligned}$$
(37)

and the total number of d electrons is

$$\begin{aligned} n_d=1+ \frac{\bigg (\frac{\varDelta _{pd}+U}{2}-\sqrt{\big (\frac{\varDelta _{pd}+U}{2}\big )^2+t_{pd}^2}\bigg )^2}{ t_{pd}^2 +\bigg (\frac{\varDelta _{pd}+U}{2}-\sqrt{\big (\frac{\varDelta _{pd}+U}{2}\big )^2+t_{pd}^2}\bigg )^2 }. \end{aligned}$$
(38)

For \(U=0\), the total number of d electrons takes the value \(n_d^0=1+x_0\), where \(x_0\) is a positive number, which can be determined by setting \(U=0\) in right term in Eq. (38); for \(U\rightarrow \infty \), we find instead \(n_d^\infty =1\).

Let us now analyze the effect of the DC correction in a mean-field treatment of the Coulomb term. In the first step we calculate the KS Hamiltonian. For simplicity we assume that this yields \(\varepsilon _d\longrightarrow \varepsilon ^{\mathrm{KS}}_d= \varepsilon _d+U_{\mathrm{KS}}\; n_d\). When \(U_{\mathrm{KS}}=U\), the correction is the Hartree term for the Hubbard model. From Eq. (38) one may see that, if \(U=0\), the exact ground state density for finite U is recovered by merely replacing \(\varepsilon _d\) with \( \varepsilon _d+U\); indeed, this \(N=3\) electron problem is an uncorrelated system in the hole representation. In the second step, we augment \(\hat{H^{\mathrm{KS}}_e}\) with \(\varDelta \hat{H^\mathrm{U^*}_e}={\hat{H}}^{U^*}_e-{\hat{H}}_{\mathrm{DC}}\), and solve the problem in the static mean-field approximation; we replace U with \(U^*\), since, in general, we do not know the exact value of U included in the Hartree term. Thus

$$\begin{aligned} \varDelta _{pd} \longrightarrow&\varDelta ^*_{pd}(n_d)=\varDelta _{pd}+Un_d+\frac{U}{2}^*(n_d-n_d^*). \end{aligned}$$
(39)

The choice of \(n_d^*\) defines the DC correction. The self-consistent condition for the number of d electrons is obtained setting \(U=0\) in Eq. (38) and replacing \(\varDelta _{pd}\) with \(\varDelta _{pd}^*\) via Eq. (39). In the FLL limit, i.e., when the double-counting energy is given by Eq. (31), for \(n_d>1\) we obtain \(n_d<n_d^*\); hence, the double-counting term tends to slightly increase the self-consistent value of \(n_d\). The opposite happens if, e.g., we set \(n_d^*=1\). The example shows that the largest effect can arise, however, from the uncertainty in the value of \(U^*\), rather than from the specific choice of \(n_d^*\).

Besides DFT+U, there are many attempt to address the problem of strong correlation by correcting the functional [28, 76,77,78,79,80,81,82,83,84]. Examples are, e.g., the self-interaction-correction approach [28, 77, 78] or the hybrid functional method [79,80,81,82,83,84]. Other approaches try instead to propose schemes alternative to the Kohn–Sham construction, such as the strictly-correlated-electron approach [85, 86] and/or focus on the total energy [87].

7 Model Hamiltonians

In this section, we will discuss the second family of approaches, the one that uses the Kohn–Sham states as the optimal one electron basis for constructing many-body models. The Hamiltonian (1), with all terms explicitly written, has the form

$$\begin{aligned} {\hat{H}}_e=-\frac{1}{2}\sum _i \nabla ^2_i - \sum _i \sum _\alpha \frac{Z_\alpha }{|\mathbf{r}_i{-}{} \mathbf{R}_\alpha |} + \sum _{i>j} \frac{1}{|\mathbf{r}_i{-}{} \mathbf{r}_j|} +\sum _{\alpha >\alpha '} \frac{Z_\alpha Z_{\alpha '}}{|\mathbf{R}_\alpha {-}{} \mathbf{R}_{\alpha '} |} , \end{aligned}$$
(40)

where \(\{\mathbf{r}_i\}\) are electron coordinates, \(\{\mathbf{R}_\alpha \}\) nuclear coordinates and \(Z_\alpha \) the nuclear charges. Using a complete one-electron basis \(\{\phi _a(\mathbf{r})\}\), where \(\{ a\}\) are quantum numbers, we can rewrite it in second quantization as

$$\begin{aligned} {\hat{H}}_e=\underbrace{-\sum _{ab} t_{ab} c_a^{\dagger }c_b}_{{\hat{H}}_0} + \underbrace{\frac{1}{2} \sum _{aa'bb'} U_{aa'bb'} \; c_a^\dagger c_{a'}^\dagger c_{b'} c_{b}}_{{\hat{H}}_U}. \end{aligned}$$

The hopping integrals are given by

$$\begin{aligned} {t}_{ab}= -\!\int \! d \mathbf{r} \; \overline{\phi _{a}}(\mathbf{r}) \Big ({-}\frac{1}{2}\nabla ^2\underbrace{-\sum _\alpha \frac{Z_\alpha }{|\mathbf{r}{-}{} \mathbf{R}_\alpha |}}_{v_{\mathrm{en}}(\mathbf{r})}\Big )\,\phi _{b}(\mathbf{r}), \end{aligned}$$

while the elements of the Coulomb tensor are

$$\begin{aligned} {U}_{aa'bb'}= \int \! d \mathbf{r}_2 \!\int \! d \mathbf{r}_2 \; \overline{\phi _{a}}(\mathbf{r}_1)\,\overline{\phi _{a'}}(\mathbf{r}_2)\;\frac{1}{|\mathbf{r}_{1}{-}{} \mathbf{r}_2|}\;\phi _{b'}(\mathbf{r}_2)\,\phi _{b}(\mathbf{r}_1). \end{aligned}$$
Fig. 5
figure 5

(Color online) NMTO Wannier functions showing orbital-order in TbMnO\(_3\), as obtained by LDA+DMFT calculations. Reprinted figure with permission from Flesch et al. [51]. Copyright (2012) by the American Physical Society

All complete one-electron bases are of course equivalent in theory. In practice, since, in the general case, we cannot solve the \(N_e\)-electron problem exactly, some bases are better than others. In this spirit, the Kohn–Sham orbitals \(\{\phi ^{\mathrm{KS}}_a(\mathbf{r})\}\) have many advantages, since they provide KS spectra reasonably in line with experiments for weakly correlated systems. Even for strongly-correlated systems, they are often sufficiently accurate for empty and fully occupied states. In order to use the KS orbitals as a basis, we first rewrite \(v_{\mathrm{en}}(\mathbf{r})\) in terms of the effective KS potential, \(v_{\mathrm{KS}}(\mathbf{r})\). The latter differs from \(v_{\mathrm{en}}(\mathbf{r})\) via the Hartree term \(v_{\mathrm{H}}(\mathbf{r})\) and the (approximate) exchange-correlation contribution, \(v_\mathrm{xc}(\mathbf{r})\)

$$\begin{aligned} {v_{\mathrm{en}}(\mathbf{r})} =v_{\mathrm{KS}}(\mathbf{r})-{v_{\mathrm{H}}(\mathbf{r})} - v_{\mathrm{xc}}(\mathbf{r}). \end{aligned}$$
(41)

We thus introduce the Kohn–Sham hopping integrals

$$\begin{aligned} {\tilde{t}}_{ab}= -\!\int \! d \mathbf{r} \; \overline{\phi ^\mathrm{KS}_{a}}(\mathbf{r}) \Big ({-}\frac{1}{2}\nabla ^2 +{v_{\mathrm{KS}}(\mathbf{r})}\Big )\phi ^{\mathrm{KS}}_{b}(\mathbf{r}). \end{aligned}$$
(42)

With this definition, the total Hamiltonian, in the KS basis, is

$$\begin{aligned} {\hat{H}}_e=\underbrace{-\sum _{ab} {\tilde{t}}_{ab}\, c_a^{\dagger }c_b}_{{\hat{H}}_0={\hat{H}}^\mathrm{KS}_{e}} +\underbrace{\frac{1}{2} \sum _{aba'b'} {\tilde{U}}_{aa'bb'} \; c_a^\dagger c_{a'}^\dagger c_{b'} c_{b} -{\hat{H}}_{\mathrm{DC}}}_{\varDelta {\hat{H}}_U}, \end{aligned}$$
(43)

where \({\hat{H}}_{\mathrm{DC}}\) is the double-counting correction. Typically \(\varDelta {\hat{H}}_U\) is a short range operator since long-range effects are already well captured by the KS potential, and strong correlation effects are a manifestation of a large local U.

The key step for studying correlation effects in materials is building, starting from the general Hamiltonian (43), minimal material-specific models containing all essentials degrees of freedom and interactions. This was made possible by the development of methods for constructing Kohn–Sham localized Wannier functions spanning specific bands [88, 89] and/or projection schemes with similar capabilities.Footnote 7 An example of a Wannier-like state is shown in Fig. 5. These approaches yield single-electron Hamiltonians which very accurately reproduce the Kohn–Sham bands in specific energy windows. A quite different story is the estimate of screened Coulomb integrals. Calculating exact screening effects is in principle as hard as solving the full many-body Hamiltonian. Commonly adopted schemes are the constrained local-density approximation (cLDA) [90, 91], linear-response methods [67, 92, 93], and the constrained random-phase approximation (cRPA) [94, 95]. They have been successfully used for building models and understanding the properties of several families of correlated materials. Still, screened Coulomb parameters are only known within relatively large error bars.

Fig. 6
figure 6

LDA+DMFT spectral function for the single band systems VOMoO\(_4\) (left) and Li\(_2\)VOSiO\(_4\) (right) at 380 K and for \(1<U<5\) eV. The linewidth increases with increasing U in steps of 1 eV. Rearranged from Ref. [53]

It is important to underline that the model construction is not a mechanical procedure. Although modern techniques provide the tools, identifying the relevant degrees of freedom and interactions remains bound to our physical understanding of the problem and system analyzed. A model is typically a work in progress, a compromise between what one would like to include and what can be included in actual calculations. It evolves with time, when new facts come to light. For most transition-metal oxides, typical low-energy starting models, augmented with Coulomb corrections, have the form \({\hat{H}}_e={\hat{H}}^{\mathrm{KS}}_e+\varDelta {\hat{H}}_U\) of a generalized Hubbard Hamiltonian, with

$$\begin{aligned} {\hat{H}}^{\mathrm{KS}}_e&=-\sum _{ii'} \sum _\sigma \sum _{m m^\prime } t_{m\sigma , m^\prime \sigma '}^{i,i'}\, c^\dagger _{im\sigma } c_{i' m^\prime \sigma '} \end{aligned}$$
(44)
$$\begin{aligned} \varDelta {\hat{H}}_U&= \frac{1}{2} \sum _i\sum _{\sigma \sigma ^\prime } \sum _{m m^\prime } \sum _{p p^\prime } U_{m p m^\prime p^\prime }\, c^\dagger _{im\sigma } c^\dagger _{ip\sigma ^\prime } c_{ip^\prime \sigma ^\prime } c_{im^\prime \sigma } -{\hat{H}}_{\mathrm{DC}}, \end{aligned}$$
(45)

with m belonging to the d shell or to a subset of crystal-field states.

Fig. 7
figure 7

Left: correlated band structure of VOMoO\(_4\) and Li\(_2\)VOSiO\(_4\) for a realistic \(U=5\) eV, calculated at \(\sim 200\) K. The dots are the poles of the Green function and yield the energy dispersion. Right: corresponding real-axis self-energy \(\varSigma (\omega )\). Reprinted figure with permission from Kiani and Pavarini [53]. Copyright (2016) by the American Physical Society

8 DMFT and DFT+DMFT

The essential step forward for the theoretical description of the Mott metal-insulator transition in the Hubbard model came with the development of the dynamical mean-field theory (DMFT) [96,97,98,99,100,101,102,103]. The central idea consists in mapping the Hubbard Hamiltonian onto an auxiliary quantum-impurity model. The latter can be solved exactly, differently from the original lattice model. A typical QIM is the Anderson Hamiltonian

$$\begin{aligned} {\hat{H}}^{\mathrm{QIM}}&= \underbrace{\sum _{\mathbf{k}\sigma }\varepsilon ^s_\mathbf{k}{\hat{n}}_{\mathbf{k}\sigma }}_{{\hat{H}}_{\mathrm{bath}}} +\underbrace{\sum _{\mathbf{k}\sigma } \left( V^s_\mathbf{k} c^\dagger _{\mathbf{k} \sigma } c_{d \sigma } +\mathrm{h.c.}\right) }_{{\hat{H}}_\mathrm{hyb}} + \underbrace{\varepsilon _d\sum _\sigma {\hat{n}}_{d\sigma }+U {\hat{n}}_{d\uparrow }{\hat{n}}_{d\downarrow }}_{{\hat{H}}_{\mathrm{imp}}}. \end{aligned}$$

It describes a single correlated d impurity in a bath of non-interacting electrons s. The Anderson model was originally introduced to describe the Kondo effect in diluted metallic alloys [104,105,106]. Solution methods, approximated and numerically exact, were therefore already available when DMFT was introduced. The main difference with respect to the single-impurity case is that in DMFT the parameters of the quantum-impurity model are not known from the start, but obtained via the self-consistency condition, requiring that the quantum impurity self-energy is as close as possible to the local self-energy of the original model. Non-local self-energy terms are neglected. DMFT is exact for \(U=0\), in the atomic limit, in the single-impurity limit, and, most remarkably, in the infinite coordination number limit [96,97,98,99]. For realistic three-dimensional lattices it is an excellent approximation.

Fig. 8
figure 8

Schematic representation of the DFT+DMFT approach. Upper boxes: model building. Central boxes: map to quantum impurity model and quantum impurity solver (QIS). Lower boxes: self-consistency check and calculation of one and two-particle local Green functions. The figure is from Ref. [107], ch. 7

The success of DMFT relies on the fact that it describes the Mott paramagnetic metal to paramagnetic insulator transition, differently from all static mean-field approaches. This is shown in a representative case in Fig. 6 for two systems well described by the single-band Hubbard model. For small U/W, where W is the band width, the spectral function is characterized by a central quasi-particle peak and two Hubbard bands. Increasing U/W the central peak becomes more narrow; at the same time the effective mass of the quasi-electrons increases, and their lifetime shrinks. Eventually the central peak vanishes and the system becomes a paramagnetic insulator. In the large U/W limit, the Mott gap in the spectral function is approximately equal to \(U-W\). The paramagnetic Mott phase is characterized by a self-energy which diverges in the gap, and not by a static site-dependent potential which has the effect of reducing the symmetry, as in static mean-field theory. This is shown, for the same systems of Fig. 6, in the right panels of Fig. 7.

In the last decades, three crucial advances made the application of DMFT to materials a reality. The first is the fact that techniques to build system-specific models became available, as discussed in the previous section. The second is the development of new flexible quantum-impurity solvers, such as continuous-time quantum Monte Carlo [108], enabling us to study always more realistic quantum impurity models. The third is the (so-far increasing) power of massively parallel supercomputers, which made actual calculations possible in practice. All this gave birth to the DFT+DMFT approach [41, 75, 101, 107, 109,110,111,112,113], described in a schematic way in Fig. 8.

Fig. 9
figure 9

Orbital polarization p (left) and occupied state \(|\theta \rangle \) (right) as a function of temperature, calculated with LDA+DMFT. These calculations were instrumental for solving the problem of the origin of orbital ordering in LaMnO\(_3\). Reprinted figure with permission from Pavarini and Koch [25]. Copyright (2010) by the American Physical Society

The initial step of a DFT+DMFT calculation consists in model building (green boxes in Fig. 8). This involves the choice of the states for which the local Coulomb interaction in Eq. (43) is explicitly taken into account. The second step consist in mapping the lattice model to a quantum impurity model (QIM). The third its solution via a quantum impurity solver. Last comes the self-consistency loop. For QMC solvers, the most flexible so far, the quantum impurity model is defined via the associated bath Green function, \({{\mathcal {G}}}^{\sigma } (i\nu _n)\); since QMC solvers all work on the imaginary axis, we express the latter as a function of the fermionic Matsubara frequency \(\nu _n\). Solving the QIM yields the impurity Green function \(G^{\sigma }_{\mathrm{QIM}}(i\nu _n)\). From the local Dyson equation for the QIM we can calculate the impurity self-energy

$$\begin{aligned} \varSigma ^{\sigma }_{\mathrm{QIM}} (i\nu _n)= \big ({{G}}^{0\sigma }_\mathrm{QIM}(i\nu _n)\big )^{-1}- \big (G^{\sigma }_{\mathrm{QIM}}(i\nu _n)\big )^{-1}. \end{aligned}$$

The self-energy of the original Hubbard model is then set equal to the impurity self-energy, so that the local Green function is given by

$$\begin{aligned} G^\sigma _{i,i}(i\nu _n)=&\frac{1}{N_\mathbf{k}} \sum _\mathbf{k} \frac{1}{ i\nu _n+\mu -\varepsilon _\mathbf{k}-\varSigma ^\sigma _{\mathrm{QIM}}(i\nu _n)}, \end{aligned}$$

where \(N_\mathbf{k}\) is the number of \(\mathbf{k}\) points. The local Dyson equation is used once more, this time to calculate a new bath Green function \({{\mathcal {G}}}^{\sigma }(i\nu _n)\), which in turn defines a new QIM. This procedure is repeated until self-consistency is reached, i.e., till a fixed point for the self-energy is found. At the end of the self-consistency loop, further quantities, such as two particle Green functions, can be calculated.

The successes of the DFT+DMFT approach in describing the properties of strongly-correlated materials are beyond doubt. The technique has been instrumental to shed light on many open problems. Thanks to it, one could explain the nature of the metal–insulator transition in representative families of systems [47, 88, 114, 115] or clarifying the origin of orbital ordering [25, 46, 50,51,52]. Other times new categories of systems were identified, such as Hund’s metals [116,117,118,119]. Along the years the DFT+DMFT method has been extended to calculate properties beyond simple spectral functions, to the study of phase transitions [102, 120,121,122], total energies and phonons [123, 124], superconductivity [125,126,127,128,129,130] or spin-orbit effects [48, 49, 131]. The type of systems to which it can be applied has been extended from transition–metal oxides to heterostructures and f electron compounds. A representative, although not exhaustive, list of success stories can be found, e.g., in Refs. [41, 101, 107, 110,111,112,113]. As an example we show in Fig. 9 the super-exchange induced-orbital order transition in LaMnO\(_3\), a calculation instrumental to clarify the origin of orbital ordering in rare-earth manganites [25].

Fig. 10
figure 10

Opening of gap in the Hartree–Fock solution of the single band Hubbard model for a square lattice at half filling. Only the nearest neighbor hopping t is taken into account. Left: band folding due to the doubling of the cell, \(U=0\). Right: gap opening at the Brillouin Zone boundary. The magnetization is defined as \(m=(n_\uparrow -n_\downarrow )/2\). From Ref. [132], ch. 3

Before concluding this section it is important to analyze the relation between the DFT+DMFT and the DFT+U method.Footnote 8 The latter can be seen a special case of DFT+DMFT, in which, instead of solving the quantum-impurity problem exactly, we solve it in the Hartree–Fock (HF) approximation. For the single-band Hubbard model at half filling, approximating the DC correction in the FLL limit, the corresponding local self-energy is

$$\begin{aligned} \varSigma _{\mathrm{HF}}^{\sigma }(i\nu _n)= U \bigg (n_{i-\sigma }-n_i+\frac{1}{2}\bigg )= U \bigg (\frac{1}{2}-n_{i\sigma }\bigg ) \end{aligned}$$
(46)

This yields the DFT+U level shift, Eq.  (27) and to the HF bands shown in Fig. 10. In DMFT, the HF self-energy is merely the large frequency limit of \(\varSigma ^\sigma (i\nu _n)\), however. We can then distinguish two cases. In the paramagnetic phase, where the DMFT self-energy is strongly frequency dependent, the difference is qualitative. This may be seen for a representative case in Fig. 11. The situation changes completely in the magnetic phases. Here the HF term is a leading contribution in the DMFT self-energy, and the difference between HF and DMFT are correspondingly much smaller. Still, although ground state and spectral function are, in this situation, reasonably well described in HF, and hence in DFT+U, response functions are not, as we will see later.

Fig. 11
figure 11

Left: LDA band structure of KCuF\(_3\) around the Fermi level. Only the \(e_g\) bands are shown. Second panel: LDA+DMFT correlated band structure in the orbitally-ordered phase [46, 50]. Dots: poles of the Green function. Third panel: LDA+HF bands. Right: DMFT self-energy matrix in the basis of natural orbitals. Full lines: real part. Dotted lines: imaginary part. Parameters: \(U=7\) eV, \(J=0.9\) eV. Rearranged from Ref. [133], ch. 9

9 Quantum-impurity solvers

A crucial aspect of the LDA+DMFT technique is the method employed to solve the quantum impurity problem, the so-called quantum impurity solver. Quantum impurity solvers can be grouped in two major categories. The first consists of numerically exact ones and the second of all the rest.

Numerically exact solvers include exact diagonalization and Lanczos [134,135,136], various renormalization group methods [6,7,8, 137,138,139,140,141,142,143,144,145,146] and techniques based on quantum Monte Carlo [55, 108, 147, 148]. Approximate solvers encompass a variety of different approaches. They range from Hubbard approximations and rotationally-invariant slave fermions [149], to hybridization-expansion techniques, e.g., the non-crossing [150,151,152,153,154], and the one-crossing approximation [155]. A simple and successful approximate solver is the iterative perturbation theory [156,157,158,159]. The main advantage of approximate techniques is that they are typically computationally not too demanding. The disadvantage is that crucial many-body effects might not be captured. Furthermore, usually, approximate methods work well only within certain limits. Adopting approximate solvers, even powerful and flexible ones, implies that ultimately results have to be checked against numerically exact methods. In this article, to avoid confusion, I therefore adopt the labels DMFT and DFT+DMFT only when I refer to calculations performed with numerically exact solvers.

Fig. 12
figure 12

Scaling of QMC solvers on a massively parallel supercomputer. Black line: Hirsch–Fye (HF) solver, 2 orbitals. The dark and light lines are continuous time hybridization-expansion (CT-HYB) quantum Monte Carlo calculations. Dark lines: Krylov CT-HYB solver with truncation of the local trace (empty symbols, K-t) and without (full symbols, K). Results are for 2 (circles) and 3 (triangles) orbitals. Light lines: segment CT-HYB solver (S), 5-band model (pentagons). All points correspond to calculations of high quality (and with comparable error bars) for the systems considered in this work. For \(T\sim 165\) K the 5-orbital segment solver is about as fast the 3-orbital Krylov with trace truncation or the 2-orbital Krylov without trace truncation, and it is remarkably faster than the 2-orbital HF. Reprinted figure with permission from Flesch et al. [26]. Copyright (2013) by the American Physical Society

It is crucial to underline that even numerically exact solvers are not universal tools. Typically, each of them works well for specific classes of problems. Nevertheless, remarkable progress has been achieved to expand their regime of applicability, combining algorithmic advances with the growing performance of massively parallel supercomputers. The most general purpose solvers are based on QMC. They include the Hirsch–Fye method [55], the interaction-expansion [108, 148] and the hybridization-expansion continuous-time QMC technique [108, 147]. The bottleneck of all QMC approaches is that the computational time can become quickly prohibitively long with increasing the number of degrees of freedom and lowering the temperature. A representative scaling plot for Hirsch–Fye and the hybridization-expansion continuous-time QMC solvers is given in Fig. 12. Furthermore, due to the infamous sign problem, depending on the case, QMC solvers can be limited in the number of orbitals/sites, the type of interactions included in the model, and the lowest temperature that can be reached. Finally, they work on the imaginary axis; comparing with experimental data requires the analytic continuation of the results, an ill-posed problem. This has boosted the search for efficient analytic continuation approaches that can work for noisy data. Among these we find the maximum-entropy technique [160,161,162,163,164], stochastic approaches of various kind [165,166,167], and the recently developed average-spectrum method [168,169,170]. Complementary to QMC solvers are exact diagonalization and Lanczos [134,135,136]. Here, the bottleneck is not the computational time but the actual size of the Hilbert space, which grows very rapidly with the number of sites/orbitals. The advantage is that these approaches provide results on the real axis directly, and yield the \(T\rightarrow 0\) limit, typically not easily accessible with QMC solvers. Other methods that can work on the real axis are renormalization group techniques, such as the numerical renormalization group (NRG) [137,138,139,140], best designed to describe physics at the Kondo energy scale; more recently, the density-matrix renormalization group (DMRG) [6,7,8] was added to the powerful tools for solving the quantum impurity problems in materials [7, 143,144,145,146].

10 Lessons from a toy model: the two-site Hubbard model

The simplest model that can be studied to explain the DMFT approach is the two-site Hubbard model

$$\begin{aligned} {\hat{H}}=\varepsilon _d \sum _{i\sigma } {\hat{n}}_{i\sigma } -t\sum _{\sigma } \left( c_{1\sigma }^\dagger c_{2\sigma } + c_{2\sigma }^\dagger c_{1\sigma }\right) +U\sum _i{\hat{n}}_{i\uparrow }{\hat{n}}_{i\downarrow }, \end{aligned}$$
(47)

with \(i=1,2\). For \(N=2\) electrons (half filling) the ground state of this Hamiltonian is the singlet

$$\begin{aligned} |G\rangle _H=\frac{a_2(t,U)}{\sqrt{2}}\left( c_{1\uparrow }^\dagger c_{2\downarrow }^\dagger -c_{1\downarrow }^\dagger c_{2\uparrow }^\dagger \right) |0\rangle +\frac{a_1(t,U)}{\sqrt{2}}\left( c_{1\uparrow }^\dagger c_{1\downarrow }^\dagger +c_{2\uparrow }^\dagger c_{2\downarrow }^\dagger \right) |0\rangle \, \end{aligned}$$
(48)

The prefactors are given by

$$\begin{aligned} a_1^2(t,U)=&\frac{1}{\varDelta (t,U)}\frac{\varDelta (t,U)-U}{2},&a_2^2(t,U)=&\frac{4t^2}{\varDelta (t,U)}\frac{2}{\varDelta (t,U)-U}, \end{aligned}$$
(49)

where

$$\begin{aligned} \varDelta (t,U)=\sqrt{U^2+16t^2}. \end{aligned}$$
(50)
Fig. 13
figure 13

Spectral function of the Hubbard dimer (H) for \(U/t=3\). Local (left), \(k=0,\pi \) (right). It is obtained from Eq. (52), setting \(\nu _n \longrightarrow \omega +i\delta \), where \(\delta \rightarrow 0\). The spectral function exhibits four peaks, corresponding to the four poles of the Green function. The associated weights are: \(w_+=(1+w(t,U))/4\) for the pole at energies \(E_G(2){-}E_-(1){-}\mu \) and the one at energy \(E_-(3){-}E_G(2){-}\mu \), and \(w_-=(1-w(t,U))/4\) for the remaining two, at energies \(E_G(2){-}E_+(1){-}\mu \) and \(E_+(3){-}E_G(2){-}\mu \). Figure rearranged from Ref. [107], ch. 7

The associated ground-state energy is

$$\begin{aligned} E_G(2)=2\varepsilon _d+\frac{1}{2}\big (U-\varDelta (t,U)\big ). \end{aligned}$$
(51)

The eigenvalues for \(N_e=1\) and \(N_e=3\) electrons are instead

$$\begin{aligned} E_\pm (1)&=\varepsilon _d\pm t\\ E_\pm (3)&=3\varepsilon _d\pm t+U \end{aligned}$$

In the \(T\rightarrow 0\) limit, the local Matsubara Green function for spin \(\sigma \) (see Fig. 13 for the associated spectral function) is then the sum of four terms

$$\begin{aligned} G^\sigma _{i,i}(i\nu _n){=}&\frac{w_+}{i\nu _n-(E_G(2)-E_-(1){-}\mu )} {+}\frac{w_-}{i\nu _n-\big (E_G(2)-E_+(1){-}\mu \big )} \nonumber \\ +\quad&\frac{w_-}{i\nu _n{-}\big (\!-E_G(2){+}E_+(3){-}\mu \big )} {+}\frac{w_+}{i\nu _n{-}\big (\!-E_G(2){+}E_-(3){-}\mu \big )}, \end{aligned}$$
(52)

where \(\nu _n=\pi (2n{+}1)/\beta \) are fermionic Matsubara frequencies, the chemical potential is \(\mu =\varepsilon _d+U/2\), and

$$\begin{aligned} w_\pm&=\frac{1\pm w(t,U)}{4} \end{aligned}$$
(53)
$$\begin{aligned} w(t,U)&=2a_1(t,U)a_2(t,U). \end{aligned}$$
(54)

This can be rewritten as the average of the Green function for the bonding (\(\mathbf{k}=0)\) and the anti-bonding (\(\mathbf{k}=\pi )\) state, i.e.,

$$\begin{aligned} G^\sigma _{i,i}(i\nu _n)=\frac{1}{2}&\Bigg ( \underbrace{ \frac{1}{i\nu _n+\mu -\varepsilon _d+t-\varSigma ^\sigma (0,i\nu _n)} }_{G^\sigma (0,i\nu _n)} {+}\underbrace{ \frac{1}{i\nu _n+\mu -\varepsilon _d-t-\varSigma ^\sigma (\pi ,i\nu _n)} }_{G^\sigma (\pi ,i\nu _n)} \Bigg ). \end{aligned}$$
(55)

The \(\mathbf{k}\)-dependent self-energy is given by

$$\begin{aligned} \varSigma ^\sigma (\mathbf{k},i\nu _n)=&\frac{U}{2}+\frac{U^2}{4}\frac{1}{i\nu _n+\mu -\varepsilon _d-\frac{U}{2}-\; e^{ik}\;3t}. \end{aligned}$$
(56)

The local self-energy is, by definition, the \(\mathbf{k}\)-average

$$\begin{aligned} \varSigma ^\sigma _{l}(i\nu _n)=\varSigma ^\sigma _{i,i}(i\nu _n)&= \frac{1}{2}\bigg (\varSigma ^\sigma (\pi ,i\nu _n){+}\varSigma ^\sigma (0,i\nu _n)\bigg ) \nonumber \\&= \frac{U}{2}{+} \frac{U^2}{4}\frac{i\nu _n+\mu -\varepsilon _d-\frac{U}{2}}{(i\nu _n+\mu -\varepsilon _d-\frac{U}{2})^2-{(3t)^2}}. \end{aligned}$$
(57)

The difference

$$\begin{aligned} \varDelta \varSigma _{l}^\sigma (i\nu _n)&= \frac{1}{2}\Big (\varSigma ^\sigma (\pi ,i\nu _n){-}\varSigma ^\sigma (0,i\nu _n)\Big ) \nonumber \\&=\frac{U^2}{4}\frac{3t}{(i\nu _n+\mu -\varepsilon _d-\frac{U}{2})^2-{(3t)^2}}, \end{aligned}$$
(58)

is the non-local part of the self-energy. The local Green function can then be expressed in the form

$$\begin{aligned} G^\sigma _{i,i}(i\nu _n)=\frac{1}{i\nu _n+\mu -\varepsilon _d-F^\sigma (i\nu _n)-\varSigma ^\sigma _l(i\nu _n)} , \end{aligned}$$
(59)

where

$$\begin{aligned} F^\sigma (i\nu _n)= \frac{\big (t+\varDelta \varSigma _l^\sigma (i\nu _n)\big )^2}{i\nu _n+\mu -\varepsilon _d-{\varSigma ^\sigma _{l}(i\nu _n)}} \end{aligned}$$
(60)

is the hybridization function. For \(U=0\) it becomes

$$\begin{aligned} {{F}}^\sigma _0(i\nu _n)= \frac{t^2}{i\nu _n}. \end{aligned}$$
(61)

The local Green function satisfies a local Dyson equation

$$\begin{aligned} \varSigma ^\sigma _{l}(i\nu _n) =\frac{1}{{\mathcal{G}^{\sigma }_{i,i}(i\nu _n)}}-\frac{1}{G^\sigma _{i,i}(i\nu _n)}, \end{aligned}$$
(62)

where

$$\begin{aligned} {\mathcal{G}^{\sigma }_{i,i}(i\nu _n)}=\frac{1}{i\nu _n+\mu -\varepsilon _d-F^\sigma (i\nu _n)}. \end{aligned}$$
(63)
Fig. 14
figure 14

Spectral functions of the Hubbard dimer (\(U/t=3\)) and the auxiliary Anderson molecule (\(\varepsilon _s=\varepsilon _d+U/2\)) in the zero temperature limit. Left panel: Hubbard dimer (H), exact solution and local self-energy approximation (\(\varDelta \varSigma _l^\sigma (\omega )=0\)). Right panels: Anderson molecule (A) vs local self-energy approximation of the Hubbard dimer. Dashed lines: poles of the Green function of the Anderson molecule (right) or Hubbard dimer (left). Figure rearranged from Ref. [107], ch. 7

So far we discussed the exact solution of the Hubbard dimer. Next, we map the Hubbard dimer into an auxiliary two-site Anderson model

$$\begin{aligned} {\hat{H}}^{\mathrm{A}}=\varepsilon _s \sum _{\sigma }{\hat{n}}_{s\sigma } -t\sum _{\sigma } \left( c_{d\sigma }^\dagger c_{s\sigma } + c_{s\sigma }^\dagger c_{d\sigma }\right) +\varepsilon _d \sum _{\sigma }{\hat{n}}_{d\sigma }+U{\hat{n}}_{d\uparrow }{\hat{n}}_{d\downarrow }, \end{aligned}$$
(64)

where d labels the correlated impurity and s the uncorrelated bath site. The parameters of this quantum-impurity model have to be chosen such that the impurity self-energy equals the local self-energy of the original model. The first constraint is that the ground state of the Anderson molecule has the same occupation numbers of the 2-site Hubbard model; at half filling, \(n_d=n_s=1\). This self-consistency condition is satisfied if \(\varepsilon _s=\mu =\varepsilon _d+U/2\). The ground state of the Anderson molecule is then

$$\begin{aligned} |G\rangle _A&=\frac{a_2(t,U/2)}{\sqrt{2}}\left( c_{d\uparrow }^\dagger c_{s\downarrow }^\dagger {-}c_{d\downarrow }^\dagger c_{s\uparrow }^\dagger \right) |0\rangle +\frac{a_1(t,U/2)}{\sqrt{2}}\left( c_{d\uparrow }^\dagger c_{d\downarrow }^\dagger {+}c_{s\uparrow }^\dagger c_{s\downarrow }^\dagger \right) |0\rangle , \end{aligned}$$
(65)

where the coefficient are the same functions defined in Eq. (49). The impurity Green function is then given by

$$\begin{aligned} G_{d,d}^\sigma (i\nu _n)=&\frac{w'_+}{i\nu _n-(E_0(2)-E_-(1)-\mu )} {+}\frac{w'_-}{i\nu _n-(E_0(2)-E_+(1)-\mu )} \nonumber \\ {+}&\quad \frac{w'_-}{i\nu _n-(E_+(3)-E_0(2)-\mu )} {+}\frac{w'_+}{i\nu _n-(E_-(3)-E_0(2)-\mu )} , \end{aligned}$$
(66)
Fig. 15
figure 15

Spectral function of the Hubbard dimer (H, left) and Anderson molecule (A, right) in the zero temperature limit for increasing U/t from 0 to 4 in equal steps. In the right panel, the solution of the Hubbard dimer in the DMFT approximation is also shown (dark blue lines). Figure rearranged from Ref. [107], ch. 7

where

$$\begin{aligned} E_0(2){-}E_\pm (1)-\mu&{=-}\bigg (E_\pm (3)-E_0(2)-\mu \bigg )\nonumber \\&{=-}\frac{1}{4}\bigg (2\varDelta (t,U/2)\pm \varDelta (t,U)\bigg ), \end{aligned}$$
(67)

and

$$\begin{aligned} w'_\pm&=\frac{1\pm w'(t,U)}{4},&\quad w'(t,U)&= \frac{1}{2} \frac{32t^2-{U^2}}{\varDelta (t,U)\varDelta (t,U/2)}. \end{aligned}$$
(68)

This can be recast into the simpler expression

$$\begin{aligned} G^{\sigma }_{d,d}(i\nu _n)=&\frac{1}{i\nu _n+\mu -\varepsilon _d- {{\mathcal {F}}}_{0}^{\sigma }(i\nu _n)-\varSigma ^\sigma _A(i\nu _n)} . \end{aligned}$$
(69)

One can verify that the impurity self-energy equals the local self-energy of the Hubbard dimer is

$$\begin{aligned} \varSigma ^\sigma _A(i\nu _n)=&\frac{U}{2}+\frac{U^2}{4} \frac{i\nu _n}{(i\nu _n)^2-{(3t)^2}}. \end{aligned}$$

The non-interacting hybridization function is given by

$$\begin{aligned} {{\mathcal {F}}}_{0}^{\sigma }(i\nu _n)=&\frac{t^2}{i\nu _n}. \end{aligned}$$

The impurity Green function thus satisfies the impurity Dyson equation

$$\begin{aligned} \varSigma _A^{\sigma }(i\nu _n)&=\frac{1}{G_{d,d}^{0\sigma }(i\nu _n)}-\frac{1}{{ G}^\sigma _{d,d}(i\nu _n)}. \end{aligned}$$
(70)

By comparing the expressions, one can see that the exact impurity Green function \(G_{d,d}^{\sigma }(i\nu _n)\) is not equal to \(G_{i,i}^{\sigma }(i\nu _n)\), the local Green function of the Hubbard dimer in the local self-energy approximation. The deviation, however, arises in very large part from the non-local self-energy, which does not appear in \(G_{d,d}^{\sigma }(i\nu _n)\). If we could neglect it (local self-energy approximation), there would still be a remaining discrepancy, but this is now merely a small correction. It arises from

$$\begin{aligned} \varDelta \mathcal{F}_l(i\nu _n)={{\mathcal {F}}}_l^{\sigma }(i\nu _n)-{{\mathcal {F}}}^{\sigma }_0(i\nu _n) = t^2p^2\bigg (-\frac{2}{i\nu _n} +\frac{1}{i\nu _n-\varepsilon _{a}}+ \frac{1}{i\nu _n+\varepsilon _{a}}\bigg ) \end{aligned}$$

where \( p^2={U^2}/{8\varepsilon _{a}^2}\) and \(\varepsilon _{a}=\sqrt{9t^2+U^2/4}.\) This can be seen in Fig. 14.

There are several important lessons to learn from this toy model. The first is that DMFT is a method for solving a many-body lattice problem.Footnote 9 It is exact in infinite dimension, and it is an excellent approximation for real correlated system. This is because in materials the co-ordination number is high, and hopping integrals do not decay fast with the distance. For a dimer, or small molecules with few correlated sites, it is instead not a good approximation, in particular at low energy. This is shown in Fig. 14. The second lesson is that, when the local-self-energy approximation is a good approximation, as it happens for correlated materials, the auxiliary Anderson model reproduces excellently the spectral function. This can be seen in Fig. 14, right panel. The third lesson is that this success is not limited to a specific range of parameters, as Fig. 15 clearly shows. This allows us to study phase transitions. Finally, leaving DMFT now aside, each of the two toy models discussed capture different aspects of the physics of an actual correlated lattice with infinite sites (thermodynamic limit). The Hubbard dimer captures the evolution of the Hubbard bands in a correlated lattice; the latter move far apart with increasing U. The Anderson molecule captures the evolution of the central “quasi-particle” peak in the spectral function a strongly-correlated lattice; the latter becomes progressively more narrow with increasing the Coulomb repulsion, and vanishes at the metal-insulator transition.

Fig. 16
figure 16

Spectral function of the Hubbard dimer (H) in the zero temperature limit for increasing V/t from 0 to \(V/t=U/t=3\) in equal steps. Left: total local spectral function. Right: k-resolved spectral function

11 Non-local Coulomb terms

The bare Coulomb interaction is strong and long ranged, as can be seen from the expression of the Hartree term. Nevertheless, in many systems, the effects of the long-range part of the Coulomb interaction are, to a large extent, well described already in the KS Hamiltonian, thanks to the effective KS potential. What is not captured are mostly the peculiar effects of the local Coulomb interaction, the Hubbard U in the example of the Hubbard dimer. The reason can be grasped reconsidering once more the Hubbard dimer, including this time, in addition to U, the inter-site Coulomb interaction, which we label with V. The Hamiltonian is

$$\begin{aligned} {\hat{H}}=&\varepsilon _d\sum _{i\sigma }{\hat{n}}_{i\sigma } -t \sum _\sigma \left( c^\dagger _{1\sigma } c_{2\sigma } +c^\dagger _{2\sigma } c_{1\sigma } \right) +U\sum _{i=1,2} {\hat{n}}_{i\uparrow } {\hat{n}}_{i\downarrow } +V\sum _{\sigma \sigma '}{\hat{n}}_{1\sigma } {\hat{n}}_{2\sigma '} \end{aligned}$$

For \(N_e=1\) and \(N_e=3\), this model has the same eigenvectors as in the \(V=0\) case. For \(N_e=2\), the form of the eigenvectors is also identical, provided that \(U\longrightarrow U-V\). Thus we can obtain the \(T\rightarrow 0\) Green function from Eq. (52) with the following transformations

$$\begin{aligned} E_\pm (N_e=1,U;V)&=E_\pm (N_e=1,U;0)\\ E_\pm (N_e=3,U;V)&=E_\pm (N_e=3,U;0)+2V\\ E_G(N_e=2,U;V)&=E_G(N_e=2,U{-}V;0)+V\\ \mu (U;V)&=\mu (U)+V \\ w(t,U;V)&=w(t,U{-}V;0) \end{aligned}$$

The result is shown in Fig. 16. Increasing V makes the spectra closer to a non-correlated system with only two central peak, but an effectively enhanced hopping integral. This example illustrates that a hypothetical system in which the Coulomb interaction strength is independent on the distance between sites (here \(U{=}V\)) is likely to be well described, at least in first approximation, via an effective non-interacting Hamiltonian. Obviously, in real materials, the effects of long-range Coulomb integrals can be more complex than in the two-site model just discussed, but the general considerations made here remain qualitatively true.Footnote 10 The smaller the difference between U and longer range terms, the less strongly-correlated the system is. V-type Coulomb terms can be taken into account in DMFT, e.g., via the E-DMFT extension [173, 174].

12 Linear response functions

Linear response functions can be obtained in DFT+DMFT in two steps. First the local susceptibility \(\chi (i\omega _m)\) is calculated solving the quantum-impurity problem. This can be done at the end of the self-consistency loop, as shown in Fig. 8. Calculating \(\chi (i\omega _m)\) is a much harder task than calculating single-particle Green functions, in particular because one needs to compute the full three-frequency tensor \(\chi (i\omega _m;i\nu _n,i\nu _{n'})\). Once this is done, the Bethe–Salpeter (BS) equation has to be solved

$$\begin{aligned} \chi (\mathbf{q};i\omega _m)&= {\chi }_0(\mathbf{q};i\omega _m) + {\chi }_0(\mathbf{q};i\omega _m) \varGamma (\mathbf{q};i\omega _m) \chi (\mathbf{q};i\omega _m) . \end{aligned}$$
(71)

The BS equation is known from standard many-body perturbation theory. In DMFT the “bubble” term \(\chi _0(\mathbf{q};i\omega _m)\) is obtained however using the DMFT Green function as the propagator. The next step consist in finding approximations for the vertex \(\varGamma (\mathbf{q};i\omega _m)\), which is unknown. In the infinite coordination limit, the vertex can however be replaced by a local quantity [102, 103, 175]

$$\begin{aligned} {\varGamma (\mathbf{q};i\omega _m)}\longrightarrow {\varGamma (i\omega _m)}. \end{aligned}$$

This leads to the modified BS equation

$$\begin{aligned} \chi (\mathbf{q};i\omega _m)&= {\chi }_0(\mathbf{q};i\omega _m) + {\chi }_0(\mathbf{q};i\omega _m) \varGamma (i\omega _m) \chi (\mathbf{q};i\omega _m) . \end{aligned}$$
(72)

By solving Eq. (72) we find, formally

$$\begin{aligned} \chi ^{-1}(\mathbf{q};i\omega _m) = \chi _0^{-1}(\mathbf{q};i\omega _m ) - \varGamma (i\omega _m). \end{aligned}$$
(73)

The local vertex can be obtained from Eq. (73), however replacing all quantities with their local counterparts

$$\begin{aligned} \varGamma (i\omega _m)= \chi _0^{-1}(i\omega _m ) -\chi ^{-1}(i\omega _m). \end{aligned}$$
(74)

Despite its apparent simplicity, solving the BS equation, even in the local vertex approximation, is a delicate task in practice. It involves the inversion of, in principle, infinite-dimension matrices. This can be seen from the analytic formula of the magnetic susceptibility of the Hubbard model in the atomic limit [53, 176, 177]

$$\begin{aligned} \chi _{zz}(i\omega _m;0)&= \sum _{\sigma \sigma ^\prime } \, \sigma \sigma ^\prime \,\chi ^{n,{n^\prime }}_{i\sigma \sigma ,i\sigma ^\prime \sigma ^\prime }(0)\nonumber \\&= M_{n'}\frac{d M_n}{dy} +M_n\frac{d M_{n'}}{dy} \nonumber \\&\quad -\beta n(y)\bigg [ \delta _{n,n'}+\delta _{n,-n'} \bigg ] \frac{d M_n}{dy} + \beta n(-y) M_n M_{n'}\nonumber \\&\quad - \frac{1}{y} \bigg \{ M_{n'} -\beta \bigg [n(y)\delta _{n,-n'}-{ n(-y)\delta _{n,n'} } \bigg ]\bigg \} M_n \; \end{aligned}$$
(75)

where

$$\begin{aligned} M_{n}= \frac{1}{i\nu _n-y}-\frac{1}{i\nu _n+y}, \end{aligned}$$
(76)

and \(n,n'\) are fermionic Matsubara frequencies. The expression (75) evidences the slow decay in the fermionic frequencies, in particular along the diagonals \({n=n'}\) and \({n=-n'}\). Only summing up all fermionic Matsubara frequencies one recovers the expression

$$\begin{aligned} \chi _{zz} (i\omega _m)=\delta _{m,0} \frac{n(-U/2)}{4k_BT} \end{aligned}$$
(77)

where n(x) is the Fermi distribution function. To complicate the matter, in actual materials the susceptibility tensor depends on the flavors (spin, site, orbitals) associated with the four (creation and annihilation) operators entering the expression of the generalized correlation function [176]. It is thus important to exploit symmetries and find compact representations, in order to reduce both truncation errors and computational time. More compact representation are based on polynomial expansions, for example Legendre polynomials [178].

In order to better explain how susceptibilities are calculated in DMFT we make further use of the toy model previously introduced, the two-site Hubbard model, and calculate the magnetic susceptibility in this case. For simplicity we perform in advance all sums over fermionic Matsubara frequencies, so that only the bosonic frequency \(\omega _m\) remains. The site-dependent components, defined in imaginary time as

$$\begin{aligned} \chi ^{ij}_{zz}(\tau )=&\langle {{\mathcal {T}}} S^i_z (\tau ) S^j_z(0) \rangle - \langle {{\mathcal {T}}} S^i_z (\tau )\rangle \langle {{\mathcal {T}}}S^j_z(0) \rangle , \end{aligned}$$
(78)

where \({{\mathcal {T}}}\) is the time-ordering operator and \(\tau \) the imaginary time. The total magnetic susceptibility is

$$\begin{aligned} \chi _{zz}(\mathbf{q}; i\omega _m)=\frac{1}{N_s} \sum _{ij} \int _0^\beta d\tau e^{i\omega _m\tau } \chi ^{ij}_{zz}(\tau ) e^{-i{\mathbf{q}\cdot (\mathbf{R}_i-\mathbf{R}_j}) } \end{aligned}$$
(79)

where \(\mathbf{R}_i\) are vectors identifying the \(N_s=2\) atomic sites; the only possible \(\mathbf{q}\) values are here 0 and \(\pi \), as we have already seen. The first step consists in calculating the impurity susceptibility for the (self-consistent) auxiliary model, the Anderson molecule. For this is necessary to consider not only the ground state, Eq. (65), but also the excited triplet states, with energy \(E_T\). They are given by

$$\begin{aligned} |1,0\rangle _A=&\frac{1}{\sqrt{2}}\bigg [c_{d\uparrow }^\dagger c_{s\downarrow }^\dagger - c_{d\downarrow }^\dagger c_{s\uparrow }^\dagger \bigg ] |0\rangle \end{aligned}$$
(80)
$$\begin{aligned} |1,1\rangle _A=&c_{d\uparrow }^\dagger c_{s\uparrow }^\dagger |0\rangle \end{aligned}$$
(81)
$$\begin{aligned} |1,-1\rangle _A=&c_{d\downarrow }^\dagger c_{s\downarrow }^\dagger |0\rangle . \end{aligned}$$
(82)

At sufficiently low temperature higher energy states can be neglected. We define the energy difference \(E_T=E_G+E_X\) with

$$\begin{aligned} E_X=-\frac{1}{4}\bigg [U-2\varDelta (t,{U}/{2})\bigg ]=\frac{U}{4}\bigg [-1+\sqrt{ 1+\frac{16\varGamma _{\mathrm{SE}}}{U} }\bigg ]. \end{aligned}$$
(83)

For sufficiently small \(\varGamma _{\mathrm{SE}}/U\) one then finds

$$\begin{aligned} \chi _{zz}(i\omega _m)&\sim \frac{1-e^{-\beta E_X}}{4Z} \bigg (\frac{1}{i\omega _m+E_X}- \frac{1}{i\omega _m-E_X}\bigg ) + \delta _{m,0} \beta \frac{e^{-\beta E_X}}{2Z} \end{aligned}$$

where \(Z\sim 1+3e^{-\beta E_X}\). At high temperature this becomes

$$\begin{aligned} \chi _{zz}(0)&\sim \frac{1}{4}\beta \delta _{m,0}, \end{aligned}$$
(84)

which is the atomic limit susceptibility. Next we have to calculate the DMFT bubble term. The latter can be written as

$$\begin{aligned} \chi _{zz}^0(\mathbf{q};i\omega _m)=-(g\mu _B)^2\frac{1}{4} \sum _{\sigma }\frac{1}{\beta } \sum _n \frac{1}{N_\mathbf{k}} \sum _\mathbf{k} G_{\mathbf{k} \sigma }(i\nu _n)G_{\mathbf{k}\sigma }(i\nu _n+i\omega _m), \end{aligned}$$
(85)

where \(G_{\mathbf{k} \sigma }(i\nu _n)\) is the DMFT Green function. Decomposing all the products of poles in sums of single poles, this can be recast into [53, 176, 177, 179]

$$\begin{aligned} \chi ^0_{zz}(i\omega _m;\mathbf{q})=(g\mu _B)^2\frac{1}{2} \frac{1}{N_\mathbf{k}}\sum _\mathbf{k} \bigg (-I^{++}_\mathbf{k,\mathbf{q}}-I^{--}_\mathbf{k,\mathbf{q}}+I^{+-}_\mathbf{k,\mathbf{q}}+I^{-+}_\mathbf{k, \mathbf{q}} \bigg ). \end{aligned}$$
(86)

The addends are defined as

$$\begin{aligned} I^{\alpha \gamma }_\mathbf{k,\mathbf{q}}=\frac{E^\alpha _\mathbf{k} E^\gamma _\mathbf{k + \mathbf{q}} }{(E^+_\mathbf{k}-E^-_\mathbf{k})(E^+_\mathbf{k+\mathbf{q}}-E^-_\mathbf{k+\mathbf{q}})} \frac{n(E^\alpha _\mathbf{k})-n(E^\gamma _\mathbf{k +\mathbf{q}})}{E^\alpha _\mathbf{k}-E^\gamma _\mathbf{k +\mathbf{q}}+i\omega _m}. \end{aligned}$$
(87)

The energies are

$$\begin{aligned} E_\mathbf{k}^\pm =\frac{1}{2}\varepsilon _\mathbf{k}\pm \frac{1}{2}\sqrt{\varepsilon _\mathbf{k}^2+U^2}, \end{aligned}$$
(88)

with \(\varepsilon _\mathbf{0}=-t\) and \(\varepsilon _{{\varvec{\pi }/a}}=+t\). Thus

$$\begin{aligned} \chi ^0_{zz}(i\omega _m;{ 0})&\approx -\frac{(g\mu _B)^2}{4} \frac{1}{N_\mathbf{k}}\sum _\mathbf{k}\frac{U^2}{\sqrt{\varepsilon _\mathbf{k}^2+{U^2}}} \frac{1}{(i\omega _m)^2-({\varepsilon _\mathbf{k}^2+U^2}) }\\ \chi ^0_{zz}(i\omega _m; {\varvec{\pi }})&\approx -\frac{(g\mu _B)^2}{4} \frac{1}{N_\mathbf{k}}\sum _\mathbf{k}\frac{\left( \varepsilon _\mathbf{k} +\sqrt{\varepsilon _\mathbf{k}^2+U^2}\right) ^{3}}{\varepsilon _\mathbf{k}^2+U^2} \frac{1}{(i\omega _m)^2 -\left( \varepsilon _\mathbf{k}+ \sqrt{\varepsilon _\mathbf{k}^2+U^2}\right) ^2}. \end{aligned}$$

For small \(\varGamma _{\mathrm{SE}}/U\) the static local term is then

$$\begin{aligned} \chi ^0_{zz}(0)&\sim \frac{(g\mu _B)^2}{4U} \bigg (1-\frac{\varGamma _\mathrm{SE}}{4U}\bigg ). \end{aligned}$$
(89)

In the last step we solve the simplified Bethe–Salpeter equation

$$\begin{aligned} \frac{1}{\chi _{zz}(i\omega _m;\mathbf{q})}=\frac{1}{\chi ^0_{zz}(i\omega _m;\mathbf{q})}-\frac{1}{\chi ^0_{zz}(i\omega _m)}+\frac{1}{\chi _{zz}(i\omega _m)}. \end{aligned}$$
(90)

In the high-temperature regime this leads for the static susceptibility to the Curie–Weiss expression

$$\begin{aligned} {\chi _{zz}(0;\mathbf{q})}=\frac{(g\mu _B)^2}{4}\frac{1}{k_BT-\varGamma _\mathrm{SE}(\mathbf{q})/4}, \end{aligned}$$
(91)

where \(\varGamma _{\mathrm{SE}}({0})=-\varGamma _{\mathrm{SE}}/2\) and \(\varGamma _{\mathrm{SE}}({\pi })=+\varGamma _{\mathrm{SE}}/2\). Few observations are in place. First, the DMFT bubble term is dominated by the energy scale U, and therefore it has poles for spin excitations at incorrect positions. Instead, the local susceptibility from the quantum impurity problem and the final \(\mathbf{q}\)-dependent susceptibility are rightly controlled by the scale \(\varGamma _\mathrm{SE}\), which is the energy scale of spin excitations. This also holds when the results are extended to the lattice [53, 176, 177]. Second, the bubble term is weakly temperature dependent. The Curie–Weiss temperature behavior arises entirely from the local susceptibility of the quantum-impurity problem. The \(\chi ^0_{zz}(i\omega _m;\mathbf{q})\) term is thus a poor approximation of the actual magnetic response function.

Fig. 17
figure 17

Spin waves in the \(t-t'\) Hubbard model at half filling, calculated with DMFT. The white lines is the standard linear spin-wave dispersion obtained from the associated Heisenberg model (large t/U limit of the Hubbard model). Reprinted figure with permission from Musshoff et al. [180]. Copyright (2021) by the American Physical Society

The example illustrates why one cannot obtain the correct response function using the Kohn–Sham eigenenergies calculated with DFT+U, even in the magnetic phase, i.e., when in the self-energiy are similar in DFT+DMFT and DFT+U. This corresponds to taking into account only the “bubble” term \(\chi ^0_{zz}(i\omega _m;\mathbf{q})\). Instead, DMFT provides a good description. It can be used to calculate the susceptibility above and below the critical temperature, as well as bosonic excitations. This may be seen, e.g., in Figs. 17 and 18, taken from Ref. [180].

Fig. 18
figure 18

The transition to an antiferromagnetic phase in the DMFT solution of the \(t-t'\) Hubbard model at half filling. It shows that DMFT approximately yields the mean-field solution of the associated Heisenberg model. Reprinted figure with permission from Musshoff et al. [180]. Copyright (2021) by the American Physical Society

13 Methods for non-local correlations

In the last decades, the LDA+DMFT approach has become the state of the art for strongly-correlated systems. Its success has been striking in describing trends, explaining experiments and highlighting new phenomena in correlated materials. It is also obvious where its main limitation lie. We have already met them when discussing the toy model, the two-site Hubbard model. Non-local correlations are not captured by a local self-energy. Still, here one should make a distinction. Long-range ordered phases themselves can be well described in single-site DMFT, simply using supercells. A different story is the description of non-local effects in phase transitions. In such a case, DMFT gives basically the static mean-field critical exponents, i.e., the Curie–Weiss law. This is shown for a representative case in Fig. 18. Straightforward extensions that go beyond DMFT are cellular dynamical mean-field theory and the dynamical cluster approximation [122]. With cluster approaches it is possible to recover the exact critical behavior, but the cluster size required might be very large, depending on the problem analyzed and the type of clusters used. This is shown in Refs. [122, 181] for the two-dimensional Hubbard model. The example chosen is an extreme case since the exact \(T_N\) is zero due to the Mermin–Wagner theorem, while at the same time the co-ordination number is relatively small. Still, even in three-dimensional materials a large cluster size might be necessary, depending on the case, for accurately reproducing the details of phase transitions. Alternatives to cluster approaches are diagrammatic methods, such as the GW+DMFT approach [182], the dynamical-vertex approximation [183, 184], or the dual-fermion technique [185,186,187].

14 Conclusion

In this short article, I introduced some of the modern methods for studying strong correlations in materials. I emphasized successes and shortcomings of the current state-of-the art approaches, and briefly discussed developments that could play an important role in the future. The article is a personal view of a field in evolution, of which snapshots have been captured. More detailed descriptions of specific aspects, solvers and applications can be found, e.g., in Refs. [41, 107, 108, 110,111,112,113, 122].

In this conclusion, I would like to underline few points. The impressive progress achieved in the last decades in the description of correlation effects in materials relies on schemes which, instead of solving the full many-body Hamiltonian, solve low-energy materials-specific models. The latter are built to capture both the essence of the mechanisms at play and the key specific characteristics which distinguish a system from the others. Furthermore, these models are solvable, if not exactly within very good approximations. This successful strategy combines the best of many-body physics and ab-initio approaches—overcoming the historical dichotomy between the two. Anderson, had once coined the ironical expression “the great solid-state physics dream machine” when referring to DFT-based first principles methods [24]. The main criticism was that, in a many-body system, it is extremely hard to predict truly novel emergent properties—and indeed it rarely happens. Methods like DFT+DMFT and its extensions, although not fully ab-initio, give us the impression that eventually reductionism could win.

Even in the best case scenario, this will always be true only in a limited sense, however. Traditional ab-initio methods show us the nature of the restriction. Using Kohn–Sham eigenvalues for interpreting experiments implies a stark assumption on the form of the model and on the type of phenomena described. State-of-the-art methods for correlated systems allow us to study a vast variety of new systems and phenomena, and have been instrumental to clarify several open problems. They are becoming progressively more powerful, slowly but steadily approaching predictive power. Still, they suffer from the same core problem—they do not exhaust all the possibilities in a many-electron system. If something truly new is found, we have to return to the basics—trying first to guess and solve the basic model, and only then make it materials-specific. In fact, despite all progress, even old fundamental problems remain unsettled, not least the nature of superconductivity in cuprates.

Finally, we have seen that together with theoretical progress and algorithmic advances, the advent of massively parallel supercomputers was crucial to study realistic models in practice. And yet, the exponential growth of the Hilbert space poses limit on what can be done with a classical computer. In the future, quantum computers hold the promise to bring a new revolution for the description of strongly correlated systems. If this promise is realized, we will be able to tackle always more complex problems. Still, numerical calculations alone will never be sufficient. A simulation that exactly reproduce the world will not give us, alone, more information then experiments themselves. The final frontier will always be the understanding, disentangling the crucial aspects from all the rest. In that respect, the many-body problem and its emergent phenomena should keep us busy for a long time to come.