Beyond Quantum Cluster Theories: Multiscale Approaches for Strongly Correlated Systems

Strongly correlated systems include some of the most technologically promising materials of our time. To harness their significant promise, understanding the fundamental mechanisms responsible for their intriguing properties is essential. Wolf et al. (2001); Morosan et al. (2012); Dagotto and Tokura (2008); Tokura and Nagaosa (2000); Dagotto (2005) This understanding remains a challenge for the condensed matter community despite several decades of intense effort. For instance, although the discovery of high temperature superconductors dates back to 1987, Bednorz and Müller (1986) the underlying superconducting mechanism remains the subject of intense research activity. Following their discovery, the Hubbard model was postulated to contain the ingredients necessary to explain the properties of high temperature superconductors and their low-energy excitations. Zhang and Rice (1988) But despite its simplicity, an exact solution of the Hubbard model beyond one dimension remains elusive. Lieb and Wu (2003, 1968) Therefore, numerical methods have played a crucial role. These methods are however constrained by the minus sign problem for Quantum Monte Carlo (QMC), or by the exponential scaling of the Hilbert space for exact diagonalization, to relatively small system sizes. Embedding schemes have emerged as an important avenue to treat the problem in the thermodynamic limit. These schemes map the lattice problem onto an impurity, for the case of dynamical mean field theory, or onto a cluster for the case of its cluster extensions, dynamical cluster approximation or cellular dynamical mean field theory, embedded into a mean field. Maier et al. (2005); Georges et al. (1996) Embedding approaches allow the exact treatment of short length scales at the level of the cluster or the impurity, and the treatment of longer length scales at the mean field level. Multiscale many body methods follow this logic to its natural next step by incorporating between the previous two length scales, an intermediate one at which correlations are treated with diagrammatic perturbation theory. Slezak et al. (2009)

In general, the difficulty in understanding correlated systems lies in the fact that there are no simple theories to explain both the weak interaction limit of the metallic state and the strong interaction limit of the Mott insulating phase.Mott (1968, 1982) The most successful theory of interacting fermions is the Fermi liquid theory. Landau (1956, 1957, 1959) The basic underlying assumption is that the interaction can be turned on adiabatically from the non-interacting free fermions limit. The consequence is that the quantum numbers of the non-interacting fermions remain unchanged. Electrons can be treated as quasi-particles in a rather stable state with a lifetime that becomes very long for those states near the Fermi level.

The Fermi liquid theory is a very efficient description of interacting fermions in a metallic phase Landau (1956, 1957, 1959). It is applicable to almost all metallic phases, except for special circumstances such as the notable exception of one-dimensional systems. The theory owes its simplicity to being an effective renormalized single particle theory. Once the system is beyond the simple single particle description, there is no universal prescription to handle the competition or cooperation between the different degrees of freedom and the interplay between the kinetic and the potential energies. Precisely for this reason, numerical methods are often inevitable for practical calculations.

Widely used mean field methods factorize the interaction terms in the Hamiltonian to reduce the problem to an effective single particle theory in a static potential. The mean field, Hartree-Fock, approximation often provides reasonable results Hartree (1928); Fock (1930a, b); Slater (1930), but its shortcomings are also obvious, in particular for intermediate interaction strengths where quantum fluctuations are large. The Hartree-Fock approximation quenches the quantum or temporal fluctuations completely. This may be a reasonable assumption if the interactions are overpowered by the kinetic energy terms. However, for many physical realizations of strongly correlated systems, perhaps the most well known one being the cuprate superconductors, the interaction is of the same order of magnitude as the bandwidth. Naively factorizing the interaction term to suppress all the quantum fluctuations is questionable at best. Indeed, there is currently no simple mean field theory that can explain most features of the cuprate superconductors. Understanding the metallic phases beyond Fermi liquid theory is key for understanding broken symmetry phases, such as d-wave superconducting pairing in the cuprates. While one can construct a phase with no explicit broken symmetry and use the mean field method to understand the effective theory, this always involves fractionalized particles and strong constraints such as those of gauge theories. Lee et al. (2006)

Beyond mean field theory, there exists a plethora of techniques based on weak coupling expansion. They are typically based on low order perturbative methods, such as second order perturbation theory, or on selecting a certain class of diagrams and summing them up to infinite order. A typical example is the random phase approximation (RPA), which selects the class of ladder diagrams and sums them up to infinite order. Bohm and Pines (1951); Pines and Bohm (1952); Bohm and Pines (1953); Gell-Mann and Brueckner (1957) A more sophisticated approach is to sum a large class of diagrams in an iterative way. For example, parquet diagrams are generated when second order diagrams are inserted iteratively into the interaction vertex. This generates a class of diagrams that can only be separated into two disconnected pieces by cutting at least two fermion lines. Landau et al. (1954a, b, c) The advantage of the parquet approach compared to second order perturbation theory is in the ability to sum up a large variety of diagrams including those at infinite order. This, in principle, allows the instability towards a broken symmetry to be captured Bickers (1998); Bickers and White (1991); Bickers and Scalapino (1993). Its main advantage over random phase approximation is in its unbiased sum of diagrams in different scattering channels to enable the study of the competition among different broken symmetries.

Instead of using a diagrammatic expansion approach, the Dynamical Mean Field Theory (DMFT) maps the strongly correlated lattice onto an impurity site embedded in a self-consistently determined effective medium. The interest in the high spatial dimension limit of strongly correlated models in the late 80’s and early 90’s led to the understanding that in this limit, strongly correlated models with local interactions can be greatly simplified. This is due to the fact that an expansion in terms of the hopping amplitude in infinite dimension leads to the vanishing of all diagrams except the local ones, and, for a translationally invariant system, the model loses all spatial dependence. This simplification led to the dynamical mean field theory.

DMFT remains the subject of active research efforts, particularly because there is no universal quantum impurity solver. Various methods have been proposed over the past few decades. These include semi-analytical methods based on perturbation theory or modified mean field theories. The more well known methods include the iterative perturbation theory and the local moment approximation. Numerical approaches include various kinds of Quantum Monte Carlo and exact diagonalization methods. Recently, density matrix renormalization group and matrix product state methods have also been explored. Quantum computing algorithms for solving the quantum impurity problem have been proposed recently. Bauer et al. (2016); Rungger et al. (2020); Keen et al. (2019); McClean et al. (2016) After all, solving even a single impurity is a non-trivial problem, as the mean field hybridization function is not given by a simple form that can allow an analytical solution.

It is worth noting that the DMFT can be viewed as a formal generalization of the coherent potential approximation (CPA) proposed by Soven in the 60’s. Soven (1967); Shiba (1971) The CPA has since been extensively used for studies of random disorder models with negligible interactions, in particular for random alloys. Various extensions of CPA have been proposed over the years. The earliest one that goes beyond the single site approximation is the molecular CPA. Hass et al. (1984); Lempert et al. (1987) It embeds the cluster into an effective medium that possesses the same structure as the cluster itself. It thus generates an effective medium different from that of the CPA. The obvious deficiency of the method is the breaking of translational invariance.

A different scheme of embedding cluster methods is the Dynamical Cluster Approximation (DCA) Hettler et al. (2000, 1998). By construction, the method naturally preserves the translational invariance of the original model by directly working with both the cluster and the effective medium in momentum space. The method has been extensively employed on the Hubbard model Gutzwiller (1963); Kanamori (1963); Hubbard (1963), Anderson model Anderson (1958), periodic Anderson model Anderson (1961), and Falicov-Kimball model Falicov and Kimball (1969). The cluster extension is not just a quantitative improvement on the DMFT. It is necessary to produce important features that are absent in the DMFT results. Perhaps the most important one is the DCA’s ability to capture nonlocal correlations such as that of d-wave pairing, which is obviously absent for approximations that do not consider spatial dependence explicitlyTsuei et al. (1994). The method has also been considered in the context of multiple scattering theory where it is re-branded as non-linear coherent potential approximationRowlands et al. (2003).

