Internal consistency of multi-tier $GW$+EDMFT – Supplementary Material

I $GW$+EDMFT formalism

I.1 $GW$ approximation

The self-energy $\Sigma$ encodes how correlations affect the propagation of electrons in a solid. To avoid problems with diverging self-energy diagrams when expanding in terms of the bare Coulomb interaction $v$, Hedin [1] introduced an expansion in powers of the screened interaction $W$. The first two terms of the self-energy expansion read:

$$\Sigma(1,2)=iG(1,2)W(1^{+},2)-\int G(1,3)G(3,4)G(4,2)W(1,4)W(3,2)d3d4+\ldots,$$ (1)

where $G$ is the interacting electron Green’s function and $(i,j)$ are space-time indices. Keeping only the leading order term $\Sigma(1,2)=iG(1,2)W(1^{+},2)$ corresponds to the well-known $GW$ approximation, which partially captures non-local correlations via polarization effects. This can be seen by noting that $W$ obeys a bosonic Dyson equation

$$W(1,2)=v(1,2)+\int v(1,3)\Pi(3,4)W(4,2)d3d4,$$ (2)

where $\Pi(1,2)$ is the electron polarization function which in the Random Phase Approximation can be obtained as

$$\Pi(1,2)=-iG(1,2)G(2,1^{+}).$$ (3)

Equations (2) and (3), plus the first order term in Eq. (1) and the corresponding Fermionic Dyson equation $G(1,2)=G_{0}(1,2)+\int G_{0}(1,3)\Sigma(3,4)G(4,2)d3d4$ (with $G_{0}$ the noninteracting Green’s function) define the $GW$ self-consistency loop. In practical first principles calculations, one replaces the exchange correlation potential in the DFT Kohn-Sham propagators by the single-shot or self-consistent (quasi-particle) $GW$ self-energy, which has been shown to yield improved descriptions of weakly and moderately correlated materials.

I.2 Extended Dynamical Mean Field Theory

Within (E)DMFT, correlations are assumed to be local (i.e. the self-energy and polarization are momentum independent) and the system is self-consistently mapped onto an impurity model with a bath Green’s function $\mathcal{G}$ (or hybridization function $\Delta$) and an effective dynamically screened bare interaction $\mathcal{U}$ [2, 3]. $\mathcal{G}$ and $\mathcal{U}$ are determined self-consistently, by imposing that the physical impurity correlation functions are equal to their local lattice counterparts: $G^{\mathrm{imp}}=\frac{1}{N_{k}}\sum_{k}G^{\mathrm{lat}}(k)$ and $W^{\mathrm{imp}}=\frac{1}{N_{k}}\sum_{k}W^{\mathrm{lat}}(k)$. The EDMFT self-consistency loop contains the following steps:

1. Make an initial guess for the impurity self-energy and polarization $\Sigma^{\mathrm{imp}}$, $\Pi^{\mathrm{imp}}$.

2. Identify the impurity self-energy and polarization with the corresponding lattice quantities (correlations are assumed to be local within EDMFT): $\Sigma^{\mathrm{lat}}=\Sigma^{\mathrm{imp}},\Pi^{\mathrm{lat}}=\Pi^{\mathrm{imp}}$.

3. Compute the local lattice Green’s function and screened interaction:

   • $G^{\mathrm{loc}}=\frac{1}{N_{k}}\sum_{k}G^{\mathrm{lat}}(k)=\frac{1}{N_{k}}\sum_{k}(G_{0}^{-1}(k)-\Sigma^{\mathrm{lat}})^{-1}$,

   • $W^{\mathrm{loc}}=\frac{1}{N_{q}}\sum_{q}W^{\mathrm{lat}}(q)=\frac{1}{N_{q}}\sum_{q}v(q)(1-\Pi^{\mathrm{lat}}v(q))^{-1}$.

4. Impose the self consistency conditions: $G^{\mathrm{imp}}=G^{\mathrm{loc}}$ and $W^{\mathrm{imp}}=W^{\mathrm{loc}}$.

