Influence of Fock exchange in combined many-body perturbation and dynamical mean field theory
Thomas Ayral Silke Biermann Philipp Werner Lewin Boehnke




Institut de Physique Thorique (IPhT), CEA, CNRS, UMR 3681, 91191 Gif-sur-Yvette, France
Department of Physics and Astronomy, Rutgers University, Piscataway, NJ 08854, USA

Centre de Physique Thorique, Ecole Polytechnique, CNRS, Universit Paris Saclay, 91128 Palaiseau, France

Department of Physics, University of Fribourg, 1700 Fribourg, Switzerland

Department of Physics, University of Fribourg, 1700 Fribourg, Switzerland


In electronic systems with long-range Coulomb interaction, the nonlocal Fock exchange term has a band-widening effect. While this effect is included in combined many-body perturbation theory and dynamical mean field theory schemes, it is not taken into account in standard extended DMFT (EDMFT) calculations.
Here, we include this instantaneous term in both approaches and investigate its effect on the phase diagram and dynamically screened interaction. 
We show that the largest deviations between previously presented EDMFT and GW+EDMFT results originate from the nonlocal Fock term, and that the quantitative differences are especially large in the strong-coupling limit.
Furthermore, we show that the charge-ordering phase diagram obtained in GW+EDMFT methods for moderate interaction values is very similar to the one predicted by dual boson methods that include the fermion-boson or four-point vertex.




  Influence of Fock exchange in combined many-body perturbation and dynamical mean field theory

    Lewin Boehnke
    
=================================================================================================




Dynamical Mean Field Theory<cit.> (DMFT) 
self-consistently maps a correlated Hubbard lattice problem with local interactions onto an effective impurity problem 
consisting of a correlated orbital hybridized with a noninteracting 
fermionic bath. If the bath is integrated out, one obtains an impurity action with retarded hoppings. 
Extended dynamical mean field theory<cit.>
(EDMFT)  
extends the DMFT idea to systems with long-range interactions. It does so by mapping a lattice problem with long-range interactions onto an effective impurity model with
self-consistently determined fermionic and bosonic baths, or, in the action formulation, an impurity model with retarded hoppings and retarded interactions. 


While EDMFT captures dynamical screening effects and charge-order instabilities, it has  been found to suffer from qualitative shortcomings in finite dimensions. 
For example, the charge susceptibility computed in EDMFT does not coincide with the derivative of the average charge with respect to a small applied field<cit.>, nor does it obey local charge conservation rules<cit.> essential for an adequate description of collective modes such as plasmons.

The EDMFT formalism has an even more basic deficiency: since it is based on a local approximation to the self-energy, it does not include even the first-order nonlocal interaction term, the Fock term. The combined GW+EDMFT 
<cit.> scheme corrects this by supplementing the local self-energy from EDMFT with the nonlocal part of the GW diagram, where G is the interacting Green's function and W the fully screened interaction. 
Indeed, the nonlocal Fock term “[Gv]^nonloc” is
included in the nonlocal “[GW]^nonloc” diagram. As described in more detail in Ref. Ayral2013 (see also the appendix of Ref. Golez2017), the GW+EDMFT method 
is formally obtained by constructing an energy functional of G and W, the Almbladh<cit.> functional Ψ, and by approximating Ψ as a sum of two terms, one containing all local diagrams (corresponding to EDMFT), the other containing the simplest nonlocal correction (corresponding to the GW approximation<cit.>). This functional construction rules out double-counting of local terms in the self-energy and polarization<cit.>.
Even though it has been introduced under the name GW+DMFT
<cit.> in the literature, we will denote this full scheme by GW+EDMFT
to emphasize that it is based on the EDMFT formalism, and to distinguish it 
from simplified implementations without two-particle self-consistency, 
which have appeared in the literature (and which we will denote in the
following as GW+DMFT).


In a recent implementation of the GW+EDMFT method, Ref. Ayral2013, and related papers<cit.>, the nonlocal Fock term was omitted.<cit.>  
Here, we explore and 
highlight the role of this term and its interplay with the local correlations. We quantify the band-widening effect of the Fock term and study the consequences of its presence or absence on various observables, and on the charge-order phase boundary. 
Our self-consistent implementation goes beyond previous studies of the effect of the Fock exchange
in realistic calculations, where it was studied systematically
within GW <cit.> and GW+DMFT,<cit.> 
albeit not in a self-consistent way.


