Multitier self-consistent $GW$+EDMFT

Fredrik Nilsson; Lewin Boehnke; Philipp Werner; Ferdi Aryasetiawan

arXiv:1706.06808·cond-mat.str-el·October 2, 2017

Multitier self-consistent $GW$+EDMFT

Fredrik Nilsson, Lewin Boehnke, Philipp Werner, Ferdi Aryasetiawan

PDF

Open Access

TL;DR

This paper introduces a parameter-free, self-consistent $GW$+EDMFT method for simulating correlated materials, effectively capturing local and nonlocal interactions and screening effects with improved accuracy over previous approaches.

Contribution

The paper develops and details a multitier $GW$+EDMFT framework that combines ab initio $G^0W^0$ calculations with a self-consistent $GW$ plus EDMFT treatment for better modeling of correlated materials.

Findings

01

The method accurately captures correlation effects in intermediate regimes.

02

It predicts a Mott transition at a specific lattice spacing.

03

Impurity interactions show expected trends with lattice spacing.

Abstract

We discuss a parameter-free and computationally efficient ab initio simulation approach for moderately and strongly correlated materials, the multitier self-consistent $G W$ +EDMFT method. This scheme treats different degrees of freedom, such as high-energy and low-energy bands, or local and nonlocal interactions, within appropriate levels of approximation, and provides a fully self-consistent description of correlation and screening effects in the solid. The ab initio input is provided by a one-shot $G^{0} W^{0}$ calculation, while the strong-correlation effects originating from narrow bands near the Fermi level are captured by a combined $G W$ plus extended dynamical mean-field (EDMFT) treatment. We present the formalism and technical details of our implementation and discuss some general properties of the effective EDMFT impurity action. In particular, we show that the retarded impurity…

Tables1

Table 1. Table 1: Relation of existing formalisms to the multitier scheme. Embedding the TIER I values into the TIER II approximation is done as described in Sec. II.5 . Formalisms without an approximation to the polarization Π Π \Pi in TIER II do not update the impurity interaction. The TIER I column lists the quantities that have to be measured in the impurity problem and enter the self-consistency on TIER II. It is implicit that the corresponding doublecountings between the tiers are removed, which makes it necessary to add them back on TIER II in the case of E xc LDA superscript subscript 𝐸 xc LDA E_{\mathrm{xc}}^{\mathrm{LDA}} and E SEx LDA superscript subscript 𝐸 SEx LDA E_{\mathrm{SEx}}^{\mathrm{LDA}} , where the full LDA/SEx band values are used in the self-consistency. ∅ \emptyset indicates that the corresponding tier is empty.

	non-sc		selfconsistent
Colloquial name	TIER III		TIER II		TIER I	References / Comments
	$Σ$	$Π$	$Σ$	$Π$	used quantities
model DMFT	—	—	local	—	$Σ$	Ref. Georges et al.,1996
model EDMFT	—	—	local	local	$Σ$ , $Π$	Ref. Chitra and Kotliar,2001
model D $Γ$ A	—	—	SDE²²2Schwinger-Dyson equation/BSE³³3Bethe-Salpeter equation	—	$Σ$ , $γ$ ⁴⁴4 $γ$ denotes the fully irreducible vertex of the impurity problem	Ref. Toschi et al., 2007 (non-self-consistent)
model $G W$ +EDMFT	—	—	$G W$ +imp	$G G$ +imp	$Σ$ , $Π$	Refs. Sun and Kotliar,2002; Ayral et al.,2013; Huang et al.,2014; Ayral et al.,2017
model TRILEX	—	—	$G W Λ$ +imp	$G G Λ$ +imp	$Σ$ , $Π$ , $Λ$ ⁵⁵5 $Λ$ denotes the three-leg vertex	Refs. Ayral and Parcollet,2015, 2016a
$G^{0} W^{0}$	$G^{0} W^{0}$	$G^{0} G^{0}$	$\emptyset$	$\emptyset$	$\emptyset$	often referred to as $G W$ . Ref. Aryasetiawan and Gunnarsson,1998
$G W$	—	—	$G W$	$G G$	$\emptyset$	rarely used. Refs. Holm and von Barth,1998; García-González and Godby,2001; Koval et al.,2014
$Q P$ sc $G W$	—	—	$G^{Q P} W^{Q P}$	$G^{Q P} G^{Q P}$	$\emptyset$	Refs. van Schilfgaarde et al.,2006; Bruneval et al.,2006; Gatti and Guzzo,2013; Koval et al.,2014
LDA+DMFT	$E_{xc}^{LDA}$	—	$E_{xc}^{LDA}$ +imp	—	$Σ$	Ref. Kotliar et al.,2006
LDA+DMFT+ $U (ω)$	$E_{xc}^{LDA}$	$G^{0} G^{0}$	$E_{xc}^{LDA}$ +imp	—	$Σ$	Ref. Werner et al.,2012; Sakuma et al.,2013; Werner et al.,2015
SEx+DMFT	$E_{xc}^{SEx}$	$G^{0} G^{0}$	$E_{xc}^{SEx}$ +imp	—	$Σ$	Ref. van Roekeghem et al., 2014 (BFA⁶⁶6BFA stands for Bose Factorization ApproachCasula et al. (2012) for impurity)
$G W$ +DMFT+ $U (ω)$	$G^{0} W^{0}$	$G^{0} G^{0}$	$G^{0} W^{0}$ +imp	—	$Σ$	Ref. Tomczak et al., 2014 (BFA for impurity)
ab initio D $Γ$ A	$E_{xc}^{LDA}$	—	$E_{xc}^{LDA}$ +SDE/BSE	—	$Σ$ , $γ$	Refs. Toschi et al., 2011; Galler et al., 2017 (non-self-consistent)
$G W$ +EDMFT	$G^{0} W^{0}$	$G^{0} G^{0}$	$G W$ +imp	$G G$ +imp	$Σ$ , $Π$	Refs. Biermann et al.,2003; Boehnke et al.,2016

Equations170

Σ (1, 2) = - \int d (34) G (1, 3^{+}) W (1, 4) Γ (3, 2, 4),

Σ (1, 2) = - \int d (34) G (1, 3^{+}) W (1, 4) Γ (3, 2, 4),

G (1, 2) = G_{0} (1, 2) + \int d (34) G_{0} (1, 3) Σ (3, 4) G (4, 2),

Γ (1, 2, 3) = δ (1 - 2) δ (2 - 3),

+ \int d (4567) \frac{\partial Σ ( 1 , 2 )}{\partial G ( 4 , 5 )} G (4, 6) G (7, 5) Γ (6, 7, 3),

Π (1, 2) = \int d (34) G (1, 3) Γ (3, 4, 2) G (4, 1^{+}),

W (1, 2) = v (1, 2) + \int d (34) v (1, 3) Π (3, 4) W (4, 2) .

Σ^{G W} = - G (1, 2) W (1, 2),

Σ^{G W} = - G (1, 2) W (1, 2),

Π (1, 2) = G (1, 2) G (2, 1^{+}) .

Π (1, 2) = G (1, 2) G (2, 1^{+}) .

Ψ_{α R} (r) = w_{i R}^{*} (r) w_{j R} (r),

Ψ_{α R} (r) = w_{i R}^{*} (r) w_{j R} (r),

G_{imp} = G_{loc},

G_{imp} = G_{loc},

W_{imp} = W_{loc} .

S =

S =

