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

TL;DR

This paper investigates how including the nonlocal Fock exchange term in combined many-body perturbation and dynamical mean field theory affects the phase diagram and interactions in electronic systems, revealing significant deviations from standard EDMFT especially in strong coupling.

Contribution

It introduces the inclusion of the nonlocal Fock exchange term into EDMFT and $GW$+EDMFT, highlighting its impact on phase diagrams and interaction effects in correlated electronic systems.

Findings

Nonlocal Fock exchange causes significant deviations in phase diagrams.

Differences between EDMFT and $GW$+EDMFT are largest in strong-coupling regimes.

$GW$+EDMFT results align closely with dual boson methods for moderate interactions.

Abstract

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.

Equations

$$\mathcal{H}=-\sum _{ij}{t}_{ij}{c}_i^{†}{c}_j+\frac{1}{2}\sum _{ij}{v}_{ij}{n}_i{n}_j-\mu \sum _i{n}_i.$$

$$t_{ij}=$$

$$v_{ij}=$$

$${\epsilon }_{\mathbf{k}}=$$

$$+2t^{\prime }(\cos (k_x+k_y)+\cos (k_x-k_y))$$

$$v_{\mathbf{q}}=U$$

$$+2V^{\prime }(\cos (q_x+q_y)+\cos (q_x-q_y)).$$

$$\Sigma (\mathbf{k},i{\omega }_n)={\Sigma }_{\mathrm{imp}}(i{\omega }_n)+{\Sigma }_{G{W}_c}^{\mathrm{nonloc}}(\mathbf{k},i{\omega }_n)+{\Sigma }_{\mathrm{F}}^{\mathrm{nonloc}}(\mathbf{k}).$$

$${\Sigma }_{G{W}_c}^{\mathrm{nonloc}}(\mathbf{k},\tau )$$

$${\Sigma }_{\mathrm{F}}^{\mathrm{nonloc}}(\mathbf{k})$$

$$G_{\mathbf{k}}(i{\omega }_n)=(i{\omega }_n+\mu -{\epsilon }_{\mathbf{k}}-\Sigma (\mathbf{k},i{\omega }_n){)}^{-1}$$

$$W_{\mathbf{q}}^{\mathrm{c}}(i{\nu }_m)\equiv \frac{v_{\mathbf{q}}}{1-v_{\mathbf{q}}{P}_{\mathbf{q}}(i{\nu }_m)}-v_{\mathbf{q}},$$

$$\tilde{t}(t^{\prime },U,V,V^{\prime })=t+\delta t(t^{\prime },U,V,V^{\prime }).$$

$${\Sigma }_{\mathrm{F}}^{\mathrm{nonloc}}{}_{ij}=$$

$$=$$

$$-G_{⟨⟨ij⟩⟩}(\tau =0^{+})V^{\prime }{\delta }_{⟨⟨ij⟩⟩},$$

$$\delta t(t^{\prime },U,V,V^{\prime })=-G_{⟨ij⟩}(\tau =0^{+})V.$$

$$\delta t^{\prime }(t^{\prime },U,V,V^{\prime })=-G_{⟨⟨ij⟩⟩}(\tau =0^{+})V^{\prime }.$$

$$H_0>H_{\mathrm{xc}}^{\mathrm{F}}>H_{\mathrm{KS}}.$$

$${\chi }_{\mathbf{q}}(i{\nu }_m)=\frac{-{\Pi }_{\mathbf{q}}(i{\nu }_m)}{1-{\Pi }_{\mathbf{q}}(i{\nu }_m)v_{\mathbf{q}}}.$$

$$V_c^{\mathrm{int}}=U/4,$$

$$\delta t=\frac{tV}{2U-V}.$$

$$\delta t_c=\frac{t(\alpha U+\beta )}{2U-\alpha U-\beta }.$$

$$W^{\eta }(\mathbf{q},i\Omega )=\frac{v^{\eta }(\mathbf{q})}{1-v^{\eta }(\mathbf{q})P^{\eta }(\mathbf{q},i\Omega )},$$

$$v^{\mathrm{ch}}(\mathbf{q})$$