The manuscript is organized as follows:
In section <ref>, we recap the GW+EDMFT equations with special emphasis on the Fock term and make general statements about the expected impact.
In section <ref>, we show explicit results for the effective renormalization of the band structure by the instantaneous Fock contribution within GW+EDMFT, followed by systematic comparisons 
with the results of simplified 
formalisms in section <ref>.
Section <ref> discusses the role of the Fock term in the Mott-insulating phase, where it stays relevant up to very large values of the on-site interaction.
Finally, in section <ref>, we compare our results with 
results obtained within 
the recent dual boson method.



§ FORMALISM


We aim at solving the extended Hubbard model on the two-dimensional square lattice 
by constructing an effective impurity problem that gives the local part of the self-energy Σ and polarization P, and a diagrammatic expansion in their nonlocal components.
The model is defined by
the Hamiltonian

    ℋ=-∑_ijt_ijc^†_ic_j+1/2∑_ij v_ijn_in_j-μ∑_in_i.

Here,
t_ij are the real-space hopping matrix elements, c_i^(†) the electronic annihilator (creator) on site i, v_ij the Coulomb interaction and n_i=c^†_ic_i.
We will
restrict ourselves to models with hoppings and interactions between nearest neighbors and next nearest neighbors only,

    t_ij=   tδ_<ij>+t'δ_<<ij>>,
    
      v_ij=   Uδ_ij+Vδ_<ij>+V'δ_<<ij>>,

where δ_ij is the usual Kronecker delta, δ_<ij>(resp. δ_<<ij>>) is 1 for i and j nearest neighbors (resp.  next-nearest neighbors, along the diagonal of the square lattice) and 0 otherwise.

This results in the Fourier transforms 

    ε_k⃗=+   2t(cos(k_x)+cos(k_y))
       +2t'(cos(k_x+k_y)+cos(k_x-k_y))

and

    v_q⃗=U   +2V(cos(q_x)+cos(q_y))
       +2V'(cos(q_x+q_y)+cos(q_x-q_y)).


The full expression for the self-energy in the GW+EDMFT approximation is

    Σ(k⃗,iω_n)=Σ_imp(iω_n)+Σ_GW_c^nonloc(k⃗,iω_n)+Σ_F^nonloc(k⃗).

The last two terms correspond to the nonlocal part of the GW self-energy. They can be expressed as a function of imaginary time τ 
 and momentum k⃗ as follows:


    Σ_GW_c^nonloc(k⃗,τ)     =     -∑_q⃗G_q⃗+⃗k⃗(τ)W^c_q⃗(τ)
           +[∑_q⃗G_q⃗+⃗k⃗(τ)W_q⃗^c(τ)]_loc,
    Σ_F^nonloc(k⃗)     =     -∑_q⃗G_q⃗+⃗k⃗()v_q⃗
           +[∑_q⃗G_q⃗+⃗k⃗()v_q⃗]_loc.

Fourier transformations between τ and fermionic [bosonic] Matsubara frequencies iω_n=i(2n+1)π/β [iν_m=i2mπ/β] are assumed where needed (β denotes the inverse temperature).
The “loc” suffix denotes a sum over the first Brillouin zone.

    G_k⃗(iω_n)=(iω_n+μ-ε_k⃗-Σ(k⃗,iω_n))^-1

is the interacting lattice Green's function, and W^c is defined as

    W_q⃗^c(iν_m)≡v_q⃗/1-v_q⃗P_q⃗(iν_m)-v_q⃗,