5. Compute the fermionic and bosonic dynamical mean fields:

   • $\mathcal{G}=({G^{\mathrm{imp}}}^{-1}+\Sigma^{\mathrm{imp}})^{-1}$,

   • $\mathcal{U}=W^{\mathrm{imp}}(1+\Pi^{\mathrm{imp}}W^{\mathrm{imp}})^{-1}$.

6. Solve the impurity problem with these dynamical mean fields, i.e., calculate the impurity Green’s function $G^{\mathrm{imp}}$ and the charge susceptibility $\chi^{\mathrm{imp}}$.

7. Update the impurity self-energy and polarization:

   • $\Sigma^{\mathrm{imp}}=\mathcal{G}^{-1}-{G^{\mathrm{imp}}}^{-1}$,

   • $\Pi^{\mathrm{imp}}=\chi^{\mathrm{imp}}(\mathcal{U}\chi^{\mathrm{imp}}-1)^{-1}$.

8. Go back to step 2 until convergence is reached.

I.3 Multi-tier $GW$+EDMFT

$GW$+EDMFT considers a self-energy and polarization with both local and nonlocal components. Due to the separability of the vertex function in the Hedin equations [2], the self-energy and polarization can be written as sums of their $GW$ parts and the vertex corrections: $\Sigma=\Sigma^{GW}+\Sigma^{\textrm{vc}}$ and $\Pi=\Pi^{GW}+\Pi^{\textrm{vc}}$. The vertex corrections are approximated by EDMFT: $\Sigma^{\textrm{vc}}\rightarrow\Sigma^{\mathrm{EMDFT}}$ and $\Pi^{\textrm{vc}}\rightarrow\Pi^{\mathrm{EMDFT}}$. Since $GW$ already contains local self-energy and polarization components, which are more accurately described by EDMFT, we subtract these components to avoid a double counting. This yields the following corrections to the momentum dependent self-energy and polarization:

$$\Sigma^{GW+\text{EDMFT}}(k)=\Sigma^{GW}(k)-\frac{1}{N_{k}}\sum_{k}\Sigma^{GW}(k)+\Sigma^{\text{EDMFT}},$$ (4)

$$\Pi^{GW+\text{EDMFT}}(q)=\Pi^{GG}(q)-\frac{1}{N_{q}}\sum_{q}\Pi^{GG}(q)+\Pi^{\text{EDMFT}},$$ (5)

where $\Sigma^{GW}$ and $\Pi^{GG}$ are defined in Sec. I.1. The self-consistency cycle is equivalent to the one described in Sec. I.2, with the replacement of $\Sigma^{\mathrm{lat}},\Pi^{\mathrm{lat}}$ by Eqs. (4) and (5).

In the multi-tier scheme, different energy windows are treated at different levels of theory. In the “full” space (Tier III), composed of $\sim 100$ bands, we consider the one-shot self-energy $\Sigma^{G^{0}W^{0}}$ and polarization $\Pi^{G^{0}G^{0}}$. The intermediate space (Tier II), spanning typically $\sim 10$ bands, is solved with self-consistent $GW$. Within the correlated space (Tier I), which can be equal to or smaller than the intermediate space, we treat the local correlations with EDMFT. The Tier I and Tier II spaces considered in the current study are illustrated in Fig. 1. The Tier III subspace corresponds to the $d$-orbital space. In practice, the basis used to describe the intermediate and correlated spaces is that of Maximally Localized Wannier Functions (MLWFs), obtained by Wannierizing an appropriate low-energy sector of the DFT bands. The one-shot $\Sigma^{G^{0}W^{0}}$ and $\Pi^{G^{0}G^{0}}$ are also rotated to this basis and define the “bare” propagators and interactions in Tier II.

Figure 1: Different model spaces used in the $GW$+EDMFT simulations of the cubic Perovskites. The most correlated $t_{2g}$ and $t_{2g}+e_{g}$ local subspaces are treated with EDMFT, while their non-local parts, as well as the lower lying Oxygen $p$-states are treated by self-consistent $GW$.

I.4 One-shot $G^{0}W^{0}$ and Constrained Random Phase Approximation

