# Causal versus local $GW$+EDMFT scheme and application to the triangular-
lattice extended Hubbard model

Jiyu Chen  Department of Physics, University of Fribourg, 1700 Fribourg,
Switzerland  Francesco Petocchi  Department of Physics, University of
Fribourg, 1700 Fribourg, Switzerland  Philipp Werner  Department of Physics,
University of Fribourg, 1700 Fribourg, Switzerland

###### Abstract

Using the triangular-lattice extended Hubbard model as a test system, we
compare $GW$+EDMFT results for the recently proposed self-consistency scheme
with causal auxiliary fields to those obtained from the standard
implementation which identifies the impurity Green’s functions with the
corresponding local lattice Green’s functions. Both for short-ranged and long-
ranged interactions we find similar results, but the causal scheme yields
slightly stronger correlation effects at half-filling. We use the two
implementations of $GW$+EDMFT to compute spectral functions and dynamically
screened interactions in the parameter regime relevant for 1$T$-TaS2. We
address the question whether or not the sample used in a recent photoemission
study [Phys. Rev. Lett. 120, 166401 (2018)] was half-filled or hole-doped.

###### pacs:

71.10.Fd

##  I Introduction

Strongly correlated electron materials exhibit remarkable phenomena such as
quasi-particles with strongly renormalized masses and correlation-induced
metal-insulator (Mott) transitions.Mott (1968); Imada et al. (1998) The
screening of the Coulomb interaction determines the effective electron-
electron interaction in solids and thus plays an important role in the
modeling of this class of materials.Hedin (1999); Werner and Casula (2016) The
dynamically screened interaction has to be computed self-consistently, because
the correlated electronic structure determines the relevant low-energy
screening modes, while the interaction determines band renormalizations,
quasi-particle lifetimes and Mott gaps. A simple model which allows to study
the interplay between correlations and screening is the single-band extended
Hubbard model with local and nonlocal (possibly long-ranged Hansmann et al.
(2013)) interactions.Ayral et al. (2012); Huang et al. (2014)

Here, we use a combination of the Extended Dynamical Mean-Field Theory (EDMFT)
method Sun and Kotliar (2002) and the $GW$ method ($GW$+EDMFT) Biermann et al.
(2003); Aryasetiawan et al. (2004); Biermann et al. (2004) to study the phase
diagram and correlation functions of the single-band Extended Hubbard model on
the two-dimensional triangular lattice. EDMFT provides a local self-energy and
polarization, and allows to capture strong correlation effects and Mott
physics, while $GW$ provides nonlocal self-energy and polarization components
which capture nonlocal correlation and screening effects. The latter have been
shown to be relevant for an accurate description of materials.Petocchi et al.
(2021) In order to treat long-ranged interactions, and compare the results for
models with short-ranged and long-ranged nonlocal interactions, we implement
the Ewald lattice summation.Ewald (1921); Harris (1998); Hansmann et al.
(2013)

One purpose of this study is to systematically test the recently proposed
manifestly causal $GW$+EDMFT scheme of Backes, Sim and Biermann.Backes et al.
(2020) The mapping of the lattice model to a single-site impurity model yields
two dynamical mean fields: a fermionic mean field (the hybridization
function), and a bosonic mean field (the effective local interaction). In the
conventional $GW$+EDMFT implementation, which accomplishes the mapping by
identifying the local lattice Green’s function and local screened interaction
with the corresponding impurity quantities, these auxiliary fields are not
necessarily causal. This is not in principle a problem, as long as all the
physical quantities (e.g. the screened interaction or the Green’s function)
are causal,Nilsson et al. (2017) but it can create numerical issues. Such
problems are avoided by the modified self-consistency scheme of Ref. Backes et
al., 2020, which by construction yields causal auxiliary fields. To clarify
the influence of the implementation on the auxiliary fields and physical
observables, we present a systematic comparison between the $GW$+EDMFT results
obtained with the conventional (“local”) and the new “causal” self-consistency
scheme for the extended Hubbard model.

Our second purpose is to provide reference data for strongly correlated
electron systems on a triangular lattice, such as 1$T$-TaS2. We present the
$GW$+EDMFT phase diagram of the triangular lattice model as a function of the
onsite interaction $U$ and nonlocal interaction $V$, and determine the
parameters appropriate for single-band extended Hubbard simulations of
1$T$-TaS2. We then use these parameters to study the spectral functions and
dynamically screened interactions of half-filled and chemically doped systems
in the temperature region corresponding to the C phase ($T\lesssim 180$ KWang
et al. (2020)), and compare to the photoemission spectrum from Ref. Ligges et
al., 2018. The latter study concluded, based on the dynamics of photo-doped
doublons, that their nominally half-filled sample may have been hole doped.
Our comparison to the occupation functions for different electron fillings
however demonstrates a good match between the experimental equilibrium
spectrum and the result for the undoped model.

The paper is organized as follows. In Sec. II.1 we describe the model and in
Sec. II.4 the $GW$+EDMFT method, as well as the two self-consistency schemes.
In Sec. III, after presenting the $GW$+EDMFT phase-diagram of the triangular
lattice extended Hubbard model (III.2), we show the comparison of the two
schemes in different relevant parameter regimes (III.3). We then present the
results representative of 1$T$-TaS2 in Sec. III.4 and compare our results to
the equilibrium spectrum from Ref. Ligges et al., 2018. Section IV contains
the conclusions.

##  II Model and method

###  II.1 Model

We consider the single-band Extended ($U$-$V$) Hubbard model on the two-
dimensional triangular lattice. The grand-canonical Hamiltonian is given by

$\displaystyle H_{\text{Hubbard}}=$
$\displaystyle-\sum_{ij\sigma}t_{ij}c^{\dagger}_{i\sigma}c_{j\sigma}+U\sum_{i}n_{i\uparrow}n_{i\downarrow}$
(1) $\displaystyle+\frac{1}{2}\sum_{ij}V_{ij}n_{i}n_{j}-\mu\sum_{i}n_{i},$

where the operators $\hat{c}_{i\sigma}$ and $\hat{c}^{\dagger}_{i\sigma}$
denote the annihilator and the creator of an electron with spin
$\sigma\in\\{\uparrow,\downarrow\\}$ at the $i$-th lattice site.
$n_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$ and
$n_{i}=n_{i\uparrow}+n_{i\downarrow}$ are the particle number operators for
site $i$. We will restrict the hoppings to nearest neighbor (NN) sites
$\langle ij\rangle$, and write $t_{ij}=t\delta_{\langle ij\rangle}$. $U$ is
the on-site repulsive interaction between two electrons with opposite spins on
the same site, while $V_{ij}$ is the inter-site repulsive interaction, for
which we choose the Coulomb form $V_{ij}\propto 1/r_{ij}$. Here, $r_{ij}$
denotes the physical distance between the sites $i$ and $j$. $\mu$ is the
chemical potential for the grand-canonical ensemble.

We will use two approximations for the nonlocal interaction, to get some
insights into the effects of the range of the interaction. On the one hand, we
truncate the interaction to the NN components, $V_{ij}=V_{01}\delta_{\langle
ij\rangle}$. On the other hand, we will treat the long-range Coulomb
interaction via the Ewald lattice summation, Hansmann et al. (2013) as
explained below. This will define the second approximation, denoted as
“Ewald.”

###  II.2 Ewald lattice summation

It is convenient to Fourier transform the Hamiltonian from real space to
momentum space ($k$-space). For our two-dimensional homogeneous lattice model,
this transformation yields
$A(\mathbf{k})=\sum_{i}A_{0i}e^{i\mathbf{k}\cdot\mathbf{r}_{0i}}$, where
$A_{ij}$ stands for the hopping $t_{ij}$ or the Coulomb repulsion $V_{ij}$.
For the triangular lattice, in the translationally invariant NN case, using
$A_{ij}\equiv A_{01}\delta_{\langle ij\rangle}$, the band dispersion becomes
$\epsilon_{\mathbf{k}}=-2t(\cos
k_{x}+\cos(-\frac{1}{2}k_{x}+\frac{\sqrt{3}}{2}k_{y})+\cos(\frac{1}{2}k_{x}+\frac{\sqrt{3}}{2}k_{y}))$
and the bare interaction $v_{\mathbf{k}}=U+V_{\mathbf{k}}$ with
$V_{\mathbf{k}}=2V_{1}(\cos
k_{x}+\cos(-\frac{1}{2}k_{x}+\frac{\sqrt{3}}{2}k_{y})+\cos(\frac{1}{2}k_{x}+\frac{\sqrt{3}}{2}k_{y}))$.
In addition to this, we also consider long-range interactions using the Ewald
lattice summation, as was previously done in Ref. Hansmann et al., 2013. More
specifically, we will assume an infinite-range Coulomb repulsion,
$V_{ij}=V_{1}/r_{ij}$, which scales as one over the distance between the
sites, and which is parametrized by the nearest-neighbor interaction $V_{1}$.
The Ewald lattice summation method Ewald (1921); Harris (1998) provides a way
of regularizing the interaction energy and avoiding problems of conditional
convergence. To this end, the long-range Coulomb repulsion is represented as
the sum of two parts, corresponding to a real-space summation and a
reciprocal-space summation.Harris (1998) First, we re-express the Coulomb
potential in terms of a Gaussian integral. Second, we divide the Gaussian
integral into two parts, by means of a parameter $\eta$ which fixes the
partitioning between short-range and long-range components. Third, the
summation in the two parts is performed in real space and reciprocal space,
respectively. Details are presented in Appendix A. The final expression for
the Ewald long-range interaction becomes

$V_{ij}=V_{1}\Bigg{\\{}\frac{1}{r_{ij}}\operatorname{erfc}\left(\frac{r_{ij}}{2\eta}\right)+\frac{1}{N_{\mathbf{k}}}\sum_{\mathbf{k}\neq
0}e^{i\mathbf{k}\cdot\mathbf{r}_{ij}}\frac{2\pi}{k}\operatorname{erfc}\left(k\eta\right)\Bigg{\\}}.$
(2)

Figure 1: Comparison of the interaction $V(q)$ in the first Brillouin zone for
three different potentials, with the same NN interaction $V_{1}=1$ and lattice
spacing $a=1$. ($\Gamma$, M, K denote the centre, the midpoint of edge and the
corner of the hexagonal first Brillouin zone, respectively.)

In Fig. 1 we plot $V(k)$ along the indicated path in the Brillouin zone (BZ).
The green curve shows the interaction obtained by truncating the Coulomb
interaction at the 32nd neighbors. The orange line plots the Ewald interaction
(2) for $\eta=16$. In our implementation, we also use a cutoff at the 32nd
neighbors in the Ewald case, which explains the negative values in some
regions of $k$ space. For comparison, we also show a Yukawa potential
$V_{\text{Yukawa}}(r)=\frac{V_{1}}{r}e^{-r/\alpha}$ for $\alpha=16$. One sees
that the regularization with the $\eta$ parameter in the Ewald summation has a
similar effect as a cutoff parameter $\alpha=\eta$ in the Yukawa case,
smoothing the potential and regularizing the divergence at the $\Gamma$ point.

###  II.3 Action

In order to compute the Green’s function and other physical quantities, it is
convenient to reformulate the model in terms of the action using the coherent-
state path-integral representation.Negele and Orland (1988) For an arbitrary
Hamiltonian $H(c^{\dagger}_{\alpha},c_{\alpha})$, where $c^{\dagger}_{\alpha}$
and $c_{\alpha}$ are fermionic creators and annihilators indexed by $\alpha$,
the grand-canonical partition function $Z=\operatorname{Tr}e^{-\beta H}$ can
be expressed in the path-integral form as
$Z=\int_{\phi(\beta)=-\phi(0)}\mathcal{D}[\phi_{\alpha}^{*}(\tau),\phi_{\alpha}(\tau)]e^{-S[\phi^{*},\phi]}$
with action

$\displaystyle S[\phi^{*},\phi]\equiv$
$\displaystyle\int_{0}^{\beta}d\tau\bigg{\\{}\sum_{\alpha}\phi_{\alpha}^{*}(\tau)(\partial_{\tau}-\mu)\phi_{\alpha}(\tau)$
(3) $\displaystyle+H(\phi_{\alpha}^{*}(\tau),\phi_{\alpha}(\tau))\bigg{\\}}.$