$$v^{\mathrm{sp}}(\mathbf{q})$$

$$U^{\mathrm{ch}}=(3\alpha -1)U,U^{\mathrm{sp}}=(\alpha -2/3)U.$$

$$\Sigma (\mathbf{k},z)=\frac{U^2}{4z}+2\delta t(\cos k_x+\cos k_y)$$

$$G(\mathbf{k},z)=\frac{1}{z-{\tilde{\epsilon }}_{\mathbf{k}}-\frac{U^2}{4z}}$$

$${\tilde{\epsilon }}_{\mathbf{k}}\equiv 2(t+\delta t)(\cos k_x+\cos k_y)$$

$$G(\mathbf{k},z)=\frac{z}{(z-z_{+}(\mathbf{k}))(z-z_{-}(\mathbf{k}))}$$

$$z_{\pm }(\mathbf{k})=\frac{{\tilde{\epsilon }}_{\mathbf{k}}\pm \sqrt{{\tilde{\epsilon }}_{\mathbf{k}}^2+U^2}}{2}.$$

$$G(\mathbf{k},z)=\frac{A_{+}(\mathbf{k})}{z-z_{+}(\mathbf{k})}+\frac{A_{-}(\mathbf{k})}{z-z_{-}(\mathbf{k})}$$

$$A_{\pm }(\mathbf{k})\equiv \frac{1}{2}\left(1\pm \frac{{\tilde{\epsilon }}_{\mathbf{k}}}{\sqrt{{\tilde{\epsilon }}_{\mathbf{k}}^2+U^2}}\right).$$

Full text

Dynamical Mean Field TheoryGeorges et al. (1996) (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 theorySachdev and Ye (1993); Sengupta and Georges (1995); Kajueter (1996); Si and Smith (1996); Parcollet and Georges (1999); Chitra and Kotliar (2000, 2001); Smith and Si (2000); Si et al. (2001) (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 fieldvan Loon et al. (2015), nor does it obey local charge conservation rulesHafermann et al. (2014) 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 Biermann et al. (2003); Sun and Kotliar (2004); Ayral et al. (2013) 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]^{\mathrm{nonloc}}$” is included in the nonlocal “$[GW]^{\mathrm{nonloc}}$” diagram. As described in more detail in Ref. Ayral et al., 2013 (see also the appendix of Ref. Golez et al., 2017), the $GW$+EDMFT method is formally obtained by constructing an energy functional of $G$ and $W$, the AlmbladhAlmbladh et al. (1999) functional $\Psi$, and by approximating $\Psi$ 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$ approximationHedin (1965)). This functional construction rules out double-counting of local terms in the self-energy and polarizationBiermann et al. (2004); Ayral et al. (2013). Even though it has been introduced under the name $GW$+DMFT Biermann et al. (2003) 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. Ayral et al., 2013, and related papersAyral et al. (2012); Huang et al. (2014), the nonlocal Fock term was omitted.Ayral et al. (2016) 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$ Miyake et al. (2013) and $GW$+DMFT,Tomczak et al. (2014) albeit not in a self-consistent way.

The manuscript is organized as follows: In section I, we recap the $GW$+EDMFT equations with special emphasis on the Fock term and make general statements about the expected impact. In section II, 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 III. Section IV 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 V, we compare our results with results obtained within the recent dual boson method.

I 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 $\Sigma$ and polarization $P$, and a diagrammatic expansion in their nonlocal components. The model is defined by the Hamiltonian

[TABLE]

Here, $t_{ij}$ are the real-space hopping matrix elements, $c_{i}^{(\dagger)}$ the electronic annihilator (creator) on site $i$, $v_{ij}$ the Coulomb interaction and $n_{i}=c^{\dagger}_{i}c_{i}$. We will restrict ourselves to models with hoppings and interactions between nearest neighbors and next nearest neighbors only,

[TABLE]

where $\delta_{ij}$ is the usual Kronecker delta, $\delta_{\left<ij\right>}$(resp. $\delta_{\left<\left<ij\right>\right>}$) 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

[TABLE]

and

[TABLE]

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