We perform a one-shot $G^{0}W^{0}$ calculation in the TIER III space encompassing $\sim 100$ bands with the SPEX code [4] and use the resulting Green’s functions as the bare propagators of the intermediate TIER II space [2].

The effective frequency- and momentum-dependent bare interaction of the intermediate space is computed within the constrained Random Phase Approximation (cRPA) [5], which excludes screening processes occurring within the target space.

II Susceptibility calculations

II.1 Density-density response function

The density-density correlation function $\chi_{n_{q}n_{-q}}(t)$, defined by [6]

$$\chi_{n_{q},n_{-q}}(t)=-i\lim_{\eta\rightarrow 0^{+}}\int_{0}^{\infty}\langle[n_{q}(t),n_{-q}]\rangle e^{i(\omega+i\eta)t}dt\,$$ (6)

probes the charge excitations of a system and in the noninteracting case is related to the polarization function

$$\Pi_{ijkl}(q,i\nu_{m})=\frac{1}{\beta}\sum_{i\omega_{n},k}G_{ik}(k,i\omega_{n})G_{lj}(k-q,i\omega_{n}-i\nu_{m})$$ (7)

via

$$\chi_{n_{q},n_{-q}}(q,\omega)=\sum_{ij}\Pi_{iijj}(q,\omega),$$ (8)

where $i,j$ are orbital indices.

The polarization function probes the density response to the total potential, $\Pi(1,2)\sim\frac{\delta n(1)}{\delta V_{\mathrm{tot}}(2)}$, where $V_{\mathrm{tot}}=V_{\mathrm{ext}}+V_{\mathrm{ind}}$ (with $V_{\mathrm{ind}}$ the Hartree potential), while the physically measurable charge susceptibility is the density response to the external potential, $\chi(1,2)\sim\frac{\delta n(1)}{\delta V_{\mathrm{ext}}(2)}$. The two are related via $\chi=\Pi+\Pi v\chi$, where $v$ is the bare Coulomb interaction.

Within the $GW$+EDMFT formalism the polarization function $\Pi$ is defined in Eq. (5) and the effective bare Coulomb interaction is taken to be the cRPA $U$, yielding

$$\chi^{GW+\mathrm{EDMFT}}(q,i\nu_{n})=\Pi(q,i\nu_{n})[1-U^{\mathrm{cRPA}}(q,i\nu_{n})\Pi(q,i\nu_{n})]^{-1}.$$ (9)

Here, the contraction of rank-4 tensor quantities is implemented as a matrix multiplication in a product basis of Wannier functions [2]. The full scalar charge susceptibility is then computed as $\chi^{GW+\mathrm{EDMFT}}_{c}(q,i\nu_{n})=\sum_{i,j}\chi^{GW+\mathrm{EDMFT}}_{ii,jj}(q,i\nu_{n})$. This function can be analytically continued with the Maximum Entropy (MaxEnt) method [7].

II.2 Dielectric function

The matrix denominator in the expression for $\chi$ corresponds to the multi-orbital dielectric function $\epsilon=1-v\Pi$. This function can be probed by Electron Energy Loss Spectroscopy (EELS), which measures a related quantity $-\mathrm{Im}\frac{1}{\epsilon}$, known as the loss function [8]. The experimentally measured dielectric function typically probes all inter- and intra-band transitions. Evaluated with Fermi’s golden rule, its imaginary part reads [9]

$$\mathrm{Im}[\epsilon(q,\omega)]\sim\frac{4\pi^{2}e^{2}}{q^{2}}\sum_{\alpha,\beta,k}|\langle\alpha,k+q|e^{iq\hat{r}}|\beta,k\rangle|^{2}(1-f(E_{\alpha,k+q}))f(E_{\beta,k})\delta(E_{\alpha,k+q}-E_{\beta,k}-\omega)$$ (10)

where $\alpha,\beta$ index the energy bands, $f(E)$ is the Fermi-Dirac distribution function and $v$ is the scalar bare Coulomb interaction $\frac{4\pi^{2}e^{2}}{q^{2}}$. With this in mind, and to facilitate experimental comparisons, we define a scalar dielectric function by summing all the density-density components $\epsilon(q,i\nu_{n})=\sum_{i,j}\epsilon^{GW+\mathrm{EDMFT}}_{ii,jj}(q,i\nu_{n})$, analogously to Sect. II.1. We find that the inverse of this expression is causal and can be analytically continued with MaxEnt.