+ \frac{1}{2} \int_{0}^{β} d τ d τ^{'} σ σ^{'} \sum ab c d \sum U_{ab d c} (τ - τ^{'})

\frac{1}{2} \int_{0}^{β} d τ d τ^{'} \times c_{aσ}^{†} (τ) c_{bσ} (τ) c_{c σ^{'}}^{†} (τ^{'}) c_{d σ^{'}} (τ^{'})

G = (Σ^{imp} + G^{imp}^{- 1})^{- 1}

G = (Σ^{imp} + G^{imp}^{- 1})^{- 1}

U = W^{imp} (\mathbbm 1 + Π^{imp} W^{imp})^{- 1} .

U = W^{imp} (\mathbbm 1 + Π^{imp} W^{imp})^{- 1} .

H_{K} =

H_{K} =

- ab \sum J (spin - flip c_{a ↑}^{†} c_{a ↓} c_{b ↓}^{†} c_{b ↑} + pair - hopping c_{b ↑}^{†} c_{b ↓}^{†} c_{a ↑} c_{a ↓}) .

H = \frac{1}{2} σ σ^{'} \sum ab c d \sum U_{ab d c} (c_{aσ}^{†} c_{bσ} c_{c σ^{'}}^{†} c_{d σ^{'}})

H = \frac{1}{2} σ σ^{'} \sum ab c d \sum U_{ab d c} (c_{aσ}^{†} c_{bσ} c_{c σ^{'}}^{†} c_{d σ^{'}})

Σ = Σ^{G W} + Σ^{vc},

Σ = Σ^{G W} + Σ^{vc},

Π = Π^{GG} + Π^{vc} .

Σ^{G W + EDMFT} (k) = Σ^{G W} (k) + Σ^{EDMFT}

Σ^{G W + EDMFT} (k) = Σ^{G W} (k) + Σ^{EDMFT}

- k \sum Σ^{G W} (k),

Π^{G W + EDMFT} (k) = Π^{GG} (k) + Π^{EDMFT}

- k \sum Π^{GG} (k) .

Z = \int D [ψ ψ^{†}] exp (- S [ψ, ψ^{†}]),

Z = \int D [ψ ψ^{†}] exp (- S [ψ, ψ^{†}]),

S [ψ, ψ^{†}] = \int d x ψ^{†} (x) (\partial_{τ} - \frac{\nabla ^{2}}{2 m} + V_{c} (x)) ψ (x)

S [ψ, ψ^{†}] = \int d x ψ^{†} (x) (\partial_{τ} - \frac{\nabla ^{2}}{2 m} + V_{c} (x)) ψ (x)

- \frac{1}{2} \int d x d x^{'} ψ^{†} (x) ψ^{†} (x^{'}) V (x - x^{'}) ψ (x^{'}) ψ (x) .

S [ϕ, ψ, ψ^{†}] = \int d x ψ^{†} (x) (\partial_{τ} - \frac{\nabla ^{2}}{2 m} + V_{H} (x)) ψ (x)

S [ϕ, ψ, ψ^{†}] = \int d x ψ^{†} (x) (\partial_{τ} - \frac{\nabla ^{2}}{2 m} + V_{H} (x)) ψ (x)

- \frac{1}{2} \int d x d x^{'} ϕ (x) V^{- 1} (x - x^{'}) ϕ (x^{'})

- i α \int d x ϕ (x) ψ^{†} (x) ψ (x) .

S [ϕ, ψ, ψ^{†}, J_{f}, J_{b}] = S [ϕ, ψ, ψ^{†}]

S [ϕ, ψ, ψ^{†}, J_{f}, J_{b}] = S [ϕ, ψ, ψ^{†}]

- \int d x d x^{'} J_{f} (x, x^{'}) ψ^{†} (x) ψ (x^{'})

- \frac{1}{2} \int d x d x^{'} J_{b} (x, x^{'}) ϕ (x) ϕ (x^{'}) .

Ω (J_{f}, J_{b}) = - ln (Z (J_{f}, J_{b})) .

Ω (J_{f}, J_{b}) = - ln (Z (J_{f}, J_{b})) .

Γ [G, W] = Ω (J_{f}, J_{b}) - J_{f} G - \frac{J _{b}}{2} W =

Γ [G, W] = Ω (J_{f}, J_{b}) - J_{f} G - \frac{J _{b}}{2} W =

Tr [ln (G)] - Tr [(G_{H}^{- 1} - G^{- 1}) G] - \frac{Tr [ ln ( W )]}{2}

+ \frac{Tr [( V ^{- 1} - W ^{- 1} ) W ]}{2} + Ψ [G, W],

Ψ (G, W) = \int d α \int d x ⟨ ϕ (x) ψ^{†} (x) ψ (x)⟩ .

Ψ (G, W) = \int d α \int d x ⟨ ϕ (x) ψ^{†} (x) ψ (x)⟩ .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Taxonomy

TopicsPhysics of Superconductivity and Magnetism · Advanced Physical and Chemical Molecular Interactions · Electronic and Structural Properties of Oxides

Full text

Multitier self-consistent GW+EDMFT

F. Nilsson

[email protected]

Department of Physics, Division of Mathematical Physics, Lund University, Professorsgatan 1, 223 63 Lund, Sweden

L. Boehnke

[email protected]

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

P. Werner

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

F. Aryasetiawan

Department of Physics, Division of Mathematical Physics, Lund University, Professorsgatan 1, 223 63 Lund, Sweden

Abstract

We discuss a parameter-free and computationally efficient ab initio simulation approach for moderately and strongly correlated materials, the multitier self-consistent $GW$ +EDMFT method. This scheme treats different degrees of freedom, such as high-energy and low-energy bands, or local and nonlocal interactions, within appropriate levels of approximation, and provides a fully self-consistent description of correlation and screening effects in the solid. The ab initio input is provided by a one-shot $G^{0}W^{0}$ calculation, while the strong-correlation effects originating from narrow bands near the Fermi level are captured by a combined $GW$ plus extended dynamical mean-field (EDMFT) treatment. We present the formalism and technical details of our implementation and discuss some general properties of the effective EDMFT impurity action. In particular, we show that the retarded impurity interactions can have non-causal features, while the physical observables, such as the screened interactions of the lattice system, remain causal. As a first application, we present ab initio simulation results for SrMoO3, which demonstrate the existence of prominent plasmon satellites in the spectral function not obtainable within LDA+DMFT, and provide further support for our recent re-interpretation of the satellite features in the related cubic perovskite SrVO3. We then turn to stretched sodium as a model system to explore the performance of the multitier self-consistent $GW$ +EDMFT method in situations with different degrees of correlation. While the results for the physical lattice spacing $a_{0}$ show that the scheme is not very accurate for electron-gas like systems, because nonlocal corrections beyond $GW$ are important, it does provide physically correct results in the intermediate correlation regime, and a Mott transition around a lattice spacing of $1.5a_{0}$ . Remarkably, even though the Wannier functions in the stretched compound are less localized, and hence the bare interaction parameters are reduced, the self-consistently computed impurity interactions show the physically expected trend of an increasing interaction strength with increasing lattice spacing.

pacs:

71.10.Fd

I introduction

Accurate descriptions of materials with strong electron-electron interactions remain one of the main challenges in modern condensed matter theory. The $GW$ approximation, proposed in 1965 by Hedin Hedin (1965) was one of the first successful attempts to go beyond density functional theory (DFT)Kohn and Sham (1965) for real materials. One of the major successes of the $GW$ approximation in its one-shot ( $G^{0}W^{0}$ ) version is that it cures the band-gap underestimation of the local density approximation (LDA) for a wide range of semi-conductors.Aryasetiawan and Gunnarsson (1998); van Schilfgaarde et al. (2006) However, even for weakly correlated metals such as sodium, the $G^{0}W^{0}$ description yields a too small band narrowing,Lyo and Plummer (1988) and the theory does not capture the Mott physics that is crucial for understanding the properties of many strongly-correlated 3 $d$ and 4 $f$ materials. Fully self-consistent $GW$ calculations are rarely performed, because they are computationally expensive and generally worsen the one-shot results, which are in better agreement with experiment. This implies that the quasiparticle Green’s function, $G^{0}$ , includes (in a somewhat uncontrolled way) vertex corrections needed in the fully self-consistent $GW$ calculations. Inspired by this observation, schemes based on quasiparticle self-consistency have been developed and shown to improve the $GW$ method for a number of different materials.van Schilfgaarde et al. (2006); Bruneval et al. (2006)

For materials with open 3 $d$ or 4 $f$ shells the valence electrons are relatively localized around the atomic sites, and these materials therefore exhibit strong electron-electron interaction effects. For this class of materials, which includes e.g. many different types of high- $T_{\mathrm{c}}$ superconductors, the combination of density functional theory (usually based on the local density approximation (LDA)) and dynamical mean-field theory (LDA+DMFT)Georges et al. (1996) has been the method of choice. DMFT provides a good description of onsite correlations but neglects the intersite correlations. Furthermore, in LDA+DMFT calculations, the local interactions are often treated as adjustable parameters, and the combination of a density-functional based and a diagrammatic scheme requires the introduction of a double-counting parameter, which is supposed to compensate the local correlation effects already contained in the LDA bandstructure. This parameter, which can have a substantial effect on the simulation results,Karolak et al. (2010) is difficult if not impossible to define in a consistent manner. For this reason, LDA+DMFT cannot provide a true ab initio description of materials.

Systematic procedures such as the constrained random phase approximation (cRPA) Aryasetiawan et al. (2004) in principle allow one to calculate the interaction parameters appropriate for LDA+DMFT type calculations, by taking into account the screening effect of the bands outside the low-energy subspace considered in the LDA+DMFT description. However, these interactions typically have substantial nonlocal contributions, and in a recent publication we showedBoehnke et al. (2016) that the resulting non-local screening has a big influence on the effective local interactions. This is even the case for materials such as SrVO3, which were considered to be strongly correlated.Pavarini et al. (2004); Lechermann et al. (2006) A proper description of the low-energy model thus requires an extended dynamical mean field theory (EDMFT) formalism,Si and Smith (1996); Smith and Si (2000); Chitra and Kotliar (2001); Sun and Kotliar (2002) in which both self-energy and polarization effects are treated in a consistent manner. A natural way to perform ab initio simulations based on EDMFT is to combine this scheme with the $GW$ method.Biermann et al. (2003) This formalism, which we will call $GW$ +EDMFT in this paper,111This formalism has been introduced in Biermann et al. (2003) as $GW$ +DMFT. Here, we use the term $GW$ +EDMFT to distinguish this approach from recent implementations that discard the screening of the impurity interaction involves a fully self-consistent calculation of the interaction parameters, and takes into account the effects of local and nonlocal screening. As a fully diagrammatic scheme, it also does not suffer from the LDA+DMFT type double counting problems, and it is a promising formalism for the nonequilibrium simulation of strongly correlated materials.Golež et al. (2017)

While the $GW$ +EDMFT method was formulated more than a decade ago,Sun and Kotliar (2002); Biermann et al. (2003) challenges associated with the numerical treatment of retarded interactions have prevented its implementation for many years.Werner and Casula (2016) Recently, some non-self-consistent,Sakuma et al. (2013) or partially self-consistent Tomczak et al. (2014); Choi et al. (2016) calculations as well as self-consistent model studiesKarlsson (2005); Ayral et al. (2013); Huang et al. (2014); van Loon et al. (2016); Stepanov et al. (2016a); Ayral et al. (2017) have been presented, but the fully self-consistent ab initio scheme has been realized so far only in Ref. Boehnke et al., 2016. An important issue with regard to self-consistency is to what extent the onsite vertex corrections provided by EDMFT counteract the detrimental effects of self-consistency in $GW$ . In this paper we describe the details of our multitier $GW$ +EDMFT implementation and test it on materials with different degrees of correlation.

The paper is organized as follows. In Sec. II we review the basic theory of the $GW$ +EDMFT method. The method is derived from the free-energy functional $\Psi$ and we show that the causality breakdown of the hybridization function that has been reported using a Baym-Kadhanoff version of $GW$ +DMFT in Ref. Lee and Haule, 2017 for the hydrogen dimer, is related to the lack of the bosonic self-consistency in the latter implementation and will likely not be present in the complete $\Psi$ -derivable theory. We also discuss the multitier approach in general terms and clarify the relation between the present approach and commonly used methods such as LDA+DMFT. We then proceed to a detailed discussion of our implementation in Sec. III. In Sec. IV we use an exactly solvable dimer model to analyze some general causality features of the effective impurity interaction and in Sec. V present and discuss the results of the full $GW$ +EDMFT calculations. First, the method is applied to the moderately correlated perovskite SrMoO3, where we find that, contrary to LDA+DMFT which cannot reproduce the experimental satellite features,Wadati et al. (2014) this material is well described by $GW$ +EDMFT. By comparing the spectra with the effective impurity interaction and fully screened interaction we can distinguish Hubbard bands from plasmonic satellites and thereby deduce the nature of the satellite features. By systematically changing the occupation we also use SrMoO3 as a model to investigate the causality of the effective impurity interaction and relate these results to our more general causality-considerations in Sec. IV.

We then focus on sodium as a model system to investigate the performance of the method in different regimes. By successively increasing the lattice constant we investigate the effects of self-consistency and long-range screening for situations with different degrees of local correlations. We show that the method performs well in the moderately to strongly correlated regime and captures the Mott-Hubbard metal to insulator transition at some critical value of the streching. Furthermore we show that the self-consistency is essential to capture the correct trend in the impurity interaction. Section VI summarizes the most important findings.

II Method

II.1 GW-approximation

By expanding the self-energy in terms of the screened Coulomb interaction Hedin derived the following set of exact coupled differential equations defining the self-energy $\Sigma$ and the Green’s function $G$ in terms of the polarization $\Pi$ and vertex function $\Gamma$ Hedin (1965),

[TABLE]

The above equations are for finite temperature with the notation $1=(\tau_{1},\mathbf{r}_{1})$ , etc., and the spin has been omitted for simplicity. By approximating the vertex function by the first term, $\Gamma(1,2,3)\approx\delta(1-2)\delta(2-3)$ , we arrive at the well known $GW$ approximation for the self-energy,

[TABLE]

and the random-phase approximation for the polarization function,

[TABLE]

With a suitable choice of basis functions the integral equations (2) and (5)-(7) can be mapped to matrix equations which can be treated in computer codes with standard linear algebra libraries.

II.2 Extended Dynamical Mean-Field Theory

By introducing a localized basis set $w_{n\mathbf{R}}(\mathbf{r})$ , where $n$ is an orbital index and $\mathbf{R}$ a site index, it is possible to separate the correlations into local (onsite) correlations and nonlocal (offsite) correlations. The key assumption in EDMFT is that all correlations are local. In practice this means that the selfenergy and polarization can be described using only the local basis functions. While the selfenergy and the Green’s functions are one-particle quantities that can be expanded directly in the local one-particle basis, the polarization, which is a two-particle quantity, requires a local product basis

[TABLE]

where, $\alpha=(i,j)$ . For realistic materials the localized basis is typically chosen as linear muffin tin orbital (LMTO) basis functionsAryasetiawan and Gunnarsson (1998) or maximally localized Wannier functions (MLWFs)Marzari and Vanderbilt (1997); Mostofi et al. (2008) that are derived from the LDA band structure.

In EDMFT the full lattice problem is mapped to an impurity problem with a dynamical bare propagator $\mathcal{G}$ and interaction $\mathcal{U}$ Werner and Casula (2016). These so-called Weiss-fields are determined self-consistently by requiring that the local part of the lattice Green’s function ( $G_{\mathrm{loc}}$ ) and fully screened interaction $W_{\mathrm{loc}}$ (as defined by the projection onto the local one- and two- particle basis) should be equal to their impurity counterparts, $G_{\mathrm{imp}}$ and $W_{\mathrm{imp}}$ , respectively,

[TABLE]

In Refs. Hügel et al., 2016; Golež et al., 2017 these self-consistency conditions were formally derived by constructing the free energy functional of the impurity $\Gamma^{\prime}$ and considering the variation of $\Gamma-\Gamma^{\prime}$ , where $\Gamma$ is the lattice system free energy functional.

The impurity action is given by

[TABLE]

and the EDMFT self-consistency cycle takes the following form:

Start with an inital guess for $\Sigma^{\mathrm{imp}}$ and $\Pi^{\mathrm{imp}}$ . 2. 2.

Use these for the local quantities $\Sigma^{\mathrm{loc}}=\Sigma^{\mathrm{imp}}$ and $\Pi^{\mathrm{loc}}=\Pi^{\mathrm{imp}}$ (EDMFT approximations) . 3. 3.

Use $\Sigma_{\mathbf{k}}=\Sigma^{\mathrm{loc}}$ and $\Pi_{\mathbf{q}}=\Pi^{\mathrm{loc}}$ . 4. 4.

Calculate $G^{\mathrm{loc}}=\sum_{\mathbf{k}}\left(G^{(0)}_{\mathbf{k}}{}^{-1}-\Sigma_{\mathbf{k}}\right)^{-1}$ and $W^{\mathrm{loc}}=\sum_{\mathbf{q}}v_{\mathbf{q}}\left(\mathbbm{1}-\Pi_{\mathbf{q}}v_{\mathbf{q}}\right)^{-1}$ , 5. 5.

Use $G^{\mathrm{imp}}=G^{\mathrm{loc}}$ and $W^{\mathrm{imp}}=W^{\mathrm{loc}}$ (EDMFT self-consistency conditions). 6. 6.

Calculate the fermionic Weiss field

[TABLE]

and the effective impurity interaction

[TABLE] 7. 7.

Numerically solve the impurity problem to obtain $G^{\mathrm{imp}}$ and the impurity charge susceptibility $\chi^{\mathrm{imp}}$ . 8. 8.

Use the current $\mathcal{G}$ and $\mathcal{U}$ to calculate $\Sigma^{\mathrm{imp}}=\mathcal{G}^{-1}-G^{\mathrm{imp}}{}^{-1}$ and $\Pi^{\mathrm{imp}}=\chi^{\mathrm{imp}}\left(\mathcal{U}\chi^{\mathrm{imp}}-\mathbbm{1}\right)^{-1}$ . The fully screened interaction $W^{\mathrm{imp}}=\mathcal{U}-\mathcal{U}\chi^{\mathrm{imp}}\mathcal{U}$ only enters the calculations through the self-consistency condition in step 5. 9. 9.

Go back to step 2.

In view of pure model calculations, but also with simulations of cubic $t_{2g}$ materials in mind, it is instructive to explicitly determine the product basis representation of the Kanamori interactionKanamori (1963)

[TABLE]

It is straightforward to show that this is the special case of the general rank-4 tensor representation of the interaction

[TABLE]

with the choice $U=U_{\mathrm{K}}$ ,

[TABLE]

(At the borders or the matrix, we show the indices of the product basis according to Eq. (8).)

We first note that (15) will have a determinant of zero and is hence not invertible. While this forbids the use the inverted forms of the bosonic Dyson equation, e.g. $\mathcal{U}^{-1}=\Pi+W_{\mathrm{loc}}^{-1}$ , this does not pose a problem because all steps of the EDMFT self-consistency cycle can be formulated without, as shown above.

II.3 $GW$ +EDMFT

II.3.1 Self-energy and polarization

While (E)DMFT treats the strong local correlations in a nonperturbative fashion to all orders, the non-local correlations are omitted. Due to the separability of the vertex function in Eq. (3) into the sum of the trivial vertex function used in $GW$ and the non-trivial vertex correction, the Hedin equations for the self-energy (Eq. (1)) and the polarization (Eq. (4)) also separate into the $GW$ contribution and the contribution from the vertex corrections:

[TABLE]

Within the combined $GW$ +EDMFT scheme $\Sigma^{\mathrm{vc}}$ and $\Pi^{\mathrm{vc}}$ are approximated by their impurity counterparts. In this scheme the double counting is well defined and is simply the local projection of the $GW$ self-energy and polarization, respectively,

[TABLE]

The self-consistency cycle is equivalent to the EDMFT self-consistency cycle (see Sec. II.2) but with $\Sigma(\mathbf{k})$ and $\Pi(\mathbf{k})$ in step 3 replaced by the corresponding quantities defined in Eqs. (18)-(19).

II.3.2 Functional derivation

In Ref. Biermann et al., 2003 it was shown that the $GW$ +EDMFT formalism can be rigorously derived as an approximation to the free-energy functional $\Psi$ ALMBLADH et al. (1999); Chitra and Kotliar (2001). In this section we briefly review the derivation following Refs. Biermann et al., 2003; Chitra and Kotliar, 2001; Ayral et al., 2013; Rohringer et al., 2017 and also make a connection with the Baym-Kadanoff formulation of $GW$ +DMFT in, e.g., Ref Lee and Haule, 2017.

The partition function in the grand canonical ensemble for interacting electrons moving in the crystal potential $V_{c}$ is given by

[TABLE]

where $S$ is the action

[TABLE]

The electron-electron interaction term is decoupled using a Hubbard-Stratonovic transformation which yields an additional bosonic field $\phi$ ,

[TABLE]

Here $V_{H}$ includes both the crystal and the Hartree potential and $\alpha$ is a coupling constant that is set to 1 for the physical case. Now we introduce additional source fields that couple to the fermionic and bosonic propagators, which yields the total action

[TABLE]

The free energy of the system in the presence of the external fields is defined as

[TABLE]

By defining the fermionic Green’s function $G(x,y)=-\langle T\psi(x)\psi^{\dagger}(y)\rangle=\frac{\partial\Omega}{\partial J_{f}}$ and the corresponding bosonic propagator $W(x,y)=\langle T\phi(x)\phi(y)\rangle=\frac{2\partial\Omega}{\partial J_{b}}$ and performing a double Legendre transform of $\Omega$ we obtain the free energy functional $\Gamma$ :

[TABLE]

where $G_{H}$ is the Hartree Green’s function of the solid and $\Psi$ contains all further contributions,

[TABLE]

Physically, $\Psi$ includes all two-particle irreducible diagrams constructed with the electron-boson vertex. Setting the source terms to zero and requiring stationarity of $\Gamma$ yields the Dyson equations

[TABLE]

from which we can identify the self-energy and the polarization as

[TABLE]

A clear advantage of the functional formalism is that the derived methods satisfy conservation laws in the Baym-Kadanoff sense. However, when the self-consistency is restricted to a subspace of the full Hilbert space, the conservation laws may not be fulfilled anymore. The $GW$ approximation corresponds to the lowest order approximation to $\Psi$ in $W$ :

[TABLE]

The EDMFT functional is defined by making a local approximation to $G$ and $W$ . It should be noted that by local in this context we mean an onsite approximation in the localized basis defined in Section II.2, and the local product basis (Eq. (8)) for the bosonic quantities. This yields the full $GW$ +EDMFT functional:

[TABLE]

where the last term constitutes the doublecounting term between the $GW$ and EDMFT-functionals. Figure 1 gives a visual representation of these terms of the diagrams.

II.3.3 Double counting

The double-counting terms for the polarization and selfenergy are obtained by evaluating the functional derivatives in Eqs. (29)-(30) for $\Psi^{DC}$ defined in Eq. (32):

[TABLE]

A summation over repeated indices is assumed and the matrix elements are taken in the localized one- and two-particle basis functions defined in Sec. II.2, which for $W$ corresponds to

[TABLE]

For the case where the orbital subspaces for the EDMFT and $GW$ calculations are the same, $\Sigma^{\mathrm{DC}}$ reduces to the local projection ( $\mathbf{k}$ -sum) of the full $GW$ selfenergy. If the orbital subspace for the EDMFT calculation is smaller than that of the $GW$ calculation the difference between $\Sigma^{\mathrm{DC}}$ and the local projection of the full $GW$ self-energy is

[TABLE]

where $r$ is the part of the full Hilbert space that is not included in the EDMFT subspace.

In Ref. Lee and Haule, 2017 a different formulation of $GW$ +DMFT based on the Baym-Kadanoff functional was suggested and evaluated for the hydrogen dimer. This formulation does not include the bosonic self-consistency of the full formulation presented above and therefore the double-counting between the $GW$ and the DMFT parts is different:

[TABLE]

where $U$ is the impurity interaction, which in this case is fixed to the local model interaction, and all sums are restricted to the DMFT orbital subspace. In Ref. Lee and Haule, 2017 it was shown that this formulation yields non-causal hybridization functions in the strongly correlated regime due to the non-causality of the difference between the local projection of the $GW$ selfenergy and the double-counting term

[TABLE]

An alternative double-counting was then introduced to replace this derived double-counting with a causal one.

The full $\Psi$ -derivable formulation of $GW$ +EDMFT on the other hand is not expected to experience the same causality problem, since in that case the double-counting that follows naturally from the derivation coincides with the ‘causal’ double-counting introduced in Ref. Lee and Haule, 2017 when the orbital subspaces for the EDMFT and $GW$ calculations are the same. Furthermore, if the EDMFT subspace is smaller than the $GW$ subspace the difference $\left(\Sigma^{\mathrm{loc}}_{GW}-\Sigma^{\mathrm{DC}}\right)$ reduces to the expression in Eq. (36), which is causal by construction.

II.4 Multitier Self-Consistent GW+EDMFT

In the present multitier $GW$ +EDMFT the full Hilbert space is divided into three subspaces and each subspace is treated with an appropriate level of approximation. The aim of our approach is to accurately describe systems that have both strong local correlation and nonlocal correlation effects, typically 3 $d$ - or $4d$ -compounds such as transition-metal oxides, transition metals, and high $T_{\mathrm{c}}$ cuprate superconductors, and to do so at a reasonable computational cost. A common feature of many of these compounds is strongly correlated partially filled 3 $d$ -states mixed with less correlated extended $s$ or $p$ states. While local vertex contributions from the DMFT-type impurity problem are needed for the $d$ states, the $s$ and $p$ states are typically well described within the $GW$ approximation or even LDA. However, due to the nontrivial mixing with the correlated $d$ states the $s$ and $p$ states need to be included in the self-consistency cycle with an appropriate double counting term. These considerations suggest adopting the following three-step procedure:

•

TIER III: Perform a one-shot $G^{0}W^{0}$ calculation in the full Hilbert space. Choose a basis on an intermediate subspace, typically 3-8 bands, and calculate the effective interaction $U(\omega)$ for this subspace using the constrained random phase approximation. $\Sigma^{G^{0}W^{0}}$ is kept in TIER III only.

•

TIER II: Within the intermediate subspace the self-energy and polarization are calculated self-consistently using a fully self-consistent $GW$ implementation

•

TIER I: Choose a correlated subspace, smaller or equal to the intermediate subspace, for which local vertex contributions are calculated using an EDMFT-type impurity problem.

Self-consistent calculations are performed in the intermediate subspace and local vertex contributions from EDMFT are included inside the correlated subspace at each step in the self-consitency loop. The multitier $GW$ +EDMFT scheme is illustrated in Fig. 2.

Within the multitier approach all double countings are well defined. The full Green’s function is given by

[TABLE]

and the corresponding equation for the bosonic propagators by

[TABLE]

The explicit $\mathbf{q}$ -dependence of the polarization and self-energy, which couples back to the impurity problem in the self-consistency cycle, makes the $GW$ +EDMFT solution dependent on the details of the bare lattice problem (e. g. the geometry of the lattice or the $\mathbf{q}$ -dependence of the bare interaction) to a much greater extent than what one typically finds in approximations like DMFT, where the observables primarily depend on integrated quantities like the bare local density of states.

II.5 Connecting the tiers

The multitier approach allows one to systematically remove the least relevant degrees of freedom from the description and to replace them by an effective medium into which the more relevant degrees of freedom are embedded. The general recipe for this procedure does not depend on the actual separation into ‘less relevant’ and ‘more relevant’ spaces, the only prerequisite is that there exists a physically motivated approximation for the matrix-elements of the self-energy and polarization outside of the ‘more relevant’ space. To make this more explicit, let the available degrees of freedom like lattice site and orbital be $S$ and $s\subset S$ denote some small subspace thereof. At this point there are no restrictions for these subspaces. They may represent the local part, a limitation to a certain number of (correlated) bands or any combination thereof.

In the following, an uppercase character represents an index in $S$ , a lowercase character an index in $s$ . The same character in upper- and lowercase represents the same index. For example $A_{CD}+B_{cd}$ would be a shorthand notation for

[TABLE]

With this notation, the usual Dyson equation can be written as

[TABLE]

with

[TABLE]

The only restriction of the partitioning of $\Sigma_{CD}$ into its two constituents is that $\Sigma^{\mathrm{d}}$ must be non-zero on $s$ only. We do not impose any conditions on $\Sigma^{\mathrm{r}}$ ; specifically, it can be zero or non-zero on $s$ . That this is an exact rewriting can be understood from the expression

[TABLE]

If the aim is to evaluate Eq. (45) on $s$ , it can be rewritten to be completely contained within that space,

[TABLE]

thus providing the bare propagator $\widetilde{G}^{0}_{ab}$ on the $s$ space, which is effectively retarded through the channels in $S$ that are not contained in $s$ ( $S\setminus s$ ). It is easy to show that this reduces to the DMFT formalism when $S$ is chosen as the whole lattice, $s$ a single site of that and $\Sigma^{\text{r}}$ to be the periodic continuation of $\Sigma^{\text{d}}$ .

Equations (43) through (48) are equally valid for $W$ , $U$ , and $\Pi$ instead of $G$ , $G^{0}$ , and $\Sigma$ , which leads to a description of the the effectively retarded $s$ -space interaction. Care has to be taken in that case of the rank 4 structure of the involved quantities, e.g., using the handling described in Sec. III.1.1.

With this formalism in place, cRPA can be identified as “ $S$ : all bands, $s$ : a limited subset of those, and $\Pi^{r}$ the $G^{0}G^{0}$ polarization bubble outside $s$ .” EDMFT, in addition to the DMFT construction discussed above, includes the analog considerations for the effectively retarded local interaction.

The general formulation of the multitier approach is a versatile tool, and the implementation with $G^{0}W^{0}$ @TIER III + $GW$ @TIER II + QMC@TIER I used in this paper and beforeBoehnke et al. (2016) harvests its full potential with the currently available tools. This is just a snapshot of the development however and as better tools become available as solvers on the respective tiers, they can be systematically incorporated.

Table 1 lists a number of formalisms that can be understood within this framework. Some complimentary approaches, like dual fermion/dual boson (DF/DB) based formalismsRubtsov et al. (2008, 2012); Stepanov et al. (2016b) and QUADRILEXAyral and Parcollet (2016b) are related, but cannot directly be cast into the multitier form. DF/DB are derived from a reformulation of the lattice action, which yields a different form of the Dyson equations with a non-local bare propagator and while QUADRILEX has a two-particle self-consistency cycle that modifies the impurity interaction, this is accomplished without introducing an explicit polarization in bosonic variables.

III Computational details

III.1 TIER III-II

In TIER III the one-shot $GW$ selfenergy as well as the effective interaction on the intermediate subspace ( $U(\omega)$ ) are calculated within cRPA using the all electron FLAPW code SPEX Friedrich et al. (2010); The FLEUR group . The intermediate subspace for the self-consistency loop on TIER II and I is defined using MLWFs as implemented in the Wannier90 library.Marzari and Vanderbilt (1997); Mostofi et al. (2008); Freimuth et al. (2008); Sakuma (2013) TIER II is treated using a custom finite-temperature Matsubara axis self-consistent $GW$ implementation on the intermediate subspace.

III.1.1 TIER II: Matrix Formulation at Finite Temperature

As single-particle basis for TIER II we use MLWFs,

[TABLE]

where $w_{n\mathbf{R}}(\mathbf{r})$ is a Wannier function centered in the unit cell specified by $\mathbf{R}$ and $n$ is an orbital index.

The two-particle quantities $W$ and $\Pi$ require a product basis. In this work we use a restricted product basis of the form

[TABLE]

where $\alpha=\{i,j\}$ . Note that in general $\langle\Psi_{\alpha\mathbf{q}}(\mathbf{r})|\Psi_{\beta\mathbf{q}}(\mathbf{r})\rangle\neq\delta_{\alpha\beta}$ . Formally, it is possible to work in this nonorthogonal basis in the following way: Define the matrix elements of the screened interaction as $W_{\alpha\beta}$ from downfolding of the two-particle object $W(\mathbf{r},\mathbf{r}^{\prime})$ by

[TABLE]

Further define the matrix elements of the polarization $\Pi_{\alpha\beta}$ implicitly by

[TABLE]

Hence, $\Pi_{\alpha\beta}$ is an object that is upfoldable to a real-space representation $\Pi(\mathbf{r},\mathbf{r}^{\prime})$ . The nonorthogonality of the product basis $\Psi_{\alpha}(\mathbf{r})$ ( $\mathbf{k}$ -index dropped for simplicity) implies that

[TABLE]

However, this does not prevent us from mapping equations of the type

[TABLE]

to matrix operations in the product basis notation by applying Eq. (51):

[TABLE]

where $d\mathbf{R}=d\mathbf{r}_{1}d\mathbf{r}_{2}d\mathbf{r}d\mathbf{r}^{\prime}$ . As a rule of thumb, convolutions alternating between upfoldable and downfolded quantities map to matrix operations. This mapping to matrix operations also implies a corresponding definition for the inversion of these quantities.

With these definitions and using the Wannier basis in Eq. (49) for the single-particle quantities and the product basis in Eq. (50) for the two-particle quantities the integral equations (2) and (5)-(7) can be mapped to the following matrix equations (summation over repeated indices assumed):

[TABLE]

In TIERs II and I it has to be ensured that a sufficiently large number of Matsubara frequencies and a fine enough imaginary time discretization are used. The substantial spectral weight at large energies that TIER II inherits in its bare propagators from TIER III requires a much larger number than for e.g. LDA+DMFT calculations. For SrMoO3 and SrVO3 discussed in later sections, we used at least 1500 Matsubara frequencies.

III.1.2 Gamma point handling

Both the bare Coulomb interaction and the screened Coulomb interaction diverge as $1/\mathbf{k}^{2}$ for $k\rightarrow 0$ . Because of this the $\Gamma$ -point needs to be handled with special care. In the $GW$ calculation using the FLAPW code SPEX in TIER III the treatment of the $\Gamma$ -point is greatly simplified by making use of the analytic forms of both the bare Coulomb interaction and the bare Green’s function. By a basis change to the Coulomb eigenbasis the divergence is restricted to a single eigenvalue and can be treated separately using so-called kp-perturbation theory.Friedrich et al. (2010) Furthermore, the value of the divergent term of the Coulomb interaction is redefined as the integral of the corresponding term in a small region around the $\Gamma$ -point, whose size is determined by the $\mathbf{k}$ -mesh:

[TABLE]

With this definition the integral of $V$ over the Brillouin zone will be correctly reproduced by the $\mathbf{k}$ -point sampling. By making use of the analytic form of the bare polarization a similar expression can be derived for $U(\omega)$ Friedrich et al. (2010):

[TABLE]

where the constant $c(\omega)$ is derived from the head element of the polarization matrix, $\langle E_{0}|\Pi|E_{0}\rangle$ , with $|E^{0}\rangle$ being the eigenfunction that corresponds to the divergent eigenvalue of the Coulomb interaction.

In TIER II we perform a self-consistent $GW$ calculation on the intermediate subspace starting from a downfolded frequency-dependent interaction $U(i\omega_{n})$ and retarded bare propagator $G^{0}(i\nu_{m})$ calculated in TIER III, and therefore we lack analytic expressions for both the Green’s function and the interaction. This forces us to resort to cruder approximations for the $\Gamma$ -point than in TIER III.

When $U(i\omega_{n})$ is projected onto the intermediate subspace the divergent contribution will no longer be contained in a single matrix element, but will in general give a contribution to all matrix elements of $U(i\omega_{n},\mathbf{k}=\Gamma)$ . Furthermore, since the intermediate subspace is much smaller than the space spanned by the complete product basis in TIER III it is not possible to isolate the divergent term by a simple rotation.

The dielectric function $\epsilon(\mathbf{k})$ , on the other hand, is a smooth function of $\mathbf{k}$ and will generally have an extremal point for $\mathbf{k}=\Gamma$ and can therefore be approximated as a constant function of $\mathbf{k}$ in a small region around the $\Gamma$ -point. Hence, if the $\mathbf{k}$ -mesh is chosen sufficiently dense, the inverse dielectric function may be approximated as

[TABLE]

where $\bar{\epsilon}^{-1}(i\omega_{n})$ is the average value of $\epsilon^{-1}(i\omega_{n},\mathbf{k})$ at the boundaries defined by $\mathbf{k}_{\mathrm{cutoff}}$ . Using this approximation for the inverse dielectric function the fully screened interaction can readily be calculated as

[TABLE]

where the notation $\bar{U}$ is used to emphasize that the divergent term of the interaction is redefined according to Eq. (61).

In Fig. 3 we compare the fully screened interaction calculated in the custom $GW$ implementation for TIER II with the results from SPEX for the benchmark material SrVO3. One can see that the local part of the screened interaction agrees very well with the results from SPEX. Also the $\Gamma$ -point contribution shows remarkable agreement in spite of the rather crude approximations in TIER II.

III.1.3 The Hartree contribution

During the self-consistent calculations the density within the intermediate subspace, and thus the Hartree energy, will change. In Eq. (40) the Hartree contribution was hidden in the self-energies. In this section we derive the Hartree contribution for a downfolded model and discuss the technical details of the current implementation.

The Hartree potential is given by

[TABLE]

We now take the matrix elements of the above expression and use

[TABLE]

Inserting the definition of the basis functions [Eq. (49)] into the above equation yields

[TABLE]

As long as the Wannier functions are sufficiently localized this term will be dominated by the $\mathbf{R}_{1}=\mathbf{R}_{2}$ and $\mathbf{R}_{3}=\mathbf{R}_{4}$ contribution, in which case the Hartree term reduces to

[TABLE]

Since a Hartree contribution is already contained in the LDA Hamiltonian in TIER III we only need to consider the correction to the Hartree term in TIER II due to the change in the density in the intermediate subspace ( $\Delta V^{\mathrm{H}}_{ij}$ ). Since the change in the density integrates to zero the divergent term arising from $\mathbf{q}=0$ in the Coulomb potential is eliminated,

[TABLE]

The additional tilde on the interaction matrix is used to emphasize that the constant (divergent) eigenfunction of the interaction has been removed before projecting onto the Wannier basis.

III.2 TIER I

The relevant low-energy local problem in the presence of effective bare propagators and interactions screened through all the degrees of freedom that were systematically removed in TIERs I and II is primarily due to the retarded effective interaction nontrivial to solve and some approximations cannot be avoided at this point. In our calculations we employ the CT-Hyb Werner et al. (2006); Werner and Millis (2007); Gull et al. (2011) quantum Monte-Carlo algorithm as implemented in Alps.Bauer et al. (2011); ALPS collaboration ; Hafermann et al. (2013) This implementation features the treatment of dynamic screening of the interaction by a numerically cheap reweighting of the Monte-Carlo configurations, as introduced in Ref Werner and Millis, 2010. This approach handles the retardation effects exactly, as long as the screening affects only the monopole part of the interaction. To prevent systematic errors from the fitting of high frequency tails we make use of the Legendre polynomial based compact representation of the Green’s function.Boehnke et al. (2011) This is of particular importance in multitier $GW$ +EDMFT, where structure can be expected up to very high frequencies. A fitting procedure is at risk of losing those features.

The ALPS solver,Hafermann et al. (2013) while efficient, restricts the impurity interactions to density-density like terms (“segment picture”) and it allows only the measurements of density-density like contributions to the impurity charge susceptibility.

The first restriction has nontrivial consequences, since the $\mathcal{U}_{abab}$ terms give rise to the $-J\delta{\sigma\sigma^{\prime}}n_{a\sigma}n_{b\sigma^{\prime}}$ contribution to (II.2) due to the equal-spin term in Eq. (14), which is of density-density form and should be retained, but also give rise to the spin-flip contribution for the opposite-spin term in (14), which can not be dealt with in the implementation of Ref. Hafermann et al., 2013. Thus, there will be a discrepancy between the $\mathcal{U}$ -tensor used in the self-consistency and the Hamiltonian that ultimately enters the impurity problem.

The second restriction, although more severe at first glance, turns out to be of minor importance. The orbital-resolved impurity charge susceptibility $\chi^{\text{imp}}$ enters the EDMFT scheme only through the formula

[TABLE]

from which it becomes apparent that limiting $\chi^{\mathrm{imp}}$ to $\chi^{\mathrm{imp}}_{iijj}$ type contributions merely means that only the $U$ and $U^{\prime}$ will be screened while $J$ remains unscreened, which is a physically reasonable approximation (it is much harder to screen $l\neq 0$ components of a charge distribution). To partially remedy this restriction, the non-density-density components of the local polarization are in our calculations taken into account on the $GG$ -bubble level.

Additionally, the solver used in the current implementation allows for an orbital diagonal hybridization only, which is however the case by construction for all systems discussed in this paper. It inhibits us however from applying the formalism to a cluster of atoms as the local $\mathcal{C}$ space on TIER I. Recent developments aim to remove some of these restrictions Steiner et al. (2015) by introducing a segment-picture-based expansion in the hybridization and Hund’s coupling parameter $J$ . This approach can efficiently handle the density-density components of retarded $U$ , $U^{\prime}$ , and $J$ interactions, while retarded spin-flip terms produce a severe sign-problem. Other choices for impurity solvers that can deal with the retarded interaction and bare propagator include the non- and one-crossing approximations (NCA/OCA)Grewe and Keiter (1981); Haule et al. (2002); Golež et al. (2015); Werner and Eckstein (2013) that are viable options in the large- $U$ limit. However, since the self-energies produced by these approximate solvers do not contain all the local diagrams that have been removed in TIER II, but only a partial summation, numerical artifacts can be expected.Golež et al. (2017)

III.3 Self-consistency cycle

In Fig. 4 the full flow of the calculations is illustrated. First a DFT calculation is performed using the FLAPW-code FLEUR.The FLEUR group Then, a low-energy model is defined using MLWFsMarzari and Vanderbilt (1997); Mostofi et al. (2008); Freimuth et al. (2008); Sakuma (2013) and the model interaction as well as the $G^{0}W^{0}$ selfenergy is computed with the SPEX codeFriedrich et al. (2010); The FLEUR group . This defines the bare propagators for the intermediate subspace [ $G_{I,\mathbf{k}}^{0}{}^{-1}$ and $U_{I,\mathbf{q}}^{-1}$ in Eqs. (40)-(41)]. With the initial assumptions

[TABLE]

the selfenergy and screened interaction are computed according to Eqs. (40)-(41). If the self-consistency conditions (Eq. 9) are not fulfilled the impurity selfenergy ( $\Sigma^{\mathrm{EDMFT}}\big{|}_{C,\mathrm{loc}}$ ) and polarization ( $\Pi^{\mathrm{EDMFT}}\big{|}_{C,\mathrm{loc}}$ ) are computed using the Alps CT-Hyb impurity solver.Werner et al. (2006); Gull et al. (2011) These values replace the initial guesses above, and yield a new Green’s function and screened interaction. The scheme is then iterated until the self-consistency conditions are fulfilled. This cycle on TIERs II and I is implemented using the TRIQS framework Parcollet et al. (2015). In the calculations presented in this paper, we did not have to employ any mixing to converge the results. While multiple solutions may exist Stan et al. (2015), there were no indications of additional unphysical solutions.

III.4 Analytic continuations

All analytic continuations were done using the maximum entropy method (MaxEnt)Bryan (1990); Gubernatis et al. (1991) as implemented in Ref. Boehnke, , except for the results in Fig. 5 which were obtained using a Padé analytic continuation Vidberg and Serene (1977), since MaxEnt cannot handle noncausal spectral functions.

IV Causality considerations

A recurring feature in $GW$ +EDMFT based model and material calculations is that the impurity interaction $\mathcal{U}(z)$ , $z\in\mathbb{C}$ can become non-causal, that is, it may feature poles in the upper half of the complex plane instead of the lower half. This leads to negative spectral weight of the bosonic modes $-\Im\mathcal{U}(\omega)$ . An example is shown in Fig. 5, where $\Im\mathcal{U}(\omega)$ is positive in contrast to $\Im W(\omega)$ , which has the expected negative peak. This behavior can also been observed in Sec. V.1 for SrMoO3, but earlier examples include the simulations of SrVO3 Boehnke et al. (2016) and model studies Ayral et al. (2013); Huang et al. (2014); van Loon et al. (2016); Stepanov et al. (2016a); Ayral et al. (2017), where this property can be seen at least indirectly. From the Hilbert transform of $\mathcal{U}(\omega)$ it follows that

[TABLE]

where $f(\omega)$ is a positive function on the interval $\omega\in[0,\infty)$ . Hence, if $\mathcal{U}(i\omega_{n+1})<\mathcal{U}(i\omega_{n})$ for any Matsubara frequency ( $i\omega_{n}$ ) the spectral weight must be noncausal in some frequency range. The converse is not necessarily true.

From the mean-field value of the (static) screened effective interaction for the $U$ - $V$ model, $U_{\mathrm{eff}}=U-zV$ ( $z$ being the coordination number), one can expect noncausal interactions to appear for peculiar choices of parameters, like an attractive nearest-neighbor interaction $V$ . Given however that the above-mentioned examples with noncausal impurity interactions are not in this regime, the question arises as to whether the noncausality is a consequence of the $GW$ approximation for the nonlocal components of the self-energy or a generic consequence of DMFT-type local approximations.

To answer this question, we consider the simple case of a dimer problem, where the definition of the effective impurity model is merely a formal step, and where we can take all nonlocal diagrams into account. Any formalism that maps to an auxiliary (impurity) problem does so by defining the self-consistency conditions. In the case of a derivation from a $\Psi$ functionalv ALMBLADH et al. (1999), this is unambiguouslyGolež et al. (2015) given by

[TABLE]

Dual boson uses the same conditions,Stepanov et al. (2016b) although originally a different procedure was proposed.Rubtsov et al. (2012) In addition to these self-consistency conditions, a local formalism needs to define an approximation to the local $G$ and $W$ through lattice properties and impurity observables. As introduced before, $GW$ +EDMFT takes the impurity self-energy and polarization and augments them by the nonlocal components within the $GW$ approximation. Dual Boson introduces a dual expansion which is then cut at some perturbation order.

In this section we instead use the exact values for the local $G$ and $W$ of a dimer, and determine the impurity bath Green’s function $\mathcal{G}$ and interaction $\mathcal{U}$ such that the solution of the impurity model reproduces these exact local $G$ and $W$ . Specifically, for the dimer with Hamiltonian

[TABLE]

we compute $G_{\mathrm{ex}}$ and $W_{\mathrm{ex}}$ using exact diagonalization.Strand Then the functions $\mathcal{G}$ and $\mathcal{U}$ (i.e. the corresponding bath parameters) are determined numerically in such a way that the solution of the impurity model reproduces the local parts $G_{\mathrm{ex,loc}}$ and $W_{\mathrm{ex,loc}}$ . This is achieved by executing a self-consistency loop with the local $G$ and $W$ fixed to their exact values (here $i$ is the iteration)

[TABLE]

Figure 5 shows the impurity interaction on the real frequency axis for the dimer parameters $U=8$ , $V=1$ , $t=5$ and $\beta=5$ . The results were broadened slightly for better visibility. The imaginary part is purely noncausal in this case, while for model and material calculations in $GW$ +EDMFT we typically see a mixed behavior of causal or noncausal spectral weight at low frequencies and causal spectral weight at high frequencies. Consistent with the strong antiscreening mode, the static impurity interaction is significantly increased relative to the bare local interaction $U$ .

Let us stress that the arguments presented here do not depend on any particular non-local approximation. The exact result for the dimer shows that given the self-consistency equations (75) and (76), the auxiliary impurity problem can have a noncausal retarded interaction. Hence, this is an intrinsic possibility of any local approximation.

V Results

V.1 The cubic perovskites SrMoO3 and SrVO3

V.1.1 Spectral functions

SrMoO3 crystallizes in a cubic perovskite structure.Macquart et al. (2010) The conduction states originating from the Mo 4 $d$ states are of $t_{2g}$ character and occupied by two electrons. In LDA the $t_{2g}$ states form well-isolated bands around the Fermi energy with a bandwidth of roughly 3.8 eV. SrMoO3 is in many respects similar to the previously studied SrVO3 Boehnke et al. (2016) but has twice the filling of the $t_{2g}$ shell and also a slightly larger bandwidth. The larger bandwidth is expected since the Mo 4 $d$ -states in SrMoO3 are less localized than the corresponding V 3 $d$ states in SrVO3.

The calculations in this section were performed for the inverse temperature $\beta=15\frac{1}{\mathrm{eV}}$ corresponding to the temperature $\approx 774$ K.

The photoemission spectra (PES) of SrMoO3 measured in Ref. Wadati et al., 2014 (Inset in Fig. 6) show a well defined quasiparticle peak together with a weak shoulder structure around $-2.5$ eV. Even though the specific heat coefficient is renormalized to approximately twice its LDA value the PES spectrum does not show any clear band narrowing compared to LDA. The hump around $-2.5$ eV is not seen in bandstructure calculations and therefore presents a clear sign of electron correlations within the $t_{2g}$ band. However, in Ref. Wadati et al., 2014 it was shown that this structure cannot be reproduced in LDA+DMFT using any realistic values of the instantaneous local Hubbard interaction $U$ . Only with an unphysically large interaction of 5 eV a similar structure appeared in the LDA+DMFT spectra, but the band narrowing in this case was much too large. This led the authors to speculate that the shoulder structure was of plasmonic rather than Hubbard band character.Wadati et al. (2014)

The $GW$ +EDMFT approach includes both the strong local correlations and the long-range screening effects within a single unified framework and can therefore capture both Mott physics and plasmonic excitations, as well as sub-plasmonic features originating from fluctuations within the correlated space, on an equal footing. Furthermore, the multitier approach used in the present work does not include any adjustable parameters and therefore provides an ideal tool to distinguish plasmonic features from Hubbard bands. When applied to SrVO3 we already demonstrated in Ref. Boehnke et al., 2016 that the inclusion of long-range correlations changes the interpretation of the side-bands; the satellites in SrVO3 should be interpreted as plasmonic satellites rather than Hubbard bands. The situation in SrMoO3 seems even more clearcut with no side-bands at all in plain DMFT.

In Fig. 6 we show the local spectral function for SrMoO3 computed with the present multitier self-consistent $GW$ +EDMFT scheme together with the experimental photoemission spectra from Ref. Wadati et al., 2014. The occupied part agrees well with the photoemission spectra and the shoulder structure around $-2.5$ eV is clearly visible. We also predict a satellite feature in the unoccupied part of the spectrum, potentially visible in inverse photoemission experiments, centered at roughly 3 eV.

The plasmonic signature for SrMoO3 in the experimental data is more pronounced than what we find in our $GW$ +EDMFT calculations. The same is true for SrVO3, where early measurements suggestedMorikawa et al. (1995) a very strong lower satellite but later measurements with higher photon energies reveal a much reduced satellite intensity compared with the quasiparticle weightSekiyama et al. (2004). However, it is only very recently,Backes et al. (2016) that the experimental ratio in spectral weight between the lower satellite and the occupied part of the quasiparticle has been reexamined and found to be closer to 1:3, which is in reasonable agreement with the $GW$ +EDMFT calculations in Ref. Boehnke et al., 2016 that suggest a ratio of $\sim$ 1:4. The difference to previous measurements has been attributed to oxygen vacancies in the system. A similar situation seems to be present in other transition metal oxides like SrTiO3Backes et al. (2016) and its interfaces.Dudy et al. (2016); Sing et al. (2017) Furthermore, extrinsic loss that can be large for plasmonic satellitesKas et al. (2016), is not included in our calculations. We propose that the difference in spectral weight between the photoemssion experiments and our ab initio results in the inset of Fig. 6 may have a similar origin.

In Fig. 6 we also show the spectral function obtained from the $G^{0}W^{0}$ approximation. The $G^{0}W^{0}$ results agree very well with the $GW$ +EDMFT results except for the position of the satellite in the unoccupied region that is overestimated within the single-shot $GW$ approximation. The quasiparticle dispersion in the $\mathbf{k}$ -resolved spectral function (Fig. 7) also agrees with the $GW$ results and shows a small band-narrowing compared to the LDA. The dispersion of the satellites follows the dispersion of the quasiparticle bands. This behavior is similar to the behavior of plasmon satellites in the $G^{0}W^{0}$ +cumulant expansion which suggests that the satellites are of plasmonic origin Caruso et al. (2015); Lischner et al. (2015); Caruso and Giustino (2015).

V.1.2 Effective interactions

To fully answer the question about the origin of the satellite features we look at the effective impurity interaction and the fully screened interaction in Figs. 8 and 9. The static value of $\mathcal{U}$ is approximately 2.8 eV, which is clearly too small to explain the satellite features as Hubbard bands, since the separation of the satellites is approximately 5.5 eV. On the other hand the positions of the satellites agree with the pole in $W$ in Fig. 9. This leads to the conclusion that the satellites are indeed of plasmonic orgin and hence SrMoO3 and SrVO3 are qualitatively similar, despite the qualitatively different predictions based on LDA+DMFT calculations. It is noteworthy that for SrMoO3 the plasma frequency is almost identical in $G^{0}W^{0}$ and $GW$ +EDMFT while for SrVO3 the plasma frequency is reduced in the full $GW$ +EDMFT compared to one-shot $G^{0}W^{0}$ .Boehnke et al. (2016) Thus, RPA based on the LDA-bandstructure works better for SrMoO3 which can be expected since this is a more extended system and hence for this compound $W^{0}$ is expected to yield a better plasma frequency.

In the following, we use SrMoO3 as an example to discuss the frequency dependence and analytic properties of the effective impurity interaction $\mathcal{U}$ . Figure 8 shows the different components of $\mathcal{U}$ as defined in Eq. (15). It is clear that the screenings of $\mathcal{U}$ and $\mathcal{U}^{\prime}$ are similar, while $J$ remains almost unscreened, which justifies the approximations in TIER I.

From Eq. (74) it directly follows that if $\mathcal{U}(i\omega_{n})$ is smaller than the static value for any Matsubara frequency $\omega_{n}$ then $\Im\mathcal{U}(\omega^{\prime})$ has to be noncausal in some finite frequency range. This is clearly the case for SrMoO3 in Fig. 8. We further note that the effective impurity interaction for SrVO3, shown in the same figure, does not only exhibit signs of a pole in the upper half plane, for SrVO3 the pole in $\mathcal{U}$ is located on the Matsubara axis between the first and second Matsubara frequencies, yielding a kink structure in $\mathcal{U}$ . This is not a problem per se, since poles in the upper half plane of $\mathcal{U}$ are to be expected (see Sec. IV) and thus might end up on the Matsubara axis for certain parameters. Yet, extra care has to be taken in this case to avoid numerical instabilities.

V.1.3 Filling dependence

To investigate the screening behavior and the noncausality further we have computed $\mathcal{U}(i\omega_{n})$ for SrMoO3 with different fillings of the $t_{2g}$ manifold; see Fig. 10. To reduce the computational time we held $\Sigma_{r}$ and $\Pi_{r}$ fixed at the experimental filling and only recalculated the model polarization and self-energies for the new fillings. This is a reasonable approximation since all effects from within the low-energy subspace have been removed in $\Sigma_{r}$ and $\Pi_{r}$ and these quantities are therefore relatively insensitive to small changes of the chemical potential. For low fillings, $n\leq 1$ , $\mathcal{U}(i\omega_{n})$ does not display any clear signs of noncausality. As the filling is increased towards the experimental filling ( $n$ =2) an antiscreening mode develops at low frequencies, which gives an increase in the static value and a negative slope of $\mathcal{U}(i\omega_{n})$ between the first and second Matsubara frequencies. The static value continues to increase and reaches its maximum near half-filling (see inset). To interpret the trend in $\mathcal{U}$ we also have to consider the fully screened interaction ( $W$ ) for different fillings (Fig. 10). First of all one can note that even if the effective impurity interaction is noncausal for certain fillings, $W$ remains causal, as expected. The static value of $W$ follows the opposite trend to that of $\mathcal{U}$ ; i.e., it reaches its minimum at half filling, reflecting the increased screening as the number of free charge carriers is increased.

At self-consistency $W_{\mathrm{loc}}$ is obtained by screening $\mathcal{U}$ by the local polarization, $W_{\mathrm{loc}}=\mathcal{U}+\mathcal{U}\Pi_{\mathrm{loc}}W_{\mathrm{loc}}$ . Thus by comparing the two panels of Fig. 10 one can deduce that the local screening is strongest close to half filling. It is also interesting to note that $\mathcal{U}$ becomes noncausal at low frequencies for the cases where the local screening is strong.

V.2 Sodium as a model system

V.2.1 Results for different lattice constants

Sodium provides an ideal playground for testing our approach. Elemental Na has the electronic configuration [Ne]3s1 and crystallizes in a bcc structure. In the solid the 3 $s$ -states hybridize with the unoccupied 3 $p$ states to form a broad conduction band. The conduction states are very delocalized and close to an electron gas model. The main features of the occupied part of the experimental spectra is a well defined quasiparticle peak, a plasmonic satellite feature around $-7$ eV, which is repeated at approximately $-14$ eV Steiner et al. (1979).

The calculations in this section were performed for the inverse temperature $\beta=10\frac{1}{\mathrm{eV}}$ corresponding to $\approx 1160$ K.

In the current work we will use sodium as a model system to scan different degrees of correlation. By successively increasing the lattice constant we increase the degree of local correlations in a controlled manner. This analysis is similar in spirit to the calculations on stretched diatomic molecules which are commonly used to benchmark quantum chemistry methods.

To faithfully reproduce the low-energy band structure we have to consider a 4-band model, consisting of the 3 $s$ and 3 $p$ states (see Fig. 11). For this material we will utilize the full strength of the multitier approach by choosing the $sp$ 4-band subspace for the self-consistent $GW$ calculation but only considering local EDMFT corrections for the $s$ -state. Hence, the intermediate subspace is spanned by the full $s$ and $p$ Wannier functions but the correlated subspace is restricted to the $s$ -like Wannier function. We consider the experimental lattice constant $a_{0}$ , as well as the artificially increased lattice constants $1.4a_{0}$ and $1.6a_{0}$ .

When the lattice constant is increased the bandwidth of the conduction band is decreased, reflecting a decreasing hopping amplitude between the neighboring sites. The model interaction, $U$ (Fig. 12), is almost static for all lattice constants, implying that the model includes all important screening channels and hence also the dominant correlation contributions to the self-energy. Contrary to what one might expect, both $U$ and the bare interaction decreases as the lattice constant increases. A decreasing bare interaction implies that the localization of the Wannier basis states around the atomic positions is weaker for the larger lattice constants. Albeit counterintuitive, this phenomenon is well known and has been investigated for model systems in Ref. Tomczak et al., 2009 and later for manganese monoxide under pressure in Ref. Tomczak et al., 2010. Since $U$ is only weakly screened the static value will follow the same trend as the bare interaction.

The effective impurity interaction $\mathcal{U}$ on the other hand follows the correct trend; i.e., it increases as the lattice constant is increased, reflecting the increasing importance of the local correlations. That $\mathcal{U}$ and $U$ follow different trends illustrates the importance of the nonlocal $s$ - $s$ screening as well as the $s$ - $p$ screening channels for this system.

We will next discuss the fully screened interactions (Fig. 13) and the spectral functions (Fig. 14) for the different cases, comparing them both to self-consistent $GW$ and the one-shot $G^{0}W^{0}$ . In the weakly correlated regime, where plasmonic physics dominates, the screened interaction already provides useful information about the spectral function:

A peak in Im $W(\omega)$ at $\omega=\omega_{\mathrm{p}}$ will give a corresponding satellite feature in the spectral function at $E\pm n\omega^{\prime}$ , where $n$ is an integer and $E$ is the energy of the quasiparticle peak. The plus (minus) sign refers to unoccupied (occupied) states. 2. 2.

For a given quasiparticle peak $G^{0}W^{0}$ will only give a single peak in the spectral function below and above the Fermi energy at too high energy. 3. 3.

The renormalization of the quasiparticle peak will have a nontrivial dependence on the frequency and weight of the peak in Im $W$ . A peak with a large weight at low frequency will generally give the largest quasiparticle renormalization. However the $\mathbf{k}$ -dependence of the self-energy also influences the quasiparticle renormalization factor.

For elemental sodium (leftmost panels of Figs. 13 and 14) we know by comparing $G^{0}W^{0}$ calculations with the experimental spectra that the plasma frequency in $W^{0}$ is relatively good. The discrepancies with the experimental spectra are mainly related to self-energy corrections that can be accounted for using the cumulant expansion. Also the quasiparticle renormalization is slightly underestimated in $G^{0}W^{0}$ Aryasetiawan et al. (1996). Self-consistent $GW$ (sc $GW$ in Fig. 13) on the other hand severely worsens the result compared to experiment. The plasmon pole in $W$ is almost completely washed out, similar to what has been found for the electron gas.Holm and von Barth (1998) The local corrections from EDMFT in the full calculations ( $GW$ +EDMFT in Fig. 13) improves the sc $GW$ results but the strength of the pole is still severly underestimated and the position of the pole is at too high energy. This yields a weak plasmon between $-20$ to $-10$ eV in the spectral function ( $a_{0}$ case in Fig. 14), in poor agreement with experiment. Also the width of the quasiparticle peak is severely overestimated in the $GW$ +EDMFT results. The reason for the poor agreement with experiment is that the nonlocal screening is too big to be accounted for by only the first bubble diagram in the expansion. Thus, to get a good description of these kinds of very weakly correlated electron gas-like metals it is necessary to include higher-order nonlocal screening beyond RPA.

As the lattice constant is increased to 1.4 $a_{0}$ the pole in $W^{0}$ is shifted to lower energies (Fig. 13). This is expected since the plasma frequency in the electron gas can be shown to be proportional to the square-root of the density and we effectively decrease the density by increasing the lattice constant. sc $GW$ still gives a very wide and featureless plasmon similar to the original lattice constant. However, $W$ for the full $GW$ +EDMFT calculation develops a well defined plasmonic pole at slightly higher energy than the pole in $W^{0}$ . Hence, for this lattice constant we enter a regime where the $GW$ +EDMFT approximation, which only takes into account the first bubble diagram in the nonlocal polarization, becomes physically reasonable. In addition to the main peak in $W$ there is an additional shoulder structure around 15 eV. This structure, which is present also in sc $GW$ and the effective impurity interaction $\mathcal{U}$ , is related to the lack of local corrections for the $s$ - $p$ screening channel. However, the structure is smaller in $GW$ +EDMFT than in sc $GW$ which implies that the local corrections for the $s$ - $s$ screening at least partially remove the unphysical high frequency structures in sc $GW$ . Ideally though, all relevant low-energy screening channels should be included in the correlated subspace and only less important screening channels should be treated in sc $GW$ .

It is interesting to note that even though the peak in $W$ is at higher frequency in $GW$ +EDMFT compared to $G^{0}W^{0}$ , the plasmon satellite in the spectral function (Fig. 14) is closer to the quasiparticle peak. Hence, the local vertex contributions to the selfenergy corrects the error in $G^{0}W^{0}$ and pulls the satellite closer to the quasiparticle. There is also an additional satellite feature in the $GW$ +EDMFT results around $-12$ eV. This feature is a combination of a repetition of the main plasmon and a second plasmonic peak due to the shoulder structure around 15 eV in $W$ .

Finally we will discuss the largest lattice constant $1.6a_{0}$ (rightmost panels of Figs. 13 and 14). In this case the conduction band is well separated from the other bands and almost of pure $s$ -character (see Fig. 11). This means that we can directly compare the bandwidth of the conduction band with the static value of the impurity interaction in Fig. 12 to estimate the degree of local correlations. Since $\mathcal{U}$ is much larger than the bandwidth we get an insulating solution with an upper and lower Hubbard band separated by approximately the static value of the interaction (Fig. 13). $G^{0}W^{0}$ clearly fails to capture the strong local correlations driving the metal to insulator transition and yields a metallic solution with a plasmonic satellite feature below the Fermi energy. Due to the gap in the spectral function there are no excitations below $\approx 3$ eV, which means that Im $W$ is zero in this frequency range. There is a peak in $W$ at around 4 eV corresponding to transitions between the lower and upper Hubbard band and an additional peak from the $s$ - $p$ screening at higher frequency. In the spectral function for the $1.6a_{0}$ case (Fig. 14) there is a weak (barely visible) satellite feature corresponding to the first peak in Im $W$ .

V.2.2 Effect of a local approximation in extended systems

In the $GW$ +EDMFT results for elemental sodium above we found that the plasmon weight in $\Im W$ was reduced substantially and the position of the pole was pushed to higher energies. In sc $GW$ we found a similar reduction, but here the pole was pushed to lower energies instead. To investigate these differences in behavior, which must originate from the local EDMFT contributions, we first consider the extreme case of approximating $\Pi_{\mathbf{q}}=\Pi_{\mathrm{loc}}$ in a simple one-shot $G_{0}W_{0}$ -calculation for sodium and compare it to the regular $G_{0}W_{0}$ result with the full $\mathbf{q}$ -dependent polarization $\Pi_{\mathbf{q}}$ (Fig. 15). This clearly exhibits the same trend as the full $GW$ +EDMFT results, albeit even more extreme: By only including a local polarization the plasmon peak in $W$ is shifted to much higher energies and is broadened substantially. Thus, in cases where the local polarization does not have more physical relevance than the non-local terms, the $GW$ +EDMFT plasma frequency might end up getting overestimated. Namely, for such a compound the nonlocal part is simply too large to be treated with only the first bubble diagram in the self-consistent expansion. Including all local diagrams but only the first-order nonlocal diagram leads to an underestimation of the $\mathbf{q}$ -dependence compared to the local contribution, which in turn shifts the position of the pole in $W$ to higher frequencies. However, as is evident from Fig. 15, the $GW$ +EDMFT result should still be much better than what one would expect to get with an EDMFT-type approximation with only a local polarization.

To understand this behavior in more detail we have to look at the specific form of $\Pi_{\mathbf{q}}$ and $\Pi_{\mathrm{loc}}$ . We do this for the non-interacting electron gas where the analytic form of the polarization is known explicitly (Figs. 16 and 17)Fetter and Walecka (2003).

The plasmon pole in $W_{\mathbf{q}}$ occurs at the zeros of the dielectric function

[TABLE]

From Fig. 16 one can see that Eq. (80) is only fulfilled at the second crossing between the line $1/V_{\mathbf{q}}$ and $\mathrm{Re}\Pi_{\mathbf{q}}(\omega)$ , since the imaginary component of the polarization is big at the first crossing. For the chosen parameters (see the caption of Fig. 16) the position of the plasmon pole in $W_{\mathbf{q}}$ will have a small dispersion between approximately 6-10 eV which yields a sharp peak in $W_{\mathrm{loc}}$ in Fig. 18. Furthermore, there will not be a well defined plasmon peak for all $\mathbf{q}$ -values. For $q/k_{F}\gtrsim 1$ the line $1/V_{\mathbf{q}}$ will not cross $\mathrm{Re}\Pi_{\mathbf{q}}(\omega)$ , which results in a relatively broad and weak plasmonic feature for these $\mathbf{q}$ -points. Another interesting observation is that the peak in $\mathrm{Im}\Pi_{\mathbf{q}}(\omega)$ is sharper and shifted to lower frequencies for the $\mathbf{q}$ -points close to the $\Gamma$ -point compared to the large- $\mathbf{q}$ components as can also be understood from Fig. 19.

To define a quantity corresponding to the local polarization for the electron gas we choose a cubic unit cell with the same volume as the bcc unit cell for sodium. We then define $\Pi_{\mathrm{loc}}$ as the $\mathbf{q}$ -sum of $\Pi_{\mathbf{q}}$ over the first Brillouin zone. It should be noted that this definition of “local polarization” differs slightly from the local polarization in the full ab initio calculations, where the local subspace is defined by MLWF:s and the screening from all bands are included. For the electron gas we use a plane-wave basis and only screening from the first (conduction) band is included in $\Pi_{\mathrm{loc}}$ . Hence, the following discussion should only be used to gain a qualitative understanding of Fig. 15 while some quantitative differences such as the exact position and weight of the plasmon pole might differ.

The local polarization (Fig. 17) contains a large contribution from the $\mathbf{q}$ -points far away from the $\Gamma$ -point and therefore the peak in $\mathrm{Im}\Pi_{\mathrm{loc}}(\omega)$ is broad and peaked at a relatively high frequency. This gives a correspondingly broad Kramers-Kronig feature in the real part at high energies. For many $\mathbf{q}$ -points there will not be a real crossing between $1/V_{\mathbf{q}}$ $\mathrm{Re}\Pi_{\mathrm{loc}}(\omega)$ , but if there is it will occur at energies larger than 10 eV. This yields a broad plasmon in $W$ at high frequency (Fig. 18), just as we observed for Na in Fig. 15.

In Fig. 19 we show the dispersion of the plasmon together with the imaginary part of the polarization for the noninteracting electron gas. The two cases shown in Fig. 16 correspond to two points in Fig. 19. The dispersion of the plasmon (solid curve in Fig. 19) is defined as the $(\mathbf{q},\omega)$ point where $|\Pi(\mathbf{q},\omega)-1/V_{\mathbf{q}}|$ takes its minimum value [which for the case $\Im\Pi_{\mathbf{q}}(\omega_{p})=0$ is given by the crossing between the line $1/V_{\mathbf{q}}$ and $\mathrm{Re}\Pi_{\mathbf{q}}(\omega)$ as discussed above]. When $\Im\Pi_{\mathbf{q}}(\omega=\omega_{p})\neq 0$ the strength of the plasmon in $W$ will be reduced and the pole is broadened, corresponding to a finite lifetime of the plasmonic mode (e.g., Landau damping).

We can also use the non-interacting electron gas to gain a qualitative understanding of the effect of increasing the lattice constant for Na. By increasing the lattice constant we effectively shrink the first Brillouin zone (1:st BZ). Therefore the maximum difference between $\Pi_{\mathbf{q}}$ for different $\mathbf{q}\in$ 1:st BZ will decrease, that is $\Pi_{\mathbf{q}}$ becomes more local. At some point the nonlocal components of $\Pi_{\mathbf{q}}$ are sufficiently small that a truncation at the first bubble diagram is appropriate. This is the point where $GW$ +EDMFT becomes justified. At a much later point the nonlocal components of $\Pi_{\mathbf{q}}$ are sufficiently small to be ignored completely, in which case an EDMFT treatment is sufficient.

VI Conclusions

While the $GW$ +EDMFT method was first proposed in 2003 the implementation of the fully self-consistent scheme for real materials was not realized until 2016 (Ref. Boehnke et al., 2016). In this paper we provided a detailed description of the self-consistent multitier $GW$ +EDMFT implementation used in Ref. Boehnke et al., 2016 and we tested the $GW$ +EDMFT method for different systems with a focus on the effect of self-consistency. We first applied the $GW$ +EDMFT method to SrMoO3, a cubic perovskite with 4 $d$ valence electrons. SrMoO3 is in many respects similar to the 3 $d$ cubic perovskite SrVO3. However, while SrVO3 has previously been thought to be well described by LDA+DMFT, with a renormalized quasiparticle peak and Hubbard side-bands,Pavarini et al. (2004); Lechermann et al. (2006); Backes et al. (2016) LDA+DMFT obviously fails to give a reasonable description of SrMoO3. The interaction needed to produce the observed satellite features yields a much too strong renormalization of the quasiparticle peak, which indicates that the satellites could be of plasmonic rather than Hubbard band character.Wadati et al. (2014) In this work we have shown that the parameter-free multitier $GW$ +EDMFT scheme is able to describe the satellite features of SrMoO3 and that the satellites are indeed of plasmonic origin. These results can also be connected to our previous investigation of SrVO3 and further support the main conclusion reached there,Boehnke et al. (2016) namely that $GW$ +EDMFT is able to describe the satellite structures in terms of plasmonic fluctuations in a moderately correlated material. We also used SrMoO3 as a model system to provide insights into the effectiveness of screening and the causality of the impurity interaction by systematically changing the occupation of the $t_{2g}$ manifold. Close to half-filling the local screening is strong, which yields a noncausal effective impurity interaction while all physical observables remain causal.

Finally we used sodium as a model system to investigate the performance of the method for different degrees of local correlations. Starting from a weakly correlated metal with the experimental lattice constant, we successively increased the lattice constant and thus the degree of local correlations. This allowed us to scan the metal to insulator transition in a controlled manner for a realistic material and to evaluate the performance of the multitier $GW$ +EDMFT approach in a wide parameter range. We showed that the method performs well in the moderately to strongly correlated regime but underperforms in the very weakly correlated regime where the nonlocal screening is comparable to the local one. In the latter case it is no longer sufficient to treat the nonlocal correlations within the sc $GW$ approximation, i.e., nonlocal diagrams beyond RPA are needed. The $G^{0}W^{0}$ approximation on the other hand, works relatively well in this regime, which implies that the partial cancellation of diagrams in $G^{0}W^{0}$ is not restricted to the local case. In the intermediate to strongly correlated regime the local vertex contributions to the self-energy remedy the well known problems with the one-shot $G^{0}W^{0}$ approximation, such as the overestimation of the satellite position in the spectral function compared to the pole in $\Im W$ . $GW$ +EDMFT also correctly captures the Mott-Hubbard metal to insulator transition. The quasiparticle bandwidth is reduced compared to sc $GW$ and, contrary to the latter scheme, where the plasmon pole in $\Im W$ is washed out, the $GW$ +EDMFT calculation yields well defined plasmonic peaks in $\Im W$ . The static model interaction, computed with the cRPA, exhibits the wrong trend for stretched sodium; it decreases as the lattice constant is increased. In this work we showed that the additional screening channels included in $GW$ +EDMFT solve this problem and provide the physically expected result, namely an impurity interaction which increases with increasing lattice constant.

On the conceptual level, a noteworthy feature in most previous self-consistent $GW$ +EDMFT calculations, both for models Ayral et al. (2013); Huang et al. (2014); van Loon et al. (2016); Stepanov et al. (2016a); Ayral et al. (2017) and real materials,Boehnke et al. (2016) is that the effective interaction for the impurity problem $\mathcal{U}$ can become noncausal. Using an exactly solvable dimer model, we showed that noncausal impurity interactions are not a specific feature of $GW$ +EDMFT, and that any method that includes dynamic long-range polarizations is bound to yield noncausal effective impurity interactions in certain parameter regimes. However, we also found that quantities that relate to physical observables, as opposed to auxiliary ones, remain causal.

The results of our study demonstrate that multitier self-consistent $GW$ +EDMFT enables the ab initio simulation of a broad range of materials in the intermediate to strongly correlated regime. The method is relatively cheap in terms of computational effort and memory requirements, and hence applicable to multiorbital systems. It does not contain any adjustable parameters, once the different subspaces, or tiers, have been fixed. In particular, the method provides a self-consistent calculation of the dynamically screened interaction parameters, and captures the effect of collective charge excitations. This distinguishes the method from standard LDA+DMFT simulations, and makes it the first true ab initio simulation method for moderately to strongly correlated materials. As a fully Green’s function based approach, multitier self-consistent $GW$ +EDMFT is also a promising method for the study of nonequilibrium phenomena in correlated solids.

Acknowledgements.

F.N. and F.A. acknowledge financial support from the Swedish Research Council (VR). L.B. is supported by the Swiss National Science Foundation through NCCR MARVEL. P.W. acknowledges support from ERC Consolidator Grant No. 724103. The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at LUNARC and at the CSCS Dora cluster provided by MARVEL. LB would like to thank Denis Golež and Thomas Ayral for valuable discussions and Hugo Strand for providing his exact diagonalization code.Strand F.N. and L.B. contributed equally to this work.

Bibliography94

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Hedin (1965) L. Hedin, Phys. Rev. 139 , A 796 (1965) . · doi ↗
2Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140 , A 1133 (1965) . · doi ↗
3Aryasetiawan and Gunnarsson (1998) F. Aryasetiawan and O. Gunnarsson, Reports on Progress in Physics 61 , 237 (1998) . · doi ↗
4van Schilfgaarde et al. (2006) M. van Schilfgaarde, T. Kotani, and S. Faleev, Phys. Rev. Lett. 96 , 226402 (2006) . · doi ↗
5Lyo and Plummer (1988) I.-W. Lyo and E. W. Plummer, Phys. Rev. Lett. 60 , 1558 (1988) . · doi ↗
6Bruneval et al. (2006) F. Bruneval, N. Vast, and L. Reining, Phys. Rev. B 74 , 045102 (2006) . · doi ↗
7Georges et al. (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68 , 13 (1996) . · doi ↗
8Karolak et al. (2010) M. Karolak, G. Ulm, T. Wehling, V. Mazurenko, A. Poteryaev, and A. Lichtenstein, Journal of Electron Spectroscopy and Related Phenomena 181 , 11 (2010) , proceedings of International Workshop on Strong Correlations and Angle-Resolved Photoemission Spectroscopy 2009. · doi ↗

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Taxonomy

Multitier self-consistent GW+EDMFT

Abstract

pacs:

I introduction

II Method

II.1 GW-approximation

II.2 Extended Dynamical Mean-Field Theory

II.3 GWGWGW+EDMFT

II.3.1 Self-energy and polarization

II.3.2 Functional derivation

II.3.3 Double counting

II.4 Multitier Self-Consistent GW+EDMFT

II.5 Connecting the tiers

III Computational details

III.1 TIER III-II

III.1.1 TIER II: Matrix Formulation at Finite Temperature

III.1.2 Gamma point handling

III.1.3 The Hartree contribution

III.2 TIER I

III.3 Self-consistency cycle

III.4 Analytic continuations

IV Causality considerations

V Results

V.1 The cubic perovskites SrMoO3 and SrVO3

V.1.1 Spectral functions

V.1.2 Effective interactions

V.1.3 Filling dependence

V.2 Sodium as a model system

V.2.1 Results for different lattice constants

V.2.2 Effect of a local approximation in extended systems

VI Conclusions

Acknowledgements.

II.3 $GW$ +EDMFT