The difficulty of solving a cluster impurity (or embedded cluster) problem scales exponentially. Roughly speaking, quantum Monte Carlo based methods scale exponentially with the number of impurity sites (cluster size), with inverse temperature, or with the interaction strength Hirsch and Fye (1986); Gull et al. (2008, 2011); Rubtsov et al. (2005). The exception is strong-coupling expansion based Monte Carlo methods, but this is usually limited to a rather small number of impurity sitesWerner et al. (2006); Werner and Millis (2006). Another class of impurity/cluster solvers is based on diagonalization of the effective finite size Hamiltonian. For these Hamiltonian-based solvers, such as exact diagonalization Liebsch et al. (2008); Sénéchal (2008); Liebsch and Ishida (2011); Gunnarsson and Schönhammer (1983); Hong and Kee (1995); Koch et al. (2008), the Hilbert space grows exponentially with the cluster size and thus both computing memory and time requirements grow at the same rate. This is also true to a large extent for another Hamiltonian-based approach, the numerical renormalization group method Bulla et al. (2008); Krishna-murthy et al. (1980a, b); Wilson (1975). In general, for practical calculations, the maximum number of impurity sites is rather modest ($\sim 10$ sites).

Over the last couple of decades various novel methods have been proposed. These include the density matrix renormalization group Fernández et al. (2014); Zhu et al. (2017), the related matrix product wavefunction Wolf et al. (2014a) and, even more recently, different forms of machine learning approaches Arsenault et al. (2014); Walker et al. (2020). These more recent methods may have potential for certain applications. For instance, they may be more efficient for calculating real time Green functions in nonequilibrium problems Wolf et al. (2014b). Approaches based on machine learning could also be more efficient in solving a large set of impurity problems, and this may be useful for applications on random systems that require averaging over random disorder realizations Terletska et al. (2018). However, none of these novel impurity solvers are suitable for the calculation of the vertex function which is essential for most methods that are built on an expansion on top of the DMFT solution. Additionally, Monte Carlo sampling of the partition function provides more flexibility for controlling the error as the impurity cluster size is increased. Also, although it has been proven that the single impurity problem does not exhibit a minus sign problem, the absence of minus sign in the Monte Carlo sampling can not be assumed for a generic impurity problem Yoo et al. (2005).

Following the logic of embedding a small system into a mean-field host, one can anticipate better accuracy in the result if an intermediate length scale is inserted between the previous two. Since short length scales are appropriately treated with exact solvers and the long length scale by a mean field, this intermediate length scale can be treated reliably with diagrammatic methods. This is the essence of multiscale many body methods for strongly correlated systems as formulated in the early 2000’s Slezak et al. (2009) and the subject of continued efforts since then Rohringer et al. (2018a).

This review focuses exclusively on the methods for studying systems in equilibrium. It is noteworthy to point that effort has been devoted to the generalization of the multiscale many body (MSMB) approaches to non-equilibrium problems Jung et al. (2012); Zhou and Guo (2019); Chen et al. (2019). We refer the reader to the original article for details Jung et al. (2012), as we do not discuss these nonequilibrium approaches in the present paper.

The rest of the review is structured as follows. In section II, we summarize the DMFT method and its connection to the Anderson impurity problem by approaching it from its ”cavity method” formulation. In section III, we discuss two cluster extensions of DMFT, dynamical cluster approximation (DCA) and cellular dynamical mean field theory (CDMFT). We proceed in section IV with a discussion of the extended dynamical mean field theory (EDMFT) that extends DMFT to the treatment of nonlocal interaction. Section V is focused on the $1/d$ expansion, a systematic expansion of DMFT with respect to the hopping amplitude. In section VI, we describe the original formulation of the multiscale many body method. The parquet formalism, which encompasses various commonly used diagrammatic approximations, and which is essential for the diagrammatic treatment of intermediate length scales, is described in section VII. Following this, we briefly describe different implementations to incorporate nonlocal corrections into the DMFT/DCA starting with the dynamical vertex approximation in section VIII, and then the dual fermion method in section IX. In section X, we present the dual boson extension for corrections to EDMFT. In section XI, we discuss efforts to incorporate nonlocal correction into the GW approximation (an approach that obtains the self-energy from the single particle Green function (G) and the screened Coulomb interaction (W)) in the form of the “triply irreducible local expansion” (TRILEX). In section XII, we discuss functional renormalization group and its usage for nonlocal corrections to DMFT. In section XIII, an important computational challenge, the numerical representation of the vertex functions in memory, is discussed. In section XIV, we discuss important physical constraints on the methods. In section XV, we summarize results obtained on different models with various implementations of the multiscale many-body approach before ending with our conclusions.

Following the discovery of high temperature superconductors, it was argued that the Hubbard model (1) captures their low energy properties. Zhang and Rice (1988) This deceptively simple model describes itinerant electrons that can hop between nearest-neighbor sites $\langle i,j\rangle$ on a lattice with a hopping integral $t_{ij}$, and are subject to a Coulomb interaction $U$ when a site is doubly occupied. The model is schematically depicted in Fig. 1 and defined as:

$\displaystyle\hat{H}=-\sum_{<i,j>,\sigma}t_{ij}c^{\dagger}_{i,\sigma}c_{j,\sigma}+\sum_{i}U\hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow},$

(1)

where $c_{i,\sigma}$ $(c_{i,\sigma}^{\dagger})$ is the destruction (creation) operator that destroys (creates) an electron with spin $\sigma$ at site $i$. $\hat{n}_{i,\sigma}=c^{\dagger}_{i,\sigma}c_{i,\sigma}$ is the number of particles of spin $\sigma$ at site $i$.

The DMFT method in itself is a generalization of the usual mean field theory. But unlike the usual mean field theory, say for the Ising model for example, the mean field here is not an order parameter. Instead it is a function of time or frequency. Thus, the approach captures all temporal fluctuations. Focusing on the translationally invariant paramagnetic phase, there is no explicit order parameter. One may also construct an explicit order parameter to represent the broken symmetry, but this is not done routinely and is not the focus of this discussion.

Within the DMFT solution of a translationally invariant system, the spatial fluctuations are completely suppressed as in a traditional mean field theory. As we will see, a major way to improve the approximation is to systematically incorporate the corrections due to the spatial dependence of the model into the DMFT solution. Also, since the self-consistency condition is valid only at the single particle level, there is no guarantee that high order Green functions, such as the different susceptibilities, are matched between the lattice and impurity models.

When the DMFT was developed in the early 90’s Metzner and Vollhardt (1989a, b); Müller-Hartmann (1989a, b); Georges and Kotliar (1992); Jarrell (1992); Georges et al. (1996); Janiś (1991); Metzner and Vollhardt (1990), solving even a single impurity problem was rather challenging. With the formulation of new numerical algorithms and advances in computing power, a single impurity problem can, in general, be numerically solved quite efficiently. While there is still not a completely satisfactory method that is accurate for a wide range of parameters, particularly for solving interacting problems with random disorder which requires a large ensemble of different impurity realizations for disorder averaging Terletska et al. (2018); Dobrosavljević et al. (2003); Ekuma et al. (2015), the single impurity problem can mostly be handled at the present time.

As a side note, it is worth mentioning that earlier work on DMFT can be traced back to the study of the transverse-field Sherrington-Kirkpatrick model.Sherrington and Kirkpatrick (1975) In that model the spatial fluctuations are completely suppressed as the model consists of spin couplings in the fully connected network with the transverse external magnetic field, thus it can be mapped to a single site problem. Bray and Moore (1980) The idea of similarly handling fermion problems only appeared in the late 80’s in studies of fermionic systems in the infinite dimensional limit by Vollhardt and Metzer. Metzner and Vollhardt (1989a, b) A series of papers by Müller-Hartmann were also influential in the development of DMFT and later generalizations to the dynamical cluster approximation. Müller-Hartmann (1989a, b) In 1991, the physics of the Hubbard model via DMFT was discussed by Georges and Kotliar. Georges and Kotliar (1992) The first numerically “exact” DMFT solution of the Hubbard model was presented by Jarrell. Jarrell (1992)

The DMFT formalism can be explained quite transparently from a path integral formulation. Consider the action of the Hubbard model on a lattice,

			$\displaystyle S=-\sum_{\mathbf{r}_{i},\mathbf{r}_{j},\sigma}\int_{0}^{\beta}\int_{0}^{\beta}d\tau_{i}d\tau_{j}\psi_{\sigma}^{*}(\mathbf{r}_{i},\tau_{i}){G}_{0}^{-1}(\mathbf{r}_{i},\tau_{i},\mathbf{r}_{j},\tau_{j})\psi_{\sigma}(\mathbf{r}_{j},\tau_{j})$		(2)
		$\displaystyle+$	$\displaystyle U\sum_{\mathbf{r}_{i}}\int^{\beta}_{0}d\tau_{i}\psi_{\uparrow}^{*}(\mathbf{r}_{i},\tau_{i})\psi_{\uparrow}(\mathbf{r}_{i},\tau_{i})\psi^{*}_{\downarrow}(\mathbf{r}_{i},\tau_{i})\psi_{\downarrow}(\mathbf{r}_{i},\tau_{i}),$