II.3 Spin response function

One can probe the system’s tendency towards magnetic order by studying its multi-orbital spin susceptibility $\chi_{s}$, which in the RPA formalism is defined as [10]

$$\chi_{s}(q,i\omega_{n})=\chi_{0}(q,i\omega_{n})[\mathbb{I}+U_{s}(i\omega_{0})\chi_{0}(q,i\omega_{n})]^{-1}.$$ (11)

Here, the uncommon plus sign in the denominator originates from the minus sign in the definition of the charge susceptibility in Eq. (6).

The local spin interaction $U_{s}$ is obtained from the decomposition of the generalized multi-orbital Hubbard interaction on a given site,

$$H_{\mathrm{int}}=\frac{1}{4}\sum_{\alpha_{1},...,\alpha_{4}}U^{(\alpha_{1},\alpha_{4}),(\alpha_{3},\alpha_{2})}c^{{\dagger}}_{\alpha_{1}}c^{{\dagger}}_{\alpha_{2}}c_{\alpha_{3}}c_{\alpha_{4}},$$ (12)

into spin and charge channels:

$$U^{(\alpha_{1}\alpha_{4})(\alpha_{3}\alpha_{2})}=-\frac{1}{2}U_{s}^{(l_{1}l_{4})(l_{3}l_{2})}\mathbf{\sigma}_{\sigma_{1},\sigma_{4}}\cdot\mathbf{\sigma}_{\sigma_{2},\sigma_{3}}+\frac{1}{2}U_{c}^{(l_{1}l_{4})(l_{3}l_{2})}\delta_{\sigma_{1},\sigma_{4}}\delta_{\sigma_{2},\sigma_{3}},$$ (13)

where $\alpha_{i}=(l_{i},\sigma_{i})$ combines orbital and spin degrees of freedom. For a Kanamori-type $t_{2g}$ Hamiltonian, the spin and charge interaction matrices read [11]

$$U^{(ll^{\prime})(nn^{\prime})}_{s,c}=\begin{cases}U&\left(l=l^{\prime}=n=n^{\prime}\right)\\ U^{\prime},\,-U^{\prime}+2J&\left(l=n\neq n^{\prime}=l^{\prime}\right)\\ J,\,2U^{\prime}-J&\left(l=l^{\prime}\neq n^{\prime}=n\right)\\ J^{\prime}&\left(l=n^{\prime}\neq n=l^{\prime}\right)\end{cases},$$ (14)

where $U$, $U^{\prime}$, $J$, $J^{\prime}$ are the conventional intra-orbital Coulomb, inter-orbital Coulomb, Hund’s coupling and inter-orbital pair hopping interactions, respectively.

Due to the multi-orbital nature of the model, we define a scalar spin susceptibility $\chi_{s}$ by summing over all the density-density components:

$$\chi_{s}=\sum_{i,j}\chi_{s}^{(ii)(jj)}.$$ (15)

We use the static value of the local bare impurity interaction $\mathcal{U}$ to approximate the values of $U$, $U^{\prime}$ and $J$ to construct $U_{s}$, assuming $J^{\prime}=J$. We note that directly inserting the obtained $U_{s}$ into Eq. (11) results in unphysical values of $\chi_{s}$ due to an eigenvalue of $-U_{s}\chi_{0}$ becoming greater than 1. To avoid this, it is necessary to rescale the spin interaction as $U_{s}\rightarrow\eta U_{s}$. A value of $\eta=0.1$ is used in all subsequent spin susceptibility plots. This renormalization is consistent with the strongly renormalized spin vertex found, for example, in the two-particle self-consistent approach [12, 13]. Similar interaction tunings have been used before in RPA calculations [11].

Figure 2: Bare ($G^{0}G^{0}$) and RPA spin susceptibilities for the three cubic perovskites.

Figure 3: Imaginary part of the local self-energy within different schemes for cubic perovskites in different model spaces.