[TABLE]

The last two terms correspond to the nonlocal part of the $GW$ self-energy. They can be expressed as a function of imaginary time $\tau$ and momentum $\mathbf{k}$ as follows:

[TABLE]

Fourier transformations between $\tau$ and fermionic [bosonic] Matsubara frequencies $i\omega_{n}=i(2n+1)\frac{\pi}{\beta}$ [$i\nu_{m}=i2m\frac{\pi}{\beta}$] are assumed where needed ($\beta$ denotes the inverse temperature). The “loc” suffix denotes a sum over the first Brillouin zone.

[TABLE]

is the interacting lattice Green’s function, and $W^{\mathrm{c}}$ is defined as

[TABLE]

with $P$ the polarization function. All results are given in units of $D=4|t|$ (which is the half bandwidth when $t^{\prime}=0$), and the momentum discretization is $N_{k}=32\times 32$ points in the first Brillouin zone, unless otherwise stated. We use the original formulation of the $GW$+EDMFT scheme,Biermann et al. (2003) corresponding – within a functional formulation – to a Hubbard-Stratonovich decoupling of the full interaction term, dubbed “HS-$UV$ decoupling” in Ref. Ayral et al., 2013. 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. (6), which is real-valued and instantaneous, renormalizes the bandwidth. It can become quite large and momentum-dependent. This is illustrated in Fig. 1 for the parameters $U=2$, $V=0.4$ (and $t^{\prime}=0$, $V^{\prime}=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. (4) by defining a $U$- and $V$-dependent hopping

[TABLE]

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

[TABLE]

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

Thus, for the case of Fig. 1, where $V^{\prime}=0$, the Fock term enters Eq. (9) as a renormalization of the nearest-neighbor hopping by

[TABLE]

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^{\prime}$ also gets renormalized according to

[TABLE]

Note that such a term breaks the particle-hole symmetry of the lattice. In the particle-hole symmetric case with $t^{\prime}=0$ and half filling, the next-nearest-neighbor occupation term vanishes: $G_{\left<\left<ij\right>\right>}({\tau\!=\!0^{+}})=0$.

II Effective bandstructure

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

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

To study a model with broken particle-hole symmetry, we introduce a nearest-neighbor hopping $t^{\prime}=t/\sqrt{2}$ and fix the filling at $\left<n\right>=0.8$ as well as the nearest-neighbor interaction $V$. We then calculate $\delta t$ as a function of the next nearest nearest-neighbor interaction $V^{\prime}$. As shown in Fig. 3, the main effect on the hopping renormalization comes from the (essentially linear) dependence on $V$. The qualitative effect of $V^{\prime}$ 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_{\mathrm{KS}}$ of Density Functional Theory (DFT),

(ii) the Hartree Hamiltonian $H_{0}=H_{\mathrm{KS}}-V[v_{\mathrm{xc}}(\mathbf{r})]$, where $v_{\mathrm{xc}}(\mathbf{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(\mathbf{r})$, while “nonlocal” denotes a function depending on two variables $f(\mathbf{r},\mathbf{r^{\prime}})$),

(iii) the nonlocal exchange Hamiltonian $H_{\mathrm{xc}}^{\mathrm{F}}=H_{0}+V[v_{\mathrm{xc}}^{\mathrm{F}}(\mathbf{r},\mathbf{r^{\prime}})]$ (where $v_{\mathrm{xc}}^{\mathrm{F}}(\mathbf{r},\mathbf{r^{\prime}})$ denotes an exchange-correlation potential including “nonlocal” Fock exchange).

Then the hierarchy of the bandwidths in a metallic system is

[TABLE]

Thus, $H_{\mathrm{xc}}^{\mathrm{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. Hirayama et al., 2015 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. Hirayama et al., 2015, 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. Hirayama et al., 2015 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. van Roekeghem and Biermann, 2014) 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 van Roekeghem et al. (2014, 2016a)). We refer the interested reader to Refs. Hirayama et al., 2015; van Roekeghem et al., 2014; van Roekeghem and Biermann, 2014 for details.

III 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. Ayral et al., 2013 (which contains only the $GW_{c}$ term, see Ref. Ayral et al., 2016), and “$GW$+EDMFT” the formula with the self-energy expression (6). For comparison, we also show results for “$Gv$+EDMFT”, a scheme where $\Sigma$ is the sum of the impurity self-energy and of the nonlocal Fock term only (the first and third terms of Eq. (6)). 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. Ayral et al., 2013, steps (5) (a) and (5) (b) of section V.

In Fig. 4, 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. Ayral et al., 2013, which compared $GW_{c}$+EDMFT to EDMFT.

In the top panels of Fig. 5, we replot panels (a) and (b) of Fig. 15 of Ref. Ayral et al., 2013, 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. 6, we plot the corresponding spectral functions (the EDMFT and $GW_{c}$+EDMFT results are identical to Fig. 2(b) of Ref. Ayral et al., 2012). 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. 7, we plot the dependence of the inverse charge susceptibility $\chi^{-1}_{\mathbf{q}=\pi,\pi}(i\nu_{m}=0)$ on the nearest-neighbor repulsion $V$, with $\chi$ defined by

[TABLE]

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. 8, where we plot the results from Fig. 5 of Ref. Ayral et al., 2013 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).

IV Inside the Mott phase

In the large-interaction regime of the phase diagram (Fig. 8), 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}\rightarrow 0$) and zero-temperature limit of the extended Hubbard model by exact Monte-Carlo simulations, as e.g. in Ref. Pawlowski, 2006, and is given by the analytical expression