where $G_{0}$ is the bare Green function, $\psi_{\sigma_{i}}^{*}(\mathbf{r}_{i},\tau_{i})$ and $\psi_{\sigma_{i}}(\mathbf{r}_{i},\tau_{i})$ are the Grassmann fields for spin $\sigma_{i}$ at location $\mathbf{r}_{i}$ and imaginary time $\tau_{i}$.

The first part of the action contains the kinetic energy as characterized by the bare Green function $G_{0}$. It is simply obtained from the bare dispersion of the considered model. The second term includes the interaction characterized by the parameter $U$, which we assume to be local. Any interaction beyond the local Hubbard term will involve further approximations in the context of the dynamical mean field theory.

The exact Green function of the above action can be completely characterized by the self-energy $\Sigma$. If we write the self-energy in the frequency-momentum space, the relation between the bare Green function and the exact Green function $G$ is given by the Dyson equation,

$\displaystyle G(\mathbf{k},\omega)=\frac{1}{G_{0}^{-1}(\mathbf{k},\omega)-\Sigma(\mathbf{k},\omega)}.$

(3)

In the simplified case of a translationally invariant system, the idea of the dynamical mean field theory is to relate the full lattice problem with spatial dependence to a single site problem. To this end, DMFT reframes (2) into an effective action for a single site with a bare Green’s function $\mathcal{G}$:

			$\displaystyle S_{eff}=-\int d\omega\sum_{\sigma}\psi_{\sigma}^{*}(\omega)\mathcal{G}^{-1}(\omega)\psi_{\sigma}(\omega)$
		$\displaystyle+$	$\displaystyle U\int_{\omega_{1}+\omega_{3}=\omega_{2}+\omega_{4}}d\omega_{1}d\omega_{2}d\omega_{3}d\omega_{4}\psi_{\uparrow}^{*}(\omega_{1})\psi_{\uparrow}(\omega_{2})\psi^{*}_{\downarrow}(\omega_{3})\psi_{\downarrow}(\omega_{4}).$

Consider the Anderson impurity model characterizing an impurity coupled to a band of conducting electrons and given by the Hamiltonian:

	$\displaystyle H_{AIM}$	$\displaystyle=$	$\displaystyle Un_{\uparrow}n_{\downarrow}-\mu\sum_{\sigma}c_{\sigma}^{\dagger}c_{\sigma}$		(5)
		$\displaystyle+$	$\displaystyle\sum_{j,\sigma}V_{j}\left(f_{j\sigma}^{\dagger}c_{\sigma}+c_{\sigma}^{\dagger}f_{j\sigma}\right)+\sum_{j\sigma}\epsilon_{j}f_{j\sigma}^{\dagger}f_{j\sigma}.$

Where $c_{\sigma}^{\dagger}$, $c_{\sigma}$ are creation and destruction operators for the impurity electrons while $f_{j\sigma}^{\dagger}$, $f_{j\sigma}$ are those of the conduction electrons. Its action is:

			$\displaystyle S_{AIM}=-\int d\omega\sum_{\sigma}\psi_{\sigma}^{*}(\omega)\mathcal{G}^{-1}_{AIM}(\omega)\psi_{\sigma}(\omega)$
		$\displaystyle+$	$\displaystyle U\int_{\omega_{1}+\omega_{3}=\omega_{2}+\omega_{4}}d\omega_{1}d\omega_{2}d\omega_{3}d\omega_{4}\psi_{\uparrow}^{*}(\omega_{1})\psi_{\uparrow}(\omega_{2})\psi^{*}_{\downarrow}(\omega_{3})\psi_{\downarrow}(\omega_{4}),$

where $\mathcal{G}_{AIM}$ is the single impurity Anderson model non-interacting Green’s function defined by:

$$\mathcal{G}^{-1}_{AIM}(i\omega_{n})=i\omega_{n}+\mu-\Delta(i\omega_{n})$$

(7)

with

$$\Delta(i\omega_{n})=\int_{-\infty}^{+\infty}d\omega\frac{1}{i\omega_{n}-\omega}\sum_{j\sigma}V_{j}^{2}\delta(\omega-\epsilon_{j}).$$

(8)

The AIM action (II) is equivalent to (II) with $\mathcal{G}$ playing the role of the non-interacting AIM Green’s function. The construction can be justified via the concept of the ”cavity method” in the infinite dimension limit whereby all degrees of freedom are integrated out except for the site labelled by the index $0$. In this limit of $d\to\infty$ for a hypercubic lattice, the hopping amplitude is rescaled as $t_{ij}=t/\sqrt{2d}$ so that the kinetic energy and the interaction energy remain of the same order. The effective action, $S_{eff}$, in this process is defined by:

	$\displaystyle\frac{1}{Z_{eff}}exp(-S_{eff}[\psi_{0,\sigma}^{*}\psi_{0,\sigma}])=$
	$\displaystyle\frac{1}{Z}\int\prod_{i\neq 0,\sigma}exp(-S[\{\psi_{i,\sigma}^{*},\psi_{i,\sigma}\}]).$		(9)

Where $Z$ and $Z_{eff}$ are the partition functions associated with $S$ and $S_{eff}$ respectively. The effective action can be written as:

			$\displaystyle S_{eff}[\psi_{0,\sigma}^{*}\psi_{0,\sigma}]=S_{0}$		(10)
		$\displaystyle+$	$\displaystyle\sum_{n=1}^{N}\sum_{i_{1},j_{1},\cdot\cdot\cdot i_{n},j_{n}}\eta^{*}_{0,i_{1}}\eta^{*}_{0,i_{2}}\cdot\cdot\cdot\eta^{*}_{0,i_{n}}\eta_{j_{1},0}\eta_{j_{2},0}\cdot\cdot\cdot\eta_{j_{n},0}$
			$\displaystyle\prod_{i,j=1,\cdot\cdot\cdot,n}\int d\tau_{i}d\tau_{j}G^{(0)}(i_{1},\tau_{i_{1}},i_{2},\tau_{i_{2}},\cdot\cdot\cdot,i_{n},\tau_{i_{n}};$
			$\displaystyle j_{1},\tau_{j_{1}},j_{2},\tau_{j_{2}},\cdot\cdot\cdot,j_{n},\tau_{j_{n}}).$

Where $S_{0}$ is the local action at site ”0”:

	$\displaystyle S_{0}$	$\displaystyle=$	$\displaystyle-\int_{0}^{\beta}\int_{0}^{\beta}d\tau d\tau^{{}^{\prime}}\sum_{\sigma}\psi_{\sigma}^{*}(\tau){G}_{0}^{-1}(\tau,\tau^{{}^{\prime}})\psi_{\sigma}(\tau^{{}^{\prime}})$		(11)
		$\displaystyle+$	$\displaystyle U\int^{\beta}_{0}d\tau\psi_{\uparrow}^{*}(\tau)\psi_{\uparrow}(\tau)\psi^{*}_{\downarrow}(\tau)\psi_{\downarrow}(\tau).$

$G^{(0)}$ is the Green’s function connecting the cavity to the impurity. $\eta_{i,0}=t_{i,0}\psi_{0,\sigma}$, with $t_{i,0}$ the hopping from site $i$ to $0$.

Only terms of order $n=2$ survive the expansion $(\ref{eq:actionEff3})$ in the $d\to\infty$ limit. Leading to:

	$\displaystyle S_{eff}[\psi_{0,\sigma}^{*}\psi_{0,\sigma}]=S_{0}$
	$\displaystyle+\sum_{i_{1},j_{1}=1}^{N}t_{0,i_{1}}t_{j_{1},0}\int d\tau d\tau^{{}^{\prime}}\psi_{0,\sigma}^{*}\psi_{0,\sigma}G(i_{1},\tau_{i_{1}},j_{1},\tau_{j_{1}}).$		(12)

Rewriting this in frequency space, gives the AIM action (II) with $\mathcal{G}_{AIM}$ replaced by $\mathcal{G}$ such that:

$\displaystyle\mathcal{G}^{-1}(i\omega_{n})=i\omega_{n}+\mu-\sum_{i,j}t_{0,i}t_{j,0}G(i\omega_{n}).$

(13)

This relation connects the impurity Green function to the lattice Green function. For the Bethe lattice, Bethe (1935) $\mathcal{G}^{-1}(i\omega_{n})=i\omega_{n}+\mu-t^{2}G(i\omega_{n})$. For a general lattice, the connection between the lattice Hubbard model and the single impurity model is established by setting the self-energy,

$\displaystyle\Sigma_{lattice}(\mathbf{k},\omega)=\Sigma_{impurity}(\omega).$

(14)

Since the self-energy of the original lattice Hubbard model has spatial dependence while that of the single impurity Anderson model does not, to construct the lattice Green function one has to rely on the coarse-graining process that assumes the self-energy of the lattice model to be the same in the entire Brillouin zone.

	$\displaystyle G_{lattice}(\mathbf{k},\omega)=\frac{1}{G_{lattice,0}^{-1}(\mathbf{k},\omega)-\Sigma(\mathbf{k},\omega)}$		(15)
	$\displaystyle\approx\frac{1}{G_{lattice,0}^{-1}(\mathbf{k},\omega)-\Sigma_{lattice}(\omega)}.$

The missing link between the Hubbard model and the Anderson model is to determine the effective bare Green function of the Anderson model. For the $D$-dimension case, this is given by:

	$\displaystyle G_{impurity}(\omega)=\sum_{\mathbf{k}}\frac{1}{(2\pi)^{D}}G_{lattice}(\mathbf{k},\omega)$		(16)
	$\displaystyle=\sum_{\mathbf{k}}\frac{1}{(2\pi)^{D}}\frac{1}{G_{lattice,0}^{-1}(\mathbf{k},\omega)-\Sigma_{lattice}(\omega)}.$

The summation over the momentum can thus be replaced by an integral over the bare density of states to simplify the calculation. The bare density of states for the hypercubic lattice at the limit of infinite spatial dimensions can be exactly calculated. Müller-Hartmann (1989a, b) We have gathered all the ingredients for the DMFT approximation and the algorithm can be summarized as in Fig. 2.

A natural and direct avenue to generalize the DMFT is to incorporate nonlocal correlations by including more than a single impurity site, i.e by formulating the theory around a cluster of multiple sites in a self-consistently determined host. This type of cluster DMFT remains an important method for the study of strongly correlated systems, as it allows, by increasing the cluster size, a systematic correction unlike perturbative expansion methods. Moreover, one can envision that a perturbative expansion on top of the cluster method would produce an even better result, since the bare effective Hamiltonian or action for the perturbative expansion, which corresponds to the DCA or CDMFT solution on the smaller system, is presumably more accurate and already includes a substantial amount of nonlocal correlations.

For the classical spin model, the first attempt of a multiple site mean field theory was the so-called Bethe-Peierls-Weiss approximation using a cluster of $(z+1)$ sites, with one site at the center surrounded by $z$ sites on the shell. Weiss (1948); Bethe (1935); Peierls (1936); Kikuchi (1951) The interaction between the center spin and its $z$ nearest neighbor spins is treated explicitly while the interaction between the remaining $z$ spins with the other spins outside their own cluster is treated by a mean field.

Another approach by Oguchi is a more direct generalisation of the mean field method. Oguchi (1955) A cluster of $N_{c}$ spins is considered. The interaction among these $N_{c}$ spins is treated explicitly, and the interaction between the spins at the edge of the cluster and spins outside of the cluster is treated by a mean field. Unlike the Bethe-Peierls-Weiss method where pairwise interactions are treated explicitly only for pairs involving the central spin, pairwise interactions among all spins of the $N_{c}$ spin cluster are treated explicitly in the Oguchi approach.

A cluster of impurities in real space is considered instead of a single one for the cellular dynamical mean field theoryKotliar et al. (2001); Lichtenstein and Katsnelson (2000). There is a technical problem with using such an approach for the paramagnetic solution of the Hubbard model as the cluster naturally breaks translational invariance. A procedure for restoring the symmetry is needed for a periodic solution. Biroli et al. (2004)

A further approach for a cluster generalization of DMFT is based on the idea of coarse-graining that is central to the dynamical cluster approximation (DCA). Aryanpour et al. (2002); Fotso et al. (2012); Hettler et al. (2000); Jarrell et al. (2001) A cluster of impurities is used, but after the local cluster is solved the lattice quantities are averaged over different patches of the first Brillouin zone for a coarse-grained quantity. The advantage is that the method is manifestly translationally invariant. This allows a perturbative expansion to be implemented on top of the DCA solution more naturally. As we will discuss, almost all perturbative expansion methods become simpler and less cumbersome to implement when formulated in momentum space.

One can consider that the DMFT impurity bare Green function is the coarse-grained lattice Green function. Since the self-energy which is obtained by solving the impurity model does not have spatial dependence, the coarse-graining procedure for DMFT is done over the entire first Brillouin zone with one single impurity site. Effectively, from a diagrammatic point of view, DMFT neglects momentum conservation at the internal vertices. The dynamical cluster approximation systematically restores momentum conservation at these internal vertices. To this end, it divides the Brillouin zone into $N_{c}$ cells with each cell (of linear size $\Delta k$) represented by a cluster momentum K in the center of the cell. The DCA then requires that momentum conservation in the internal vertices be respected for momentum transfers between cells (momentum transfers larger than $\Delta k$), but neglected for momentum transfers within a cell (less than $\Delta k$). In this way momentum conservation is fully recovered in the limit of $N_{c}\to\infty$, while the DMFT result is obtained for $N_{c}=1$.

The DCA coarse-graining process is illustrated by Fig.(3) for $N_{c}=8$. In DCA, the self-energy is no longer momentum-independent. Rather, we have for a lattice momentum $\mathbf{k}$, $\Sigma(\mathbf{k})=\Sigma(\mathbf{K})$. Where $\mathbf{K}$ is the momentum at the center of the cell containing $\mathbf{k}$. The impurity, or in this case the cluster, Green’s function is here related to the lattice Green’s function by:

	$\displaystyle G_{impurity}(\mathbf{K},\omega)\approx\sum_{\mathbf{k}}\frac{1}{(2\pi)^{D}}G_{lattice}(\mathbf{k},\omega)$		(17)
	$\displaystyle=\sum_{\mathbf{k}}\frac{1}{(2\pi)^{D}}\frac{1}{G_{0,lattice}^{-1}(\mathbf{k},\omega)-\Sigma_{lattice}(\mathbf{K},\omega)},$

where the summation over $\mathbf{k}$ is restricted to the patch corresponding to the impurity/cluster site with momentum $\mathbf{K}$. The impurity solver is now a cluster solver for the momentum-dependent self-energy $\Sigma(\mathbf{K})$.

The algorithm for DCA is summarized in Fig. (4).

The earliest attempt to include the effect of nonlocal interactions was motivated by the competition between the local interaction and the RKKY interaction in heavy fermion materials. Si and Smith (1996) The basic idea was to include the density-density and/or spin-spin interaction term and scale the interaction strength with respect to the spatial dimensionality so that its fluctuations are non-zero in the high dimension limit. Smith and Si (2000a); Si and Smith (1996); Smith, J. L. and Si, Q. (1999); Chitra and Kotliar (2000); Smith and Si (2000b); Sun and Kotliar (2002) This procedure leads to an impurity embedded in a self-consistent fermionic bath and, at the same time, a self-consistent bosonic bath due to the nonlocal interactions.

