# Beyond Quantum Cluster Theories: Multiscale Approaches for Strongly
Correlated Systems

Herbert F. Fotso  Department of Physics, University at Albany (SUNY), Albany,
New York 12222, USA  Ka-Ming Tam  Department of Physics and Astronomy,
Louisiana State University, Baton Rouge, LA 70803, USA  Center for Computation
and Technology, Louisiana State University, Baton Rouge, LA 70803, USA  Juana
Moreno  Department of Physics and Astronomy, Louisiana State University, Baton
Rouge, LA 70803, USA  Center for Computation and Technology, Louisiana State
University, Baton Rouge, LA 70803, USA

###### Abstract

The degrees of freedom that confer to strongly correlated systems their many
intriguing properties also render them fairly intractable through typical
perturbative treatments. For this reason, the mechanisms responsible for these
technologically promising properties remain mostly elusive. Computational
approaches have played a major role in efforts to fill this void. In
particular, dynamical mean field theory (DMFT) and its cluster extension, the
dynamical cluster approximation (DCA) have allowed significant progress.
However, despite all the insightful results of these embedding schemes,
computational constraints, such as the minus sign problem in Quantum Monte
Carlo (QMC), and the exponential growth of the Hilbert space in exact
diagonalization (ED) methods, still limit the length scale within which
correlations can be treated exactly in the formalism. A recent advance to
overcome these difficulties is the development of multiscale many body
approaches whereby this challenge is addressed by introducing an intermediate
length scale between the short length scale where correlations are treated
exactly using a cluster solver such QMC or ED, and the long length scale where
correlations are treated in a mean field manner. At this intermediate length
scale correlations can be treated perturbatively. This is the essence of
multiscale many-body methods. We will review various implementations of these
multiscale many-body approaches, the results they have produced, and the
outstanding challenges that should be addressed for further advances.

##  I Introduction

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) However, this 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 for their low-energy description. Zhang and Rice (1988)
But despite its simplicity, an exact solution of the Hubbard model remains
elusive beyond one dimension. 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, 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 to 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 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 kinetic
energy and the potential energy due to couplings. 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 interacting 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 the 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 is 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 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-consistent 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 hoping amplitude in infinite dimension
leads to the vanishing of all diagrams except the local ones, and the model
loses all spatial dependence. This simplification led to the dynamical mean
field theory.

The 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 which allows an
analytical solution.

It is worth noting that the DMFT can be viewed as a formal generalization of
the coherent potential approximation (CPA) which was 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, but 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 the 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 other various novel methods have been
proposed. These include the density matrix renormalization group Fernández
_et al._ (2014); Zhu  _et al._ (2017), 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, it has
been proven that the single impurity problem does not exhibit a minus sign
problem. However, 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 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._ (2018)

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 the DMFT, the dynamical cluster
approximation (DCA) and the cellular dynamical mean field theory (CDMFT). We
proceed in section IV with a discussion of the extended dynamical mean field
theory 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 that encompasses
various commonly used diagrammatic approximations, and that 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 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 methods before ending with our
conclusions.

##  II Dynamical Mean Field Theory

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 on a
lattice that can hop between nearest-neighbor sites with a hopping integral
$t$, 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}=-t\sum_{<i,j>,\sigma}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$.

Figure 1: Schematic illustration of the Hubbard model on a square lattice.
Electrons are allowed to hop between sites and the interaction is local.

Within DMFT the Hubbard model is reduced to a single impurity Anderson model
with the hybridization determined by a self-consistent condition that sets the
Green function of the impurity problem equal to that of the Hubbard model in
the limit of infinite spatial dimension. 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$ $\displaystyle=$
$\displaystyle-\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 Grassmaann 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 model being considered. 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)

The idea of the dynamical mean field theory is to relate the full model with
spatial dependence to one without spatial dependence. Stripping off the
spatial dependence of the above action, we obtain:

$\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}})$
(4) $\displaystyle+$ $\displaystyle
U\int^{\beta}_{0}d\tau\psi_{\uparrow}^{*}(\tau)\psi_{\uparrow}(\tau)\psi^{*}_{\downarrow}(\tau)\psi_{\downarrow}(\tau).$

This is exactly the action of the single impurity Anderson model with the
hybridization given by $G_{0}$. The next assumption of the DMFT is that the
optimal approximation of the full Hubbard model is given by the single
impurity Anderson model when the self-energy of both models is the same. This
can be justified via the concept of the “cavity method” whereby all degrees of
freedom are integrated out except for the site labelled by the index $0$.
Then, the effective action, $S_{eff}$, becomes:

$\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}\\}]).$ (5)

For a system with $N+1$ sites, the effective action can be formally written as

$\displaystyle S_{eff}[\psi_{0,\sigma}^{*}\psi_{0,\sigma}]=S_{0}$ (6)
$\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(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 $\eta_{i,0}=t_{i,0}\psi_{0,\sigma}$, with $t_{i,0}$ the hopping from
site $i$ to $0$, and $S_{0}$ is the local action as in Eq. 4.

Consider the expansion of $S_{eff}$ with respect to the hopping from the 0-th
site to the rest of the lattice, by truncating terms higher than the second
order we arrive 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}}).$
(7)

Rewriting the effective action in frequency space, we get:

$\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}),$

where the effective action is equivalent to that of the single impurity
Anderson model and $\mathcal{G}$ acts as an effective (mean) field,

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

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).$
(10)

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)}$
(11)
$\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 as

$\displaystyle
G_{impurity,0}(\omega)=\sum_{\mathbf{k}}\frac{1}{(2\pi)^{D}}G_{lattice}(\mathbf{k},\omega)$
(12)
$\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 state to simplify the calculation. The bare density of state
for the hypercubic lattice at the limit of infinite spatial dimensions can be
exactly calculated. Müller-Hartmann (1989a, b) The DMFT algorithm can be
summarized as in Fig. 2

Figure 2: Algorithm for the DMFT. Step 1.The impurity solver provides the
impurity Green function, $G_{imp}$. Step 2. The Dyson equation is used to
extract the self-energy $\Sigma$. Step 3. The self-energy is coarse-grained
over the entire first Brillouin zone to obtain the lattice Green function.
Step 4. The bath Green function for the impurity problem is calculated and fed
into the impurity solver. Step 5. Repeat steps 1 to 4 until convergence is
obtained for the self-energy or Green function.

In brief the goal of the DMFT is to set a scheme to achieve a solution for the
self-consistency equation which equate the self-energy of the lattice and
impurity models. The method in itself is a generalization of the usual mean
field theory, but unlike the usual mean field theory, say for the Ising model,
the mean field here is not an order parameter. Instead it is a function of the
time or frequency. Indeed, there is no explicit order parameter, as no broken
symmetry is considered. The ‘order parameter’ here addresses the dynamics of
the system.

Within DMFT the spatial fluctuations are completely suppressed as in a
traditional mean field theory. As we will see, a major improvement of 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 various
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 which is accurate for a wide
range of parameters, particularly for solving interacting problems with random
disorder 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 already completely suppressed, and the quantum fluctuations
of the transverse field can be addressed through a quantum impurity problem.
Bray and Moore (1980) The idea of 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)

##  III Cluster Route for extending the DMFT

A natural and direct avenue to generalize the DMFT is to consider the spatial
fluctuation by including more than a single impurity site. This type of
cluster DMFT remains an important method for the study of strongly correlated
systems, as it allows a systematic correction with respect to the increase of
the cluster size unlike perturbative expansion methods. Moreover, one can
envision that a perturbative expansion on top of the cluster method would
produce an even better result as the bare effective Hamiltonian or action is
presumably more accurate and already includes a substantial amount of spatial
fluctuations.

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. A cluster of impurities in the 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)

Another approach for a cluster generalization of the 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) Instead of working in real space, the cluster
of impurities is used, with impurities representing different patches of the
first Brillouin zone. The lattice quantities are averaged over a patch for a
coarse-grained quantity. The advantage is that the method is manifestly
translational 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, however the self-energy which is obtained by
solving the impurity model does not have spatial dependence. Therefore, the
coarse-graining for DMFT is done over the entire first Brillouin zone with one
single impurity. In the dynamical cluster approximation, coarse-graining is
performed on a cluster that spans the entire first Brillouin zone and is made
of $N_{c}$ impurities. Each of the $N_{c}$ impurities represents an equal
patch of the Brillouin zone.

$\displaystyle
G_{0,impurity}(\mathbf{K},\omega)\approx\sum_{\mathbf{k}}\frac{1}{(2\pi)^{D}}G_{lattice}(\mathbf{k},\omega)$
(13)
$\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 sites with momentum $\mathbf{K}$.

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

Figure 3: Algorithm for the DCA. 1. The impurity solver provides the impurity
Green function, $G_{imp}(K)$. 2. The Dyson equation is used to extract the
self-energy $\Sigma(K)$. 3. The self-energy is coarse-grained over the patches
that fill the first Brillouin zone to obtain the lattice Green function,
$G(K)$ 4\. The bath Green function, $\mathcal{G}(K)$, for the impurity problem
is feed into the impurity solver. 5. Repeat from step 1 until convergence is
obtained for the self-energy or Green function.

##  IV Extended Dynamical Mean Field Theory