Here, $\phi^{*}_{\alpha}$ and $\phi_{\alpha}$ denote the corresponding anti-
periodic Grassmann numbers satisfying
$\phi_{\alpha}(\tau+\beta)=-\phi_{\alpha}(\tau)$. Specifically, for the
Hamiltonian of the extended Hubbard model defined in Eq. (1), the action is
$S[c^{*},c]$
$=\int_{0}^{\beta}d\tau\\{\sum_{i,\sigma}c_{i\sigma}^{*}(\tau)(\partial_{\tau}-\mu)c_{i\sigma}(\tau)-\sum_{ij\sigma}t_{ij}c_{i\sigma}^{*}(\tau)c_{j\sigma}(\tau)$
$+U\sum_{i}n_{i\uparrow}(\tau)n_{i\downarrow}(\tau)+\sum_{i<j}V_{ij}n_{i}(\tau)n_{j}(\tau)\\}.$
By applying the identity $n_{i}n_{i}=n_{i}+2n_{i\uparrow}n_{i\downarrow}$ and
introducing $v_{ij}=U\delta_{ij}+V_{ij}$ and $\tilde{\mu}=\mu+U/2$, we can
write it in the more compact form

$\displaystyle S\left[c^{*},c\right]=$
$\displaystyle\int_{0}^{\beta}d\tau\bigg{\\{}\sum_{ij\sigma}c_{i\sigma}^{*}(\tau)\left[\left(\partial_{\tau}-\tilde{\mu}\right)\delta_{ij}-t_{ij}\right]c_{j\sigma}(\tau)$
$\displaystyle+\frac{1}{2}\sum_{ij}v_{ij}n_{i}(\tau)n_{j}(\tau)\bigg{\\}}.$
(4)

###  II.4 $GW$+EDMFT

####  II.4.1 General remarks

Within EDMFT the lattice model is mapped to an auxiliary single-site impurity
model with two dynamical mean fields: the fermionic bath Green’s function
$\mathcal{G}^{-1}=i\omega_{n}+\mu-\Delta$, which contains information on how
electrons are connected to the bath through the hybridization function
$\Delta$, and the bosonic field $\mathcal{U}$, which is a frequency-dependent
effective impurity interaction.Sun and Kotliar (2002); Ayral et al. (2013)
This scheme identifies the local projection of the lattice self-energy
$\Sigma$ and polarization $\Pi$ with the corresponding fields provided by an
effective impurity model, which is solved in a self-consistent manner.

By adding the nonlocal self-energy and polarization components obtained with
the $GW$ methodHedin (1965) to the EDMFT local solution, the
$GW$+EDMFTBiermann et al. (2003) scheme is able to simultaneously treat
screening effects induced by local and nonlocal charge fluctuations and strong
local correlations. Even if short-range correlations and frustration effects
produced by the lattice geometry are not fully captured, $GW$+EDMFT in
realistic materials contexts has been shown to yield accurate results for the
effective interaction strengths in solids.Petocchi et al. (2020a); Petocchi et
al. (2021)

The standard (“local”) self-consistency scheme of $GW$+DMFT fixes the
auxiliary fields by identifying the local self-energy and local polarization
with the impurity counterparts.Biermann et al. (2003); Ayral et al. (2013)
This can lead to noncausal dynamical mean fields, which is by itself not a
problem (since the physical observables such as the Green’s function and
screened interaction are causal),Nilsson et al. (2017) but it can potentially
lead to numerical issues, such as sign problems in Monte-Carlo based impurity
solvers, or incompatibilities with standard analytical continuation procedures
that assume causality.

In Ref. Backes et al., 2020, a modified $GW$+EDMFT self-consistency scheme has
been proposed, which by construction produces causal auxiliary fields. We
sketch below the local and the new causal $GW$+EDMFT self-consistency loops,
while a more detailed derivation of the causal self-consistent equations can
be found in Appendix B.

####  II.4.2 EDMFT

Even though the causal self-consistency loop differs from the local one only
in the presence of nonlocal self-energy and polarization components, we will
lay out the equations in this section in a way that is compatible with the new
scheme of Backes et al.Backes et al. (2020) The EDMFT self-consistent
equations map the lattice action defined in Eq. (II.3) to the impurity action

$\displaystyle S^{\text{EDMFT}}=$
$\displaystyle-\int_{0}^{\beta}\sum_{\sigma}c^{\dagger}_{\sigma}(\tau)\mathcal{G}^{-1}(\tau-\tau^{\prime})c_{\sigma}(\tau^{\prime})$
$\displaystyle+\frac{1}{2}\int_{0}^{\beta}d\tau
d\tau^{\prime}n(\tau)\mathcal{U}(\tau-\tau^{\prime})n(\tau^{\prime}),$ (5)

which contains the (fermionic) Weiss field $\mathcal{G}$ and the (bosonic)
effective interaction $\mathcal{U}$ (see Fig. 2).

This mapping involves approximations in the case of finite-dimensional
lattices,Georges et al. (1996) but it becomes exact in the large-connectivity
limit,Si and Smith (1996) i.e. for coordination number $z\to\infty$, if one
properly renormalizes the hoppings ($t\to t/\sqrt{z}$) and interactions ($V\to
V/z$). As noted in Ref. Sun and Kotliar, 2002, there are in principle two ways
to formulate the self-consistency loop, either using the charge density $n$
and the density-density correlation function $\chi$, or (after a Hubbard-
Stratonovich transformation of the lattice action Eq. (II.3)) using the
auxiliary bosonic field $\phi$ and bosonic Green’s function $W$. For the local
EDMFT self-consistent equations the two formulations are equivalent,Sun and
Kotliar (2002) while this equivalence may not hold any more for the causal
EDMFT equations. In this work, following Ref. Backes et al., 2020, we use the
formulation involving $n$ and $\chi$. The result of the generalized self-
consistency equations (see Appendix B for a detailed derivation) can be
expressed in terms of corrections to the impurity Weiss fields:Backes et al.
(2020)

$\displaystyle\mathcal{G}^{-1}\left(i\omega_{n}\right)$
$\displaystyle=\quad\langle
G\rangle^{-1}_{k}+\langle\Sigma\rangle_{k}+\Delta_{\text{cor}},$ (6)
$\displaystyle\mathcal{U}(i\nu_{n})$ $\displaystyle=\left[\langle
W\rangle^{-1}_{k}+\langle\Pi\rangle_{k}\right]^{-1}+\mathcal{U}_{\text{cor}},$

where $G$ and $W$ denote the lattice Green’s function and screened
interaction, respectively, while $\Sigma$ is the self-energy and $\Pi$ the
polarization. $\langle\rangle_{k}$ denotes an average over momentum in the BZ.
The correction terms of the casual self-consistency scheme are given byBackes
et al. (2020)

$\displaystyle\Delta_{\text{cor}}$ $\displaystyle=-\langle\Sigma
G\Sigma\rangle_{k}+\langle\Sigma G\rangle_{k}\langle G\rangle_{k}^{-1}\langle
G\Sigma\rangle_{k}-2\langle\Sigma\rangle_{k}+\langle\Sigma G\rangle_{k}\langle
G\rangle_{k}^{-1}+\langle G\rangle_{k}^{-1}\langle G\Sigma\rangle_{k},$ (7)
$\displaystyle\mathcal{U}_{\text{cor}}$ $\displaystyle=-\langle\Pi
W\rangle_{k}[\langle\Pi\rangle_{k}+\langle\Pi W\Pi\rangle_{k}]^{-1}\langle
W\Pi\rangle_{k}+\langle\Pi\rangle_{k}\langle
W\rangle_{k}[\langle\Pi\rangle_{k}+\langle\Pi\rangle_{k}\langle
W\rangle_{k}\langle\Pi\rangle_{k}]^{-1}\langle
W\rangle_{k}\langle\Pi\rangle_{k}.$

Note that momentum-independent $\Sigma$ and $\Pi$ yield by construction
vanishing corrections $\Delta_{\text{cor}}=0$ and
$\mathcal{U}_{\text{cor}}=0$. Therefore, the EDMFT self-consistency scheme,
which assumes $\Sigma=\Sigma_{\text{imp}}$, $\Pi=\Pi_{\text{imp}}$ is not
altered.Ayral et al. (2013) The correction terms however have a nontrivial
effect if $\Sigma$ and $\Pi$ are $k$-dependent.

Figure 2: Left: The extended Hubbard model on the triangular lattice, with NN
hopping $t$, on-site interaction $U$ and nonlocal interactions $V$. Right:
EDMFT maps the lattice model to an impurity model defined by a hybridization
function $\Delta$ (or fermionic Weiss field
$\mathcal{G}^{-1}=i\omega_{n}+\mu-\Delta$) and a bosonic field (retarded
interaction) $\mathcal{U}$.

####  II.4.3 $GW$ method

In this subsection, we discuss the $GW$ method,Aryasetiawan and Gunnarsson
(1998); Sun and Kotliar (2002); Biermann et al. (2003) which allows to
construct a $k$-dependent self-energy $\Sigma^{GW}(k,i\omega)$ and
polarization $\Pi^{GG}(k,i\omega)$ by truncating Hedin’s vertex function to
its first order.Hedin (1965); Pavarini et al. (2011) This is a weak-coupling
approach that can be expected to work for weak or moderate interactions, or
(in combination with EDMFT) for the less dominant nonlocal components in more
strongly interacting systems. Hedin derived a set of exact coupled
differential equations which link the Green’s function $G$, self-energy
$\Sigma$, vertex function $\Gamma$, screened interaction $W$, and polarization
$\Pi$:Hedin (1965)