In DMFT, the hopping is considered as a function of spatial dimension $d$, $t_{ij}\sim 1/\sqrt{d}^{|i-j|}$ in the infinite dimension limit. The EDMFT includes the nonlocal interaction in a similar manner: $V_{ij}\sim 1/\sqrt{d}^{|i-j|}$. Smith and Si (2000a); Si and Smith (1996); Smith, J. L. and Si, Q. (1999) Here, we only consider the density-density interaction. A generic two-body interaction can include spin-spin, correlated hopping, and pair hopping. For example, if the nonlocal interaction is density-density interaction, such as in the extended Hubbard model, the effective action for the impurity problem acquires an extra term including a retarded density-density coupling mediated by the charge susceptibility. Si and Smith (1996); Chitra and Kotliar (2000)

	$\displaystyle S_{0}$	$\displaystyle=$	$\displaystyle\int d\tau d\tau^{{}^{\prime}}\sum_{\sigma}\psi_{\sigma}^{*}(\tau){G}_{0}^{-1}(\tau,\tau^{{}^{\prime}})\psi_{\sigma}(\tau^{{}^{\prime}})$		(18)
		$\displaystyle-$	$\displaystyle U\int d\tau n_{\uparrow}(\tau)n_{\downarrow}(\tau)$
		$\displaystyle+$	$\displaystyle\int d\tau d\tau^{{}^{\prime}}\phi^{*}(\tau){D}_{0}^{-1}(\tau,\tau^{{}^{\prime}})\phi(\tau^{{}^{\prime}})$
		$\displaystyle+$	$\displaystyle\int d\tau\phi(\tau)n(\tau).$

The effective action is equivalent to that of the single impurity model where $\mathcal{G}$ and $\mathcal{D}$ act respectively as effective fields for the fermionic bath and for the bosonic bath in the Anderson impurity model. They can be obtained as usual via the Dyson equations. For the Bethe lattice, Chitra and Kotliar (2000)

$\displaystyle\mathcal{G}^{-1}(i\omega_{n})=i\omega_{n}+\mu-\sum_{i,j}t_{0,i}t_{j,0}G(i\omega_{n}),$

(19)

and

$\displaystyle\mathcal{D}^{-1}(i\omega_{n})=\sum_{i,j}V_{0,i}^{-1}V_{j,0}^{-1}D(i\omega_{n}).$

(20)

Since the retarded coupling can be understood in term of a bosonic field coupled to the impurity charge density, EDMFT remedies the limitation of DMFT where the interaction is strictly restricted to the local on-site Hubbard interaction. We note that the decoupling of the interaction by Hubbard-Stratonovich fields is not unique: one can decouple the coupling in different ways. Ayral et al. (2013) The decoupling presented above only decouples the density-density interaction, the $V$ term in the extended Hubbard model.

The algorithm is thus similar to that of the DMFT, except that at each iteration both the mean fields for the fermionic and for the bosonic baths have to be updated. The algorithm can be summarized as in Fig. 5

One of the earliest attempts to consider the effect of finite dimension and move away from the infinite dimension limit is the systematic expansion with respect to the hopping amplitude. This approach has been studied in other contexts, especially in the closely related CPA for disordered systems. Gebhard (1990); Vlaming and Vollhardt (1992); Uhrig (1996); Janiš and Vollhardt (2001); van Dongen (1994) Since the hopping term carries the factor of $1/\sqrt{d}$, an expansion with respect to the hopping regains the dimensional dependence at the expense of including multiple sites. Work by Schiller and Ingersent studied the case of two impurities for the Falicov-Kimball model. Schiller and Ingersent (1995) The approach has also been applied to the half-filled Hubbard model. Zaránd et al. (2000)

The $1/d$ corrections can be recovered by considering a multi-impurity problem. For a single impurity ($n=1$) or a two-impurity $(n=2)$ problem,

	$\displaystyle S^{(n)}$		$\displaystyle=\sum_{\sigma}\sum_{\alpha,\beta=1}^{n}\int d\tau\int d\tau^{\prime}{\psi}^{*}_{\alpha\sigma}(\tau)({\cal G}^{(n)}_{\alpha,\beta})^{-1}(\tau-\tau^{\prime})\psi_{\beta\sigma}(\tau^{\prime})$		(21)
			$\displaystyle+\sum_{\alpha=1}^{n}U\int d\tau n_{\alpha\uparrow}(\tau)n_{\alpha\downarrow}(\tau)\;,$

where $\alpha$ and $\beta$ label the sites for $n=2$. Schiller and Ingersent (1995); Zaránd et al. (2000) The mean fields ${\cal G}^{(1)}$ and ${\cal G}^{(2)}$ are chosen in such a way that the impurity Green functions $G^{(1)}$ and $G^{(2)}$ coincide with the full on-site and nearest neighbor lattice propagators, $G^{\rm lattice}_{00}$ and $G^{\rm lattice}_{01}$:

$$G^{(1)}=G^{(2)}_{11}=G^{\rm lattice}_{00}\;\phantom{n},\phantom{nnn}G^{(2)}_{12}=G^{\rm lattice}_{01}\;.$$

(22)

For skeleton diagrams (diagrams without self-energy dressing or vertex correction) of order ${\cal O}(1/d)$, the impurity self energies $\Sigma^{(1)}$ and $\Sigma^{(2)}_{\alpha\beta}$ and the diagonal and off-diagonal lattice self energies, $\Sigma^{\rm lattice}_{0}$ and $\Sigma^{\rm lattice}_{1}$ are related by Schiller and Ingersent (1995); Zaránd et al. (2000)

	$\displaystyle\Sigma^{\rm lattice}_{0}$	$\displaystyle=$	$\displaystyle\Sigma^{(1)}+2d(\Sigma^{(2)}_{11}-\Sigma^{(1)})\;,$		(23)
	$\displaystyle\Sigma^{\rm lattice}_{1}$	$\displaystyle=$	$\displaystyle\Sigma^{(2)}_{12}\;.$		(24)

The lattice Green function can be obtained from the lattice self-energy:

$$G_{lm}^{\rm lattice}(i\omega)={1\over 1+\sqrt{d}\;\Sigma^{\rm lattice}_{1}}G^{0}_{lm}\left({i\omega-\Sigma^{\rm lattice}_{0}(i\omega)\over 1+\sqrt{d}\;\Sigma^{\rm lattice}_{1}(i\omega)}\right)\;,$$

(25)

where $G^{0}_{lm}(z)$ and $G_{lm}^{\rm lattice}(z)$ denote the unperturbed and dressed lattice propagators between sites $l$ and $m$, respectively. Schiller and Ingersent (1995); Zaránd et al. (2000)

Among the shortcomings of the $1/d$ expansion are its limitations in the description of long range fluctuations. Perhaps more importantly, the truncation at finite order may lead to non-analytic properties of some dynamical quantities. Pruschke et al. (2001) The method is an example of the nested cluster scheme and is also related to the recently proposed self-energy embedding theory. Biroli et al. (2004); Kananenka et al. (2015); Zgid and Gull (2017)

Given the challenge of solving a large cluster by numerical methods and the lack of accurate analytical approaches, perturbation theory is a possible route for the exploration of physics beyond local theories. Initial efforts aimed to use perturbative methods as solvers for the cluster impurity problem. Notably, the fluctuation exchange method was used as an impurity cluster solver. Aryanpour et al. (2003) The fluctuation-exchange or FLEX is one of the simplest methods to incorporate correlations among different channels. It is, therefore, a conceptually appealing approach for systems in which particle-hole or particle-particle vertex fluctuations are not small. If correlations of the two particle fluctuations are ignored, one obtains the second order perturbation theory for which the irreducible vertex is replaced by the bare vertex. FLEX allows vertex contributions from different channels to be correlated, thus the name of fluctuation exchange approximation.

The general scheme of the multiscale many body method as envisioned by Jarrell and collaborators is to construct a theory which allows for the treatment of different length scales by different approaches. Jarrell et al. (2007); Slezak et al. (2009) The short length scale is addressed by some highly accurate numerical approach, such as quantum Monte Carlo. The long length scale is treated at a mean field level, while the intermediate length scale is treated by some form of perturbative technique.

An early proposition was to supplement the Quantum Monte Carlo calculation on small cluster sizes with FLEX on larger clusters via a self-energy self-consistency scheme between the two methods. J.P.Hague et al. (2004) It was generalized in 2006 to address short and long length scales within the DCA formalism while the intermediate length scale would be treated by the parquet formalism. Slezak et al. (2009) In this approach, the connection between the two methods is established by the fully irreducible vertex from the DCA calculation on the cluster. The promise of the approach stems from the better scaling of the parquet formalism compared to that of QMC as a cluster solver.

The general construction of the multiscale many body method, as summarized in Fig. 7, is obviously rather general. There is in fact plenty of freedom in the choice of a solver for the intermediate length scale. Indeed, this area of research has been the subject of significant activity over the past decade or so. Nevertheless, most of the developed methods are based on some simplification of the original proposal, either in picking a certain subclass of diagrams from the parquet formalism or in simplifying the solution by numerical techniques.