The earliest attempt to include the effect of nonlocal interaction was
motivated by the competition between 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 interaction.

In the 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 $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 of 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}})$
(14) $\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 and
$\mathcal{G}$ and $\mathcal{D}$ act as effective fields for the fermionic bath
and 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})$
(15)

and

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

The retarded coupling, can be understood in term of a bosonic field coupled to
the impurity charge density. The EDMFT remedies the limitation of DMFT in
which 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 field for the fermionic and for the bosonic baths have
to be updated. The algorithm can be summarized as in Fig. 4

Figure 4: Algorithm for the EDMFT. 1. The impurity solver provides the
impurity Green function for the electrons, $G_{imp}$ and the impurity Green
function for the Bosonic field, $D_{imp}$. 2. Dyson equation is used to
extract the self-energy $\Sigma$ and the polarization, $\Pi$. 3. The self-
energy and the polarization is coarse-grained over the entire first Brillouin
zone to obtain the lattice Green function. 4. The bath Green function for the
electrons and the Bosonic field of the impurity problem is feed into the
impurity solver. 5. Repeat from step 1 until convergence is obtained for the
self-energy and polarization or Green functions.

##  V $1/d$ Expansion

One of the earliest attempts to consider the effect of finite dimension 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) In such an expansion, the hopping term carries the factor of
$1/\sqrt{d}$. Thus, an expansion with respect to the hopping regains the
dimensional dependence at the expense of multiple sites. Work by Schiller and
Ingersent has studied the case of two impurities for the Falicov-Kimball
model. Schiller and Ingersent (1995) The approach has also been applied on the
half-filled Hubbard model. Zaránd  _et al._ (2000)

The $1/d$ corrections can be recovered by considering the 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),-1}_{\alpha,\beta}(\tau-\tau^{\prime})\psi_{\beta\sigma}(\tau^{\prime})$
(17) $\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}\;.$
(18)

For skeleton diagrams of the order of ${\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)})\;,$ (19)
$\displaystyle\Sigma^{\rm lattice}_{1}$ $\displaystyle=$
$\displaystyle\Sigma^{(2)}_{12}\;.$ (20)

The lattice Green function can be obtained with 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)\;,$ (21)

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 limitations of the $1/d$ expansion is its difficulty to capture long
range fluctuations. Perhaps more importantly, the truncation at finite order
may lead to non-analytic properties of dynamical quantities. Pruschke  _et
al._ (2001) The method is an example of the nested cluster scheme and it also
related to the recently proposed self-energy embedding theory. Biroli  _et
al._ (2004); Kananenka  _et al._ (2015); Zgid and Gull (2017)

##  VI Multiscale many body

Given the challenge of solving a large cluster by numerical simulations and
the lack of accurate analytical methods, perturbation theory is a possible
route for the exploration of physics beyond local theories. Initial efforts
aimed to use perturbative methods for solving 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 method to allow the correlation among different channels. It is,
therefore, conceptually appealing approach for systems in which particle-hole
or particle-particle two 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. The FLEX allows the vertex contributions from different
channels to be correlated, thus the name of fluctuation exchange
approximation.

The general scheme of the multiscale many body 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 in a mean field, while the intermediate length scale is
treated by some form of a 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 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 the QMC as a cluster
solver.

The general construction of the multiscale many body method, as summarized in
Fig.6, 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.

Figure 5: Conceptual idea of the multiscale many body approach. The original
lattice model is mapped onto a small cluster (yellow-orange-red colors),
embedded in a larger cluster (rainbow colors) embedded in a mean field (orange
color). The information passed between the clusters and the effective medium
is composed of irreducible quantities like the self-energy and the fully
irreducible vertex function  Figure 6: Algorithm for the multiscale many body
method. 1. The fully irreducible vertex ($\Lambda$) and optionally the self-
energy ($\Sigma$) are obtained from the DMFT or DCA solver. 2. The parquet
method or some other perturbative method is used to calculate the full vertex
to restore spatial dependence or enhance spatial resolutions in the full
vertex ($F$). 3. The full vertex is fed into Schwinger Dyson equation and
Dyson equation to update the self-energy and the Green function respectively.
4. The Green function can be coarse-grained to be fed back to the DMFT or DCA
solver as the bath Green function for self-consistent calculations.

##  VII Diagrammatic methods and the Parquet Formalism

A major route for perturbative expansions around the DMFT solution is based on
the parquet diagrams. This route encapsulates many of the approximations that
have been proposed 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, at the lowest orders, 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 adapt the notation $F^{ph}(2-4)_{1,3}$ for the numerical
implementation of the single band Hubbard model. Other representations are
available. 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 method.

The vertices can be categorized by extending the notion of diagram
reducibility to the two-particle level as illustrated by Fig.7. 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 other 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 is reducible in the same channel is known as the irreducible vertex
$\Gamma$. The subset of irreducible vertices that is irreducible in any
channel is called the fully irreducible vertex $\Lambda$.

Figure 7: Illustration of the reducibility of a diagram in the particle-hole
horizontal channel. (a) The diagram is reducible in the particle-hole channel
in the sense that the diagram can be separated by cutting two horizontal Green
function lines. (b) Examples of diagrams which are irreducible with respect to
the particle-hole horizontal channel. The one on the left hand side is
reducible in the particle-hole vertical channel. The one on the right hand
side is reducible in the particle-particle hole channel. (c) Examples of fully
irreducible diagrams. They cannot be separated into two parts by cutting to
Green functions lines in any one of the three channels. From RefYang  _et al._
(2009).

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 are not going
to explore this direction in detail here. The above decomposition is the most
natural one in the sense that the method can be easily understood in term of
an iteration process. The higher order diagrams are all generated by
iteratively replacing the vertex function at a lower order approximation.

Furthermore, since 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},$
(22) $\displaystyle\Gamma_{m}$ $\displaystyle=$
$\displaystyle\Gamma^{PH}_{\uparrow\uparrow;\uparrow\uparrow}-\Gamma^{PH}_{\uparrow\uparrow;\downarrow\downarrow},$
(23) $\displaystyle\Gamma_{s}$ $\displaystyle=$
$\displaystyle\Gamma^{PP}_{\uparrow\downarrow;\uparrow\downarrow}-\Gamma^{PP}_{\uparrow\downarrow;\downarrow\uparrow},$
(24) $\displaystyle\Gamma_{t}$ $\displaystyle=$
$\displaystyle\Gamma^{PP}_{\uparrow\downarrow;\uparrow\downarrow}+\Gamma^{PP}_{\uparrow\downarrow;\downarrow\uparrow},$
(25)

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}}){}$
(26)
$\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),$ (27)

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}},$
(28)
$F_{r^{\prime}}(Q)_{P,P^{\prime}}=\Gamma_{r^{\prime}}(Q)_{P,P^{\prime}}+\Psi_{r^{\prime}}(Q)_{P,P^{\prime}},$
(29)

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\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;$
(30)
$\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\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;$
(31)
$\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}$ 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.
28 and 29) the crossing symmetries (that are a consequence of the Pauli
exclusion principle) 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}$
(36)
$\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}}.$

These relations can allow us to restore the crossing symmetry for the full
vertices without heavy computational overhead. Additionally, the parquet
formalism enables the identification of separate channel contributions in
discerning the origin of instabilities identified in the susceptibility. This
is seen in the expression of the full vertex of Eqs.(VII, VII, VII, VII) that
include cross-channel terms.Yang  _et al._ (2011a)

Disentangling the nature of diverging correlations: A key feature of the
parquet equations is that they express the contribution to a given scattering
channel of processes from other channels. This is essential in understanding
the mechanisms for key processes and phase transitions…

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 can be neatly summarized as in the Fig. 8.

  * •

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 G-based 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 and a second order perturbation.

  * •

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)

Figure 8: Hierarchy of approximations encompassed by the parquet
approximation. See the text for details. S.D. and B.S. stand for Schwinger-
Dyson equation and Bethe-Salpeter equation respectively.

The parquet formalism dates back to the 50’s, but as can be seen from the
analytical forms of the governing equations, a general analytical solution is
rather hard to track. Pomeranchuk  _et al._ (1956); De Dominicis and Martin
(1964) Various approximations to simplify the equations have been proposed to
solve a set of problems ranging from the Anderson impurity to nuclear
structure over the years. 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)

Perhaps the more 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 antiferromagnetic instability has a strong divergence
for momentum transfer $(\pi,\pi)$. For this reason, the vertex is represented
by numbers spanning a wide range with many orders of magnitude. This clearly
is a recipe for possible instabilities in the iterative solution for the
vertex functions at low temperate and on the verge of instability towards long
range order. 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 algorithm for the numerical solution of the parquet formalism is
summarized in Fig.9

Figure 9:  Algorithm for the parquet method. 1. Read in the fully irreducible
vertex, $\Lambda$. 2. Solve the Bethe Salpeter equation for irreducible
vertex, $\Gamma$. 3. Solve the parquet equations for the full vertex, $F$. 4.
Use the Dyson Schwinger equation to solver for the self-energy. 5. The Green
function is updated by the Schwinger equation. Repeat from the step 2 until
convergence in the full vertex and self-energy is obtained.