[TABLE]

where 4 corresponds to the number of nearest neighbors. This line is plotted as a dashed grey line in Fig. 8. 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. 8, 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.Pawlowski (2006) Since the quantum tunneling (hopping) has a physical effect similar to temperature in classical systems,Biermann et al. (1998) 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. (16)) imposes that the charge-ordering transition should occur for negative values of $v_{\mathbf{q}}$, since the polarization $\Pi_{\mathbf{q}}$ is always negative (for the parameters studied here). For the square lattice, this implies $V_{c}>\frac{U}{4}$, which is the classical energy estimate (Eq. (17)). 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 $\delta t$ (defined in Eq. (13)) coming from the Fock term. Approximating the self-energy as the sum of the atomic limit (in the spirit of the Hubbard-I approximationHubbard (1963)) and of the Fock self-energy, as described in more detail in Appendix A, we obtain

[TABLE]

Thus $\delta 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. 8, we can parametrize the phase boundary in the large-$U$ limit as a constant slope, i.e $V_{c}=\alpha U+\beta$. [Within the $U$-range that we can simulate (we performed measurements up to $U=8.0$), we can estimate $\alpha=1.25$ and $\beta=-2.5$.] With this parametrization, we obtain

[TABLE]

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

V 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,Rubtsov et al. (2012) which, in its full-fledged implementation,van Loon et al. (2014); Stepanov et al. (2016a) computes the susceptibilityHafermann et al. (2014); van Loon et al. (2015) 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 ApproximationToschi et al. (2007); Katanin et al. (2009); Schäfer et al. (2013); Valli et al. (2015); Schäfer et al. (2015, 2016) (which was recently shown to be a simplified version of QUADRILEX, a method consisting in an atomic approximation of the four-particle irreducible functionalAyral and Parcollet (2016a)), can nonetheless only be obtained at a considerable computational expense and require a proper parametrization and treatment of their asymptotic behavior.Kuneš (2011); Boehnke et al. (2011); Li et al. (2016); Wentzell et al. (2016) Consequently, it is not possible to use them routinely in multi-orbital calculations (see however Boehnke et al. (2014); Galler et al. (2017)) 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 methodAyral and Parcollet (2015, 2016b) and simplified dual boson schemes such as DB+$GW$ or DB+$GW\gamma$.Stepanov et al. (2016b) 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 $\gamma$ 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.Tomczak et al. (2012); Hansmann et al. (2013, 2016); Boehnke et al. (2016)