A major route for perturbative expansions around the DMFT solution is based on the parquet diagrams. This approach encapsulates many of the approximations that have emerged in this field, including the dynamical vertex approximation. For this reason, we review the parquet method in rather self-contained detail in this section.

Standard diagrammatic perturbative expansions attempt to describe all the scattering processes as single or two-particle Feynman diagrams. In the single-particle formulation the self-energy describes the many-body processes that renormalize the motion of a particle in the interacting background of all the other particles. In the two-particle context, one is able to probe the interactions between particles using the so-called vertex functions, which are matrices describing two particle scattering processes. For example, the reducible (full) two-particle vertex $F^{ph}(12;34)$ describes the scattering amplitude of a particle-hole pair from its initial state $\left|3,4\right>$ into the final state $\left|1,2\right>$. Here, $i=1,2,3,4$ represents a set of indices which combines the momentum $\mathbf{k}_{i}$, the Matsubara frequency $i\omega_{n_{i}}$ and, if needed, the spin $\sigma_{i}$ and band index $m_{i}$. Since the total momentum and energy of the vertex are conserved, it is convenient to adopt the notation $F^{ph}(2-4)_{1,3}$ for the numerical implementation on the single band Hubbard model. Other representations are also possible. Karrasch et al. (2008)

In general, depending on how particles or holes are involved in the scattering processes, one can define three different two-particle scattering channels. These are the particle-hole (p-h) horizontal channel, the p-h vertical channel and the particle-particle (p-p) channel. The parquet formalism is, in essence, a method for summing up diagrams that characterize scattering processes at the two-particle level. From another perspective, the diagrams are generated by inserting the one loop, second order, diagrams repeatedly into itself. Without channel mixing, this is equivalent to the random phase approximation. With the mixing of three channels, this becomes the parquet formalism.

The vertices can be categorized by extending the notion of diagram reducibility to the two-particle level as illustrated by Fig. 8. At the one particle level, a diagram is said to be reducible if it can be split in two disconnected parts by breaking a single Green’s function line. A two-particle diagram will be said to be irreducible if it can not be separated in two disconnected parts by breaking two Green function lines in the same channel. Bickers (1998) It will be said to be fully irreducible if it can not be separated in two disconnected parts by breaking two-Green’s function lines in any channel. In the single particle formalism, the Green function is related to the self-energy containing all single-particle irreducible diagrams by the Dyson equation. A connection is made between the single-particle and the two-particle diagrams by the Schwinger-Dyson equation that connects the self-energy to the (full/reducible) vertex $F$ containing all allowed diagrams in a given channel. The subset of all two-level diagrams in the full vertex that are irreducible in the same channel is known as the irreducible vertex $\Gamma$. The subset of irreducible vertices that are irreducible in any channel is called the fully irreducible vertex $\Lambda$.

It is worth noting that the above idea for the decomposition is not unique. One can devise other possible decompositions. For example, a recent attempt is to decompose diagrams in terms of the fermion-boson vertex. We will not explore this direction in detail in this review. The above decomposition is the most natural one in the sense that the method can be easily understood in terms of an iterative process. The higher order diagrams are all generated by iteratively replacing the vertex function at a lower order approximation.

Furthermore, we are mostly interested in models that preserve the $SU(2)$ spin rotation symmetry. Since this symmetry is always obeyed for the two-dimensional calculations at non-zero temperature, it is convenient to preserve this symmetry. This is accomplished by decomposing the vertices in the so-called spin-diagonalized representation. In this representation, the spin degrees of freedom decompose the particle-hole channel into the density and the magnetic channels, and the particle-particle channel into the spin singlet and the spin triplet channels which we denote as $d$-channel, $m$-channel, $s$-channel, and $t$-channel respectively. Bickers (1998) They are defined for the irreducible vertex as follows,

	$\displaystyle\Gamma_{d}$	$\displaystyle=$	$\displaystyle\Gamma^{PH}_{\uparrow\uparrow;\uparrow\uparrow}+\Gamma^{PH}_{\uparrow\uparrow;\downarrow\downarrow},$		(26)
	$\displaystyle\Gamma_{m}$	$\displaystyle=$	$\displaystyle\Gamma^{PH}_{\uparrow\uparrow;\uparrow\uparrow}-\Gamma^{PH}_{\uparrow\uparrow;\downarrow\downarrow},$		(27)
	$\displaystyle\Gamma_{s}$	$\displaystyle=$	$\displaystyle\Gamma^{PP}_{\uparrow\downarrow;\uparrow\downarrow}-\Gamma^{PP}_{\uparrow\downarrow;\downarrow\uparrow},$		(28)
	$\displaystyle\Gamma_{t}$	$\displaystyle=$	$\displaystyle\Gamma^{PP}_{\uparrow\downarrow;\uparrow\downarrow}+\Gamma^{PP}_{\uparrow\downarrow;\downarrow\uparrow},$		(29)

and similarly for $F$ and $\Lambda$.

The formalism is completed by equations that connect the different types of vertices. The full vertex is related to the irreducible vertex by the Bethe-Salpeter equation and the irreducible vertex is in turn related to the fully irreducible vertex by the parquet equation. We reproduce the full set of equations for the parquet formulation in the spin diagonalized representation in the following.

The Schwinger-Dyson equation that connects the vertex to the self-energy is

	$\displaystyle\Sigma(P)$	$\displaystyle=$	$\displaystyle-\frac{UT^{2}}{4N}\sum_{P^{\prime},Q}\{G(P^{\prime})G(P^{\prime}+Q)G(P-Q)(F_{d}(Q)_{P-Q,P^{\prime}}-F_{m}(Q)_{P-Q,P^{\prime}}){}$		(30)
			$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;+G(-P^{\prime})G(P^{\prime}+Q)G(-P+Q)(F_{s}(Q)_{P-Q,P^{\prime}}+F_{t}(Q)_{P-Q,P^{\prime}})\},$

where $G$ is the single-particle Green function, which itself can be calculated from the self-energy using the Dyson equation,

$\displaystyle G^{-1}(P)$

$\displaystyle=$

$\displaystyle G_{0}^{-1}(P)\;-\;\Sigma(P),$

(31)

where $G_{0}$ is the bare Green function. Here, the indices $P$, $P^{\prime}$ and $Q$ combine momentum ${\bf k}$ and Matsubara frequency $i\omega_{n}$, i.e. $P=({\bf k},i\omega_{n})$.

The reducible and the irreducible vertices in a given channel are related by the Bethe-Salpeter equation,

$$F_{r}(Q)_{P,P^{\prime}}=\Gamma_{r}(Q)_{P,P^{\prime}}+\Phi_{r}(Q)_{P,P^{\prime}},$$

(32)

$$F_{r^{\prime}}(Q)_{P,P^{\prime}}=\Gamma_{r^{\prime}}(Q)_{P,P^{\prime}}+\Psi_{r^{\prime}}(Q)_{P,P^{\prime}},$$

(33)

where $r=d\;\mbox{or}\;m$ for the density and magnetic channels and $r^{\prime}=s\;\mbox{or}\;t$ for the spin singlet and spin triplet channels. The vertex ladders are defined as

	$\displaystyle\Phi_{r}(Q)_{P,P^{\prime}}\equiv\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;$		(34)
	$\displaystyle\sum_{P^{\prime\prime}}F_{r}(Q)_{P,P^{\prime\prime}}\chi_{0}^{ph}(Q)_{P^{\prime\prime}}\Gamma_{r}(Q)_{P^{\prime\prime},P^{\prime}},$
	$\displaystyle\Psi_{r^{\prime}}(Q)_{P,P^{\prime}}\equiv\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;$		(35)
	$\displaystyle\sum_{P^{\prime\prime}}F_{r^{\prime}}(Q)_{P,P^{\prime\prime}}\chi_{0}^{pp}(Q)_{P^{\prime\prime}}\Gamma_{r^{\prime}}(Q)_{P^{\prime\prime},P^{\prime}},$

where $\chi_{0}$, the bare susceptibility, is the product of two single-particle Green functions.