##  VIII Dynamical vertex approximation (D$\Gamma$A)

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 computation costs, Held (2014) besides the two-dimensional Hubbard
model, it has been applied to several problems including three-dimensional
Hubbard model, 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 according 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 picking up 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 dominant.
For most interesting regimes, such as the d-wave pairing in the Hubbard model,
presumably all channels could contribute and one has to resort to solving the
full set of equations of the parquet formalism.

The algorithm is summarized in Fig. 10

Figure 10:  Algorithm for the dynamical vertex approximation. 1. The fully
irreducible vertex is obtained from the DMFT solver. 2. Perturbative methods
or the parquet method is used to obtain the self-energy. 3. The Dyson equation
is solved for the Green function. The results can be read out after this step.
For self-consistent calculation, 4. The Green function can be coarse-grained
to obtain the bath Green function for the impurity problem and the procedure
is repeated from step 1 until convergence in the Green function is obtained.

##  IX Dual Fermions

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],$ (40)

where $mu$ is the chemical potential, $h_{k}$ is the lattice bare Green
function, $c^{*}_{\omega,k,\sigma}(c_{\omega,k,\sigma})$ are the Grassmann
variables corresponding to creation and annihilation operators, and $S_{loc}$
is the local part of the action, the lattice problem can be rephrased into
that of a set of impurities and an additional term to account for 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}.$
(41)

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]$ (42) $\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]$ (43) $\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]$
(44) $\displaystyle=$ $\displaystyle
Z_{imp}\exp\left(-\sum_{\omega,\sigma}f^{*}_{\omega,\sigma}g^{-1}_{\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,...$.

An explicit expression for the dual potential is obtained by expanding both
sides of this equation and comparing the resulting expressions order by order.
The dual potential to lowest order reads

$\displaystyle V[f^{*},f]$
$\displaystyle=\frac{1}{4}\sum_{\omega\omega^{\prime}\nu}\sum_{\sigma_{1},\sigma_{2},\sigma_{3},\sigma_{4}}\gamma_{\sigma_{1},\sigma_{2},\sigma_{3},\sigma_{4}}(i\omega,i\omega^{\prime},i\nu)$
(45) $\displaystyle\times
f_{\omega+\nu,\sigma_{1}}^{*}f_{\omega,\sigma_{2}}f_{\omega^{\prime},\sigma_{3}}^{*}f_{\omega^{\prime}+\nu,\sigma_{4}}+\cdots.$

Thus, nonlocal correlations are addressed by solving the many-body problem
with bare green’s function $g$ and interaction potential $V$. The lattice
fermions can finally be integrated out to produce an action that only depends
on the dual fermions:

$S_{D}[f^{*},f]=-\sum_{\omega,k,\sigma}f^{*}_{\omega,k,\sigma}G^{d,0}_{\sigma}(i\omega,k)f_{\omega,k,\sigma}+\sum_{i}V[f^{*}_{i},f_{i}].$
(46)

With the bare dual fermion Green’s function defined by:

$G^{d,0}_{\sigma}(i\omega,k)=-\frac{g^{2}(i\omega)}{[g(i\omega)+\left(\Delta(i\omega)-h_{k}\right)^{-1}]}.$
(47)

The action of Eq.(46) is the tool for accounting for nonlocal correlations. It
can be treated using diagrammatic perturbation theory. In this context, the
interaction potential is usually truncated to the 4-point vertex. For most
practical calculations, higher order terms of the dual potential are
truncated, Gukelberger  _et al._ (2017) though they may have non-negligible
effect. Ribic  _et al._ (2017); Katanin (2013) However, the formalism is
shown, by construction, to be convergent both in the strong coupling and in
the weak coupling regime. The process for solving the formalism follows a
typical diagrammatic procedure. The impurity Green function and the vertex are
obtained from the DMFT calculation. The impurity Green function is used to
evaluate the Dual fermions non-interacting Green function. The vertex and the
non-interacting Green’s function are then used for a self-consistent
diagrammatic solution with a given subset of all the allowed diagrams. This
solution produces the dressed dual fermion Green’s function. The dual Green
function can subsequently be used to evaluate the lattice Green function.
Alternatively, phase transitions can be studied directly using the dual
fermion diagrams since the instability identified from the dual fermions
diagrams is found to be equivalent to that of the lattice Green function. In
general, upon solving the dual fermions diagrams, a new expression for the
impurity self-energy or of the impurity hybridization can be extracted and fed
back into the impurity and the entire procedure repeated iteratively until
convergence. The algorithm of dual fermion method is summarized in Fig.11.

Figure 11: Algorithm for the dual fermion method. 1. The DMFT solver provides
the Green function and the full (reducible) vertex functions. 2. The Green
function and the irreducible vertex functions from the DMFT solver are used to
construct the dual fermion bare Green function and the dual potential. 3. A
perturbative method or the parquet method can be used to calculate the self-
energy and the Green function of the dual fermions. 4. Inverse transform of
the Green function from dual fermions to physical fermionic degree of freedom.
The results can be read out after this step. For self-consistent calculation.
5. The Green function can be coarse-grained to obtain the bath Green function
of the DMFT impurity problem.

The dual fermion formalism, as discussed above, was initially introduced to
add nonlocal correlations to the DMFT result. It was subsequently extended to
DCA.Yang  _et al._ (2011b); Iskakov  _et al._ (2016a) In this context, it
serves as a way to address intermediate length scales beyond the short ones
that are treated by the cluster solver of the DCA formalism.

Because the dual fermions concept is rather general, it has been the subject
of numerous developments. Hafermann  _et al._ (2009a, b); Iskakov  _et al._
(2016b); Krivenko  _et al._ (2010); Astretsov  _et al._ (2020); Hafermann  _et
al._ (2008) A notable generalization is the treatment of disorder with the
formalism. Terletska  _et al._ (2013); Haase  _et al._ (2017) The effect of
disorder in correlated systems has long been an important outstanding problem
in condensed matter physics, in particular for the two dimensional case. While
experimentally available systems such as semi-conductors usually involve long
range Coulomb interactions, the problem of a perhaps simpler Anderson-Hubbard
model which have both interaction and disorder at local sites still represents
an outstanding challenge which has attracted a lot of attention.

Another important development is the the treatment of nonlocal interactions,
such as that of the extended Hubbard model with nearest neighbor interaction.
The inclusion of nonlocal interaction opens up the possibility for interesting
physics such as charge density wave and perhaps more exotic bond order wave.
The idea of dual fermion has also been generalized for the extended dynamical
mean field theory (EDMFT) to the dual bosons theory as we discuss in the
following section. Rubtsov  _et al._ (2012)

The dual fermion method has also been used to generalize the real space DMFT
to the real space dual fermion method which allows the study of systems with
open boundary condition. Takemori  _et al._ (2016)

##  X Dual Bosons, extension of EDMFT

One can apply ideas similar to those of the dual fermions for DMFT to EDMFT.
Rubtsov  _et al._ (2012) As discussed in section IV, the nonlocal part of the
EDMFT is not limited to the hopping term, but also includes the interaction.
Si and Smith (1996) For this problem, the action can be written as:

$\displaystyle
S\left[c^{*},c\right]=\sum_{\omega,k,\sigma}c^{*}_{\omega,k,\sigma}\left[i\omega+\mu-
h_{k}\right]c_{\omega,k,\sigma}$
$\displaystyle+\frac{1}{2}\sum_{\omega,k}V(k)n_{\omega,-k}n_{\omega,k}+\sum_{i}S_{local}[c^{*},c],$
(48)

The next step is to introduce a Hubbard-Stratonovich transformation as in the
dual fermion formulation for the nonlocal bilinear term for the kinetic
energy. In addition, the nonlocal density-density potential energy term is
decoupled by another Hubbard-Stratonovich transformation for the bosonic
charge density. In keeping with the dual fermion method parallel, the original
fermion degrees of freedom, and the bosonic charge density can formally be
integrated out exactly, yielding a new effective potential given by the
vertices of the impurity model. The action becomes

$\displaystyle
S_{D}[f^{*},f]=-\sum_{\omega,k,\sigma}f^{*}_{\omega,k,\sigma}G^{d,0}_{\sigma}(i\omega,k)f_{\omega,k,\sigma}$
$\displaystyle-\sum_{\omega,k}\phi^{*}_{-\omega,-k}D^{d,0}(i\omega,k)\phi_{\omega,k}+\sum_{i}V[f^{*}_{i},f_{i},\phi_{i}].$
(49)

Similar to the dual fermions, the potential of the dual variables can be
expanded and truncated at finite order via the vertex functions of the
impurity problem. Perturbative methods can be used to solve the effective
problem expressed in term of the dual variables. The method has been applied
to problems with nonlocal interaction, such as the extended Hubbard model and
models with long range Coulomb coupling. Van Loon  _et al._ (2014); Stepanov
_et al._ (2016); van Loon  _et al._ (2016); Vandelli  _et al._ (2020); Peters
_et al._ (2019); Van Loon  _et al._ (2014)

The algorithm of dual bosons method is summarized in Fig. 11

Figure 12: Algorithm for the dual boson method. This algorithm is parallel
that of the dual fermion method, except that both the bosonic Green function
and the dual potential with the bosonic degree of freedom are needed. 1\. The
DMFT solver provides the Green function and the irreducible vertex functions.
2. The Green function and the reducible vertex functions for both fermion-
fermion vertex ($F^{f-f}_{imp}$) and fermion-boson vertex ($F^{f-b}_{imp}$)
from the DMFT solver are used to construct the Green functions ($G_{dual}$ and
$D_{dual}$) and the dual potential ($V_{dual}$). 3. Perturbative methods can
be used to calculate the self-energy and the Green function of the dual
fermions and dual bosons. 4. Inverse transformation of the Green function from
dual fermions and dual bosons to physical fermionic and bosonic degrees of
freedom. The results can be read out after this step. For the full self-
consistent calculation, 5. The Green function can be coarse-grained to obtain
the bath Green function of the DMFT impurity problem.

##  XI GW approximation with three-particle irreducible vertex

A conventional scheme to address quantum fluctuations beyond Hartree-Fock is
the so-called GW approximation. It was developed for solutions of the electron
gas problem in the 50’s. Due to its simplicity, it has been widely adapted in
density functional theory calculations and it can be derived from many-body
perturbation theory. The form of the self-energy in the GWA is kept as that of
the Hartree-Fock, but the interaction, originally Coulomb interaction, is
dynamically screened. Quinn and Ferrell (1958); DuBois (1959a, b); Hedin
(1965); Aryasetiawan and Gunnarsson (1998)

Recent proposals have utilized the properties of the GW approximation and
generalized it for not only the self-energy, but also the polarizability of
the two-particle vertex functions. The three points vertex is calculated via
the effective impurity problem similar to that of the EDMFT. It was introduced
by Aryal and Parcollet as Triply Irreducible Local Expansion (TRILEX). Ayral
and Parcollet (2015); Vučičević  _et al._ (2017) The method can be formally
derived from a functional of the vertex given by three-particle irreducible
diagrams de Dominicis and Martin (1964); De Dominicis and Martin (1964).

Figure 13: Diagrams for generating momentum dependence for the self-energy of
the single particle function and the polarization of the two particle
function.

The key step for incorporating nonlocal fluctuations is to consider the
polarization by the local boson-fermion vertex from the numerical solution of
the impurity problem, see Fig.13 for the diagrams. The rest of the algorithm
is parallel to that of the EDMFT. The updated self-energy and the
polarizability are coarse-grained and fed back into the EDMFT effective
impurity problem.

The algorithm for the TRILEX method is summarized in Fig.11.

Figure 14: Algorithm for the TRILEX method. 1. The DMFT solver provides the
boson-fermion vertex. 2. The Hatree diagrams of the GW approximation is
calculated with the vertex from the DMFT for both the self-energy and the
polarization to generate spatial dependence. 3. The Dyson equation is employed
to calculate the Green function for both the fermions ($G$) and the bosons
($D$). 4. The coarse-grained Green functions can be fed back to the effective
DMFT impurity problem. Repeat from the step 1 until convergence is attained.

##  XII Functional Renormalization Group

The functional renormalization group for fermionic systems evolved from
Luttinger liquid calculations for the one dimensional system. Voit (1995);
Tomonaga (1950); Luttinger (1963); Giamarchi (2003) The full vertex in this
context is usually denoted as “g”, for this reason, the method is often
referred as g-ology. Menyhárd and Sólyom (1973); Sólyom (1973) From the
renormalization group perspective, the frequency dependence of the vertex
function is not the most relevant contribution, thus for systems that don’t
have explicit frequency dependence in the bare vertex, such as electron-phonon
coupling, the frequency or time dependence is often neglected. Shankar (1994)

The generalization of g-ology beyond one dimension gained increased attention
after the work by Shankar. Shankar (1991, 1994) Traditional implementation of
Wilsonian renormalization group only considers the flow of a handful of
coupling constants, which is physically the case for the one dimensional
systems as there are only two Fermi points instead of a surface. The
functional formalism for all coupling constants are considered when the method
is generalized to high spatial dimensions.

The idea of performing renormalization group on a system with extended Fermi
surface was suggested by Anderson. Anderson (1984) Benfatto and Gallavotti
Benfatto and Gallavotti (1990), and Feldman and Trubowitz Feldman and
Trubowitz (1990, 1991) studied the stability of the Fermi liquid against
perturbations. Effective theories based on the idea of renormalization group
for Fermi liquid were derived by Polchinski, Polchinski (1992) and that for
the superconductor by Weinberg. Weinberg (1994)

Most studies in the late 90’s and early 00’s focused on the two dimensional
Hubbard model. Zanchi and Schulz (2000, 1998); Honerkamp  _et al._ (2001);
Halboth and Metzner (2000a, b); Kampf and Katanin (2003); Tsai and Marston
(2001) The two important developments paved the way for using FRG to improve
DMFT solutions. First, originally, the conventional renormalization group uses
the energy cutoff. Shankar (1994) This is clearly not a unique choice to
characterize the approach to divergences and other cutoffs such as
temperature, interaction, and even the hybridization cutoff for the impurity
problem have been proposed and implemented. Honerkamp  _et al._ (2004);
Honerkamp and Salmhofer (2001); Kinza and Honerkamp (2013) Second is the
implementation of the frequency dependent vertex. This is largely motivated by
the interest in studying electron-phonon coupling. Honerkamp and Salmhofer
(2005); Tsai  _et al._ (2005) Since it is not the most relevant contribution
driving the instability, the frequency dependence even at the single particle
level is not considered, that is, the self-energy is not renormalized. The
situation changes when explicitly retarded interaction is considered in the
models with electron-phonon coupling and Anderson impurity model. Tam  _et
al._ (2007a, b); Sédéki  _et al._ (2000); Karrasch  _et al._ (2008)

It was realized early in the study of the frequency dependent vertex that the
method can be used to study mean field fluctuations. Attempts have been made
to use it for the study of Gaussian fluctuations of the slave boson mean field
solution for the t-J model. Kotliar and Liu (1988) Progress in this research
direction has been limited because the conventional assumption that the
bandwidth should be larger than the interaction is not met for a wide range of
parameters of the slave boson t-J model solution.

Related ideas have been revived in the past few years to consider the mean
field fluctuations. Instead of a static mean field solution, the solution of
DMFT is considered. Taranto  _et al._ (2014) The technique developed as a
result can be directly applied to build up the spatial dependence from the
DMFT or DCA solution. The main idea is that the bare vertex is replaced by the
full vertex from the DMFT or DCA, and the self-energy is set to that of the
DMFT or DCA for the initial conditions of the renormalization flows.

From a general perspective, the renormalization group can be viewed as summing
up diagrams that are often referred to as forming the leading divergence.
Zheleznyak  _et al._ (1997a) In particular if second order perturbation is
used in the renormalization group calculation, the diagrams generated have the
same topology as those of the parquet method discussed above.

There is a subtlety though, the integration over the internal frequency and
momentum are not the same as those of the parquet method. The difference
depends on the different cutoff schemes being used, but the summation is never
performed over the full range of frequency and momentum, instead it is done on
a shell of energy range, and iteratively approaches the desired low energy or
low temperature manifold.

It is long believed that the leading divergence of the 1-loop FRG and the
parquet should be equivalent. Diekmann and Jakobs (2020); Sólyom (1973);
Zheleznyak  _et al._ (1997b) A recent proposal of multi-loop flow equation for
the four-point vertex framework showed that the FRG flow consisting of
successive one-loop calculations is equivalent to a solution of the parquet
equations. This further supports the idea of FRG as a possible alternative to
solving the parquet equations.Kugler and von Delft (2018a, b)

Given the same topology of the diagrams, one might expect that the leading
divergence among these two methods should be the same even outside of the
multi-loop setting. Therefore, for the purpose of looking for the instability
from metal to ordered phase, these two approaches should be expected to give
the same result. On the other hand, for the metal phase with no proximity to
an instability, the results are not naively equivalent.

The clear advantage of functional renormalization group is the ease of
numerical calculations. Unlike, the parquet or even the simplified dynamical
vertex approximation, the functional renormalization group works directly on
the full reducible vertex. The numerical solution does not involve solving
self-consistent equations. It is given by the flow of the full vertex as the
cutoff is lowered to the desired energy or temperature. The formulation is
represented in term of ordinary differential equations.

The main challenge of the parquet method is the instability of the numerical
simulation. Even if we assume that the solution exists, a robust method for
attaining the solution is highly non-trivial. This is at the end of the day
the main issue with methods that seek self-consistent solutions. Unlike the
dynamical mean field theory in which only single particle quantities are
involved, the two-particle methods, full parquet or simplified form, involve
solving for the vertex through self-consistent equations. Divergence in the
vertex function is expected to occur as the temperature is lowered. Therefore,
equations with variables spanning a large range of values with many orders of
magnitude have to be solved. This is clearly a non-trivial task from the
perspective of numerical simulation. Although one can solve the equations at
high temperature, it is not always clear how far the solution can be pushed
down in temperature. Contrary to this, the functional renormalization group
method sums up the same set of diagrams without the necessity of solving
relevant equations self-consistently. The divergence is approached step by
step rather than via shooting as is done in the self-consistent solution.

While the advantage of functional renormalization group from the point of view
of numerical simulation is clear. The justification of its usage for nonlocal
corrections with the vertex function from the DMFT is not obvious. The
conventional wisdom of renormalization group calculations is that the
bandwidth should be large and the interaction is a small parameter. The full
vertex from the DMFT is not necessarily small compared to the effective
bandwidth. Therefore, the conventional wisdom of justifying the low order
expansion is not incontrovertibly fulfilled. Moreover, these couplings
generate self-energy corrections, which have been shown to be important for
studying systems with retardation effects, leading to the renormalization of
the Fermi velocity and quasiparticle lifetime. Thus, these couplings have to
be kept even though they are often ignored in the non-retarded system. The
effective system being solved is retarded although the original Hubbard model
is not.

To leading second order expansion, the renormalization group equations for the
scale $(\Lambda)$ dependent full vertex function,
$F_{\Lambda}(k_{1},k_{2},k_{3})$, Zanchi and Schulz (1998, 2000) and the self-
energy, $\Sigma_{\Lambda}(k)$ for a spin rotational invariant two-body
interacting system, are given by:

$\displaystyle\partial_{\Lambda}F_{\Lambda}(k_{1},k_{2},k_{3})=$
$\displaystyle-$ $\displaystyle\int
dp\partial_{\Lambda}[G_{\Lambda}(p)G_{\Lambda}(k)]F_{\Lambda}(k_{1},k_{2},k)F_{\Lambda}(p,k,k_{3})$
(50) $\displaystyle-$ $\displaystyle\int
dp\partial_{\Lambda}[G_{\Lambda}(p)G_{\Lambda}(q_{1})]F_{\Lambda}(p,k_{2},q_{1})g_{\Lambda}(k_{1},q_{1},k_{3})$
$\displaystyle-$
$\displaystyle\int\\!dp\partial_{\Lambda}[G_{\Lambda}(p)G_{\Lambda}(q_{2})][-2g_{\Lambda}(k_{1},p,q_{2})F_{\Lambda}(q_{2},k_{2},k_{3})$
$\displaystyle+$ $\displaystyle
F_{\Lambda}(p,k_{1},q_{2})F_{\Lambda}(q_{2},k_{2},k_{3})+F_{\Lambda}(k_{1},p,q_{2})F_{\Lambda}(k_{2},q_{2},k_{3})],$
$\displaystyle\partial_{\Lambda}\Sigma_{\Lambda}(k)\\!=\\!-\\!\int\\!dp\partial_{\Lambda}[G_{\Lambda}(p)][2F_{\Lambda}(p,k,k)\\!-\\!F_{\Lambda}(k,p,k)],$
(51)

where $k=k_{1}+k_{2}-p$, $q_{1}=p+k_{3}-k_{1}$, $q_{2}=p+k_{3}-k_{2}$, $\int
dp=\int d\mathbf{p}\sum_{\omega}1/(2\pi\beta)$, and $G_{\Lambda}$ is the self-
energy corrected propagator at cutoff $\Lambda$.

Figure 15: RG equations for the full vertex. There is a pair of red arrows
which represent the derivative of the product of propagators with respect to
the RG scale $\Lambda$.

The RG equation can be presented in term of diagrams. For a non-retarded
system, the low energy instability can be obtained from the renormalization
flow of the couplings, and different phases can be identified by the fixed
points corresponding to the relevant spin and charge modes. This procedure may
not be appropriate for systems with vertex coupling. One can explicitly
construct the flows of the susceptibilities of different order parameters. For
example the pairing susceptibility, $\chi^{\delta}(\mathbf{k},\omega)$, is
defined by:

$\displaystyle\chi^{\delta}_{\Lambda}(0,0)\\!\\!=\\!\\!\int\\!\\!\int
dp_{1}dp_{2}\langle
c_{p_{1},\downarrow}c_{-p_{1},\uparrow}c_{-p_{2},\uparrow}^{\dagger}c_{p_{2},\downarrow}^{\dagger}\rangle_{\Lambda}.$
(52)

The RG equations are:

$\displaystyle\partial_{\Lambda}\chi^{\delta}_{\Lambda}(0,0)=\int
dp\partial_{\Lambda}[G_{\Lambda}(p)G_{\Lambda}(-p)](Z^{\delta}_{\Lambda}(p))^{2},$
(53) $\displaystyle\partial_{\Lambda}Z^{\delta}_{\Lambda}(p)=-\int
dp^{\prime}\partial_{\Lambda}[G_{\Lambda}(p^{\prime})G_{\Lambda}(-p^{\prime})]$
(54) $\displaystyle
Z^{\delta}_{\Lambda}(p^{\prime})g_{\Lambda}(p^{\prime},-p^{\prime},-p,p).$

The function $Z_{\Lambda}^{\delta}(p)$ is the effective vertex in the
definition for the susceptibility $\chi_{\Lambda}^{\delta}$. The RG equations
for susceptibilities are solved with initial condition
$\chi^{\delta}_{\Lambda=\Lambda_{0}}(0,0)=0$. The dominant instability in the
ground state is given by the most divergent susceptibility by solving the
renormalization equations numerically. Similar equations can be derived for
other susceptibilities.

For the FRG boosted DMFT approach, DMF2RG, Taranto  _et al._ (2014) the
initial condition for the full vertex functions and the self-energy are both
given by the DMFT. Taranto  _et al._ (2014) The scale dependent bare
propagator is defined as an interpolation between the DMFT propagator and the
bare lattice propagator as

$\displaystyle G_{\Lambda}^{0}(\mathbf{k},i\omega)^{-1}=\Lambda
G_{imp}^{0}(i\omega)^{-1}$ (55)
$\displaystyle+(1-\Lambda)G_{lattice}^{0}(\mathbf{k},i\omega)^{-1}.$

Since one can treat the functional renormalization group as a vehicle for
summing diagrams, it can in principle be applied on any effective fermionic
interacting system. For example, it has recently been used as a solver for the
dual boson method. Katanin (2019) FRG on auxiliary fermion or dual fermion has
also been considered. Wentzell  _et al._ (2015); Katanin (2015)

The DMF2RG algorithm can be summarized as in the Fig. 16

Figure 16: Algorithm for the functional renormalization group boosted DMFT. 1.
The full reducible vertex and the Green function is obtained from the DMFT
solver. 2. They are fed into the FRG method as the initial conditions. 3.
Solve the flow equations for the FRG to obtain the Green function, full
vertex, and susceptibilities.

##  XIII Numerical methods to represent vertex functions

The various methods described above in principle provide nonlocal corrections
to the conventional dynamical mean field theory or additional nonlocal
corrections in the case of the dynamical cluster approximation or cellular
dynamical mean field theory. The ultimate goal is to attain a better
approximation to the exact solution. The physically most interesting regimes
are often in the intermediate coupling away from the non-trivial limits that
could serve as a good basis for perturbative treatments. After all, all these
methods are based on some truncation with respect to the interaction or more
complicated object as that of the “dual” variables approach. Moreover, the
actual numerical implementation is sometimes rather challenging. Unlike the
DMFT in which only single particle quantities are involved in the self-
consistent equations, the storage requirement of two particle vertex function
is increased by two powers of the space-time grid size. Simply storing those
vertices is in itself a rather difficult ask. Of course, that largely depends
on the interaction and the temperature range. One can naively expect that the
intermediate interaction at low temperature would require more fine resolution
in the space-time grid. For high temperature and at very weak or very strong
coupling, a rather sparse grid for discretizing space-time functions could be
sufficient.

The above storage problem can be resolved to an extent in today’s machines
with tens of thousands of computing nodes. How to efficiently manipulate an
object with such large memory requirement, and involving all-to-all memory
swaps, is an active research problem in computer science. Wagle  _et al._
(2018); Tam  _et al._ (2013) Altogether, with the improvement of computer
hardware and better implementations, the storage problem can be mitigated.

The main issue that is ultimately common to most of these methods is the
existence of a solution. There is in general no guarantee that a self-
consistent solution exists. While this is also true in the single particle
self-consistent method, there is no guarantee that there is a a uniform
convergence to the solution. The situation is more acute in the theories that
require two particle self-consistency. The vertex function is a measure of the
instability towards an ordered phase, therefore it should become singular as
the instability is approached. One can often work in the range of interactions
and temperatures where singularities are far from being reached. However, this
may altogether defeat the purpose of the methods as the regimes near the
instabilities are usually those where corrections to the dynamical mean field
theory are sought.

Even if the solution does exist, the method of searching for the solution may
not be capable of providing a stable solution. This problem was discussed in
the early days when the full parquet formalism for the single impurity
Anderson model was studied. It is fair to say that a universal solution has
yet to be found.

For the DMFT, besides the simple iteration method, the more powerful Broyden
method which utilizes the gradient has been implemented and is occasionally
used in cases where convergence is difficult to achieve. Žitko (2009) For the
two-particle theories, given the complexity of the equations being solved,
more sophisticated methods have been attempted. One of then is the homotopy
method, in which a known convergence solution is relaxed to hopefully lead to
a solution for another temperature or interaction strength. Practically, the
methods are not easy to apply and convergence is not guaranteed.

The above difficulties in solving the two particle vertex function self-
consistently may give an advantage to the FRG based method in which no self-
consistent solution is sought for the method. The solution is obtained not by
solving non-linear integral self-consistent equations, but rather by solving
differential equations with initial conditions using numerical methods. This
allows more flexibility in the numerical solution.

For substantial progress to be accomplished, it is essential for numerical
methods to perform well while requiring reasonable computational resources.
Although storing the full vertex is a daunting task, the amount of information
it actually contains is, in practice, not very large. At least for the weak
coupling case, the vertex functions contain very small entropy in the sense
that it can be compressed numerically to a large extent and still retains most
of the information. It has been suggested for a long time that, a possible
route to represent the vertex function is in the spectral representation.
Shvaika (2016, 2006) A clear advantage in the spectral representation is that
the high frequency information is built in the representation. The spectral
representation has been further explored in recent studies. Shinaoka  _et al._
(2018, 2017)

Various efforts have been devoted to understanding the frequency structure of
the vertex function and these may help with new ideas on approximation schemes
for the vertex function. Kuneš (2011); Katanin (2020); Li  _et al._ (2016a);
Thunström  _et al._ (2018); Wentzell  _et al._ (2020); Rohringer  _et al._
(2012); Chalupa  _et al._ (2018) The latest proposal is to use the tensor
network representation. Shinaoka  _et al._ (2020).

An intuitive scheme is to consider an inhomogeneous grid to represent
momentum-frequency space indices. Generically, the frequency or momentum
dependence is described by an interpolation scheme. This can be justified
specifically for the frequency indices because the low frequency information
should be more important, contains most of the information and thus requires
higher resolution. Moreover the high frequency contribution can be well fitted
by certain functions for convenient storage. Generally, such methods which are
based on interpolations, can be seem as approximating the vertex as follows:

$F(\omega_{k},\omega_{k^{{}^{\prime}}},\omega_{q})=S(\tilde{F}(\tilde{\omega}_{k},\tilde{\omega}_{k^{{}^{\prime}}},\tilde{\omega}_{q}),\omega_{k},\omega_{k^{{}^{\prime}}},\omega_{q}).$
(56)

$F$ is the vertex in the uniform frequency grid $\omega$ and $\tilde{F}$ is
the actual data stored in the grid of some arbitrary basis in
$\tilde{\omega}$. $S$ is the interpolating function or the basis function that
maps $\tilde{F}$ to $F$.

Note that there is no inverse for such an interpolation or basis expansion.
Similar ideas can also be extended to the space grid. This has been studied
recently in the context of FRG and also the parquet method. Eckhardt  _et al._
(2020); Schober  _et al._ (2018); Honerkamp (2018); Lichtenstein  _et al._
(2017); Eckhardt  _et al._ (2018) An extreme case of only retaining the
transfer frequency and momentum has been proposed as a simplification of the
parquet equations for analytical solutions for single impurity problems.
Augustinský and Janiš (2011); Janiš and Augustinský (2007, 2008)

The lack of an inverse transformation to the original vertex means that the
self-consistent equations are, in principle being altered by the
transformation. The representation is however a controlled approximation in
the sense that the larger number of basis functions or the larger number of
grid points can reduce the acquired error. On the other hand, strictly
speaking, the crossing symmetries of the vertex, that are a manifestation of
the Pauli exclusion principle, are also broken by such interpolation or basis
expansion schemes. These symmetries are one of the key features of the parquet
algorithm.

A scheme that has been studied in the context of FRG is to factorize the
vertex function approximately. Taranto  _et al._ (2014) In the representation
of the vertex function with two fermionic frequencies and one bosonic transfer
frequency, the frequency transfer has the most dominant contribution. Thus,
one can argue physically that the frequency or momentum transfer part of the
vertex function can be factorized. This drastically reduces the computational
effort, in particular the storage requirement. A less drastic approximation is
the so-called two-level approximation.Li  _et al._ (2016a)

Another recent proposal is to consider the momentum dependence by expansion in
terms of, say, spherical harmonics. This can help reduce the storage for the
spatial dimension and thus allow studies of larger cluster sizes. Eckhardt
_et al._ (2020)

It is important to point out out that the corrections on top of the DMFT
solution, while they are indeed quite meaningful, are often not large at all,
particularly for weak or very strong interactions. It is, in this sense,
appropriate to inquire what are the effects of numerical errors on a method
that would, in principle, be capable of significant corrections to DMFT. Might
the approximate numerical solutions overwhelm the expected solution? Of
course, one could expect that corrections in the intermediate regime are
larger, but the convergence problem in this regime may prevent a satisfactory
solution.

Overall, a number of challenges still remain for general numerical solutions.
Breakthroughs in the techniques, in the representation of the vertex function
and in the solvers, are needed for general converging solutions over a wider
range of parameters.

##  XIV Exact Constraints

###  XIV.1 Conserving approximation

A milestone in the theory of condensed matter physics is the introduction of
the so-called conserving approximation by Baym and Kadanoff. Baym and Kadanoff
(1961); Baym (1962) The idea of constructing approximation based on a
functional of the Green function was pioneered by Luttinger and Ward,
$\Phi(G)$. Luttinger and Ward (1960) It was shown that the functional
derivative with respect to the Green function can be used to define the self-
energy $\delta\Sigma=tr\Sigma\delta G$. Putting it together with the Dyson
equation $G^{-1}=G_{0}^{-1}-\sigma$, the correlation functions generated by
the functional obey conservation laws and thus the Ward identity.

Conventional methods that are conserving approximations include self-
consistent solutions of the Hatree-Fock, T-matrix approximation, and GW
approximations. Another well known method which is also a conserving
approximation is the fluctuation exchange (FLEX) approximation.

Dynamical mean field theory and its cluster extensions, CDMFT and DCA, are
both conserving approximations although the Ward identity on the lattice is
strictly speaking not satisfied due to the coarse-graining procedure. Though
“conserving approximation” does not necessary imply “better approximation”,
conserving approximations may be more important for broken symmetry cases,
that are not routinely studied by the DMFT.

It has been argued that both the dual fermion and the dual boson approximation
gives conserving approximations for ladder summation. Rubtsov  _et al._ (2012)
That should also be the case for the dynamical vertex approximation, at least
for the ladder approximation.

###  XIV.2 Mermin-Wagner-Coleman theorem

The theorem forbids broken continuous symmetry for a system in less two
dimensions at a finite temperature. Mermin and Wagner (1966); Hohenberg
(1967); Coleman (1973) A particularly interesting case is that of the spin
rotation symmetry in the two dimensional Hubbard model that should not be
broken and thus no true phase transition should occur.

We are not aware of a proof that any of the methods presented in this review
satisfies this constraint. There is plenty of evidence that the DMFT and its
cluster extensions do show finite transition temperatures through measurements
of the susceptibilities.

A recent study by the parquet method seems to provide numerical evidence that
Mermin-Wagner theorem could be satisfied. Eckhardt  _et al._ (2020) The dual
boson theory has also been argued to fulfill the Mermin-Wagner theorem by
enforcing the Hubbard model double occupancy to be the same as that of the
impurity model which is solved numerically in an exact manner. Peters  _et
al._ (2019)

###  XIV.3 Causality

DMFT, CDMFT, and DCA can be shown to be causal. We are not aware of proof of
causality for dual fermion, dual boson, dynamical vertex approximation,
parquet, and FRG Hettler  _et al._ (1998, 2000). However, in practice, that
does not seem to be a serious problem.

###  XIV.4 Crossing symmetry, Pauli principle

An important feature of the parquet method is that the crossing symmetry is
fulfilled. The full FRG also satisfies the crossing symmetries. However, in
practical implementations, the crossing symmetries may sometimes be broken due
to the approximations made on the vertex functions.

##  XV Applications

In this section we discuss the major applications of the multiscale many body
methods discussed above. The list of applications here is not intended to be
exhaustive, in particular we focus exclusively on studies of lattice models.
Efforts to combine density functional theory with these post-DMFT methods are
the subject of current active research. Readers interested in these efforts
will find them discussed in references [Lichtenstein, 2013; Boehnke and
Lechermann, 2012; Lechermann  _et al._ , 2017; Galler  _et al._ , 2019, 2018,
2017; Tomczak  _et al._ , 2017; Held  _et al._ , 2011; Schüler  _et al._ ,
2017; Biermann and Lichtenstein, 2017; Toschi  _et al._ , 2011].

###  XV.1 Dynamical vertex approximation

Most of the applications are focused in particular on the Hubbard model in two
dimensions. The one dimensional Hubbard model has also been studied by the
full parquet dynamical vertex approximation. Valli  _et al._ (2015) For the
three dimensional Hubbard model, it has been found that the antiferromagnetic
phase develops incommensurate magnetic ordering as the doping increases.
Schäfer  _et al._ (2017) Interesting results on the critical exponent have
also been obtained. Rohringer  _et al._ (2011)

The general feature of including more spatial fluctuations is the suppression
of the tendency towards ordering. The dynamical vertex approximation has been
shown in both the two and three dimensions to reduce the transition
temperatures of antiferromagnetic ordering in the half-filled Hubbard model.
As we discussed in the previous section, neither DMFT nor dynamical vertex
approximation fulfill the Mermin-Wagner theorem, thus broken symmetry is
allowed at a finite temperature even for the two-dimensional case. We
summarize the major results from dynamical vertex approximation in the
following.

####  XV.1.1 Hubbard Model

Valli et al. used the parquet dynamical vertex approximation to study the
electronic self-energies and the spectral properties of the finite-size one-
dimensional Hubbard models with periodic boundary conditions for which the
Fermi liquid theory is invalid and that should present a rather challenging
case for any perturbative expansion on top of the DMFT solution.Valli  _et
al._ (2015) They suggested that for a non-degenerate bare dispersion, the
parquet dynamical vertex approximation quantitatively reproduces the exact
many-body solution of the system. Given that the system should be a Luttinger
liquid which cannot be adiabatically tuned into a Fermi liquid, this is a very
encouraging result.

Scháfer et al. have studied the two-dimensional Hubbard model on a square
lattice. Schäfer  _et al._ (2016) They defined two transition lines in the
phase diagram: one for the gap cutting across the nodal direction and the
other for a gap throughout the Fermi surface. The self-energy data shows that
the evolution between the two regimes occurs in a gradual way, not through a
phase transition, and also that at low enough temperatures the whole Fermi
surface is always gapped.

Schäfer showed that the electron self-energy is well separable into a local
dynamical part and static nonlocal contributions for the three dimensional
Hubbard model. Schäfer  _et al._ (2015a) The quasiparticle weight remains
essentially momentum independent for different fillings, including in the
presence of overall large nonlocal corrections to the self-energy.

Pudleiner at al. computed the self-energy for the half-filled Hubbard model on
a square lattice using lattice quantum Monte Carlo simulations and the
dynamical vertex approximation. Pudleiner  _et al._ (2016) The self-energy is
strongly momentum-dependent, but it can be parametrized via the noninteracting
energy-momentum dispersion $\epsilon_{k}$, except for pseudogap features right
at the Fermi edge.

Schäfer et al. studied the two dimensional Hubbard model by combining
dynamical vertex approximation, lattice quantum Monte Carlo, and variational
cluster approximation. Schäfer  _et al._ (2015b) They demonstrated that
scattering at long-range fluctuations due to paramagnons opens a spectral gap
at weak-to-intermediate couplings, irrespective of the preformed localized or
short-ranged magnetic moments. They argued that the two-dimensional Hubbard
model has a paramagnetic phase which is insulating at low enough temperatures
for any finite interaction and no Mott-Hubbard transition is observed.

Schäfer et al. found that the antiferromagnetic phase transition of the
Hubbard model in three dimensions is in contradiction with the conventional
Hertz-Millis-Moriya theory. Schäfer  _et al._ (2017) They argued that the
quantum critical behavior is driven by the Kohn anomalies of the Fermi
surface, even when electronic correlations become strong.

Rohringer at al. studied the three dimensional half-filled Hubard model.
Rohringer  _et al._ (2011) They found the Neel temperature is lowered from
that of the DMFT as expected. More interestingly, they found the critical
exponents to be the same as those of the three dimensional Heisenberg
antiferromagnet in contrast to mean field exponent. This demonstrates that
non-mean-field behavior can indeed be obtained by these systematic nonlocal
corrections.

Rohringer and Toschi studied several spectral and thermodynamic properties of
the Hubbard model in two and three dimensions. Rohringer and Toschi (2016)
Specifically, by evaluating the electronic scattering rate and the
quasiparticle mass renormalization in the low energy regime, they
characterized the gradual deterioration by nonlocal correlations of Fermi
liquid physics as a function of the interaction strength. They found that the
kinetic energy either increases or decreases compared to that of the DMFT
depending on the interaction strength being weak or strong respectively. They
argued that these results correspond to the evolution of the ground state from
a nesting-driven (Slater) to a superexchange-driven (Heisenberg)
antiferromagnet.

####  XV.1.2 Attractive Hubbard Model

Lorenzo Del Re et al. studied the attractive Hubbard model in three dimensions
for the pairing or charge density wave ordering. Del Re  _et al._ (2019) They
found that the fitted critical exponents from the ladder D$\Gamma$A results
were larger not only than the DMFT ones, but also larger than the exact ones
belonging to the corresponding universality class.

####  XV.1.3 Periodic Anderson Model

Schäfer et al. studied the phase diagram and quantum critical region of the
periodic Anderson model. Schäfer  _et al._ (2019) They found a phase
transition between a zero-temperature antiferromagnetic insulator and a Kondo
insulator. In the quantum critical region, they determined a critical exponent
$\gamma=2$ for the antiferromagnetic susceptibility and it becomes $\gamma=1$
at high temperature.

###  XV.2 Dual Fermion

The dual Fermion method has been used to investigate models from the Falicov-
Kimball to the Hubbard model in both two and three dimensions. The Anderson
model for random disorder has also been studied and so far, this is the only
post-DMFT method for the study of disorder. Perhaps a more interesting study
is that of the Anderson-Hubbard model, which investigates the long standing
problem of the competition between Mott insulator and Anderson insulator. We
summarize the major results from dual fermion method in the following.

####  XV.2.1 Hubbard Model

H. Hafermann et al. used the ladder diagrams for the dual fermion and found
that the critical Néel temperature of the mean-field solution is suppressed in
the ladder approximation of two-dimensional Hubbard. Hafermann  _et al._
(2009c)

A. N. Rubtsov at al. found that the antiferromagnetic pseudogap, Fermi-arcs
formation, and non-Fermi-liquid effects due to the Van Hove singularity are
correctly reproduced by the lowest-order diagrams for the two-dimensional
Hubbard model. Rubtsov  _et al._ (2009)

Otsuki et al. obtained the phase diagram for the two-dimensional Hubbard
model. Otsuki  _et al._ (2014) Phase separation is obtained in the low-doping
regime around the Mott insulator.

Astretsov et al. mapped out the phase diagram of the 2D Hubbard model as a
function of temperature and doping. They identified an antiferromagnetic
region at low doping and a superconducting dome at higher doping. Astretsov
_et al._ (2020) Their results support the role of the van Hove singularity as
an important ingredient for the high value of $T_{c}$ at optimal doping. At
small doping, the destruction of antiferromagnetism is accompanied by an
increase of charge fluctuations supporting the scenario of a phase-separated
state driven by quantum critical fluctuations.

Tanaka studied the square-lattice Hubbard model at half-filling using the
ladder dual fermion approximation. He found that the almost simultaneous
creation of the pseudogap and the loss of the Fermi liquid feature is
consistent with what is expected in the Slater regime. Although the pseudogap
still appears in the quasi-particle-like single peak for $U\leq 4.0t$, the
Fermi-liquid feature is partially lost on the Fermi surface already at higher
temperatures as expected in the Mott-Heisenberg regime, where local spins are
preformed at high temperatures. A sharp crossover from a pseudogap phase to a
Mott insulator at finite $U\simeq 4.7t$ was found to occur below the
temperature of the pseudogap formation.

van Loon et al. applied the dual fermion approach with a second-order
approximation to the self-energy for the Mott transition in the two-
dimensional Hubbard model. van Loon  _et al._ (2018a) A strong reduction of
the critical interaction and inversion of the slope of the transition lines
with respect to single-site dynamical mean-field theory was observed.

Katanin et al. showed the suppression of the quasiparticle weight in three-
dimensional Hubbard model. Katanin  _et al._ (2009) With the additional
correction in the susceptibility to fulfill the Mermin-Wagner theorem, Moriya
(2012) they also found a dramatically stronger impact of spin fluctuations in
two dimensions where the pseudogap is formed at low enough temperatures. They
proposed that the origin of the pseudogap at weak-to-intermediate coupling is
in the splitting of the quasiparticle peak.

Hirschmeier et al. studied the three dimensional Hubbard model. They reported
that in the weak-coupling regime, spin-flip excitations across the Fermi
surface are important. The strong-coupling regime is described by Heisenberg
physics. Hirschmeier  _et al._ (2015) In the intermediate-interaction, aspects
of both local and nonlocal correlations appear. They also found that the
critical exponents of the transition in the strong-coupling regime is
consistent with Heisenberg physics down to an interaction of $U=10t$. Again
the non mean-filed exponents is an interesting finding.

Antipov et al. demonstrated that diagrammatic multiscale methods anchored
around local approximations are indeed capable of capturing the non-mean-field
nature of the critical point of lattice models.Antipov  _et al._ (2014) This
is an interesting result as the mean field theory describes the longest length
scale in the problem.

van Loon et al. studied the two-dimensional square-lattice for small to
moderate interaction strengths. van Loon  _et al._ (2018b) The nonlocal
correlations beyond dynamical mean-field theory induce a pseudogap in the
density of states. The upper bounds on the crossover temperature are found to
be significantly lower than previously reported dynamical vertex approximation
results at U=t.

####  XV.2.2 Hubbard Model on Triangular Lattice

Yudin at al. studied the Hubbard model on a triangular lattice. Yudin  _et
al._ (2014) They showed that the band flattening is driven by correlations and
is well pronounced even at sufficiently high temperatures, of the order of
0.1–0.2 times the hopping parameter.

Lee et al. studied the Hubbard model on the triangular lattice at half
filling. They determined the metal-insulator transition and the hysteresis
associated with a first-order transition in the double-occupancy and nearest-
neighbor spin-correlation functions as functions of temperature. Lee  _et al._
(2008) An enhancement of antiferromagnetic correlations and evidence for
magnetically ordered phases by calculating the spin susceptibility were found.

Antipov et al. studied the half-filled Hubbard model on an isotropic
triangular lattice with a spin polarized extension of the dual fermion
approach. Antipov  _et al._ (2011) They found that the dual fermion
corrections drastically decrease the energy of a spin liquid state while
leaving the non-collinear magnetic states almost non-affected. This makes the
spin liquid become a preferable state in a certain interval of interaction
strength of the order of the magnitude of the bandwidth.

Li et al, studied both the half-filled and the doped Hubbard model on a
triangular lattice. Li  _et al._ (2016b) The phase diagrams are obtained.

####  XV.2.3 Hubbard Model on Honeycomb Lattice

Hirschmeier et al. studied the Hubbard model on the honeycomb lattice in the
vicinity of the quantum critical point by means of a multiband formulation of
the dual fermion approach. They found that the critical interaction strength
of the quantum phase transition from a paramagnetic semimetal to an
antiferromagnetic insulator is in good agreement with other numerical methods.
Hirschmeier  _et al._ (2018) They also argued that the Hubbard model on the
honeycomb lattice behaves like a quantum nonlinear $\sigma$ model, while
displaying signs of non-Fermi-liquid behavior.

####  XV.2.4 Falicov-Kimball Model

Astleithner at al. studied the Falicov-Kimball model. Using the full parquet
dynamical vertex approximation, they argued that weak localization corrections
in the particle-particle channel are not the dominant vertex corrections to
the optical conductivity. Astleithner  _et al._ (2020)

####  XV.2.5 Kondo Lattice Model

Otsuki studied the Kondo lattice model to explore possible superconductivity
emerging from the critical antiferromagnetic fluctuations. Otsuki (2015) The
d-wave pairing is found to be the leading instability only in the weak-
coupling regime. As the coupling is increased, a change of the pairing
symmetry into a p-wave spin-singlet pairing was found.

####  XV.2.6 s-d Exchange Model

Sweep et al. studied the critical values of the s-d exchange coupling
constant. Sweep  _et al._ (2013) They reported a difference between the DMFT
and dual fermion results that is more than a factor of two for the square
lattice and spin one-half localized electrons.

####  XV.2.7 Anderson Disorder Model

Terletska et al. generalized the he dual fermion approach to disordered
systems using the replica method. Terletska  _et al._ (2013); Yang  _et al._
(2013) The developed method utilizes the exact mapping to the dual fermion
variables, and includes inter-site scattering via diagrammatic perturbation
theory in the dual variables. nonlocal effects missed in the CPA are captured.

####  XV.2.8 Anderson Falicov-Kimball Model

Yang at al. generalized the dual-fermion formalism for disordered fermionic
systems to include the effect of interactions. The phase diagram for the two
dimensional Anderson-Falicov-Kimball model was obtained. Yang  _et al._ (2014)

####  XV.2.9 Anderson-Hubbard Model

Haase et al. studied the three-dimensional Anderson Hubbard model. They report
that the dual-fermion approach leads to quantitative as well as qualitative
improvement of the dynamical mean-field results, and it allows one to
calculate the hysteresis in the double occupancy in three dimensions, taking
into account nonlocal correlations. Haase  _et al._ (2017)

Otsuki studied the Kondo lattice model. Otsuki (2015) and found that different
superconductivity pairing symmetries emerge from the critical
antiferromagnetic fluctuations. The d-wave pairing is found to be the leading
instability only in the weak-coupling regime. As the coupling is increased, a
change of the pairing symmetry into a p-wave spin-singlet pairing is observed.
The competing superconductivities are ascribed to a crossover between small
and large Fermi surfaces, which occurs with the formation of heavy
quasiparticles.

###  XV.3 Dual Boson

The dual boson method has been applied on models with nonlocal interaction.
These include the truncated long range Coulomb coupling, the nearest neighbor
interaction in the extended Hubbard model, and the anisotropic dipolar
coupling for cold atoms. A key difference is that those nonlocal density-
density type couplings can lead to the competition between charge fluctuations
and spin fluctuations which results in charge density wave ordering and
possibly bond order wave ordering. Nakamura (2000)

####  XV.3.1 Extended Hubbard Model

Vandelli et al. proposed to use quantum Monte Carlo to sample the diagrams
from the dual boson theory for the extended Hubbard model. Vandelli  _et al._
(2020) They proposed that the single-particle Green function allows one to
estimate the transition point to the charge density wave phase.

####  XV.3.2 Hubbard model with Dipolar Coupling

van Loon et al. studied the Hubbard model with long-range dipole-dipole
interactions. van Loon  _et al._ (2015) This is an interesting model in the
context of experiments with cold atoms on optical lattice. Lu  _et al._ (2012)
Besides the stripe phase and the checkerboard phase, based on their dual boson
calculation they suggest that there is a novel phase with “ultralong-rang”
density correlations at distances of tens of lattice sites.

####  XV.3.3 Hubbard model with Coulomb coupling

Hafermann et al. and van Loon et al. studied the polarization for the two
dimensional Hubbard model with the long range Coulomb coupling. Van Loon  _et
al._ (2014); Hafermann  _et al._ (2014) They found that plasmon spectra are
qualitatively different from those of the random-phase approximation: they
exhibit a spectral density transfer and a renormalized dispersion with
enhanced deviation from the canonical $\sqrt{q}$ behavior.

###  XV.4 TRILEX

Applications of TRILEX are mostly on the two dimensional Hubbard model.

####  XV.4.1 Hubbard Model

Aryal et al. found that the local vertex acquires a strong frequency
dependence, driving the system to a Mott transition for the half-filled
Hubbard model on a square lattice. Ayral and Parcollet (2015) At low enough
temperatures the momentum-dependence of the self-energy is enhanced due to
large spin fluctuations. Upon doping, they find a Fermi arc in the spectral
function.

Vučičević et al. studied the dependence of the superconducting temperature on
the bare dispersion at weak coupling, which shows a clear link between strong
antiferromagnetic correlations and the onset of superconductivity. Vučičević
_et al._ (2017) combination of hopping amplitudes particularly favorable to
superconductivity at intermediate doping.

###  XV.5 Parquet

The parquet formalism is quite demanding from the point of view of
computational resources and, for this reason, it has not been as extensively
pursued as other diagrammatic methods despite its promise. We highlight two
earlier studies using this method. Recent studies include applying it in the
dual fermion method and dynamical vertex approximation or simplifying the
calculation by different representations of the vertex functions as discussed
in previous subsections.

####  XV.5.1 Hubbard model

S. X. Yang et al. solved the full parquet equation for the half-filled 4x4
cluster. Yang  _et al._ (2009) Results were compared to those obtained from
Determinant Quantum Monte Carlo (DQMC), Fluctuation Exchange (FLEX), and self-
consistent second-order approximation methods. This comparison shows a
satisfactory agreement with DQMC and a significant improvement over the FLEX
or the self-consistent second-order approximation.

Kusunose solved parquet equations for both the impurity Anderson model and the
Hubbard model on a square lattice mainly for the particle-hole symmetric or
the half-filled case. Kusunose (2010) He argued that in both models the vertex
renormalization in the spin channel eliminates magnetic instabilities of mean-
field theory to ensure the Mermin-Wagner theorem. The parquet method gives the
same critical exponents as the self-consistent renormalization Moriya theory
in the quantum critical region.

Pudleiner et al. studied the Pariser-Parr-Pople (PPP) model or Hubbard model
with nonlocal interaction for the conjugated $\pi$ bonds in benzene. Pudleiner
_et al._ (2019) They found that quasiparticle renormalization is much weaker
in the PPP than in the Hubbard model, but the static part of the self-energy
enhances the band gap of the PPP model. In additional, the vertex corrections
to the optical conductivity are much more important in the PPP model.

###  XV.6 FRG

Functional renormalization group has been studied extensively for over two
decades. Recent applications on using it to improve the DMFT solution has so
far been mostly limited to the Hubbard model.

####  XV.6.1 Hubbard model

Tranto et al. first proposed to use FRG to expand the DMFT solution. Taranto
_et al._ (2014) They studied the half-filled square lattice Hubbard model and
found the method provides more prominent momentum dependence than the
conventional FRG method.

Vilardi et al. studied the doped two dimensional Hubbard model. Vilardi  _et
al._ (2019) They found strong antiferromagnetic correlations from half-filling
to 18% hole doping at low temperature and a sizable d-wave pairing interaction
driven by magnetic correlations at the edge of the antiferromagnetic regime.

##  XVI Conclusion

We have reviewed multiscale many body numerical methods to address strongly
correlated systems by appropriately treating the short length scale, the long
length scale and the intermediate length scale. The different methods
implemented to date have produced promising results despite being hindered by
the a variety of numerical challenges. Since short length scales are treated
exactly, diagrammatic methods arise as the suitable method to deal with the
intermediate length scales with a different precision level by systematically