In Fig. 9, we compare the phase diagram for model (1) obtained from various simplified variantsStepanov et al. (2016b) 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\rightarrow 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. Stepanov et al., 2016b, 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$.111Note that the dual-boson variants and HS-$V$ calculations of Ref. Stepanov et al., 2016b have been executed as a single lattice self-consistency iteration on top of the converged EDMFT solution. Interestingly, a very recent cluster-EDMFT study Terletska et al. (2017) 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 IV), while the (non-self-consistent) DB-$GW$ and DB-$GW\gamma$ 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. Stepanov et al., 2016b, 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. Ayral et al., 2013). 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. Ayral et al., 2013.

Recent works have indeed confirmed the deficiency of “HS-$V$”. For instance, all the “EDMFT+$GW$”-related variants shown in Fig. 7 of Ref. Stepanov et al., 2016b, 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. 9 (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. Stepanov et al., 2016b showed that it is formally similar to a HS-$V$ calculation. More recently, Ref. Terletska et al., 2017 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. 9, 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. 9 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 $\mathcal{U}(\omega=0)$ decreases to zero,Ayral et al. (2013) 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. Ayral and Parcollet, 2015, 2016b for implementation details. The only difference with respect to the application to the Hubbard model is that Eq. (41) of Ref. Ayral and Parcollet, 2016b must be modified to also describe nonlocal interactions, which means that Eq. (61b) of that publication becomes

[TABLE]

where $\eta$ denotes the charge ($\mathrm{ch}$) or spin ($\mathrm{sp}$) channel and

[TABLE]

and the bare on-site interactions in the charge and spin channels are parametrized, in the so-called Heisenberg decoupling,Ayral and Parcollet (2016b) by a parameter $\alpha$:

[TABLE]

In Fig. 10, 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^{\mathrm{ch}}$ (and correspondingly $U^{\mathrm{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 modelAyral et al. (2017) (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 $\alpha$) can be used as a proxy for the quantitative robustness of the present $GW$+EDMFT results.

VI 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.van Roekeghem et al. (2014); van Roekeghem and Biermann (2014); van Roekeghem et al. (2016b) 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.Aryasetiawan et al. (2004); Biermann (2014); Werner and Casula (2016)

Performing a self-consistent calculation of the screening by these higher-energy states is however computationally expensive, even within a multi-tier approach,Boehnke et al. (2016) 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.

Acknowledgements.

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. 9, as well as H. Terletska for providing us the DCA data for Fig. 9. 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 toolboxParcollet et al. (2015) and on the ALPS libraries.Bauer et al. (2011)

Appendix A 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 $\frac{U^{2}}{4z}$ self-energy locally, but goes beyond it by taking also the instantaneous nonlocal Fock contribution into account:

[TABLE]

The corresponding Green’s function reads:

[TABLE]

with $\tilde{\varepsilon}_{\boldsymbol{k}}$ denoting the effective dispersion including the Fock term:

[TABLE]

Thus, we can write:

[TABLE]

with

[TABLE]

As expected, in the atomic limit ($\tilde{\varepsilon}_{\boldsymbol{k}}\rightarrow 0$), the function has two peaks at $\pm U/2$, corresponding to the two Hubbard bands.

We can decompose the expression of Eq. (25) as:

[TABLE]

with

[TABLE]

$A_{+}$ and $A_{-}$ are the weights of the upper and lower Hubbard bands, respectively. Using Eq. (27), one writes the spectral function as:

[TABLE]

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):

[TABLE]

In the second equality, we have again used the fact that $U$ is large enough (to neglect $\tilde{\varepsilon}_{\boldsymbol{k}}^{2}$ in the square root).

On the other hand, the occupancy is related to $G(\boldsymbol{k},\tau=0^{+})$ in the following way:

[TABLE]

Lastly, $\delta t$ (in Eq. (24)) is also known, in the Fock approximation, as a function of $G_{ij}(\tau=0^{+})$ (see Eq. (13)):

[TABLE]

Putting (29-30-31) together and Fourier transforming, one gets

[TABLE]

Solving for $G_{\langle ij\rangle}(\tau=0^{+})$, one gets:

[TABLE]

and

[TABLE]

This expression is used in Section IV.