The parquet equations in the spin diagonalized representation are

	$\displaystyle\Gamma_{d}(Q)_{P{P^{\prime}}}$	$\displaystyle=$	$\displaystyle\Lambda_{d}(Q)_{P{P^{\prime}}}-{\frac{1}{2}}\Phi_{d}({P^{\prime}}-P)_{P,P+Q}-{\frac{3}{2}}\Phi_{m}({P^{\prime}}-P)_{P,P+Q}$
			$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+\;{\frac{1}{2}}\Psi_{s}(P+{P^{\prime}}+Q)_{-P-Q,-P}+{\frac{3}{2}}\Psi_{t}(P+{P^{\prime}}+Q)_{-P-Q,-P},$

	$\displaystyle\Gamma_{m}(Q)_{P{P^{\prime}}}$	$\displaystyle=$	$\displaystyle\Lambda_{m}(Q)_{P{P^{\prime}}}-{\frac{1}{2}}\Phi_{d}({P^{\prime}}-P)_{P,P+Q}+{\frac{1}{2}}\Phi_{m}({P^{\prime}}-P)_{P,P+Q}$
			$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;-\;{\frac{1}{2}}\Psi_{s}(P+{P^{\prime}}+Q)_{-P-Q,-P}+{\frac{1}{2}}\Psi_{t}(P+{P^{\prime}}+Q)_{-P-Q,-P},$

	$\displaystyle\Gamma_{s}(Q)_{P{P^{\prime}}}$	$\displaystyle=$	$\displaystyle\Lambda_{s}(Q)_{P{P^{\prime}}}+{\frac{1}{2}}\Phi_{d}({P^{\prime}}-P)_{-{P^{\prime}},P+Q}-{\frac{3}{2}}\Phi_{m}({P^{\prime}}-P)_{-{P^{\prime}},P+Q}$
			$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+\;{\frac{1}{2}}\Phi_{d}(P+{P^{\prime}}+Q)_{-{P^{\prime}},-P}-{\frac{3}{2}}\Phi_{m}(P+{P^{\prime}}+Q)_{-{P^{\prime}},-P},$

	$\displaystyle\Gamma_{t}(Q)_{P{P^{\prime}}}$	$\displaystyle=$	$\displaystyle\Lambda_{t}(Q)_{P{P^{\prime}}}+{\frac{1}{2}}\Phi_{d}({P^{\prime}}-P)_{-{P^{\prime}},P+Q}+{\frac{1}{2}}\Phi_{m}({P^{\prime}}-P)_{-{P^{\prime}},P+Q}$
			$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;-\;{\frac{1}{2}}\Phi_{d}(P+{P^{\prime}}+Q)_{-{P^{\prime}},-P}-{\frac{1}{2}}\Phi_{m}(P+{P^{\prime}}+Q)_{-{P^{\prime}},-P}.$

It is important to note that if we substitute the irreducible vertices $\Gamma$ (Eqs. VII, VII, VII, and VII) into the Bethe-Salpeter equation (Eqs. 32 and 33) the crossing symmetries (symmetry relations of the vertex that are a consequence of the Pauli exclusion principle for identical fermionic particles) in the full vertex $F$ is automatically satisfied regardless of the numerical values of the vertex ladders $\Phi$ and $\Psi$, assuming the fully irreducible vertices, $\Lambda$, obey the crossing symmetries. We write all the full vertices explicitly in the following, using only the vertex ladders, $\Phi$, $\Psi$, and the fully irreducible vertices, $\Lambda$.

	$\displaystyle F_{d}(Q)_{P,P^{\prime}}=\Lambda_{d}(Q)_{P{P^{\prime}}}-{\frac{1}{2}}\Phi_{d}({P^{\prime}}-P)_{P,P+Q}-{\frac{3}{2}}\Phi_{m}({P^{\prime}}-P)_{P,P+Q}$		(40)
	$\displaystyle+{\frac{1}{2}}\Psi_{s}(P+{P^{\prime}}+Q)_{-P-Q,-P}+{\frac{3}{2}}\Psi_{t}(P+{P^{\prime}}+Q)_{-P-Q,-P}+\Phi_{d}(Q)_{P,P^{\prime}};$

	$\displaystyle F_{m}(Q)_{P,P^{\prime}}$	$\displaystyle=$	$\displaystyle\Lambda_{m}(Q)_{P{P^{\prime}}}-{\frac{1}{2}}\Phi_{d}({P^{\prime}}-P)_{P,P+Q}+{\frac{1}{2}}\Phi_{m}({P^{\prime}}-P)_{P,P+Q}$
		$\displaystyle-$	$\displaystyle{\frac{1}{2}}\Psi_{s}(P+{P^{\prime}}+Q)_{-P-Q,-P}+{\frac{1}{2}}\Psi_{t}(P+{P^{\prime}}+Q)_{-P-Q,-P}+\Phi_{m}(Q)_{P,P^{\prime}};$

	$\displaystyle F_{s}(Q)_{P,P^{\prime}}$	$\displaystyle=$	$\displaystyle\Lambda_{s}(Q)_{P{P^{\prime}}}+{\frac{1}{2}}\Phi_{d}({P^{\prime}}-P)_{-{P^{\prime}},P+Q}-{\frac{3}{2}}\Phi_{m}({P^{\prime}}-P)_{-{P^{\prime}},P+Q}$
		$\displaystyle+$	$\displaystyle{\frac{1}{2}}\Phi_{d}(P+{P^{\prime}}+Q)_{-{P^{\prime}},-P}-{\frac{3}{2}}\Phi_{m}(P+{P^{\prime}}+Q)_{-{P^{\prime}},-P}+\Psi_{s}(Q)_{P,P^{\prime}};$

	$\displaystyle F_{t}(Q)_{P,P^{\prime}}$	$\displaystyle=$	$\displaystyle\Lambda_{t}(Q)_{P{P^{\prime}}}+{\frac{1}{2}}\Phi_{d}({P^{\prime}}-P)_{-{P^{\prime}},P+Q}+{\frac{1}{2}}\Phi_{m}({P^{\prime}}-P)_{-{P^{\prime}},P+Q}$
		$\displaystyle-$	$\displaystyle{\frac{1}{2}}\Phi_{d}(P+{P^{\prime}}+Q)_{-{P^{\prime}},-P}-{\frac{1}{2}}\Phi_{m}(P+{P^{\prime}}+Q)_{-{P^{\prime}},-P}+\Psi_{t}(Q)_{P,P^{\prime}}.$

The full parquet formalism encompasses a variety of approximations that are widely used in condensed matter physics and materials science. The hierarchy of these different approximations is neatly summarized in Fig. 9.

• •

  Hartree-Fock and second order perturbation theory: At the highest level, we might make the approximation on the two-particle Green function (analogous to the conventional single particle Hartree-Fock (HF) perturbation) such that the four-point correlation function can be factorized as a product of two two-point correlation functions. It is equivalent to ignoring the contribution from the full vertex functions. From the Schwinger-Dyson equation, the self-energy has two contributions a Hartree-Fock term and a second order perturbation theory term.

• •

  Self-consistent second order perturbation theory: Substituting the bare vertex for the full vertex in the Schwinger-Dyson equation and solving for the self-energy self-consistently results in the self-consistent second order perturbation theory.

• •

  Random-Phase Approximation (RPA): The irreducible vertex in the longitudinal charge channel is approximated by the bare Coulomb interaction. And then the Bethe-Salpeter equation is used to sum up over all the ring-type diagrams. Bohm and Pines (1951); Pines and Bohm (1952, 1952)

• •

  T-Matrix Approximation (TMA): Similar to RPA, the irreducible vertex in the transverse particle-hole channel or particle-particle channel, instead of longitudinal particle-hole channel, is approximated by the bare Coulomb interaction. And then the Bethe-Salpeter equation is used to sum all the ladder-type (instead of ring-type in RPA) diagrams. Bethe and Goldstone (1957)

• •

  Fluctuation Exchange Approximation (FLEX): A combination of RPA and T-matrix approximation, such that the fluctuations in different channels are treated equally. Bickers et al. (1989)