$\left\\{\begin{aligned} \Sigma(1,2)&=-\int
d(34)G\left(1,3^{+}\right)W(1,4)\Gamma(3,2,4),\\\ G(1,2)&=G_{0}(1,2)+\int
d(34)G_{0}(1,3)\Sigma(3,4)G(4,2),\\\ \Gamma(1,2,3)&=\delta(1-2)\delta(2-3)\\\
&+\int d(4567)\frac{\partial\Sigma(1,2)}{\partial
G(4,5)}G(4,6)G(7,5)\Gamma(6,7,3),\\\ \Pi(1,2)&=\int
d(34)G(1,3)\Gamma(3,4,2)G\left(4,1^{+}\right),\\\ W(1,2)&=v(1,2)+\int
d(34)v(1,3)\Pi(3,4)W(4,2),\end{aligned}\right.$ (8)

Here, the numbers stand for $i\equiv(\tau_{i},\mathbf{r}_{i})$ and the spin
indices are omitted. The second and the fifth equations are the fermionic and
bosonic lattice Dyson equations. The $GW$ approximation neglects the vertex
corrections in the third equation by setting
$\Gamma(1,2,3)\approx\delta(1-2)\delta(2-3)$, even though this violates the
Pauli principle. Pavarini et al. (2011) With this approximation the first and
the forth equations simplify to

$\left\\{\begin{aligned} \Sigma(1,2)&=-G(1,2^{+})W(1,2),\\\
\Pi(1,2)&=G(1,2)G(2,1^{+}),\end{aligned}\right.$ (9)

or in the notation with site indices and imaginary time,

$\left\\{\begin{aligned} \Sigma^{GW}_{ij}(\tau)&=-G_{ij}(\tau)W_{ij}(\tau),\\\
\Pi^{GG}_{ij}(\tau)&=2G_{ij}(\tau)G_{ji}(-\tau).\end{aligned}\right.$ (10)

Here the factor of two comes from the spin summation. In momentum space, we
have

$\displaystyle\Sigma^{GW}(k,\tau)$ $\displaystyle=-\frac{1}{S_{\text{BZ}}}\int
dpG(p,\tau)W(k-p,\tau)=-G*W,$ (11) $\displaystyle\Pi^{GG}(k,\tau)$
$\displaystyle=\frac{2}{S_{\text{BZ}}}\int
d{p}G({p},\tau)G({p}-{k},-\tau)=2G*G,$

where “$*$” denotes a momentum convolution, and $S_{\text{BZ}}$ the integral
over the first BZ.

In 1960 Luttinger and Ward Luttinger and Ward (1960) proposed an explicit rule
to construct a functional $\Psi[G]$, from which the self-energy $\Sigma$ can
be obtained by a functional derivative. Later, AlmbladhAlmbladh et al. (1999)
extended this “LW functional” to $\Psi[G,W]$ in a formalism which includes the
polarization $\Pi$ and the self-energy, and which yields the following
variational relations,

$\left(\frac{\delta\Psi[G,W]}{\delta
G}\right)_{W}=\Sigma,\quad\left(\frac{\delta\Psi[G,W]}{\delta
W}\right)_{G}=-\frac{1}{2}\Pi.$ (12)

Looking at Eq. (9), one finds that the $GW$ approximation corresponds to the
functional $\Psi^{GW}[G,W]=-\frac{1}{2}\operatorname{Tr}[GWG]$, which is the
lowest order approximation in $W$.Nilsson et al. (2017)

####  II.4.4 $GW$+EDMFT schemes

As mentioned before, EDMFT and $GW$ are complementary methods, in the sense
that the former allows to capture strong local correlations and the latter
weak nonlocal correlations. It is thus natural to split the functional
$\Psi[G,W]$ into two terms, such that EDMFT contributes to its local part and
the $GW$ scheme contributes to its non-local part,

$\Psi\approx\Psi^{\mathrm{EDMFT}}\left[G_{ii},W_{ii}\right]+\Psi_{\text{nonloc
}}^{GW}\left[G_{ij},W_{ij}\right].$ (13)

The local contribution of the $GW$ approximation is subtracted to avoid double
counting,

$\Psi_{\text{nonloc
}}^{GW}\left[G_{ij},W_{ij}\right]=\Psi^{GW}\left[G_{ij},W_{ij}\right]-\Psi^{GW}\left[G_{ii},W_{ii}\right].$
(14)

Applying the relations $\Sigma=\frac{\delta\Psi}{\delta G}$ and
$\Pi=-2\frac{\delta\Psi}{\delta W}$, we then obtain

$\left\\{\begin{aligned}
\Sigma_{ij}=\Sigma^{\operatorname{EDMFT}}_{ii}\delta_{ij}+\Sigma^{GW}_{ij}(1-\delta_{ij}),\\\
\Pi_{ij}=\Pi^{\operatorname{EDMFT}}_{ii}\delta_{ij}+\Pi^{GG}_{ij}(1-\delta_{ij}).\end{aligned}\right.$
(15)

Substituting this non-local self-energy and polarization into the self-
consistency equations (6) and (7), we can solve the lattice model in a self-
consistent way. Ignoring the correction terms $\Delta_{\text{cor}}$ and
$\mathcal{U}_{\text{cor}}$ in Eq. (6) yields the standard (local) $GW$+EDMFT
scheme, while including these terms corresponds to the causal $GW$+EDMFT
scheme proposed in Ref. Backes et al., 2020.

In the following, we summarize the two $GW$+EDMFT self-consistency loops, see
also Fig. 3. The terms in brackets are added in the causal variant.

  1. 1.

Start from some initial guess for the lattice self-energy
$\Sigma(k,i\omega_{n})$ and the lattice polarization $\Pi(k,i\nu_{n})$.

  2. 2.

With the lattice Dyson equation, calculate the lattice Green’s function $G$
and screened interaction $W$,

$\displaystyle
G^{-1}(k,i\omega_{n})=G_{0}^{-1}(k,i\omega_{n})-\Sigma(k,i\omega_{n}),$ (16)
$\displaystyle W^{-1}(k,i\nu_{n})=W_{0}^{-1}(k,i\nu_{n})-\Pi(k,i\nu_{n}),$

where the non-interacting lattice Green’s function $G_{0}$ and screened
interaction $W_{0}$ are given by

$\displaystyle G_{0}^{-1}(k,i\omega_{n})=i\omega_{n}+\mu-\epsilon_{k},$ (17)
$\displaystyle W_{0}^{-1}(k,i\nu_{n})=v^{-1}_{k}.$

  3. 3.

Use the (causal) EDMFT self-consistent equations in Eq. (6) to calculate the
dynamical mean fields $\mathcal{G}(i\omega_{n})$ and $\mathcal{U}(i\nu_{n})$
of the impurity model,

$\displaystyle\mathcal{G}^{-1}\left(i\omega_{n}\right)$
$\displaystyle=\quad\langle
G\rangle^{-1}+\langle\Sigma\rangle\,\,\,\,(+\,\,\Delta_{\text{cor}}),$
$\displaystyle\mathcal{U}(i\nu_{n})$ $\displaystyle=\left[\langle
W\rangle^{-1}+\langle\Pi\rangle\right]^{-1}\,\,(+\,\,\mathcal{U}_{\text{cor}}).$

  4. 4.

Solve the impurity model defined by
$\mathcal{G}(i\omega_{n}),\mathcal{U}(i,\nu_{n})$ numerically with the
continuous-time quantum Monte Carlo (CTQMC) method.Werner and Millis (2006,
2010) The solver produces the impurity Green’s function
$G_{\text{imp}}(i\omega_{n})$ and density-density correlation function
$\chi_{\text{imp}}(i\nu_{n})$. The bosonic Green’s function can be calculated
with the relation

$W_{\text{imp}}(i\nu)=\mathcal{U}(i\nu_{n})-\mathcal{U}(i\nu_{n})\chi_{\text{imp}}(i\nu_{n})\mathcal{U}(i\nu_{n}).$
(18)

Figure 3: Flow-chart of the $GW$+EDMFT schemes. The red parts in brackets are
added in the causal version of $GW$+EDMFT, while in the absence of these
terms, the scheme corresponds to the standard local version of $GW$+EDMFT.

  5. 5.

Extract the impurity self-energy and polarization from the impurity Dyson
equations,

$\displaystyle\Sigma_{\text{imp}}(i\omega_{n})=\mathcal{G}^{-1}(i\omega_{n})-G_{\text{imp}}^{-1}(i\omega_{n}),$
(19)
$\displaystyle\Pi_{\text{imp}}(i\nu_{n})=\mathcal{U}^{-1}(i\nu_{n})-W_{\text{imp}}^{-1}(i\nu_{n}).$

  6. 6.

Calculate the $GW$ self-energy and polarization,

$\displaystyle\Sigma^{GW}(k,\tau)=$
$\displaystyle-\frac{1}{N_{p}}\sum_{p}G(p,\tau)W^{c}(k-p,\tau)$ (20)
$\displaystyle-\frac{1}{N_{p}}\sum_{p}G(p,0^{+})v(k-p),$
$\displaystyle\Pi^{GG}(k,\tau)=$ $\displaystyle
2\frac{1}{N_{p}}\sum_{p}G(p,\tau)G(p-k,-\tau),$

and extract their non-local parts,

$\displaystyle\Sigma^{GW}_{\text{nonloc}}(k,\tau)$
$\displaystyle=\Sigma^{GW}(k,\tau)-\frac{1}{N_{k}}\sum_{k}\Sigma^{GW}(k,\tau),$
(21) $\displaystyle\Pi^{GG}_{\text{nonloc}}(k,\tau)$
$\displaystyle=\Pi^{GW}(k,\tau)-\frac{1}{N_{k}}\sum_{k}\Pi^{GG}(k,\tau).$

Here, $\frac{1}{N_{k}}\sum_{k}\ldots$ denotes the numerical average over the
first BZ and $W^{c}(k,\tau)=W(k,\tau)-v_{k}\delta(\tau).$

  7. 7.

Obtain the $GW$+EDMFT self-energy and polarization by combining the non-local
$GW$ contribution and the local EDMFT contribution,

$\displaystyle\Sigma^{GW+\text{EDMFT}}(k,i\omega_{n})$
$\displaystyle=\Sigma_{\text{imp}}(i\omega_{n})+\Sigma^{GW}_{\text{nonloc}}(k,i\omega_{n}),$
(22) $\displaystyle\Pi^{GW+\text{EDMFT}}(k,i\nu_{n})$
$\displaystyle=\Pi_{\text{imp}}(i\omega_{n})+\Pi^{GG}_{\text{nonloc}}(k,i\omega_{n}).$

With this, we have updated the momentum-dependent lattice self-energy
$\Sigma(k,i\omega)$ and polarization $\Pi(k,i\nu_{n})$ and can go back to step
2 for the next loop. Upon convergence, we obtain the fermionic and bosonic
Green’s functions in the imaginary frequency (or time) domain. The complete
$GW$+EDMFT scheme is illustrated in Fig. 3.

##  III Results

###  III.1 General remarks

In this section, we use the two self-consistency schemes to study the
$GW$+EDMFT solution of the triangular-lattice extended Hubbard model with
repulsive interactions $U$, $V>0$. First, we map out the phase diagram for
hopping $t=-0.02$ eV, corresponding to a bandwidth $W=|{9t}|=0.18$ eV. This
choice is motivated by the width of the Hubbard bands in the C phase of
1$T$-TaS2.Ligges et al. (2018); Wu et al. (2021) Next, we compare the results
from the causal and local self-consistency scheme in different regions of the
phase diagram. Finally, we identify the interaction parameters appropriate for
1$T$-TaS2 and study the doping dependence of the screening properties and the
electronic structure. Unless otherwise stated, we run the simulations for
inverse temperature $\beta=500$ eV-1, corresponding to $T=23$ K. In the
following, we measure energies in units of eV and omit the units for
convenience. The lattice constant is set to $a=1$.

###  III.2 Phase diagram

Figure 4: $GW$+DMFT phase diagrams of the half-filled single-band extended
Hubbard model on the triangular lattice. Top: Phase diagram with only local
and NN interactions, for an initial metallic (triangles) and Mott insulating
(stars) solution in the self-consistency loop. The black marks indicate the
$U$ and $V$ parameters which will be further studied below. The black line
corresponds to the fixed ratio $V/U=0.19$. Bottom: Phase diagram for the model
with long-range interaction treated by Ewald summation (triangles, starting
from a Mott insulating solution), and for comparison the result for the NN
interaction case (stars). The red and blue arrows indicate the shift of the
phase boundary as a result of the longer-ranged interaction. (Local self-
consistency scheme, $\beta=500$, $t=-0.02$.)

We first apply the local $GW$+EDMFT scheme to map out the phase diagram of the
half-filled extended Hubbard model on the triangular lattice.

In the weak coupling regime, $V/t\ll 1$ and $U/t\ll 1$, electrons can hop
easily between neighboring sites, which results in a metallic phase (M). For
dominant local interactions, $U/t\gg 1$ and $V/U\ll 1$, a Mott insulating
state (MI) induced by strong local correlation is expected. This state, which
minimizes the density of doublons, and thus the potential energy, is unstable
to antiferromagnetic order on a bipartite lattice, if spin symmetry breaking
is allowed. Similarly, in the $V/U\gg 1$ limit, the electrons on a bipartite
lattice would arrange themselves in an ordered pattern, with one sublattice
almost doubly occupied and the other sublattice almost empty, corresponding to
a commensurate charge-ordered (CO) phase. Since the triangular lattice of our
study is not bipartite, a suppression of the aforementioned ordering
instabilities due to lattice frustration can be expected, although this effect
is not fully captured by $GW$+EDMFT, which yields a result representative of
the bipartite situation. In the following, we will restrict our calculations
to paramagnetic solutions, without imposing any constraint on the charge
ordering tendencies for $V/U\gg 1$.

An effect which is captured by $GW$+EDMFT regardless of the lattice geometry
is that nonzero inter-site interactions $V$ produce a screening of the on-site
$U$, through charge density fluctuations on other sites. We thus expect that
the metallic phase is stabilized by increasing $V$.

Figure 5: From left to right column: real part of the impurity interaction
Re$\mathcal{U}(i\nu_{n})$, imaginary part of the hybridization function
Im$\Delta(i\omega_{n})$ and spectral function $-\text{Im}\Delta(\omega)$
obtained by Padé analytical continuation. The interaction parameters
(triangles in the top panel of Fig. 4) are chosen so as to highlight the
effect of an increase in the non-local interaction in the metallic phase. The
data are for the model with NN nonlocal interactions, treated with the local
and causal self-consistency scheme.

In Fig. 4, we present the phase diagrams of the triangular lattice model as a
function of the local and non-local interaction strength in the $U$-$V$ domain
for two kinds of non-local interactions: The top panel shows the result for
the model with $V$ truncated to the nearest neighbor (NN) lattice sites. Here
we indicate two boundaries for the MI phase, obtained by starting the self-
consistency loop from a metallic (line with triangles) and Mott insulating
(line with stars) solution, respectively. The bottom panel shows the phase
diagram for the model with long-range interactions treated with the Ewald
lattice summation (proportional to the result shown in Fig. 1). The boundary
between the Mott insulator and metal is determined from the low-frequency
behavior of Im$G(i\omega_{n})$,Sun and Kotliar (2002) while the boundary to
the charge-ordered (CO) phase is deduced from the divergence of the local
charge susceptibility $\chi_{\text{loc}}(i\nu=0)$.Huang et al. (2014) Beyond
this phase boundary we do not attempt to stabilize the CO phase, which might
anyhow be an artifact of the DMFT-based treatment. In the $V=0$ limit, i.e.
the conventional Hubbard model, the Mott transition occurs at $U/W=1.25$. As
$V$ is increased, the metallic phase is indeed stabilized and the phase
boundaries for the two models start to deviate. In the bottom panel of Fig. 4,
looking at the phase boundary to the CO phase, we note that including the
long-range interactions suppresses the CO phase a lot in the large-$V$ and
small/intermediate-$U$ region (indicated by the red arrow) and less so in the
large-$U$ region (indicated by the blue arrow). This shows that within the
$GW$+EDMFT description, the long-range interaction introduces additional
frustration, which makes it difficult to stabilize the ordered phase. A
similar effect was also observed in the case of the bipartite square
lattice.Huang et al. (2014)

###  III.3 Comparison of the two self-consistency schemes

To compare the auxiliary fields and physical quantities obtained with the two
self-consistency schemes, we simulate the half-filled extended Hubbard with
different $U$ and $V$ values representative of the metal and Mott insulator
phase, see the black markers (crosses and triangles) in the top panel of Fig.
4. The black line indicates the ratio $U/V=0.19$ which, as explained later,
corresponds to the ratio estimated for 1$T$-TaS2. Here, we focus on the
results obtained with the NN nonlocal interaction.

Figure 6: From left to right column: real part of the impurity interaction
Re$\mathcal{U}(i\nu_{n})$, imaginary part of hybridization function
Im$\Delta(i\omega_{n})$ and spectral function $-\text{Im}\Delta(\omega)$ for
increasing onsite interaction $U$ and fixed $V/U$ ratio (crosses in the top
panel of Fig. 4) with NN nonlocal interactions.

In Fig. 5 and Fig. 6, we show the evolution of the bosonic and fermionic Weiss
fields as a function of the interaction parameters of the model. The
corresponding points in the phase diagram in the top panel of Fig. 4 are
indicated by triangles for Fig. 5 and by crosses for Fig. 6. In the latter
case the first two solutions are in the metallic region, and the last one is
Mott insulating. While the quantities defined on the Matsubara axis are a
direct output of the $GW$+EDMFT simulation, the hybridization function on the
real axis, $-\text{Im}\Delta(\omega)$ has been obtained by Padé analytical
continuation.Vidberg and Serene (1977)

A noncausal impurity interaction $\mathcal{U}(i\nu_{n})$ is characterized by
the appearance of a negative slope on the Matsubara-frequency axis.Nilsson et
al. (2017) As we can see in Fig. 5 and Fig. 6, all the local effective
interactions obtained with the local self-consistency scheme are noncausal
(blue lines). In contrast, the impurity interactions obtained with the causal
scheme (orange lines) feature a positive slope at all Matsubara frequencies,
and hence are causal. This confirms that the modified self-consistency scheme
proposed in Ref. Backes et al., 2020 indeed eliminates noncausal features in
the impurity interaction.

Figure 7: From left to right column: fully screened local interaction
Re$W_{\text{loc}}(i\nu_{n})$, local self-energy
Im$\Sigma_{\text{loc}}(i\omega)$ and local fermionic spectral function
$A_{\text{loc}}(\omega)$ obtained by the local and causal $GW$+EDMFT scheme
for the model with fixed $U=0.2$ and increasing NN inter-site interactions.

Figure 8: From left to right column: fully screened local interaction
Re$W_{\text{loc}}(i\nu_{n})$, local self-energy
Im$\Sigma_{\text{loc}}(i\omega)$ and local fermionic spectral function
$A_{\text{loc}}(\omega)$ obtained for increasing onsite interaction $U$ and
fixed $V/U$ ratio.

The effective interaction $\mathcal{U}$ from the causal scheme is consistently
smaller than in the case of the local scheme, which suggests weaker
correlation effects. On the other hand, the hybridization function from the
causal scheme on the Matsubara-frequency axis is also smaller than the
corresponding one obtained from the local scheme, which implies a reduced
effective hopping and stronger correlations. The net effect of these two
opposing trends is difficult to guess, but we expect that the interaction
strength resulting from the two schemes will be similar. On the real axis
$-\text{Im}\Delta(\omega)$ does not display any prominent noncausal structure
in both schemes. This is consistent with our experience from numerous previous
$GW$+EDMFT studies,Boehnke et al. (2016); Nilsson et al. (2017); Petocchi et
al. (2020a); Petocchi et al. (2021) which showed that non-causalities in
$\Delta$ are generically not an issue in the application of the local scheme.
What the spectral functions confirm is the systematically weaker hybridization
in the solution obtained from the causal self-consistency loop.

The corrections on the Weiss fields in Eqs. (7) do not affect the causality of
the physical quantities $W_{\text{loc}}(i\nu_{n})$ and
$\Sigma_{\text{loc}}(i\omega_{n})$ that we report in Fig. 7 and Fig. 8 for the
same points in the phase diagram as considered above. The application of the
causal self-consistency scheme results, for all the $U$-$V$ points we studied
at half-filling, in larger screened interactions and self-energies, indicating
an enhanced correlation strength. The same effect can be noticed in the local
spectral functions
$A_{\text{loc}}(\omega)=-\frac{1}{\pi}\text{Im}G_{\text{loc}}(\omega)$, where
one finds that the quasi-particle peak near $\omega\approx 0$ in the metallic
solutions is narrower for the causal scheme than for the local scheme, while
the gap in the Mott insulating solution is larger.

By comparing the first column of Fig. 7 and the middle row of Fig. 8 (with the
same on-site interaction $U=0.2$) one finds that an increase of $V$ in the
metallic phase systematically enhances the screening of the local interaction
$W$, and the same holds for the impurity interaction $\mathcal{U}$.
Considering the rows in Fig. 8, where we increase $U$ while fixing the $V/U$
ratio, one can see that the relative change between the bare and static
screened interaction is reduced with increasing correlation strength (again
the same holds for $\mathcal{U}$), and that the screening effect is strongly
reduced after the transition from the metallic to the Mott insulating phase
(note the different scales).

These findings are consistent with the results reported for a dimer by Backes
et al.,Backes et al. (2020) although less pronounced. As noted above, the net
increase in the interaction strength, which we observe in the causal scheme at
half-filling, is the result of a decrease in the local effective interaction
$\mathcal{U}(i\nu_{n})$, which is over-compensated by an even stronger
reduction of the hybridization strength. As will be shown below, the opposite
trend can be found in doped systems, so that it is difficult to make a general
statement about stronger or weaker interactions resulting from one or the
other scheme. Overall, we find that the physical results ($G$, $W$, and
corresponding spectra) are rather similar for the two schemes, which is
comforting in the sense that the general approach appears to be robust. In
fact, several recent applications of the (local) $GW$+EDMFT scheme to
correlated materials provided meaningful estimates of the correlation
strength.Nilsson et al. (2017); Petocchi et al. (2020a, b); Petocchi et al.
(2021)

###  III.4 Doping and temperature dependence in a model for 1$T$-TaS2

####  III.4.1 Choice of parameters

To connect our model calculations to experimental results for the layered
transition-metal dichalcogenide $1T$-TaS2,Cho et al. (2016); Ligges et al.
(2018); Wu et al. (2021) we first estimate the realistic ratio between the on-
site interaction $U$ and the NN interaction $V$ in the half-filled, Mott
insulating system with poor screening. At low temperature ($T\leq 180$ K),
$1T$-TaS2 is in a commensurate charge density wave (C-CDW) state,
characterized by a David-star-like lattice distortion,Cho et al. (2016); Wang
et al. (2020); Wu et al. (2021) which we illustrate in Fig. 9. This structural
distortion modifies the electronic band structure, which at low energy can be
represented by a single-band model on a triangular lattice. In this model,
each lattice site corresponds to a 13-atom cluster.

Figure 9: The C-CDW pattern in the low-temperature phase of 1$T$-TaS2,
sketched based on the experimental results in Refs. Cho et al., 2016; Wang et
al., 2020; Wu et al., 2021. Purple dots are the Ta atoms, distorted in the
direction of the purple arrows, and the red lines indicate the David-star-like
super-cells, each including 13 lattice sites. The red disk shows the area with
high charge density seen by STM. $R=\sqrt{13}$ is the distance between two
neighboring super-cells and $r$ is the size of the disk.

To determine the ratio $V/U$, we make the simple assumption that the electrons
are distributed uniformly inside a disk representing the David-star-like
molecular orbital. From the scanning tunneling microscope (STM) image,Cho et
al. (2016) one can roughly estimate the size $r\approx 0.4$ of this disk (in
units of lattice spacing), while $R=\sqrt{13}$ is the distance between the
molecular orbitals (in units of the lattice spacing of the undistorted
lattice). The interactions $U$ and $V$ are then proportional to

$\displaystyle
U\sim\int_{D_{1}}\int_{D_{1}}d\vec{r}_{1}d\vec{r}_{1}^{\prime}\frac{1}{|\vec{r}_{1}-\vec{r}_{1}^{\prime}|},$
(23) $\displaystyle
V\sim\int_{D_{1}}\int_{D_{2}}d\vec{r}_{1}d\vec{r}_{2}\frac{1}{|\vec{r}_{1}-\vec{r}_{2}|},$

where $\vec{r}_{1}$ and $\vec{r}^{\prime}_{1}$ are in the same disk ($D_{1}$),
while $\vec{r}_{1}$ and $\vec{r}_{2}$ are in neighboring disks ($D_{1}$ and
$D_{2}$). The thus estimated ratio is $V/U\sim 0.19$.

To estimate the values of $U$ and $t$, we compare the separation between the
Hubbard bands and the width of the Hubbard bands in $A_{\text{loc}}(\omega)$
to the results from photoemission and STM studies. The on-site interaction $U$
essentially determines the gap between the upper Hubbard band (UHB) and the
lower Hubbard band (LHB), which is about 0.4 eV, Cho et al. (2016); Wu et al.
(2021) while the width of the Hubbard bands is approximately 0.2 eV.Cho et al.
(2016); Wu et al. (2021) A good match with experiment is obtained for $U=0.42$
and $t=-0.02$ (the bandwidth is $W=9|t|$ for the triangular lattice). Given
the onsite interaction of $U=0.42$ and the ratio $V/U=0.19$ we choose the NN
interaction as $V=V_{1}=0.08$. This “realistic” parameter set for 1$T$-TaS2 is
indicated by the black point in Fig. 4. We note that it is not close to the CO
phase boundary.

As in the previous sections, we set the inverse temperature to $\beta=500$
eV-1 ($T=23$ K), unless otherwise stated, and use eV as the unit of energy.

Figure 10: Local self-energy $\Sigma_{\text{loc}}$, fermionic spectral
function $A_{\text{loc}}(\omega)$, and distribution function
$A_{\text{loc}}(\omega)F(\omega)$ (from left to right). The doping is
increased row-wise from the undoped configuration to $5\%$ and $10\%$ hole
concentration. The interaction parameters are $U=0.42$ and $V=0.08$.

####  III.4.2 Doping dependence

In a recent time-resolved photoemission study of 1$T$-TaS2, the authors
concluded that their nominally half-filled system may in fact be substantially
hole-doped.Ligges et al. (2018) The indirect evidence for this was the
surprisingly fast recombination of photo-doped doublons and holons, which in a
pure Mott state with the given ratio between gap size and bandwidth should be
longer-lived.Eckstein and Werner (2011) This motivates us to compute the local
spectral function for the realistic parameters and different hole dopings. In
Figs. 10, the three rows correspond to half-filling, 5% hole doping ($n=0.95$)
and 10% hole doping ($n=0.9$). The half-filled system is Mott insulating,
while the doped systems are metallic.

In Fig. 10, the gap between the upper Hubbard band (UHB) and lower Hubbard
band (LHB) is $0.4$ eV, and the width of the Hubbard bands is approximately
$0.2$ eV, consistent with the experimental results. When the system is 5% or
10% hole doped, a prominent quasi-particle peak appears at the edge of the
lower Hubbard band, which is shifted to $\omega=0$, while the peaks of the LHB
and UHB are shifted to $-0.1$ eV and $0.3$ eV, respectively.

A quasi-particle peak was not reported in Ref. Ligges et al., 2018 in the
initial equilibrium spectrum. The rightmost column of Fig. 10 shows the
occupation function, i.e. the spectral function multiplied with the Fermi
function $F(\omega)=1/(e^{\beta\omega}+1)$, on a logarithmic scale. This is
the quantity measured in photo-emission experiments. The results for the
undoped and doped systems differ substantially: near the Fermi level, the
occupation function of the doped systems exhibit an almost linear decay (in
the log plot) over three orders of magnitude, in an energy interval of width
$\Delta\omega\approx 0.02$, which is determined by the temperature. In the
undoped case, on the other hand, the first three orders of magnitude decrease
from the peak value of the occupation is controlled by the shape of the lower
Hubbard band, rather than by the Fermi distribution function, since in this
case the chemical potential is in the gapped region.

We plot in Fig. 10 the results from three different calculations: blue lines
are for the local scheme with NN interactions, orange lines for the local
scheme with long-ranged interactions, and green lines for the causal scheme
with NN interactions. As far as the spectral function and occupation function
are concerned, the differences between the models and schemes are minor and on
a scale comparable to the intrinsic uncertainties of the MaxEnt approach.

Figure 11: The local fully screened interaction $W_{\text{loc}}$ in the
Matsubara frequency domain and in the real frequency domain after Padé
analytical continuation. In the top panels, the model is half-filled ($n=1$).
The middle panels are for $5\%$ hole doping ($n=0.95$), and the bottom panels
for $10\%$ hole doping ($n=0.9$). The interaction parameters are $U=0.42$,
$V=0.08$.

Figure 11 shows how the screening properties are modified by the hole doping.
The left panels plot the results for $W_{\text{loc}}$ on the Matsubara
frequency axis, and the middle and right panels the real and imaginary parts
after maximum entropy analytical continuation. Again, we show results for the
local scheme with NN interactions (blue), the local scheme with long-ranged
interactions (orange) and the causal scheme with NN interactions (green). By
comparing the static values, one notices that while the local screened
interaction in the half-filled insulator is slightly larger in the causal
scheme than in the local scheme, the opposite is true in the hole-doped
systems. Hence, the effect of the correction terms $\Delta_{\text{cor}}$ and
$\mathcal{U}_{\text{cor}}$ (Eq. (7)) on the correlation strength depends on
doping, and more generally on the choice of parameters. Overall, however, we
find a good agreement of the physical observables obtained by the two self-
consistency schemes and also between the models with NN and long-ranged
interactions.

The imaginary part of the real-frequency $W_{\text{loc}}(\omega)$ reveals the
relevant screening modes in the system. For the undoped (half-filled) system
in the first row, the dominant screening mode, both in the fully screened
local interaction Im$W_{\text{loc}}(\omega)$ and the retarded impurity
interaction Im$\mathcal{U}(\omega)$ (not shown) approximately matches the 0.4
eV energy separation between the LHB and UHB. This suggests that screening in
the Mott insulator is associated with charge excitations across the Mott
gap.Huang et al. (2014) The screening effect is however small, as one can
deduce from the small reduction of the real part, compared to the bare $U$.
For the doped systems, the peak appears around $\leavevmode\nobreak\ 0.1$ eV,
which approximately corresponds to the gap between the LHB and the quasi-
particle peak in the spectral function. Hence, these types of charge
excitations contribute primarily to the additional screening in the doped Mott
system. Because of this additional screening, the static value of
$W_{\text{loc}}$ is now substantially smaller than the bare $U$, and in the
case of $n=0.9$ close to zero (in $\mathcal{U}$, the reduction is only about
15-20%).

The effect of doping on the imaginary part of the local self-energy is plotted
in the left panels of Fig. 10. These results confirm the transition from Mott
insulating to metallic behavior, and the previous observation that in the
doped systems (especially for 5% doping), the results from the causal scheme
are less correlated. They also confirm that the differences between the local
and causal scheme are rather small, as far as the correlation strength is
concerned, and the same is true for the difference between the models with NN
and long-ranged interactions.

####  III.4.3 Temperature dependence

In Fig. 12 we show how a temperature increase changes the local spectral
function $A_{\text{loc}}(\omega)$ and the electron distribution function
$A_{\text{loc}}(\omega)F(\omega)$ for different hole concentrations. The three
rows are for half-filling, 5% hole-doping, and 10% hole-doping, respectively,
while $U=0.42$, $V=0.08$ are fixed. In each panel, we plot the results for
$T=23$ K, $30$ K and $160$ K. Up to $160$ K, i.e. inside the C-CDW phase,
temperature does not significantly alter the spectra nor any of the physical
fermionic and bosonic fields. Also the gap size is merely affected by thermal
broadening. The main effect of increasing temperature is a partial melting of
the quasi-particle peaks in the doped systems. In the occupation functions
(right panels), the higher temperature leads to a slower decay of the
occupation near the Fermi edge in the doped system, because of the broader
$F(\omega)$. In the Mott insulator, at the highest temperature, one finds a
more prominent thermal (doublon) population of the upper Hubbard band, and a
small shift of the whole spectrum relative to the chemical potential.

Figure 12: Spectral function $A_{\text{loc}}(\omega)$ and occupation function
$A_{\text{loc}}(\omega)F(\omega)$ for the half-filled, 5% hole-doped and 10%
hole-doped system, respectively, for the indicated values of temperature
($U=0.42$, $V=0.08$).

####  III.4.4 Comparison with experiment

In Fig. 13 we compare the distribution functions for the half-filled, 2.5% and
5% hole-doped system to the equilibrium photoemission spectrum reported in
Fig. 2a of Ref. Ligges et al., 2018. To enable a direct comparison, we use the
same temperature, $T=30$ K, as in the experiment. The red dashed line plots
the raw experimental data, while the purple curve in the top panel plots the
experimental result shifted by $\Delta\omega=0.065$ eV on the energy axis.
Experimentally, it is difficult to determine the chemical potential in an
insulating system, and it is thus more meaningful in the gapped case to match
the experimental curve with the upper edge of the lower Hubbard band in the
simulation data. By doing so we notice that the measured distribution function
reproduces very well the shape of the Hubbard band near the gap edge, and that
the decay over three orders of magnitude is much slower than in a metallic
system (middle and bottom panel), where it is controlled by the Fermi function
cutoff, even for a small hole doping concentration. In addition, while the
metallic quasi-particle peak in the distribution function and on a log-scale
plot is not very prominent, no hint of such a feature is evident in the
experimental data. These results suggest that the 1$T$-TaS2 sample used in the
experiment was not significantly hole-doped and that the doublon population
dynamics measured in Ref. Ligges et al., 2018 might be controlled by phenomena
that their modeling did not capture.

Figure 13: Occupation function $A_{\text{loc}}(\omega)F(\omega)$ for the half-
filled, 2.5% and 5% hole-doped system (from top to bottom). The dashed red
line corresponds to the equilibrium data plotted in Fig. 2(a) of Ref. Ligges
et al., 2018. To match the experiment, the $GW$+EDMFT calculations are carried
out here for $T=30$ K ($U=0.42$, $V=0.08$). The purple line is the
experimental distribution shifted to match the calculated results (the
chemical potential in the experiment is not precisely known).

Apart from the slope near the chemical potential, there are two notable
differences between the measured and calculated distributions. (i) The
measured spectrum shows two plateau-like structures, with the lower one
extending down to energies of approximately $-0.6$ eV. This latter feature
comes from lower-lying bands, including sulphur bands, and are thus absent in
our simulations. (ii) Even the higher-energy hump, which can be interpreted as
the lower Hubbard band, is flatter and broader than in the simulations. This
appears to be a consequence of the three-dimensional nature of 1$T$-TaS2,
which consists of a stacking of triangular-lattice-type layers. The hopping in
the stacking direction is actually large, compared to the in-plane hopping,
Lee et al. (2019); Pasquier and Yazyev (2021) so that one may expect a
broader, more 1-dimensional DOS than what comes out of our triangular-lattice
simulation. The latter is in rough agreement with the DOS measured in STM
experiments,Cho et al. (2016); Wu et al. (2021) but these experiments are less
sensitive to the stacking direction than photoemission measurements.

##  IV Conclusions

In this work, we studied the phase diagram and correlation functions of the
single-band extended Hubbard model on the two-dimensional triangular lattice
using the standard local implementation of the self-consistency
equationsBiermann et al. (2003); Ayral et al. (2013) and a recently proposed
causal variant.Backes et al. (2020) The Ewald lattice summation method has
also been implemented in order to investigate the role of longer-ranged
nonlocal interactions. Our test calculations in different regions of the phase
diagram confirmed that the modified $GW$+EDMFT self-consistency loop of
Backes, Sun and BiermannBackes et al. (2020) removes the non-causal features
in the effective impurity interaction while the hybridization function turns
out to be causal in both schemes (for the model parameters considered). The
causal variant results in slightly stronger correlations effects in the half-
filled system, while in hole-doped systems, the correlations – as measured by
the imaginary part of the self-energy – can be weakened. We showed that this
is the result of two opposing effects: while the effective interaction
$\mathcal{U}$ produced by the causal scheme tends to be smaller than in the
case of the local implementation, the hybridization function $\Delta$ is also
smaller, which suppresses the kinetic energy. Overall, this results in
physical observables (Green’s functions and screened interactions) which do
not significantly differ between the two schemes. This is comforting in the
sense that recent ab-initio studies based on the local implementation of
$GW$+DMFT have produced results which are in good agreement with
experiments.Petocchi et al. (2020a); Petocchi et al. (2021)

While noncausal structures in auxiliary quantities like $\mathcal{U}$ and
$\Delta$ are, conceptually, not a problem,Nilsson et al. (2017) the causal
dynamical mean fields produced by the scheme of Backes et al. can be
advantageous from a numerical point of view. In particular, we expect that the
real-time propagation in nonequilibrium implementations of $GW$+EDMFTGolež et
al. (2017, 2019) becomes more stable. Causality also enables the use of
maximum entropy analytical continuation for the analysis of the spectral
content of the auxiliary fields.

We presented simulation results for parameter values appropriate for the
layered dichalchogenide 1$T$-TaS2. In the low-temperature C-phase, this
material is a polaronic Mott insulator that can be described by an effective
single-band Hubbard model on the triangular lattice.Wilson et al. (1975);
Fazekas and Tosatti (1979); Sipos et al. (2008) Based on a rough estimate of
the size of the molecular orbitals, we determined the ratio $V/U=0.19$ and, by
comparing the width and separation of the Hubbard bands measured by
photoemission and STM to the $GW$+EDMFT spectra, fixed the values of $U$ and
$t$. This yields the parameters $U=0.42$ eV, $V=0.08$ eV, and $t=-0.02$ eV for
our low energy description. Models with NN and longer-ranged nonlocal
interactions treated with the Ewald summation were found to produce very
similar results. Here we should note that the bandwidth for a single layer of
TaS2, as estimated by ab-initio calculations, is much smaller than 0.2 eV.
Pasquier and Yazyev (2021) This means that the above-estimated parameters
implicitly take into account the three-dimensional nature of 1$T$-TaS2, which
has larger hoppings in the stacking direction than in the in-plane direction.
The triangular-lattice simulations do however not explicitly capture the
anisotropy of the lattice and this may be the reason why the Hubbard bands in
the simulations are more peaked than in the (bulk sensitive) photoemission
experiments.

Motivated by the results of Ref. Ligges et al., 2018, which suggested a
significant hole-doping in a nominally undoped sample of 1$T$-TaS2, we
investigated the hole-doping effect on the spectral function, self-energy,
screened interaction and effective local interaction of the model with
realistic parameters. These results show that in the absence of disorder or
other effects not captured by $GW$+EDMFT, a hole doping of a few percent
should result in a prominent quasi-particle peak at the edge of the lower
Hubbard band at $T=30$ K. This peak was not seen in the photo-emission data.
Also, the direct comparison of the experimental spectrum to the simulated
occupation functions for half-filling, 2.5% and 5% hole doping showed a very
good agreement with the undoped result (near the gap edge, and after an
appropriate energy shift of the experimental data), while there were
significant differences to the distribution functions of the hole doped
systems. In particular, the latter exhibit a much more rapid decay of the
occupation, which is controlled by the temperature, rather than by the shape
of the lower Hubbard band. These findings suggest that the 1$T$-TaS2 crystal
used in the experiments of Ref. Ligges et al., 2018 was not significantly
hole-doped, as was concluded based on the short doublon life-time. An improved
nonequilibrium DMFT description of 1$T$-TaS2 (compared to Ref. Ligges et al.,
2018) might be obtained by employing more reliable impurity solvers and by
considering photo-excitations from deeper-lying bands. It would also be
worthwhile to study the stacking effect in nonequilibrium simulations, as well
as different cooling mechanisms.Murakami et al. (2015); Werner et al. (2019)

We also considered the temperature effect in equilibrium systems, but apart
from the weight of the quasi-particle peak, the physical properties
(correlation strength, screening) are not much temperature dependent in the
range ($T\lesssim 160$ K) relevant for the C-phase of 1$T$-TaS2.

###### Acknowledgements.

We thank S. Biermann for helfpul discussions, and U. Bovensiepen for providing
the raw data for the equilibrium spectrum in Fig. 2a of Ref. Ligges et al.,
2018. This work was supported by the Swiss National Science Foundation via
NCCR Marvel and Grant No. 200021-196966. The calculations have been performed
on the Beo05 cluster at the University of Fribourg.

##  Appendix A Ewald summation

We consider a Coulomb repulsion $V(\mathbf{r})=\frac{1}{r}$ in a two-
dimensional system. Applying the Fourier transformation, one finds

$V(\mathbf{k})=\int
d\mathbf{r}e^{-i\mathbf{k}\cdot\mathbf{r}}\frac{1}{r}=\frac{2\pi}{k}.$ (24)

With the help of the the identity

$\frac{1}{k}=\frac{2}{\sqrt{\pi}}\int_{0}^{\infty}dte^{-(kt)^{2}}$ (25)

one can express $V(\mathbf{k})$ as an integral. We separate the integral into
the short-range contribution from $[0,\eta]$ and the long-range contribution
from $[\eta,\infty]$, which yields

$\displaystyle
V^{S}(\mathbf{k})=\int_{0}^{\eta}dt\frac{2\pi}{k}=4\sqrt{\pi}\int_{0}^{\eta}dte^{-(kt)^{2}},$
(26) $\displaystyle
V^{L}(\mathbf{k})=\int_{\eta}^{\infty}dt\frac{2\pi}{k}=4\sqrt{\pi}\int_{\eta}^{\infty}dte^{-(kt)^{2}}.$

As shown in the main text and in Fig. 1, the parameter $\eta$ plays a similar
role as the $\alpha$ parameter in the Yukawa potential.

Now we treat the two parts separately. For the short range part, the inverse
Fourier transformation to real space gives

$\displaystyle V^{S}(\mathbf{r})$
$\displaystyle=\frac{1}{N_{k}}\sum_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{r}}4\sqrt{\pi}\int_{0}^{\eta}dte^{-(kt)^{2}}$
(27)
$\displaystyle=\frac{1}{\sqrt{\pi}}\int_{0}^{\eta}dt\frac{1}{t^{2}}e^{-\frac{r^{2}}{4t^{2}}},$

which, using the variable change $x\equiv\frac{r}{2t}$, becomes

$\displaystyle V^{S}(\mathbf{r})$
$\displaystyle=\frac{1}{r}\text{erfc}\Big{(}\frac{r}{2\eta}\Big{)}.$ (28)

Here, $\text{erfc}(x)$ is the complementary error function, defined as
$\text{erfc}(x)=\frac{2}{\sqrt{\pi}}\int_{x}^{\infty}dxe^{-x^{2}}$. For the
long-range part, using the variable change $x\equiv kt$, one finds

$\displaystyle V^{L}(\mathbf{r})$
$\displaystyle=\frac{1}{N_{\mathbf{k}}}\sum_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{r}}4\sqrt{\pi}\int_{\eta}^{\infty}dte^{-(kt)^{2}}$
(29)
$\displaystyle=\frac{1}{N_{\mathbf{k}}}\sum_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{r}}\frac{2\pi}{k}\text{erfc}(k\eta).$

Due to charge neutrality, the $\mathbf{k}=0$ term, corresponding to the
average potential, should be removed from the sum. Adding the expressions for
the long-range and short-range potentials in Eq. (28) and Eq. (29) one obtains
Eq. (2) in the main text.

##  Appendix B Derivation of the (causal) EDMFT self-consistency equations

In this appendix, we derive, following Ref. Backes et al., 2020, the causal
EDMFT self-consistency equations (Eqs. (6) and (7)) using a generalized cavity
method.Georges et al. (1996) The basic idea of the cavity method is to
separate a lattice into two parts, a single site denoted as “site 0” and the
remaining part of the lattice, i.e. the lattice with cavity, denoted by a
superscript “(0)”. After separating the action of the lattice model into terms
related to site 0, the rest, and terms connecting the two parts, we integrate
out all the degrees of freedom from the lattice with cavity, which leaves only
those of the single site, denoted by the subscript “0”. Starting from the
lattice action, one thus obtains an effective action $S_{\text{eff}}$ of an
impurity model. Formally, one can express the relationship between the lattice
action and the action of the effective single-site model (impurity model) as

$\frac{1}{Z}\int\prod_{i\neq
0,\sigma}\mathcal{D}[c^{*}_{i\sigma},c_{i\sigma}]e^{-S[c_{i\sigma}^{*},c_{i\sigma}]}=\frac{1}{Z_{\text{eff
}}}e^{-S_{\text{eff }}\left[c_{0\sigma}^{*},c_{0\sigma}\right]}.$ (30)

In the action of the full lattice model, the terms contributed by the degrees
of freedom with no connection to the impurity site represent the action of the
lattice with cavity, $S^{(0)}$. The degrees of freedom coupled to the impurity
site are viewed as external sources, as the generators of the cavity Green’s
functions. Together with terms involving only the impurity site, contained in
$S_{0}$, we can thus split the full lattice action defined in Eq. (II.3) into
the three parts $S=S_{0}+S^{(0)}+\Delta S$, with the explicit form

$\displaystyle S_{0}=$
$\displaystyle\int_{0}^{\beta}d\tau\left\\{\sum_{\sigma}c_{0\sigma}^{*}(\tau)\left(\partial_{\tau}-\mu\right)c_{0\sigma}(\tau)+\frac{1}{2}n_{0}v_{00}n_{0}\right\\},$
(31) $\displaystyle S^{(0)}=$
$\displaystyle\int_{0}^{\beta}d\tau\left\\{\sum_{ij\neq
0,\sigma}c_{i\sigma}^{*}\left(\partial_{\tau}-\mu-
t_{ij}\right)c_{j\sigma}+\frac{1}{2}\sum_{ij,\neq
0}n_{i}v_{ij}n_{j}\right\\},$ $\displaystyle\Delta S=$
$\displaystyle\int_{0}^{\beta}d\tau\left\\{-\sum_{i\neq
0,\sigma}t_{i0}\left(c_{0\sigma}^{*}c_{i\sigma}+c_{i\sigma}^{*}c_{0\sigma}\right)+\sum_{i\neq
0}n_{i}v_{i0}n_{0}\right\\}.$

Defining $\eta_{i}\equiv t_{i0}c_{0\sigma}$ and $j_{i}\equiv v_{i0}n_{0}$ we
can rewrite the last term as

$\Delta
S[\eta_{i},\eta^{*}_{i},j_{i}]=\int^{\beta}_{0}d\tau\left\\{-\sum_{i\neq
0,\sigma}(\eta^{*}_{i}c_{i\sigma}+c^{*}_{i\sigma}\eta_{i})+\sum_{i\neq
0}j_{i}n_{i}\right\\}.$ (32)

Now, the effective action $S_{\operatorname{eff}}$ defined in Eq. (30) can be
explicitly written as

$\displaystyle\frac{e^{-S_{\operatorname{eff}}}}{Z_{\operatorname{eff}}}=\frac{e^{-S_{0}}}{Z}\int\mathcal{D}_{i\neq
0}[c^{*}_{i\sigma},c_{i\sigma}]e^{-(S^{(0)}+\Delta S)}$ (33)
$\displaystyle\Longrightarrow
S_{\operatorname{eff}}=S_{0}-\Omega+\text{const},$ (34)

where $\Omega[c^{*}_{0\sigma},c_{0\sigma}]=\ln\int\mathcal{D}_{i\neq
0}\left[c_{i}^{*},c_{i}\right]e^{-\left(S^{(0)}+\Delta S\right)}$ is the
generating functional of the cavity lattice Green’s function,Georges et al.
(1996)

$\displaystyle G_{i_{1}\ldots i_{n}j_{n}\ldots
j_{1}}^{(0)}\left(\tau_{1}\ldots\tau_{n},\tau_{1}^{\prime}\ldots\tau_{n}^{\prime}\right)$
(35) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
=(-1)^{n}\frac{\delta^{2n}\Omega}{\delta\eta_{i_{1}}^{*}\left(\tau_{1}\right)\ldots\delta\eta_{i_{n}}^{*}\left(\tau_{n}\right)\delta\eta_{j_{n}}\left(\tau_{n}^{\prime}\right)\ldots\delta\eta_{j_{1}}\left(\tau_{1}^{\prime}\right)},$
$\displaystyle\chi_{i_{1}\ldots i_{n}j_{n}\ldots
j_{1}}^{(0)}\left(\tau_{1}\ldots\tau_{n},\tau_{1}^{\prime}\ldots\tau_{n}^{\prime}\right)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ =\frac{\delta^{2n}\Omega}{\delta
j_{i_{1}}\left(\tau_{1}\right)\ldots\delta
j_{i_{n}}\left(\tau_{n}\right)\delta
j_{j_{n}}\left(\tau_{n}^{\prime}\right)\ldots\delta
j_{j_{1}}\left(\tau_{1}^{\prime}\right)}.$

Thus, we obtain the series expansion of the effective action
$S_{\operatorname{eff}}$,

$\displaystyle S_{\operatorname{eff}}$
$\displaystyle\left[\eta_{i}^{*},\eta_{i},j_{i}\right]=$ (36)
$\displaystyle-\sum_{n=1}^{\infty}\sum_{i_{1}\ldots j_{n}}\int d\tau_{1}\ldots
d\tau_{n}^{\prime}\eta_{i_{1}}^{*}\left(\tau_{1}\right)\ldots\eta_{j_{n}}\left(\tau_{n}^{\prime}\right)$
$\displaystyle\leavevmode\nobreak\ \times(-1)^{n}G_{i_{1}\ldots
i_{n}j_{n}\ldots j_{1}}^{(0)}\left(\tau_{1}\ldots\tau_{n}^{\prime}\right)$
$\displaystyle-\sum_{n=1}^{\infty}\sum_{i_{1}\ldots j_{n}}\int d\tau_{1}\ldots
d\tau_{n}^{\prime}j_{i_{1}}\left(\tau_{1}\right)\ldots
j_{j_{n}}\left(\tau_{n}^{\prime}\right)$ $\displaystyle\leavevmode\nobreak\
\times\chi_{i_{1}\ldots i_{n}j_{n}\ldots
j_{1}}^{(0)}\left(\tau_{1}\ldots\tau_{n}^{\prime}\right)+S_{0}+\text{const}.$

Thus far, every step is exact. Now, we truncate the series expansion of
$\Omega$, which is only a good approximation in a system with large
connectivity. Assuming that each lattice site has $z\to\infty$ neighbors, we
have to rescale $t\to t/\sqrt{z}$ and $V\to V/z$ in order to keep the balance
between kinetic and potential energy.111From the local perspective, only the
case when one election hops in while one hops out with amplitude
$t_{j0}t_{0i}$ contributes to the local Hamiltonian. Thus, $(zt)^{2}$ has to
be an $\mathcal{O}(1)$ quantity as $z$ goes to infinity.Müller-Hartmann (1989)
Given this rescaling, it has been proven Metzner and Vollhardt (1989); Georges
et al. (1996) that the large connectivity limit will dramatically simplify the
expansion of the effective action, because the $n^{\text{th}}$ terms scale as
$(1/z)^{n-1}$ and only the $n=1$ terms survive,

$\displaystyle S_{\operatorname{eff}}=$ $\displaystyle\int d\tau
d\tau^{\prime}\left\\{c^{*}_{i\sigma}(\tau)t_{0i}t_{j0}G^{(0)}_{ij}(\tau-\tau^{\prime})c_{j\sigma}(\tau^{\prime})\right.$
(37)
$\displaystyle\left.-n_{0}(\tau)v_{0i}v_{j0}\chi^{(0)}_{ij}(\tau-\tau^{\prime})n_{0}(\tau^{\prime})\right\\}+S_{0}+\mathcal{O}\Big{(}\frac{1}{z}\Big{)}$
$\displaystyle=$ $\displaystyle-\int^{\beta}_{0}d\tau
d\tau^{\prime}\sum_{\sigma}c_{\sigma}^{*}(\tau)\mathcal{G}(\tau-\tau^{\prime})c_{\sigma}(\tau^{\prime})$
$\displaystyle+\frac{1}{2}\int^{\beta}_{0}d\tau
d\tau^{\prime}n(\tau)\mathcal{U}(\tau-\tau^{\prime})n(\tau^{\prime})+\mathcal{O}\Big{(}\frac{1}{z}\Big{)},$

where the Weiss fields $\mathcal{G}$ and $\mathcal{U}$ are defined as,

$\displaystyle\mathcal{G}^{-1}\left(\tau_{1}-\tau_{2}\right)=$ $\displaystyle\
\delta(\tau_{1}-\tau_{2})(-\partial_{\tau_{1}}+\mu)$ (38)
$\displaystyle-\sum_{ij\neq
0}t_{0i}t_{j0}G_{ij}^{(0)}\left(\tau_{1}-\tau_{2}\right),$
$\displaystyle\mathcal{U}\left(\tau_{1}-\tau_{2}\right)=$ $\displaystyle\
\delta(\tau_{1}-\tau_{2})v_{00}$ $\displaystyle-\sum_{ij\neq
0}v_{0i}v_{j0}\chi_{ij}^{(0)}\left(\tau_{1}-\tau_{2}\right).$

Here, $G^{(0)}$ and $\chi^{(0)}$ denote the Green’s function and charge
susceptibility of the lattice with cavity,

$\displaystyle{{G}^{(0)}_{ij\sigma}}\left(\tau_{1}-\tau_{2}\right)$
$\displaystyle\equiv-\langle\mathcal{T}_{\tau}c_{i\sigma}(\tau_{1})c^{*}_{j\sigma}(\tau_{2})\rangle^{(0)},$
(39) $\displaystyle\chi_{ij}^{(0)}\left(\tau_{1}-\tau_{2}\right)$
$\displaystyle\equiv\langle n_{i}(\tau_{1})n_{j}(\tau_{2})\rangle^{(0)}.$

In the Matsubara frequency $i\omega_{n}$ $(i\nu_{n})$ domain, the Weiss fields
thus become

$\displaystyle\mathcal{G}^{-1}\left(i\omega_{n}\right)$
$\displaystyle=i\omega_{n}+\mu-\sum_{ij}t_{0i}t_{j0}G_{ij}^{(0)}\left(i\omega_{n}\right),$
(40) $\displaystyle\mathcal{U}\left(i\nu_{n}\right)$
$\displaystyle=v_{00}-\sum_{ij}v_{0i}v_{j0}\chi_{ij}^{(0)}\left(i\nu_{n}\right).$

In some previous works,Sun and Kotliar (2002); Biermann et al. (2003); Ayral
et al. (2013); Huang et al. (2014) an alternative form of the interaction part
of the action has been considered by applying the Hubbard-Stratonovich
transformationHubbard (1959) to the interaction terms in Eq. (II.3). This
allows to decouple the density-density interactions and to connect the
formalism to the LW functional $\Psi$. The resulting alternative form of the
action is

$\displaystyle
S\left[c^{*},c,\phi\right]=\\!\int_{0}^{\beta}\\!\\!d\tau\left\\{\\!-\sum_{ij\sigma}\\!c_{i\sigma}^{*}(\tau)[(-\partial_{\tau}+\tilde{\mu})\delta_{ij}-t_{ij}]c_{j\sigma}(\tau)\\!\right\\}$
(41)
$\displaystyle+\int_{0}^{\beta}d\tau\left\\{\frac{1}{2}\sum_{ij}\phi_{i}(\tau)\left[v^{-1}\right]_{ij}\phi_{j}(\tau)+i\sum_{i}\phi_{i}(\tau)n_{i}(\tau)\right\\},$

with $\phi$ the bosonic Hubbard-Stratonovich field and $\tilde{\mu}=\mu+U/2$.
Following the same steps as before, one obtains

$\displaystyle\mathcal{G}^{-1}\left(i\omega_{n}\right)$
$\displaystyle=i\omega_{n}+\mu-\sum_{ij}t_{0i}t_{j0}G_{ij}^{(0)}\left(i\omega_{n}\right),$
(42) $\displaystyle\mathcal{U}^{-1}\left(i\nu_{n}\right)$
$\displaystyle=v^{-1}_{00}-\sum_{ij}v^{-1}_{0i}v^{-1}_{j0}W_{ij}^{(0)}\left(i\nu_{n}\right),$

where $W_{ij}^{(0)}=\langle\phi_{i}(0)\phi_{j}(\tau)\rangle^{(0)}$ is the
bosonic Green’s function of the lattice with cavity. Eqs. (40) and (42) show
two different definitions of the Weiss field $\mathcal{U}$. As we will see,
they become equivalent if one uses the standard (local) self-consistency
scheme, while in the causal scheme, they are no longer obviously equivalent.

The second approximation (for lattices with finite coordination number)
concerns the relation between the correlation functions $G^{(0)}_{ij}$,
$\chi^{(0)}_{ij}$ and $W^{(0)}_{ij}$ of the lattice with cavity, and the full
lattice correlation functions $G_{ij}$, $\chi_{ij}$ and $W_{ij}$. For a
general lattice, the relations are

$\displaystyle G_{ij}^{(0)}$ $\displaystyle=G_{ij}-G_{i0}G_{0j}/G_{00},$ (43)
$\displaystyle\chi_{ij}^{(0)}$
$\displaystyle=\chi_{ij}-\chi_{i0}\chi_{0j}/\chi_{00},$ $\displaystyle
W_{ij}^{(0)}$ $\displaystyle=W_{ij}-W_{i0}W_{0j}/W_{00},$

which become exact in the large connectivity limit. The proof is done by
expressing the Green’s functions in terms of an inverse Hamiltonian matrix in
the basis of lattice sites, as described in Sec. III.C of Ref. Georges et al.,
1996, and derived originally by Hubbard.Hubbard and Flowers (1964) We can
understand the first relation by counting the paths contributing to the
Green’s function. In the large connectivity limit, the difference between
$G^{(0)}_{ij}$ and $G_{ij}$ lies in the contribution from the paths connecting
the sites $i$ and $j$ through the site $0$ only once. Symbolically, this
contribution is proportional to $G_{i0}G_{0j}$ but has to be normalized by
$G_{00}$, the probability of leaving and returning to site 0, in order to
cancel paths passing through site 0 more than once. The treatment of the
susceptibility $\chi$ and the screened interaction $W$ is analogous.

Now, inserting Eq. (43) into Eq. (40) or Eq. (42), we find

$\displaystyle\mathcal{G}^{-1}\left(i\omega_{n}\right)=$ $\displaystyle\
i\omega_{n}+\mu-\sum_{ij}t_{0i}t_{j0}G_{ij}\left(i\omega_{n}\right)$ (44)
$\displaystyle+\sum_{ij}t_{0i}G_{i0}\left(i\omega_{n}\right)G_{0j}\left(i\omega_{n}\right)t_{j0}/G_{00}\left(i\omega_{n}\right),$
$\displaystyle\mathcal{U}\left(i\nu_{n}\right)=$ $\displaystyle\
v_{00}-\sum_{ij}v_{0i}v_{j0}\chi_{ij}\left(i\nu_{n}\right)$
$\displaystyle+\sum_{ij}v_{0i}\chi_{i0}\left(i\nu_{n}\right)\chi_{0j}\left(i\nu_{n}\right)v_{j0}/\chi_{00}\left(i\nu_{n}\right),$

and in the alternative form,

$\displaystyle\mathcal{U}^{-1}\left(i\nu_{n}\right)=$ $\displaystyle\
v^{-1}_{00}-\sum_{ij}v^{-1}_{0i}v^{-1}_{j0}W_{ij}\left(i\nu_{n}\right)$ (45)
$\displaystyle+\sum_{ij}v^{-1}_{i0}W_{0i}\left(i\nu_{n}\right)W_{0j}\left(i\nu_{n}\right)v^{-1}_{j0}/W_{00}\left(i\nu_{n}\right).$

Fourier transforming to momentum space and introducing the notation
$\langle\ldots\rangle_{k}\equiv\frac{1}{N_{k}}\sum_{k}$, one findsBackes et
al. (2020)

$\displaystyle\mathcal{G}^{-1}\left(i\omega_{n}\right)=$ $\displaystyle\
i\omega_{n}+\mu-\left(\langle\epsilon G\epsilon\rangle_{k}-\langle\epsilon
G\rangle_{k}\langle G\rangle_{k}^{-1}\langle G\epsilon\rangle_{k}\right),$
(46) $\displaystyle\mathcal{U}(i\nu_{n})=$ $\displaystyle\ \langle
v\rangle_{k}-\left(\langle v\chi v\rangle_{k}-\langle
v\chi\rangle_{k}\langle\chi\rangle_{k}^{-1}\langle\chi v\rangle_{k}\right),$

or alternatively,

$\displaystyle\mathcal{U}^{-1}\left(i\omega_{n}\right)=$ $\displaystyle\
\langle v^{-1}\rangle_{k}-(\langle v^{-1}Wv^{-1}\rangle_{k}$ (47)
$\displaystyle\quad\quad\quad\quad\quad-\langle v^{-1}W\rangle_{k}\langle
W\rangle_{k}^{-1}\langle Wv^{-1}\rangle_{k}).$

We can further replace the band dispersion $\epsilon_{k}$ used in Eq. (46) by
the lattice Dyson equation
$\epsilon_{k}\to-G^{-1}(k,i\omega_{n})+i\omega_{n}+\mu-\Sigma(k,i\omega_{n})$.
In the case of the screened interaction, we replace the Coulomb interaction
$v_{k}$ used in Eq. (47) by the bosonic Dyson equation $v_{k}\to
W^{-1}(k,i\omega_{n})+\Pi(k,i\omega_{n})$ and replace the charge
susceptibility $\chi$ by the screened interaction $W$ using the identity
$W=v-v\chi v$. After some manipulations, we obtain the generalized EDMFT self-
consistency equationsBackes et al. (2020)

$\displaystyle\mathcal{G}^{-1}\left(i\omega_{n}\right)$
$\displaystyle=\quad\langle
G\rangle^{-1}+\langle\Sigma\rangle+\Delta_{\text{cor}},$ (48)
$\displaystyle\mathcal{U}(i\nu_{n})$ $\displaystyle=\left[\langle
W\rangle^{-1}+\langle\Pi\rangle\right]^{-1}+\mathcal{U}_{\text{cor}},$

and the alternative form

$\displaystyle\mathcal{U}^{-1}(i\nu_{n})$ $\displaystyle=\left[\langle
W\rangle^{-1}+\langle\Pi\rangle\right]+{\tilde{\mathcal{U}}}^{-1}_{\text{cor}},$
(49)

where the correction terms are given by

$\displaystyle\Delta_{\text{cor}}$ $\displaystyle=-\langle\Sigma
G\Sigma\rangle_{k}+\langle\Sigma G\rangle_{k}\langle G\rangle_{k}^{-1}\langle
G\Sigma\rangle_{k}-2\langle\Sigma\rangle_{k}+\langle\Sigma G\rangle_{k}\langle
G\rangle_{k}^{-1}+\langle G\rangle_{k}^{-1}\langle G\Sigma\rangle_{k},$ (50)
$\displaystyle\mathcal{U}_{\text{cor}}$ $\displaystyle=-\langle\Pi
W\rangle_{k}[\langle\Pi\rangle_{k}+\langle\Pi W\Pi\rangle_{k}]^{-1}\langle
W\Pi\rangle_{k}+\langle\Pi\rangle_{k}\langle
W\rangle_{k}[\langle\Pi\rangle_{k}+\langle\Pi\rangle_{k}\langle
W\rangle_{k}\langle\Pi\rangle_{k}]^{-1}\langle
W\rangle_{k}\langle\Pi\rangle_{k},$
$\displaystyle{\tilde{\mathcal{U}}}^{-1}_{\text{cor}}$
$\displaystyle=-\langle\Pi W\Pi\rangle_{k}+\langle\Pi W\rangle_{k}\langle
W\rangle^{-1}_{k}\langle W\Pi\rangle_{k}-2\langle\Pi\rangle_{k}+\langle\Pi
W\rangle_{k}\langle W\rangle_{k}^{-1}+\langle W\rangle_{k}^{-1}\langle
W\Pi\rangle_{k}.$

In the standard DMFT Georges et al. (1996) or EDMFT Sun and Kotliar (2002)
treatment, the self-energy $\Sigma$ and the polarization $\Pi$ are momentum
independent and calculated by the impurity solver:
$\Sigma=\Sigma_{\text{imp}}$ and $\Pi=\Pi_{\text{imp}}$. In this case, we have
$\langle G\Sigma\rangle=\langle G\rangle\Sigma_{\text{imp}}$ and $\langle
W\Pi\rangle=\langle W\rangle\Pi_{\text{imp}}$ and the correction terms defined
by Eq. (50) vanish. In this case the generalized schemes reduce to the
standard (local) self-consistency loop and there is no difference between the
two alternative procedures (Eqs. (48) and (49)). In the $GW$+EDMFT scheme,
however, one combines the local impurity contributions
$(\Sigma_{\text{imp}},\Pi_{\text{imp}})$ with the nonlocal parts of the $GW$
contributions $\Sigma^{GW}(k)$ and $\Pi^{GW}(k)$. In this case, the correction
terms defined in Eq. (50) will not cancel.

The standard (local) treatment of $GW$+(E)DMFT used in previous works Sun and
Kotliar (2002); Biermann et al. (2003); Ayral et al. (2013); Nilsson et al.
(2017) considers only the local part of the $k$-dependent self-energy and
polarization in the self-consistency equations (Eqs. (48) and (49)) which fix
the dynamical mean fields. This amounts to neglecting the correction terms.
While the local procedure becomes exact in the infinite-connectivity limit,
and thus is consistent with the first two approximations, it can produce non-
causal Weiss fields (especially $\mathcal{U}$), as discussed in Refs. Nilsson
et al., 2017; Backes et al., 2020 and the main text.

To summarize, for the EDMFT procedure or the standard (local) $GW$+EDMFT
procedure, the correction terms defined by Eq. (50) cancel, and one obtains
the local EDMFT-type self-consistency equations,

$\displaystyle\mathcal{G}^{-1}=G_{\mathrm{loc}}^{-1}+\Sigma_{\mathrm{loc}},$
(51) $\displaystyle\mathcal{U}^{-1}=W_{\mathrm{loc}}^{-1}+\Pi_{\mathrm{loc}}.$

The $GW$+EDMFT variant with the correction terms in the self-consistency
equations, Eqs. (48) or (49), will be referred to as the causal self-
consistency scheme, since as shown in Ref. Backes et al., 2020, it results in
causal dynamical mean fields. This causal scheme is relevant only in the case
of $k$-dependent self-energies and polarizations, i.e. if EDFMT is combined
with the $GW$ approximation (or some other scheme which can provide non-local
components).

In this work, we only implement and test the first of the two alternative
forms, i.e. the causal scheme defined by Eq. (48).

## References

  * Mott (1968) N. F. Mott, Rev. Mod. Phys. 40, 677 (1968), URL https://link.aps.org/doi/10.1103/RevModPhys.40.677.
  * Imada et al. (1998) M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998), URL https://link.aps.org/doi/10.1103/RevModPhys.70.1039.
  * Hedin (1999) L. Hedin, J. Phys.: Condens. Matter 11, R489 (1999), URL https://iopscience.iop.org/article/10.1088/0953-8984/11/42/201.
  * Werner and Casula (2016) P. Werner and M. Casula, J. Phys.: Condens. Matter 28, 383001 (2016), eprint 1602.00584, URL http://arxiv.org/abs/1602.00584.
  * Hansmann et al. (2013) P. Hansmann, T. Ayral, L. Vaugier, P. Werner, and S. Biermann, Phys. Rev. Lett. 110, 166401 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.110.166401.
  * Ayral et al. (2012) T. Ayral, P. Werner, and S. Biermann, Phys. Rev. Lett. 109, 226401 (2012), URL https://link.aps.org/doi/10.1103/PhysRevLett.109.226401.
  * Huang et al. (2014) L. Huang, T. Ayral, S. Biermann, and P. Werner, Phys. Rev. B 90, 195114 (2014), URL https://link.aps.org/doi/10.1103/PhysRevB.90.195114.
  * Sun and Kotliar (2002) P. Sun and G. Kotliar, Phys. Rev. B 66, 085120 (2002), URL https://link.aps.org/doi/10.1103/PhysRevB.66.085120.
  * Biermann et al. (2003) S. Biermann, F. Aryasetiawan, and A. Georges, Phys. Rev. Lett. 90, 086402 (2003), URL https://link.aps.org/doi/10.1103/PhysRevLett.90.086402.
  * Aryasetiawan et al. (2004) F. Aryasetiawan, S. Biermann, and A. Georges, arXiv:cond-mat/0401626 (2004), eprint cond-mat/0401626, URL http://arxiv.org/abs/cond-mat/0401626.
  * Biermann et al. (2004) S. Biermann, F. Aryasetiawan, and A. Georges, arXiv:cond-mat/0401653 (2004), eprint cond-mat/0401653, URL http://arxiv.org/abs/cond-mat/0401653.
  * Petocchi et al. (2021) F. Petocchi, V. Christiansson, and P. Werner, arXiv:2106.03689 [cond-mat] (2021), eprint 2106.03689, URL http://arxiv.org/abs/2106.03689.
  * Ewald (1921) P. P. Ewald, Annalen der Physik 369, 253 (1921), URL https://onlinelibrary.wiley.com/doi/abs/10.1002/andp.19213690304.
  * Harris (1998) F. E. Harris, International Journal of Quantum Chemistry 68, 385 (1998), URL https://onlinelibrary.wiley.com/doi/abs/10.1002/%28SICI%291097-461X%281998%2968%3A6%3C385%3A%3AAID-QUA2%3E3.0.CO%3B2-R.
  * Backes et al. (2020) S. Backes, J.-H. Sim, and S. Biermann, arXiv:2011.05311 [cond-mat] (2020), eprint 2011.05311, URL http://arxiv.org/abs/2011.05311.
  * Nilsson et al. (2017) F. Nilsson, L. Boehnke, P. Werner, and F. Aryasetiawan, Phys. Rev. Materials 1, 043803 (2017), URL https://link.aps.org/doi/10.1103/PhysRevMaterials.1.043803.
  * Wang et al. (2020) Y. D. Wang, W. L. Yao, Z. M. Xin, T. T. Han, Z. G. Wang, L. Chen, C. Cai, Y. Li, and Y. Zhang, Nat Commun 11, 4215 (2020), URL https://www.nature.com/articles/s41467-020-18040-4.
  * Ligges et al. (2018) M. Ligges, I. Avigo, D. Golež, H. U. R. Strand, Y. Beyazit, K. Hanff, F. Diekmann, L. Stojchevska, M. Kalläne, P. Zhou, et al., Phys. Rev. Lett. 120, 166401 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.166401.
  * Negele and Orland (1988) J. W. Negele and H. Orland, _Quantum Many-Particle Systems_ (Basic Books, 1988), ISBN 978-0-201-12593-1.
  * Ayral et al. (2013) T. Ayral, S. Biermann, and P. Werner, Phys. Rev. B 87, 125149 (2013), URL https://link.aps.org/doi/10.1103/PhysRevB.87.125149.
  * Hedin (1965) L. Hedin, Phys. Rev. 139, A796 (1965), URL https://link.aps.org/doi/10.1103/PhysRev.139.A796.
  * Petocchi et al. (2020a) F. Petocchi, F. Nilsson, F. Aryasetiawan, and P. Werner, Phys. Rev. Research 2, 013191 (2020a), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.013191.
  * Georges et al. (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996), URL https://link.aps.org/doi/10.1103/RevModPhys.68.13.
  * Si and Smith (1996) Q. Si and J. L. Smith, Phys. Rev. Lett. 77, 3391 (1996), URL https://link.aps.org/doi/10.1103/PhysRevLett.77.3391.
  * Aryasetiawan and Gunnarsson (1998) F. Aryasetiawan and O. Gunnarsson, Rep. Prog. Phys. 61, 237 (1998), eprint cond-mat/9712013, URL http://arxiv.org/abs/cond-mat/9712013.
  * Pavarini et al. (2011) E. Pavarini, E. Koch, D. Vollhardt, and A. Lichtenstein, eds., _The LDA+DMFT Approach to Strongly Correlated Materials_ , vol. 1 of _Schriften Des Forschungszentrums Jülich Reihe Modeling and Simulation_ (Forschungszentrum Jülich, Jülich, 2011), ISBN 978-3-89336-734-4, URL http://www.cond-mat.de/events/correl11.
  * Luttinger and Ward (1960) J. M. Luttinger and J. C. Ward, Phys. Rev. 118, 1417 (1960), URL https://link.aps.org/doi/10.1103/PhysRev.118.1417.
  * Almbladh et al. (1999) C.-O. Almbladh, U. V. Barth, and R. V. Leeuwen, Int. J. Mod. Phys. B 13, 535 (1999), URL https://www.worldscientific.com/doi/abs/10.1142/S0217979299000436.
  * Werner and Millis (2006) P. Werner and A. J. Millis, Phys. Rev. B 74, 155107 (2006), URL https://link.aps.org/doi/10.1103/PhysRevB.74.155107.
  * Werner and Millis (2010) P. Werner and A. J. Millis, Phys. Rev. Lett. 104, 146401 (2010), URL https://link.aps.org/doi/10.1103/PhysRevLett.104.146401.
  * Wu et al. (2021) Z. Wu, K. Bu, W. Zhang, Y. Fei, Y. Zheng, J. Gao, X. Luo, Z. Liu, Y.-p. Sun, and Y. Yin, arXiv:2105.08663 [cond-mat] (2021), eprint 2105.08663, URL http://arxiv.org/abs/2105.08663.
  * Vidberg and Serene (1977) H. J. Vidberg and J. W. Serene, J Low Temp Phys 29, 179 (1977), URL https://doi.org/10.1007/BF00655090.
  * Boehnke et al. (2016) L. Boehnke, F. Nilsson, F. Aryasetiawan, and P. Werner, Phys. Rev. B 94, 201106 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.94.201106.
  * Petocchi et al. (2020b) F. Petocchi, V. Christiansson, F. Nilsson, F. Aryasetiawan, and P. Werner, Phys. Rev. X 10, 041047 (2020b), URL https://link.aps.org/doi/10.1103/PhysRevX.10.041047.
  * Cho et al. (2016) D. Cho, S. Cheon, K.-S. Kim, S.-H. Lee, Y.-H. Cho, S.-W. Cheong, and H. W. Yeom, Nat Commun 7, 10453 (2016), URL https://www.nature.com/articles/ncomms10453.
  * Eckstein and Werner (2011) M. Eckstein and P. Werner, Phys. Rev. B 84, 035122 (2011), URL https://link.aps.org/doi/10.1103/PhysRevB.84.035122.
  * Lee et al. (2019) S.-H. Lee, J. S. Goh, and D. Cho, Phys. Rev. Lett. 122, 106404 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.122.106404.
  * Pasquier and Yazyev (2021) D. Pasquier and O. V. Yazyev, arXiv:2108.11277 [cond-mat] (2021), eprint 2108.11277, URL http://arxiv.org/abs/2108.11277.
  * Golež et al. (2017) D. Golež, L. Boehnke, H. U. R. Strand, M. Eckstein, and P. Werner, Phys. Rev. Lett. 118, 246402 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.118.246402.
  * Golež et al. (2019) D. Golež, M. Eckstein, and P. Werner, Phys. Rev. B 100, 235117 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.235117.
  * Wilson et al. (1975) J. Wilson, F. Di Salvo, and S. Mahajan, Advances in Physics 24, 117 (1975), URL https://doi.org/10.1080/00018737500101391.
  * Fazekas and Tosatti (1979) P. Fazekas and E. Tosatti, Philosophical Magazine B 39, 229 (1979), URL https://doi.org/10.1080/13642817908245359.
  * Sipos et al. (2008) B. Sipos, A. F. Kusmartseva, A. Akrap, H. Berger, L. Forró, and E. Tutiš, Nature Mater 7, 960 (2008), URL http://www.nature.com/articles/nmat2318.
  * Murakami et al. (2015) Y. Murakami, P. Werner, N. Tsuji, and H. Aoki, Phys. Rev. B 91, 045128 (2015), URL https://link.aps.org/doi/10.1103/PhysRevB.91.045128.
  * Werner et al. (2019) P. Werner, M. Eckstein, M. Müller, and G. Refael, Nat Commun 10, 5556 (2019), URL https://www.nature.com/articles/s41467-019-13557-9.
  * Metzner and Vollhardt (1989) W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, 324 (1989), URL https://link.aps.org/doi/10.1103/PhysRevLett.62.324.
  * Hubbard (1959) J. Hubbard, Phys. Rev. Lett. 3, 77 (1959), URL https://link.aps.org/doi/10.1103/PhysRevLett.3.77.
  * Hubbard and Flowers (1964) J. Hubbard and B. H. Flowers, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 281, 401 (1964), URL https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1964.0190.
  * Müller-Hartmann (1989) E. Müller-Hartmann, Z. Physik B - Condensed Matter 74, 507 (1989), URL https://doi.org/10.1007/BF01311397.