with P the polarization function.
All results are given in units of D=4|t| (which is the half bandwidth when t'=0), and the momentum discretization is N_k=32×32
points in the first Brillouin zone, 
unless otherwise stated. 
We use the original formulation of the GW+EDMFT scheme,<cit.> corresponding
– within a functional formulation – to a Hubbard-Stratonovich decoupling
of the full interaction term, dubbed “HS-UV decoupling” in
Ref. Ayral2013.
As argued there, this choice has the advantage that it treats local and 
nonlocal interactions on the same footing.




The nonlocal Fock term of Eq. (<ref>), which is real-valued and instantaneous, renormalizes the bandwidth. It can become quite large and momentum-dependent.
This is illustrated in Fig. <ref> for the parameters U=2, V=0.4 (and t'=0, V'=0, which is assumed in the following if not explicitly stated otherwise).
The figure also indicates that for the case of nearest-neighbor hopping and interaction, the nonlocal Fock term can be exactly absorbed into the bare dispersion Eq. (<ref>) by defining a U- and V-dependent hopping

    t̃(t',U,V,V')=t+δ t(t',U,V,V').

This can be understood by looking at the real-space representation of the Fock term, Eq. (<ref>),

    Σ_F^nonloc_ij
    =     -G_ij()v_ij + G_ii()v_iiδ_ij
    
      =     -G_<ij>()Vδ_<ij>
       -G_<<ij>>()V'δ_<<ij>>,

where the notation <ij> [<<ij>>] denotes
a restriction to nearest-neighbor [next-nearest-neighbor] interactions.

Thus, for the case of Fig. <ref>, where V'=0, the Fock term enters Eq. (<ref>) as a renormalization of the nearest-neighbor hopping by

    δ t(t',U,V,V')=-G_<ij>()V .

The Green's function term, which is closely related to the occupation number, has an implicit dependence on all parameters of the lattice problem.
In the presence of a next-nearest-neighbor interaction, t' also gets renormalized according to

    δ t'(t',U,V,V')=-G_<<ij>>()V' .

Note that such a term breaks the particle-hole symmetry of the lattice. In the particle-hole symmetric case with t'=0 and half filling, the next-nearest-neighbor occupation term vanishes: G_<<ij>>()=0.



§ EFFECTIVE BANDSTRUCTURE





In the simplest case of nearest-neighbor interaction V, the nearest-neighbor hopping renormalization δ t determines the band widening (see Eq. (<ref>)). Since the half bandwidth is D=4|t|, the widening will be δ D=4δ t.
Figure <ref> illustrates this effect throughout the homogeneous part of the phase diagram for the particle-hole symmetric (t'=0) half filled case.
The most obvious feature is the increase with V, that is expected from Eq. (<ref>) and the decrease close to the Mott-insulating phase.
δ t nonetheless remains significant even at very high values of U, a property that will be further investigated in Section <ref>.

Away from half filling, where the Mott-insulating phase does not exist, the corresponding suppression of δ t disappears, but otherwise the dependence on V and U is very similar to the half-filled case, see bottom panel of Fig. <ref>.


To study a model with broken particle-hole symmetry, we introduce a nearest-neighbor hopping t'=t/√(2) and fix the filling at <n>=0.8 as well as the nearest-neighbor interaction V. We then calculate δ t as a function of the next nearest nearest-neighbor interaction V'. As shown in Fig. <ref>, 
the main effect on the hopping renormalization comes from the (essentially linear) dependence on V.
The qualitative effect of V' is to slightly reduce the renormalization.

In order to make the connection to realistic electronic structure
calculations, we note as a side remark that there
the situation is slightly more subtle. The band-widening effect is indeed relative to the reference point. Let us consider three reference Hamiltonians, 

(i) the Kohn-Sham Hamiltonian H_KS
of Density Functional Theory (DFT),

(ii) the Hartree Hamiltonian H_0 = H_KS - V[v_xc(r⃗)], where v_xc(r⃗) denotes the Kohn-Sham exchange-correlation potential, which is local in the electronic structure sense
(that is, “local” denotes a function depending only on one space variable f(r⃗), while “nonlocal” denotes a function depending on two variables f(r⃗,r⃗'⃗)),

(iii) the nonlocal exchange Hamiltonian H_xc^F = H_0 + V[v_xc^F(r⃗,r⃗'⃗)] (where v_xc^F(r⃗,r⃗'⃗) denotes an exchange-correlation potential including “nonlocal” Fock exchange).

Then the hierarchy of the bandwidths in a metallic system is

    H_0 > H_xc^F > H_KS.


Thus, H_xc^F indeed widens the band with respect to density functional theory calculations. The question of the relative bandwidth changes
thus implies a question on the starting band structure.
We refer the interested reader to Ref. Hirayama2015 
for a systematic construction of explicit low-energy many-body
Hamiltonians.

Here, we only comment on the specific point of the Hartree and Fock
terms, in order to put our work on the extended Hubbard Hamiltonian
into perspective with respect to realistic electronic structure calculations.
Indeed, as argued in Ref. Hirayama2015, in realistic 
electronic structure calculations one needs to avoid double counting 
of interactions at the one- and two-particle level. 
Let us consider first the case of the Hartree terms: standard electronic
structure techniques (e.g. a DFT calculation) produce a band structure
including the Hartree contribution. This one-body potential contribution
is then already part of the effective hopping parameter determined
from this band structure.
Ref. Hirayama2015 explains how to avoid double counting 
by including – at the level of the many-body calculation – only terms 
beyond Hartree. Here, we do not need to address this point in detail, since
a Hartree term included in the model calculation would cancel out
with the corresponding shift of the chemical potential, since the
particle number is eventually determining the energetic level of the
single-orbital included in the present model.

Let us now move to the analogous question for the Fock term: one may examine the relevance of excluding it at the level of the many-body calculation, and keeping it at the level of the electronic structure calculation instead.
The answer is based on several elements: The first point to note
is that standard DFT calculations do treat exchange in a local
approximation (where “local” here means again “local in the electronic 
structure sense”, see above), which
relies on an error cancellation effect with part of the correlation
contribution (see e.g. Ref. VanRoekeghem2014a) and is not relevant here. 
The next question is therefore: why not start from a
Hartree-Fock calculation in the continuum in the full energy range
of the Coulomb Hamiltonian? Such a treatment would neglect the 
crucial screening of the bare interaction by high-energy degrees of
freedom (typically, matrix elements of the bare Coulomb interaction
in the relevant Wannier functions are of the order of several tens
of electron volts, while the effective Hubbard interactions are usually a few electron volts).
Therefore, what is relevant here is indeed the exchange term calculated
using the effective bare interaction of the low-energy Hamiltonian.
For realistic electronic structure calculations, this interaction
should correspond to a partially screened interaction, where screening
by high-energy degrees of freedom is taken into account (as done e.g.
in the screened exchange + DMFT scheme 
<cit.>).
We refer the interested reader to Refs. Hirayama2015, VanRoekeghem2014, VanRoekeghem2014a for details.




§ SIMPLIFIED VARIANTS OF GW+EDMFT



In the following, we study the effect of this V-dependent bandwidth renormalization on local observables as well as the critical value of the nearest-neighbor repulsion for the transition into the charge-ordered phase.
We will call “GW_c+EDMFT” the formula implemented in Ref. Ayral2013 (which contains only the GW_c
term, see Ref. Ayral2016a), and “GW+EDMFT” the formula with the self-energy expression
(<ref>). For comparison, we also show results for “Gv+EDMFT”,
 a scheme where Σ is the sum of the impurity self-energy
and of the nonlocal Fock term only (the first and third terms of Eq. (<ref>)). In all three schemes, the polarization is
the sum of the impurity polarization with the nonlocal part of the
interacting bubble, as described in Ref. Ayral2013,
steps (5) (a) and (5) (b) of section V.






In Fig. <ref>, we plot the self-energy and polarization
obtained from the three schemes at different momenta, and we compare
the results to the local EDMFT self-energy and
polarization. All these results are for half-filling. One can observe the following
trends:

(i) While the imaginary part of the self-energy in GW_c+EDMFT 
is larger than in EDMFT, the opposite is true for GW+EDMFT, i.e.,
the GW+EDMFT self-energy is less correlated than the EDMFT
self-energy. 

(ii) The GW+EDMFT result is more strongly correlated than Gv+EDMFT. 

(iii) At small Matsubara frequencies, the polarization is overall larger
in the GW/Gv+EDMFT method than in the GW_c+EDMFT. 

The trend in the self-energy (i.e. less correlated in GW+EDMFT than
EDMFT) can be understood easily from the broadening effect of the
nonlocal Fock term on the band: when the (effective) bandwidth gets
larger, so does the polarization P, and hence screening effects
are more important, interactions are more screened and as a result,
the imaginary part of the Matsubara self-energy is smaller in absolute
magnitude.






Less trivial is the comparison between GW+EDMFT and Gv+EDMFT. Here, the band-widening effect
is included in both calculations, and it turns out that the additional nonlocal GW contributions to the self-energy lead to stronger correlations.
This is consistent with the conclusions of  Ref. Ayral2013,
which compared GW_c+EDMFT to EDMFT. 

In the top panels of Fig. <ref>, we replot panels (a) and (b) of Fig. 15 of Ref. Ayral2013, and we show, in the bottom panels, the same observables for U=3, V=1. 
One sees that the deviation of the GW+EDMFT and Gv+EDMFT results from 
the EDMFT result is very small for U=2, V=0.4, but sizable for larger interaction values
(U=3, V=1).




In Fig. <ref>, we 
plot the corresponding spectral functions (the EDMFT and GW_c+EDMFT results are identical to Fig. 2(b) of Ref. Ayral2012).
As a logical consequence of the above observations, the GW+EDMFT and
Gv+EDMFT spectra are close to each other and slightly less correlated
than the EDMFT spectrum, in the sense that the integrated weight of
the quasiparticle peak is larger in those methods.

We next consider the phase diagram in the U-V plane.
In Fig. <ref>, we plot the dependence
of the inverse charge susceptibility χ^-1_q⃗=π,π(iν_m=0) on the nearest-neighbor repulsion V, with χ defined by

    χ_q⃗(iν_m)=-Π_q⃗(iν_m)/1-Π_q⃗(iν_m)v_q⃗.

When the inverse susceptibility vanishes, the
charge susceptibility diverges, signaling a transition to a charge-ordered
phase with a checkerboard pattern. The corresponding phase diagram is shown in Fig. <ref>, where we plot the results from Fig. 5 of Ref. Ayral2013 together with the phase boundaries for the GW+EDMFT and Gv+EDMFT methods.

At low and intermediate U, GW+EDMFT and GW_c+EDMFT yield quantitatively similar critical nonlocal interactions V_c for the transition to the charge-ordered phase over a wide range of the local interactions. More importantly, they capture the expected GW behavior at low U that EDMFT misses due to its local self-energy.
In the strong-coupling limit, the value of V_c is substantially reduced (middle and bottom panels) when going from EDMFT to GW+EDMFT or even only
Gv+EDMFT (GW_c+EDMFT is very close to EDMFT).



§ INSIDE THE MOTT PHASE


In the large-interaction regime of the phase diagram (Fig. <ref>), the nonlocal Fock term has a significant effect. The schemes which lack this instantaneous contribution, EDMFT and GW_c+EDMFT, yield a larger and steeper phase boundary than the schemes that take the Fock term into account (Gv+EDMFT and GW+EDMFT).

The exact phase boundary in the Mott phase is difficult to predict a priori. 
It can be computed in the classical (t_ij→ 0) and zero-temperature limit of the extended Hubbard model by exact Monte-Carlo simulations, as e.g. in Ref. Pawlowski2006, and is given by the analytical expression 

    V_c^int=U/4,
 
where 4 corresponds to the number of nearest neighbors.
This line is plotted as a dashed grey line in Fig. <ref>. This result can be obtained by a simple comparison between the interaction energies of the Mott-insulating phase and of the checkerboard phase. In the full-fledged model, finite temperature and quantum tunneling have to be taken into account. In the low-temperature regime (T=0.01) of Fig. <ref>, the deviation between the classical solution and the solution to the full quantum problem comes mostly from the quantum tunneling kinetic term. 

To guess the influence of the quantum tunneling term, one may observe that the effect of temperature in the classical problem is to enhance the value of V_c, i.e to disfavor the charge-ordered phase over the Mott phase.<cit.> Since the quantum tunneling (hopping) has a physical effect similar to temperature in classical systems,<cit.>  
namely to delocalize the particles 
one may speculate that it will also lead to a higher V_c in the quantum case. 

In fact, this feature is present by construction in the EDMFT and GW+EDMFT schemes. The denominator in the susceptibility (Eq. (<ref>)) imposes that the charge-ordering transition should occur for negative values of v_q⃗, since the polarization Π_q⃗ is always negative (for the parameters studied here). For the square lattice, this implies V_c>U/4, which is the classical energy estimate (Eq. (<ref>)). Therefore, by construction, in GW+EDMFT schemes, introducing hopping on the lattice will always favor the disordered phase. 
This is indeed what is seen in all variants. We also observe that the method including most diagrams, GW+EDMFT, has a phase boundary which is much closer to the classical limit than the comparatively cruder EDMFT approximation.

In order to gain a better qualitative understanding of the large-U behavior, we have performed an analytical self-consistent estimation of the value of the band-widening effect δ t (defined in Eq. (<ref>)) coming from the Fock term. Approximating the self-energy as the sum of the atomic limit (in the spirit of the Hubbard-I approximation<cit.>) and of the Fock self-energy, as described in more detail in Appendix <ref>, we obtain

    δ t=tV/2U-V.

Thus δ t may become arbitrarily
large if the nonlocal interaction coefficient exceeds twice the value of the local
one. However, even disregarding the fact that the generic case is
certainly the opposite one (local interactions in general exceed nonlocal
ones), one should be aware of the fact that in that case the Hubbard
I approximation, which is justified in the strong coupling limit,
would no longer be appropriate.

By inspecting the phase diagram in Fig. <ref>, we can parametrize the phase boundary in the large-U limit as a constant slope, i.e   V_c=α U+β. [Within the U-range that we can simulate (we performed measurements up to U=8.0), we can estimate α=1.25 and β=-2.5.] With this parametrization,
we obtain 

    δ t_c=t(α U+β)/2U-α U-β.

Hence, the bandwidth 
renormalization (proportional to δ t) 
stays relevant in the vicinity of the charge-ordering transition 
even at large U.




§ BEYOND GW+EDMFT: COMPARISON TO DUAL BOSONS AND TRILEX



As mentioned in the introduction, the EDMFT formalism suffers from certain conceptual problems, such as the lack of thermodynamic consistency and an unreliable description of collective modes. These shortcomings are alleviated in the recently developed dual boson (DB) method,<cit.> which, in its full-fledged implementation,<cit.> computes the susceptibility<cit.> after resumming an infinite number of ladder diagrams built from local impurity four-leg vertices. 

These four-leg vertices, which are also central to the Dynamical Vertex Approximation<cit.> (which was recently shown to be a simplified version of QUADRILEX, a method consisting in an atomic approximation of the four-particle irreducible functional<cit.>), 
can nonetheless only be obtained at a considerable computational expense and require a proper parametrization and treatment of their asymptotic behavior.<cit.> Consequently, it is not
possible to use them routinely in
multi-orbital calculations (see however <cit.>) 
and lightweight improvements on EDMFT, especially with realistic applications in mind, are desirable. Recent attempts to forgo the computation of four-leg vertices include the TRILEX method<cit.> and simplified dual boson schemes such as DB+GW or DB+GWγ.<cit.> Whether they retain the abovementioned conserving properties, however, is yet unclear.


In fact, the results obtained in the simplified “dual” approaches that include at least the electron-boson vertex γ are similar to those obtained by the GW+EDMFT method, which is conceptually and practically simpler than dual methods, and has hence already been applied to realistic materials in a number of works.<cit.> 


In Fig. <ref>, we compare the phase diagram for model (<ref>) obtained from various simplified variants<cit.> of the dual boson scheme, and compare it to GW_c+EDMFT, GW+EDMFT, and the GW approximation. We restrict this comparison to values of U below the Mott transition, for lack of available dual-boson results in the Mott-insulating phase. 


Let us start with the small-U limit.
For U<1.5, the phase boundaries obtained in all the GW+EDMFT as well as GW alone are almost indistinguishable.
The GW transition is a straight line for all shown values, the variants with an impurity polarization have varying degrees of up-curvature, with the GW+EDMFT line in between the GW_c+EDMFT and the Gv+EDMFT line.
For Gv+EDMFT it is notable that the U→ 0 limit does not reproduce the GW value.  

The dual boson lines start with a similar upwards trend, with the exception of the DB-GW line, which follows essentially the HS-V variant of 
GW+EDMFT 
(more properly denoted as GD + second order perturbation theory (SOPT) + EDMFT)
as shown in Ref. Stepanov2016, 
and discussed in detail at the end of this section. 
Yet, the dual-boson variants start out with a lower slope, indicating stronger ordering tendencies already for the lowest values of U, while all GW+EDMFT variants follow the slope of the `weak-coupling' GW boundary in the vicinity of U=0.[Note that   the dual-boson variants and HS-V calculations of Ref. Stepanov2016 have been executed as a single lattice self-consistency iteration on top of the converged EDMFT solution.] 
Interestingly, a very recent cluster-EDMFT study <cit.> 
 reports a similarly reduced slope in the weak-U regime. 
 We also note that only the full DB critical line is above the U/4 line (dashed grey line, discussed in section <ref>), while the (non-self-consistent) DB-GW and DB-GWγ results are not (the latter only slightly so). Further comparisons of DB with GW+EDMFT (in the HS-V decoupling) can be found in Ref. Stepanov2016, see Fig. 8.

For U>2, the phase boundary for GW+EDMFT is below the dual-boson phase boundary, while GW is even a bit lower. 
For GW_c+EDMFT, the U=2.5 point already falls in the Mott-insulating phase and is not shown in the comparison. The too low Mott-transition line for GW_c+EDMFT is not surprising, since it lacks the band-widening of the nonlocal Fock term. 

Two possible decouplings were previously discussed in the literature, “HS-UV” 
(giving rise to GW+EDMFT) and “HS-V” (
resulting in a combined “GD+SOPT+EDMFT” scheme, 
where D is the screened non-local interaction, 
see Ref. Ayral2013).
We emphasize that, contrary to the HS-UV variant and e.g. the random-phase approximation (RPA), the HS-V variants do not resum, in the nonlocal part of the self-energy, the local (U) and nonlocal (V) parts of the interaction to the same order. (They are resummed, respectively, to second and infinite order.) This arbitrary inconsistency raises questions concerning the soundness of the HS-V scheme, as already pointed out in Ref. Ayral2013.

Recent works have indeed confirmed the deficiency of “HS-V”.
For instance, all the 
“EDMFT+GW”-related variants shown in Fig. 7 of Ref. Stepanov2016, which were obtained in a HS-V flavor, yield phase boundaries which are much lower than either DB or the GW+EDMFT results in Fig. <ref> (which correspond to the HS-UV variant of the decoupling of the interaction). The only additional outlier is the simplest DB type of approximation, the DB-GW variant, for which Ref. Stepanov2016 showed that it is formally similar to a HS-V calculation.  
More recently, Ref. Terletska2016 used cluster dynamical mean field theory to study the extended Hubbard model, which allows a control on errors by increasing the size of the cluster (but neglecting inter-cluster interactions, which in the case of EDMFT and GW+EDMFT are treated via the retarded impurity interactions). These cluster results were shown to be in poor agreement with “HS-V”, but very close to the full-fledged DB method. In the lower panel of Fig. <ref>, we show that the GW+EDMFT (HS-UV) method yields a critical V_c in agreement (with a 20% accuracy or better) with the cluster results, a remarkable result in view of the reduced numerical cost of this method compared to cluster DMFT. Comparisons for the larger U values in Fig. <ref> would be of great interest. 


We end this section by examining two further questions, namely the influence of spin fluctuations and the impact of local vertex corrections. One can expect that neither are important for the charge-ordering instability under study, since (i) this is an instability in the charge channel, not the spin channel, and (ii) as V increases towards charge ordering, the effective static interaction 𝒰(ω=0) decreases to zero,<cit.> making the system behave more and more like a weakly-correlated metal, where vertex corrections are expected to be small.

In all previous implementations of the GW+EDMFT method, the interaction was formally decoupled in the charge-channel only, neglecting the possible influence of spin fluctuations. Furthermore, in GW+EDMFT, the influence of the local vertex on the nonlocal self-energy is included only through the nonlocal Green's function. 
In the TRILEX approximation, both charge and spin fluctuations are taken into account, as well as local vertex corrections to the nonlocal self-energy.

We can thus answer both questions of interest by implementing the TRILEX method for the extended Hubbard model. We refer the reader to Refs. Ayral2015,Ayral2015c for implementation details. The only difference with respect to the application to the Hubbard model is that Eq. (41) of Ref. Ayral2015c must be modified to also describe nonlocal interactions, which means that Eq. (61b) of that publication becomes

    W^η(𝐪,iΩ) = v^η(𝐪)/1 - v^η(𝐪)P^η(𝐪, iΩ) ,

 where η denotes the charge (ch) or spin (sp) channel and

    v^ch(𝐪)   =   U^ch + 2 V(cos(q_x)+cos(q_y)),
    
      v^sp(𝐪)   =   U^sp,

and the bare on-site interactions in the charge and spin channels are parametrized, in the so-called Heisenberg decoupling,<cit.> by a parameter α:


    U^ch=(3α-1)U,  U^sp=(α-2/3)U.


In Fig. <ref>, we show TRILEX results for two characteristic points of the phase diagram, namely U=1 (characteristic of the metallic phase), and U=3 (characteristic of the Mott phase). 
First, we observe that the critical V_c (computed by looking for a vanishing inverse static susceptibility, shown in the left panels) is quite close to that of GW+EDMFT, justifying our a priori intuition. This agreement is quite remarkable, since GW+EDMFT has only charge fluctuations, while TRILEX has both charge and spin fluctuations.
Second, V_c only mildly depends on the ratio of the charge to spin fluctuations, as can be seen in the right panels, where quite large variations of U^ch (and correspondingly U^sp) lead to comparatively small variations in V_c. 

Taking inspiration from the comparison of the cluster extension of TRILEX with exact benchmark results for the two-dimensional Hubbard model<cit.> (there, one observes that whenever the TRILEX solution is close to the exact solution, the dependence on the decoupling is weak), this stability (compared to charge-only GW+EDMFT, and with respect to α) can be used as a proxy for the quantitative robustness of the present GW+EDMFT results.
 



§ CONCLUSION

In conclusion, we have shown that the nonlocal Fock term
has a significant influence on the description of the charge fluctuations
in the GW+EDMFT method, especially in the strong-coupling limit. By effectively enhancing the bandwidth, it
lowers the critical value of the nonlocal interaction for the charge-ordering transition.
We have also shown that the differences between the EDMFT and GW+EDMFT phase
diagrams are to a large extent a consequence of the nonlocal Fock
term, which is not included in EDMFT. 

Another interesting result is that the simple extension from EDMFT to a Gv+EDMFT
formalism yields results similar to the full-fledged GW+EDMFT method. This suggests the possibility of studying complex multiband materials, where a full GW+EDMFT computation would be too costly, using techniques in the spirit of the recent Screened exchange + dynamical DMFT (SEx+DMFT) 
method.<cit.> 
In realistic materials, the simple single-band description is not
sufficient, and substantial screening effects resulting from the presence
of higher energy degrees of freedom must be taken into 
account.<cit.> 

Performing a self-consistent calculation of the screening by these higher-energy states is however computationally expensive, even within a multi-tier approach,<cit.>  where the updates are restricted to an intermediate energy window. A scheme which combines a properly renormalized bandstructure with a self-consistent treatment of screening effects within the low-energy subspace may provide a good basis for 
tractable, but still accurate first principles
electronic structure methods for correlated electron materials.




We acknowledge useful discussions with Y. Nomura, A. I. Lichtenstein, E. Stepanov
and E. van Loon. We thank A. Huber and A. I. Lichtenstein for providing us the dual-boson data for Fig. <ref>, as well as H. Terletska for providing us the DCA data for Fig. <ref>. 
This research was supported by SNSF through NCCR Marvel, 
IDRIS/GENCI Orsay (project number t2016091393),
ECOS-Sud MinCYT under project number A13E04, and the 
European Research Council (project number 617196). Part of the implementation is based on the TRIQS toolbox<cit.> and on 
the ALPS libraries.<cit.> 



§ ESTIMATION OF THE BANDWIDTH WIDENING WITH COMBINED HUBBARD-I AND FOCK


In this appendix we discuss a simple Hubbard I (plus Fock) type treatment of the U-V model.
These arguments are not meant to be exact or comprehensive, most notably we ignore the effect of the nearest-neighbor interaction on the local self-energy and any nonlocal screening, but they provide useful insights into the nontrivial nature of the large-U and large-V limit. 

We start with an approximation to the self-energy which follows the
spirit of the Hubbard-I approximation by taking the atomic U^2/4z self-energy
locally, but goes beyond it by taking also the instantaneous nonlocal Fock contribution into account:

    Σ(k,z)=U^2/4z+2δ t(cos k_x+cos k_y)


The corresponding Green's function reads:

    G(k,z)=1/z-ε̃_k-U^2/4z


with ε̃_k denoting the effective
dispersion including the Fock term:

    ε̃_k≡2(t+δ t)(cos k_x+cos k_y)


Thus, we can write:


    G(k,z)=z/(z-z_+(k))(z-z_-(k))


with

    z_±(k)=ε̃_k±√(ε̃_k^2+U^2)/2.

As expected, in the atomic limit (ε̃_k→0),
the function has two peaks at ± U/2, corresponding to the two
Hubbard bands.

We can decompose the expression of Eq. (<ref>) as:


    G(k,z)=A_+(k)/z-z_+(k)+A_-(k)/z-z_-(k)


with

    A_±(k)≡1/2(1±ε̃_k/√(ε̃_k^2+U^2)).

A_+ and A_- are the weights of the upper and lower Hubbard
bands, respectively. Using Eq. (<ref>), one writes the
spectral function as:

    A(k,ω)=A_+(k)πδ(ω-z_+(k))+A_-(k)πδ(ω-z_-(k))


Under the assumption that the Hubbard bands are well separated (U
large enough), only the lower Hubbard band contributes to the occupancy
(at T=0 for simplicity):

    n_k=∫_-∞^0dω/πA(k,ω)≈ A_-(k)≈1/2(1-ε̃_k/U)


In the second equality, we have again used the fact that U is large
enough (to neglect ε̃_k^2 in
the square root).

On the other hand, the occupancy is related to G(k,τ=0^+)
in the following way:

    n_k=1+G_k(τ=0^+)


Lastly, δ t (in Eq. (<ref>)) is also known,
in the Fock approximation, as a function of G_ij(τ=0^+) (see Eq. (<ref>)):

    δ t=-VG_⟨ ij⟩(τ=0^+)


Putting (<ref>-<ref>-<ref>)
together and Fourier transforming, one gets


    G_⟨ ij⟩(τ=0^+)=1/2(-t-VG_⟨ ij⟩(τ=0^+)/U).

Solving for G_⟨ ij⟩(τ=0^+), one gets:


    G_⟨ ij⟩(τ=0^+)=-t/2U-V


and 

    δ t=tV/2U-V.



This expression is used in Section <ref>.


apsrev4-1