The parquet formalism dates back to the 50’s, but as can be seen from the analytical form of the governing equations, a general analytical solution is rather hard to track. Pomeranchuk et al. (1956); De Dominicis and Martin (1964) Over the years various approximations to simplify the equations have been proposed to solve a set of problems ranging from the Anderson impurity to nuclear structure. Roulet et al. (1969); Yakovenko (1993); Brazovskiĭ (1972); Kleinert and Schlegel (1995); Chen and Bickers (1992); Janiš and Augustinský (2007, 2008); Augustinský and Janiš (2011); Janiš (2001, 2009); Bickers (1998); Bickers and White (1991); Hess et al. (1996); Luo and Bickers (1993); Janiš (1999); Kusunose (2010); Jackson et al. (1985); Jackson and Smith (1987); Pfitzner and Wölfle (1987); Weiner (1970, 1971); Yeo and Moore (1996a, b, 2001); Yeo et al. (2006); Shishanin and Ziyatdinov (2003); Aref’eva and Zubarev (1996); Bergli and Hjorth-Jensen (2011); Janiš et al. (2019) In addition, a numerical solution is computationally demanding. This difficulty is in general common to theories that involve the two-particle vertex functions. The computational difficulty arises mostly from the memory requirements to store the vertex functions as they are four-legged objects unlike the single particle quantities that only have two legs. If momentum and energy conservation are implemented, the single particle quantities scale linearly with the size of the space-time grid, on the other hand the vertex function scales as the third power of the size of the space-time grid. This challenge is not insurmountable and may be overcome with appropriate parallelization. In fact, the advent of petascale computing enabled the first solution for the two-dimensional problem.Yang et al. (2009); Tam et al. (2013)

The algorithm for the numerical solution of the parquet formalism is summarized in Fig. 10.

Perhaps the most challenging problem from the numerical point of view is the difference in the nature of the single particle and the two-particle functions. The single particle quantities in the Matsubara frequency space do not diverge. On the other hand, the two particle vertex has strong divergences when the metallic phase is unstable at the momentum and energy corresponding to an instability. For instance, in the Hubbard model, the particle-hole vertex at the verge of the antiferromagnetic instability has a strong divergence for momentum transfer $(\pi,\pi)$. For this reason, the vertex is represented by numbers spanning a wide range of values over many orders of magnitude. This clearly is a recipe for possible numerical instabilities in the iterative solution for the vertex functions at low temperature and on the verge of long range order phase transitions. As we will discuss in the rest of this review, various methods have been proposed to improve the stability of the numerical solutions. These include simplifying the equations, or abandoning the self-consistent approach for the vertex functions.

The original dynamical vertex approximation (D$\Gamma$A) is a simplification of the multiscale many-body method with a pragmatic mindset that limits the calculation to the local fully irreducible vertex using perturbation theory. Held et al. (2008); Toschi et al. (2007) As discussed in section VII above, a full parquet solution with space-time resolution is very challenging. A natural scheme to sidestep the difficulty is, in the spirit of DMFT, to consider the vertex to be only time or frequency dependent. The local fully irreducible vertex from a DMFT calculation can then be used as input for the parquet formalism. The method subsequently follows the procedure discussed above for the parquet method, therefore, this aspect will not be repeated in this section.

Since this approach is relatively transparent and numerically practical with modest computational costs, Held (2014) besides the two-dimensional Hubbard model, it has been applied to several problems including the three-dimensional Hubbard model, the attractive Hubbard model, Del Re et al. (2019) nanoscopic quantum junction systems, Valli et al. (2010) and it has also been recently combined with ab initio calculations. Galler et al. (2017) A more elaborate calculation based on parquet method has been performed for the one-dimensional Hubbard model Schäfer et al. (2017).

There are a few variations of the dynamical vertex approximation. The full parquet formalism can be solved self-consistently with the fully irreducible vertex from the DMFT solution used as input. One can also consider self-consistency at the level of both the parquet equations and the DMFT equations: this involves finding the fully irreducible vertex from the DMFT solution, then using this to solve the parquet equations. The parquet formalism provides both the full vertex and the dressed Green function. The dressed Green function can in turn be treated as input for the DMFT equation to obtain self-consistency in both the DMFT loop and the parquet equations loop.

Initial applications on the half-filled Hubbard model motivated a further simplification by decoupling the particle-hole channels from the particle-particle channels. This simplification can be justified by the physics of the systems of interest. For example, in the Hubbard model near half-filling, the density wave is driven by the nesting in the particle-hole channels and one can argue for choosing the particle-hole ladder summations. On the other hand, when the system is driven by s-wave pairing such as that in the attractive Hubbard model, one can keep only the particle-particle channel. There is no systematic universal argument on which channel should be dominant. For most interesting regimes, such as the d-wave pairing in the Hubbard model, presumably all channels could contribute and one may have no choice but to try to tackle the full set of equations of the parquet formalism.

The algorithm is summarized in Fig. 11.

Another path to the multiscale treatment of correlations in a fermionic system is that of the dual fermions approach that was built on previous analogous methods for bosonic systems. This approach systematically incorporates nonlocal correlations into the DMFT solution. The method is distinguished from others in that it maps a strongly correlated fermionic lattice onto weakly correlated delocalized fermions. This allows a perturbative treatment of nonlocal correlations using some subsets of allowed diagrams to produce satisfactory corrections on top of the short-length scale correlations that are addressed by an exact solver.

The dual fermion formalism is an extension of the theory by Sarker for strongly correlated system. Sarker (1988) He proposed a strong coupling expansion of the solution from the atomic limit that predates widespread usage of DMFT. Similar ideas have also been proposed for the study of one dimensional system. Pairault et al. (1998) In this theory, the Hubbard model is mapped onto another interacting fermionic model in which the multi-particle hopping-exchange processes appear explicitly. The formulation is equivalent to the dual fermion formalism as currently known. Rubtsov et al. (2008); Brener et al. (2008); Li et al. (2008)

Starting from the action of itinerant electrons on a lattice that can be written as:

$$S\left[c^{*},c\right]=\sum_{\omega,k,\sigma}c^{*}_{\omega,k,\sigma}\left[i\omega+\mu-h_{k}\right]c_{\omega,k,\sigma}+\sum_{i}S_{loc}[c^{*},c],$$

(44)

where $\mu$ is the chemical potential, $h_{k}$ is the hoping term in momentum space, $c^{*}_{\omega,k,\sigma}(c_{\omega,k,\sigma})$ are the Grassmann variables corresponding to the creation (annihilation) operator, and $S_{loc}$ is the local part of the action. The lattice problem can be reframed into that of a set of impurities and an additional term to account for the remaining contributions:

$$S[c^{*},c]=\sum_{i}S_{\textrm{imp}}[c_{i}^{*},c_{i}]-\sum_{\omega,k,\sigma}c_{\omega,k,\sigma}^{*}[\Delta(i\omega)-h_{k}]c_{\omega,k,\sigma}.$$

(45)

Revisiting the expression of the partition function, a Hubbard-Stratonovich transformation can be applied on the second term, introducing new fermionic degrees of freedom. The action can then be expressed as:

	$\displaystyle S[c^{*},c;f^{*},f]$	$\displaystyle=$	$\displaystyle\sum_{i}S_{restr,i}[c^{*},c;f^{*},f]$		(46)
		$\displaystyle+$	$\displaystyle\sum_{\omega,k,\sigma}\frac{f_{\omega,k,\sigma}^{*}\,f_{\omega,k,\sigma}}{g^{2}(i\omega)[\Delta(i\omega)-h_{k}]}.$

Where $g$ is the single particle DMFT Green function and $S_{restr,i}$ is the action restricted to site $i$, and is defined by:

	$\displaystyle S_{restr,i}[c^{*},c;f^{*},f]$	$\displaystyle=$	$\displaystyle S_{\textrm{imp}}[c^{*},c]$		(47)
		$\displaystyle+$	$\displaystyle\sum_{\omega,\sigma}[f_{\omega,\sigma}^{*}g^{-1}(i\omega)c_{\omega,\sigma}+h.c.].$

The lattice fermionic degrees of freedom can then be integrated out of the action restricted to the respective sites following:

		$\displaystyle\int$	$\displaystyle\exp\left(-S_{restr,i}\left[c^{*},c,f^{*},f\right]\right)D\left[c^{*},c\right]$		(48)
		$\displaystyle=$	$\displaystyle Z_{imp}\exp\left(-\sum_{\omega,\sigma}f^{*}_{\omega,\sigma}g^{-1}(i\omega)f_{\omega,\sigma}+V\left[f^{*}_{i},f_{i}\right]\right).$

This last expression introduces the dual potential $V[f^{*}_{i},f_{i}]$ in terms of the new fermionic degrees of freedom. It is shown to include all $2n$-vertices for the impurity with $n=2,3,4,...$.