Merging GW with DMFT and non-local correlations beyond

We review recent developments in electronic structure calculations that go beyond state-of-the-art methods such as density functional theory (DFT) and dynamical mean field theory (DMFT). Specifically, we discuss the following methods: GW as implemented in the Vienna {\it ab initio} simulation package (VASP) with the self energy on the imaginary frequency axis, GW+DMFT, and ab initio dynamical vertex approximation (D$\Gamma$A). The latter includes the physics of GW, DMFT and non-local correlations beyond, and allows for calculating (quantum) critical exponents. We present results obtained by the three methods with a focus on the benchmark material SrVO$_3$.


Introduction
The calculation of materials with predictive power is arguably the biggest challenge of condensed matter theory. In the 20th century we have seen the breakthrough of density functional theory (DFT) [1,2] (for reviews see Refs. [3,4]) which allows for the reliable calculation of many materials and their properties. This is quite surprising considering the fact that the approximations employed to the exchange and correlation potential, such as the local density approximation (LDA) or the generalized gradient approximation (GGA), are rather crude. Despite the success of DFT for many materials, there are entire classes of systems for which it does not work properly. This happens, e.g., for materials, in which exchange or correlation effects are large. Hence the silver bullet of method development is to find better potentials or to improve upon exchange and correlations by many-body methods [5].
Materials in which the exchange part is particularly important are, e.g., semiconductors. Here, DFT within LDA or GGA predicts consistently too small band gaps. This can be overcome by hybrid functionals [6,7,8,9] that mix part of the exact exchange to the exchange correlation functional. The amount of exact exchange that is required for an accurate modeling is, however, non-universal, i.e., material-dependent. For instance, in metals the long-range exchange is screened by long-range charge fluctuations [10]. An accurate many-body framework to capture the system-dependent screening is Hedin's GW approach [11] which calculates the screened-exchange self a e-mail: held@ifp.tuwien.ac.at arXiv:1703.08446v1 [cond-mat.str-el] 24 Mar 2017 In GW the self energy Σ is given by the interacting Green function G (black straight line) times the screened interaction W (red wiggled line) from coordinate/site R i to R j . Second line: The screened interaction W in turn is given by the bare interaction (here denoted as V ) and the screening in the random phase approximation (RPA). This RPA screening is generated by the last term which yields a ladder in terms of V and bubbles consisting of two Green functions. Third line: In DMFT the self energy is given by the local contribution of all Feynman diagrams with the local interaction U always on the same site R i . Fourth line: In AbinitioDΓ A, we take as the irreducible vertex Γ the bare non-local Coulomb interaction V q and the local vertex Γ loc which depends on orbitals (l, m ...) and frequencies (ν, ν , ω) but not momenta (k,k ,q); Γ loc also includes the local Coulomb interaction U (adapted from Ref. [16]). energy from the Green function G times the screened exchange W , see Fig. 1 for the corresponding Feynman diagram. Most GW results have been obtained using a DFT-derived Green function G 0 and an interaction W 0 that has been screened by the Lindhard function computed with G 0 . Only recently self-consistent GW calculations that use an approximate hermitianized form of the self energy, as proposed by van Schilfgaarde and Kotani [12,13], became available. In Section 2 we discuss the GW method and the calculation of the full frequency-dependence of the self energy, which is needed for spectral functions and for a self-consistency beyond the van Schilfgaarde-Kotani approximation. We detail in particular the advantages of our new imaginary-frequency implementation of GW within the Vienna ab initio simulation package (VASP) [14,15]. Both methods, hybrid functionals and GW, can lead to semiconductor band gaps in far better agreement with experiment [17,6,9], with the GW self energy acting as a "scissors operator" [17]. Beyond that, GW also describes quasiparticle renormalizations, finite life times, and improves on the total energies of, e.g., defects [18,10]. While hybrid functionals are one-electron-like by construction, also the GW-at least in all common implementations (see however the recent Refs [19,20])-is based on a Green function G 0 that is always related to a single Slater determinant. Excluding any multi-reference character in G 0 , the G 0 W 0 approach is thus not capable to treat systems in which fluctuations are strong.
Materials in which the correlation part is particularly important are, among others, transition metal oxides and heavy fermion compounds with partially filled d and f shells, respectively [21]. For treating electronic correlations in such materials, dynamical mean field theory (DMFT) [22,23] (for a review see Ref. [24]) and its merger with DFT [25,26] (for reviews see Refs. [27,28]) has been a big leap forward. DMFT takes into account a major part of the electronic correlations: the "local" ones that are confined to a single atomic site. Fig. 1 (bottom) shows the corresponding Feynman diagrams. This way, among others, quasiparticle renormalizations including kinks [29,30,31,32], Hubbard side bands, metal-insulator transitions, and magnetism can be described much more accurately than with one-particle methods, and finite temperature properties become accessible as well. Early successes of DFT+DMFT include the calculation of the Mott-Hubbard transition in V 2 O 3 [33,34,35,36], magnetism in Fe and Ni [25], and the α-γ transition in Ce [37,38]. More recently, it has also been applied to oxide heterostructures [39,40,41], surfaces [42], nanoclusters [43] and oxygen vacancies [44].
The major remaining shortcomings of DFT+DMFT are (i) the sand in the clockwork when interfacing a density functional theory with a Feynman diagrammatic approach and (ii) that only local correlations are taken into account in DMFT. Regarding (i), let us in particular mention the double counting: It is unclear which part of the DMFT correlations are taken into account already on the DFT side, and so different double counting schemes have been proposed. Most commonly used is the fully localized limit [45]. The double counting issue is particularly pronounced in so-called "d + p" DFT+DMFT calculations that in, say, oxides, include both, the transition metal d-as well as the oxygen p-orbitals and can lead to largely different results [46,47,48].
This conceptual problem can be overcome by substituting DFT by GW in the so-called GW+DMFT approach [49,50], which merges two many-body Feynman diagrammatic approaches so that one can precisely identify which diagrammatic contribution is counted twice. GW+DMFT also provides for a better treatment of the exchange contribution. This is not only of advantage for correlated semiconductors such as Ga 1−x Mn x As, and ligand-states in, e.g., transition metal oxides [51], but also for a quantitative description of effective masses of correlated electrons [52]. We discuss the GW+DMFT approach and present results in Section 3; for a more detailed introduction we refer the reader to Refs. [53,54,55,56]. One should note, however, that the treatment of non-local correlations is very limited in GW+DMFT as only charge fluctuations and only the particle-hole channel are included in GW. Moreover they are treated only in weak coupling perturbation theory, i.e., by building the particle-hole ladder only in terms of the bare Coulomb interaction V , see Fig. 1 (middle).
There are essentially two routes that deal with non-local correlations while keeping the local DMFT correlations at the same time: cluster [57,58,59] and diagrammatic extensions [60,61,62,63,64,65,66] of DMFT. The former have been successfully applied to the two dimensional Hubbard model and helped establishing the presence of superconductivity in this model. However, due to numerical restrictions, realistic multi-orbital calculations are only possibly for a handful of sites, restricting the cluster extensions essentially to nearest neighbor correlations (for a review see Ref. [67]).
Diagrammatic extensions of DMFT on the other hand can treat short-and longrange correlations on an equal footing, which allowed, among others, the calculation of critical exponents [68,69,70,71] and revealed the absence of a metal-insulator transition in the two-dimensional Hubbard model on a square lattice [72]. In Section 4 we discuss the first of these diagrammatic extensions, the dynamical vertex approximation (DΓ A) [60,73] and its extension to ab initio calculations. For the latter, AbinitioDΓ A [74,16], we take as the vertex (irreducible in the particle-hole channel) the bare non-local Coulomb interaction as well as the local Coulomb interaction and all local vertex diagrams, see Fig. 8. From this unifying framework, we naturally generate all (local and non-local) GW diagrams, all local DMFT diagrams, as well as non-local diagrams beyond. The latter include, e.g., spin fluctuations which are important in the vicinity of phase transitions, for magnons and pseudogap physics. For a pedagogical introduction see Ref. [75], and Ref. [76] for an elaborate presentation. 2 Hedin's GW method: The new VASP implementation

Method
Hedin's method is in principle an exact approach to describe many-body interactions [11,77,78,18,79]. However, in practice for computational reasons, virtually all implementations of this method are limited to the so-called GW approximation. This greatly simplifies the calculations, but also makes important approximations 1 ; the considered Feynman diagrams are shown in the top panel of Fig. 1.
The new aspect of the present VASP implementation [15] is that it is tuned for massively parallel computers and that it works in imaginary time and frequency as opposed to the earlier VASP implementation that worked along the real frequency axis [14,80] and necessitated very fine frequency grids. In the following, we will give a brief outline of the computational steps of the present code, highlighting why it is particularly convenient for a combination with DMFT. We follow previous publications but emphasize simplicity and conciseness by dropping for instance the Brillouin zone index as well as the PAW formalism [15]. The first step in a GW calculation is to determine the DFT one-electron orbitals ψ i and one-electron energies i . From the DFT orbitals the one-electron Green function follows: Generalization to finite temperature is straightforward and involves restriction of the time to −β ≤ τ ≤ β, where β is the inverse temperature, and introduction of Fermi occupancy factors n i = 1/(exp(( i − µ)β) + 1) and (1 − n i ) in the first and second equation, respectively. It is then easy to show that the function observes the antiperiodicity for Fermionic Green functions G(r, r , iτ ) = −G(r, r , i(τ − β)).
As typically done in plane wave codes, all functions are expanded in a plane wave basis and fast Fourier transformed (FFT) to real space only when this is required: Here N r is the total number of real-space grid points. Since the plane-wave basis can be chosen to be significantly smaller than the number of real space grid points [81], the plane wave expansion typically reduces the storage demand by a factor 6-8 for orbitals (one position index), and a factor 6 2 − 8 2 for Green functions (two position indices). The first crucial approximation of the GW method is that the irreducible polarizability is approximated by the independent particle polarizability (RPA). Assuming a factor 2 for spin-degenerate systems, we get that is, vertex corrections of the form P = 2GGΓ are neglected. This approximation neglects important many body effects, for instance excitonic effects [80,10] that are captured by particle-hole ladder diagrams. It has been shown that these terms become important when selfconsistent calculations are performed [80]. From the irreducible polarizability the screened interaction (see Fig. 1 middle) can be determined by: Here, V is the Coulomb kernel, and integration over repeated spatial coordinates (s and s ) is assumed. For reasons of computational efficiency, the calculation is more conveniently done in reciprocal space, where the Coulomb kernel is diagonal [15]. The Dyson-like equation for the screened interaction needs to be solved in frequency space iω. This obviously requires one to perform a Fourier transformation of the independent particle polarizability from imaginary time [compare Eq. (5)] to imaginary frequency. In previous (imaginary time) GW codes [82,83] this was a fairly cumbersome operation involving fitting, a fast Fourier transformation, and some analytic continuation at very large frequencies and times. Using a mathematical rigorous treatment, Kaltak et al. determined imaginary time and frequency grids [84] that have a number of favorable properties. (i) The grids are non-uniformly spaced. This allows to simultaneously and accurately describe intra-band transitions at very small energies (meV), as well as high energy excitations into continuum like states (up to several 100 eV). (ii) The time and frequency grids are individually optimized to allow accurate calculations of the correlation energy in second order. Convergence of the correlation energy is exponential in the number of time or frequency points, with 20 points yielding µeV convergence even for metals. (iii) The grids are dual to each other: if a Bosonic function is known at a grid of N ω frequency points ω k , k = 1, ..., N ω , the numerical error in the Bosonic function is minimal at a set of corresponding N τ = N ω imaginary time points τ j , j = 1, ..., N τ . (iv) Related to point (iii), a numerical discrete Fourier transformation exists to transform any function from imaginary time to imaginary frequency (and vice versa): This is a numerical approximation to the Fourier transformation from time to frequency : The corresponding matrix of coefficients, e.g., γ kj cos(ω k τ j ), are precalculated and stored.
The evaluation of the self energy is most conveniently done in imaginary time with the bare Coulomb kernel V subtracted before the Fourier transformation of W and the contribution GV calculated analytically. As for the polarizability, also Eq. (9) neglects vertex corrections (Σ = −GW Γ ). For non-correlated semiconductors, the vertex contributions are only of the order of 0.2 eV for states close to the Fermi-level but can reach up to 1 eV for localized d-orbitals [85].
With the evaluation of Σ c , a single shot G 0 W 0 calculation is finished, so it is worthwhile to recapitulate what can be done with the yet calculated quantities. It is straightforward to express the self energy in any basis, for instance, a set of localized Wannier functions and to export it to a DMFT solver. The advantages over a conventional GW implementation are numerous. (i) First, many GW codes avoid calculating the full frequency dependency of the self energy, and instead evaluate the self energy only at a few points close to the DFT one-electron energies Σ( DFT i ). In the present code, this is no longer necessary and one obtains the self energy at all imaginary time points by Eq. (9). There is a (small) price to pay, though: to obtain physically measurable quantities, the self energy needs to be continued to the real axis, for which continued fractions are used [15]. However, since DMFT solvers usually work in imaginary time, the interface between VASP and DMFT is simple and requires only an interpolation from the few available imaginary frequency points {iω i } to a denser Matsubara grid. (ii) Each of the individual compute steps scales (at worst) cubic in the number of grid points or plane waves and linear in the number of k-points, as opposed to conventional GW codes, which scale quartic in the number of basis functions and quadratic in the number of k-points. The favorable scaling is straightforward to see: the calculation of the polarizability [Eq. (5)] and self energy [Eq. (9)] are clearly quadratic in the number of grid points, however, cubically scaling rank one updates of matrices and matrix multiplications are required to calculate the Green function [Eq. (1)] and the screened potential [Eq. (6)]. This favorable scaling combined with the efficient frequency grids allowed for efficient calculations of the random phase approximation (RPA) of the correlation energy for isolated defects in huge supercells containing several hundred atoms [86]. (iii) The constrained RPA (cRPA) [87] is simple and straightforward to implement in the present code. One only needs to remove the polarizability P t (iτ ) = 2G t (iτ )G t (−iτ ) of some target, say t 2g , orbitals from the total polarizability P = P r + P t to obtain an effective screened interaction U : The polarizability P r (iω) then captures all screening effects, except for the one inside the target space, which will be treated in the DMFT solver. The full frequencydependent U (iω) can be calculated with very little extra cost, and after transformation from the plane wave basis to a localized target space, it can be directly imported into a DMFT continuous time quantum Monte Carlo solver. The advantages of the imaginary time and imaginary frequency representation are more obvious if self-consistency is considered. Once the GW self energy is known, the Green function can be updated by where H HF = −∇ 2 /2 + V ion + V H + V x is the Hartree-Fock Hamiltonian consisting of the kinetic energy term, the ionic V ion , Hartree V H and exact exchange potential V x .
To obtain a converging Fourier transformation when transforming to the imaginary time, the Hartree-Fock Green function (second term) needs to be subtracted and added back in imaginary time This closes the cycle and allows to continue with a re-evaluation of the independent particle-hole polarizability in Eq. (5). Clearly, it is also possible to add any local self energy in Eq. (11) (terms in square brackets) and, thus, seamlessly incorporate DMFT results. Likewise, the irreducible polarization propagator can incorporate local effects beyond the independent particle-hole approximation, if the DMFT code provides the required information (P → P GW + P imp − P GW imp , compare next section). This opens the route towards a concise implementation of GW+DMFT, as discussed in the next section (see Fig. 5). A closure of the self-consistency cycle is already possible in the present code, although some intricacies for metallic systems still need to be solved, including an approximate inclusion of Drude-like metallic screening, and an efficient update of the chemical potential, which is important to achieve robust convergence in the self-consistency cycle for metals.

Results
The GW method has now been used for almost five decades. However, despite its undisputed improvements compared to DFT, results vary significantly between different codes. Errors are actually particularly large for transition metal compounds placing a serious question mark on any quantitative predictions. Specifically in oxides, d-binding energies can vary by up to 1 eV using different codes and implementations [88]. This is clearly unacceptable, if one aims to merge GW with more accurate methods such as DMFT. One major problem of the GW method is that the convergence with respect to the basis set size is extremely slow [89,90]. Specifically, for the projector augmented wave method, as used in VASP, the partial waves, which are supposed to form a sufficiently complete basis in the vicinity of the atoms, need to be chosen such that basis set convergence can be attained. The slow convergence has been rigorously discussed by Klimes, Kaltak and Kresse in Ref. [88]. Although that paper also establishes a suitable benchmark for solid state systems, we are not aware that other comparable reference numbers have yet been published for solids. Then, how can one ascertain that the numbers predicted with VASP are accurate and reproduce the infinite basis set limit? Fortunately, the new VASP GW code allows us to address this issue. Since it is efficient for large unit cells and large basis sets, it is possible to compare the results for molecules with atomic codes that use Gaussian type orbitals (GTOs). GTOs have been used for 50 years in quantum chemistry and have matured to a point where convergence for excited state calculations can be obtained fairly easily, although careful basis set extrapolations are as important as for plane waves. Fig. 2 shows the difference between the basis set extrapolated GTO results and the VASP PAW results for the ionization potential of 100 closed shell molecules. The mean deviation between both codes is only 60 meV [91], and large outliers are practically absent. We note that the deviations between other plane wave codes and GTOs are on average twice as large, but can even reach 200 meV on average. The other important point is the large difference between theory and experiment highlighting how limited the precision of G 0 W 0 is even for simple weakly correlated systems such as small molecules. This clearly underlines the need to go beyond the random phase approximation and single shot G 0 W 0 calculations. As an illustrative example for solid state calculations, we show results for SrVO 3 . Fig. 3 shows our calculated on-site dynamical screened intra-orbital interaction U (iω), inter-orbital U (iω), and Hund's coupling J(iω) of SrVO 3 using the cRPA and V-t 2glike maximally localized Wannier functions. In imaginary frequency, U , U and J are rather smooth functions, so that it is possible to interpolate them from the optimized frequency grid to Matsubara frequencies by a Padé interpolation [92] (see the red solid lines in Fig. 3). This makes it possible to transfer them to a dynamical impurity solver. In the static limit (ω = 0), U , U , and J are calculated to be 3.38, 2.42, and 0.44 eV, respectively, agreeing perfectly with the ones directly obtained from the conventional implementation working on the real frequency axis. Moreover, they are in nice agreement with the published values [93,94]. In the high-frequency limit (ω → ∞), U , U , and J approach the unscreened (bare) counterparts (16.29, 15.07, and 0.55 eV). Figure 4(a) shows the single-shot G 0 W 0 momentum resolved spectral function. Compared to DFT, the t 2g bandwidth is reduced by 20 % in G 0 W 0 . In the G 0 W 0 approximation, spectral weight is transferred to satellites. This is much more clearly seen in the local, momentum-integrated, spectral function as shown in Fig 4(b). The plasmon satellite of the t 2g quasi-particle band at ∼3 eV arising from the t 2g contribution to the fully screened interaction at the plasmon frequency [51,95] is well reproduced. Further, a plasmon peak deriving from transitions outside the t 2g subspace is seen at ∼15 eV [93,51]. We note, however, that in our calculations repeated plasmon peaks at higher frequencies are absent. This is a well known issue of the G 0 W 0 approximation [96].
3 Screened exchange and local quantum fluctuations: GW+DMFT

Method
The key advantage of the GW approach discussed above is its treatment of dynamical screening: While standard electronic structure methodologies-such as Hartree-Fock or DFT-work with the bare Coulomb interaction V , GW explicitly incorporates the polarizability of the electronic system. As a consequence, the repulsion between electrons becomes reduced and retarded. The resulting screened-exchange self energy yields a much improved description of, e.g., sp-semiconductors gaps, whereas the retardation effects account for spectral weight transfers to (plasmon) satellite features, and finite lifetimes of electronic excitations.
However, as already mentioned in the Introduction, the perturbative GW approach is insufficient for strongly correlated materials. In fact, it fails to account for their strong mass renormalizations, Hubbard satellites, and local moments physics [21]. Our recent understanding of strong electron-electron correlations was indeed propelled by the advent of a non-perturbative technique: the dynamical mean field theory [24]. The latter maps the lattice problem onto the self-consistent solution of an Anderson impurity model, and the lattice self energy is identified with the single-site (i.e., local) self energy of the impurity [23]. This mapping becomes exact in the limit of infinite lattice coordination [22]. By construction, DMFT accounts only for correlations from on-site interactions, yet it includes-as depicted in Fig. 1-all Feynman diagrams built from the Hubbard U and Hund J interactions and the local impurity propagator. To set up a realistic DMFT calculation, the one-particle part of the Hamiltonian is taken from DFT (whence the name DFT+DMFT [25,26]) and the screened interaction parameters-U and J-can be computed from techniques such as constrained DFT [97], or, better, the constrained random phase approximation [87,98,99]. However, as mentioned in the Introduction, it is not separable how the Hubbard U already contributes to the DFT band-structure. so that there is the problem of "doublecounting" correlations when adding the DMFT self energy.
From this brief summary it is apparent that GW and DMFT are very complementary techniques: GW has no restriction on the range of the interaction or the self energy and therefore excels for sp-systems. In DMFT the interaction and the self energy are by necessity localized on an atomic site, yet their non-perturbativeness allows for a reliable description of the Kondo and Mott physics realized in many dor f -electron materials. At the same time, GW and DMFT share a common (diagrammatic) language. Therewith, both methods can profit from each other: As the RPA technique is integral part of the GW, it can provide the DMFT with a Hubbard U computed from first principles [see the preceding section and Eq. (10)]. In return, DMFT susceptibilities and self energies can add local vertex corrections to all orders to Hedin's equations for the polarization and self energy, see Eqs. (5) and (9), respectively. Contrary to DFT+DMFT, any double-counting in this combination of screening and correlations can be avoided, since a clear-cut separation is possible on the diagrammatic level.
This outlines the GW+DMFT method proposed in Ref. [50]. By elegantly combining the best of both worlds-screened exchange and local quantum fluctuations-GW+DMFT has the potential to vastly extend the realm of quantitative and predictive many-body electronic structure theory. Let us give specific examples: In many materials the separation between correlated d or f -states and ligand sp-orbitals is often severely underestimated within DFT and DFT+DMFT [100,101,102,46]. Yet, optical transitions between these states can actually be relevant for technological applications in, e.g., intelligent window coatings [103], or eco-friendly rare-earth-based pigments [104]. Calculating the red colour of CeSF indeed required incorporating a GW correction into DFT+DMFT [104]. Non-local (inter-site) self energiesà la GW were also shown to be crucial in oxides [105,102], intermetallics [106], and ironpnictides and chalcogenides [107], and in particular for explaining the non-magnetic nature of BaCo 2 As 2 [108]. Moreover, important effects of dynamical screening were found, among others, in oxides [109,51,110], pnictides [111,108] and cuprates [112]. On the other side, local vertex corrections in susceptibilities beyond RPA where shown to be crucial in, both, Hubbard models [113,114] and realistic materials, e.g., regarding the dynamical structure factor in iron-pnictides [115,116], and the absence of ferromagnetism in stoichiometric FeAl [117].
After this general rationale, we will now discuss GW+DMFT in some more detail (see Refs. [54,53,55,56] for longer reviews). The workflow of the approach is depicted in Fig. 5: On the left-in blue-is the DMFT [118,119,49] cycle with an additional self consistency for the two-particle interaction: It is required that the local screened interaction W loc equals the screened impurity interaction W imp = U + U P imp W imp = U − U χU , where χ = T n(τ )n(0) is the impurity densitydensity correlation function and U the local interaction (containing, e.g., Hubbard and Hund terms). Owing to the dynamical nature of screening, these interactions are in particular frequency-dependent, i.e., U → U (ω). At least for density-density type of terms, solving an Anderson impurity model with dynamical interactions is Fig. 5: The GW+DMFT approach. The DMFT sub-cycle is indicated in blue, the GW procedure in green, and shared quantities in grey boxes.
easily possible with quantum Monte Carlo techniques, both approximately [109] and numerically exactly [109,120,121]. On the right-in green-are Hedin's equations for the polarization P GW and the self energy Σ GW in the GW approximation. Neglecting vertex-corrections here boils it down to the RPA for P GW , Eq. (5), and the first order GW expression for Σ, Eq. (9). DMFT and GW intersect at two junctures-marked in grey-once on the twoparticle level in the polarization, and once on the one-particle/spectral level in the self energy. Both times, non-perturbative, yet local contributions from DMFT, P imp and Σ imp , are added to the GW contributions, P GW and Σ GW . Since the latter already contain some of the local diagrams of the former, these terms have to be subtracted. Indeed, at each iteration, we need to remove all contributions from the polarization and the self energy that arise when computing the impurity analogues of P and Σ on the GW level. In case of the polarization, this is achieved by subtracting a local RPA polarization P GW imp = 2G imp G imp obtained from a convolution of two impurity Green functions [122]: P = P imp + P GW − P GW imp . For the self energy, we need to subtract a term Σ GW imp that is computed as the first order contribution in an interaction W GW imp that derives from screening the impurity interaction-the Hubbard U -with the above polarization P GW imp , i.e., W GW imp = U −1 − P GW imp −1 . 2 Due to these junctures, there is an outer self-consistency that allows for a feedback of local and non-local many-body effects onto the GW and DMFT cycles, respectively. Typically such a calculation is initialized with a Green function Green function, and a guess for the self energy Σ (here including the Fock term). In the first iteration, Σ is usually replaced by the DFT exchange-correlation potential V xc , i.e., G = G DF T (called G 0 in the Introduction section).

Results
Combining two methods that have evolved and matured independently over decades into large software packages is an intricate endeavor. Therefore, the full scheme, as shown in Fig. 5, has been realized first for one-band calculations [49,113,114,123]. For realistic multi-band systems, the first implementation-Tomczak et al. [51]-resorts to simplifications, namely (1) omitting global self-consistency, i.e., performing only one-shot GW calculations starting with G = G DF T , (2) fixing the double-counting polarization P GW imp to the (dynamical!) cRPA result and approximating P imp ≈ P GW imp , (3) approximating the double-counting self energy by the local projection of the GW self energy: Σ GW imp ≈ k Σ GW , and (4) solving the DMFT impurity with dynamical U (ω) within the approximative Bose factor Ansatz [109]. In other early works, additional approximations were made: Taranto et al. used a static Hubbard U and circumvented computing a fully frequency-dependent Σ GW [124] (see also below), and Sakuma et al. combined DMFT and GW self energies from independent calculations [125].
Applied to the prototypical correlated metal SrVO 3 , GW+DMFT revealed important new insights [51,52]: The additional ingredients-the momentum-dependent self energy Σ(k, ω) and retardation effects in the Hubbard interaction U (ω)-are found to compete. The dynamics in the interaction describes, among others, spectral weight transfers to plasmon satellites (at ∼ 15eV for SrVO 3 [93]). These high energy excitations account for an additional (i.e., beyond Hubbard-model physics) reduction of the low-energy quasi-particle weight and, correspondingly, to a narrowing of the bandwidth [109,110] (by a factor Z B ∼ 0.7 in SrVO 3 [126]). Therefore, DMFT calculations that use only static interactions, have to employ a larger Hubbard U = 4 − 5.5eV [127,29,100,101] than the static limit U (ω = 0) ≈ 3.5eV [93,98,51] of the cRPA to account for the same mass enhancement; such a larger interaction is actually obtained in constrained LDA [128]. The non-local exchange self energy on the other hand widens the low-energy dispersion [107,51,129,130]. Correspondingly, effective masses of quasi-particles are reduced. With respect to the LDA reference, effective masses are given by the ratio of the LDA and GW+DMFT group velocities: Here, the denominator is related to the quasi-particle weight Z k = [1−∂ ω ReΣ(k, ω)] −1 ω=0 . In DMFT, where the self energy is local, m * /m LDA = 1/Z holds. In GW+DMFT, the extra term involving the momentum derivative of the self energy substantially counteracts the mass enhancement generated by the dynamical correlations [107,52]. Altogether this yields a similar effective mass as in the previous DFT+DMFT calculations (that use U > U (ω = 0)), but the low-energy spectral weight is different, and can be measured by transport or optics.
In Fig. 6 we compare Matsubara self energies for the t 2g orbitals of SrVO 3 obtained with our new implementation that combines the GW-code of VASP detailed in Section 1 with the w2dynamics DMFT code [131,132]. 3 We employ the same approximations (1)-(3) as in Ref. [51]. However, instead of (4) the approximative Bose-factor Ansatz [109], we use a numerical exact continuous-time quantum Monte Carlo algorithm for retarded density-density interactions [120]. Our reference is a standard DFT+DMFT calculation that uses a static Hubbard U (ω = 0) and Hund's J as provided by the cRPA (see Fig. 3). From the low-energy slope of ImΣ(iω n ) we extract a quasi-particle weight Z = 0.6. Turning on the retardation in the interaction, i.e. solving DFT+DMFT with the dynamical cRPA U (ω) adds substantial renormalizations of plasmonic origin; Z decreases to 0.3. Moreover, since the dynamical interaction recovers at high frequencies the unscreened Coulomb interaction, ReU (ω → ∞) ≈ 16eV, also the self energy Σ lives on a much larger energy scale [109] than in the standard, static DFT+DMFT case. Adding the non-local GW self energy decreases effective masses, i.e. the ratio of U over bandwidth diminishes and so does the strength of correlations: In our GW+DMFT the local quasi-particle weight is Z = 0.63, which is even slightly larger than within static DFT+DMFT. We find qualitative agreement with previous GW+DMFT results from Ref. [52]. 4 GW+DMFT spectra for the t 2g -orbitals of SrVO 3 from Ref. [51] are shown in Fig. 7 in comparison with angle-resolved photoemission spectroscopy (ARPES) results. The calculation agrees well with the experimental data. Differences to previous DFT+DMFT calculations (see, e.g., Refs. [127,29,100,101]) are however most pronounced for unoccupied states [52], that are inaccessible to ARPES experiments. The effects are in line with the above discussion: (i) W.r.t. DFT+DMFT the low-energy bandwidth is enhanced by non-local self energy contributions (e.g., the unoccupied d xy , d xz -bands at the X point move up from 0.6eV [29] to ∼ 1eV). (ii) Using an ab initio screened interaction U (ω) (instead of a larger static U adjusted so as to reproduce the experimental mass enhancement), the upper Hubbard band is placed at much lower energy (e.g., 1.2(1.9)eV instead of 2.2(2.85)eV [29] at the Γ (X) point for the d xy , d xz -components). Indeed, the upper Hubbard band merges with the quasiparticle peak in momentum-integrated spectra. This reduced importance of Hubbard physics has recently been confirmed by partially self-consistent GW+DMFT calculations [95], and is compatible with recent inverse ARPES experiment. 5 Besides the shown low-energy dispersion, also the position of ligand states, in particular the O-2p and Sr-4d improve substantially in GW+DMFT. Since the discussed GW+DMFT results are not globally self-consistent, the ligand states are at the same position as in G 0 W 0 calculations [51,124,52,133]. Fig. 7: Comparison of GW+DMFT spectra from Ref. [51] (bottom) with angle resolved experimental photoemission spectra (top) from Refs. [134,135] for the Γ and X-point.
In the setup used for these results (one-shot G 0 W 0 ), the only modification of the DMFT code comes from adding the frequency-and momentum dependent GW self energy contributions into the DMFT one-particle self-consistency. Yet, since the GW self energy is a large and unhandy object, it stands to reason to approximate it to alleviate memory consumption. Also, as mentioned in the preceding section, some GW codes cannot provide the self energy on a continuous frequency mesh. A strategy to simplify the influence of the GW self energy was pioneered by Taranto et al. in Ref. [124]: By evaluating the self energy at the Kohn-Sham energies and hermitianizing it along the lines of Ref. [12], it is possible (using an approximate double-counting correction) to cast Σ GW into Hamiltonian form. This idea was further developed into the more rigorous quasi-particle self-consistent (QS)GW+DMFT approach proposed by Tomczak in Ref. [136]-previously alluded to in Ref. [107]-in which the double-counting correction can be performed exactly, and furthermore, a global selfconsistency is performed on the QSGW [12] level. Flavors of QSGW+DMFT have subsequently been applied to cuprates and nickel oxide [137], nickel and iron [138], as well as insightful model systems [139,140]. Recently, Boehnke et al. [95] have pioneered a setup in which GW+DMFT self-consistency is performed beyond the quasi-particle approximation, yet only within a low-energy subspace, in this case the t 2g -orbitals of SrVO 3 . As anticipated in earlier model calculations for the extended Hubbard model [113,114], the self-consistent local interactions are smaller in this setup than the initial cRPA values: compared to previous non-self-consistent works [51,52] the strength of correlations is reduced and mass renormalizations in SrVO 3 become more plasmonic in origin [95].
This concludes the description of the current state-of-the-art in GW+DMFT calculations. While there is a panoply of materials to which the described methodology can be applied with great benefit, let us point out two challenges for future developments: (i) Besides influencing the low-energy dispersion, the GW self energy also effects higher lying states, e.g., the O-2p and the Sr-4d orbitals in SrVO 3 . Indeed, the O-2p orbitals are off by 1.5 − 2eV within DFT, and are pushed towards their experimental position by the GW self energy [51,52,133]. This will reduce their contribution to screening, causing an increase in the Hubbard U for the t 2g -orbitals, as indeed found when performing cRPA on top of QSGW. This effect of ligand states will counteract the reduction of correlations seen in calculations in which selfconsistency is limited to the t 2g -orbitals [95]. Hence, a GW+DMFT implementation that includes ligand states in the self-consistency is eagerly awaited. (ii) Contrary to one-shot G 0 W 0 , self-consistent GW+DMFT is a conserving theory on the oneparticle level: the GW+DMFT self energy is derivable from an approximation to a free-energy functional [141,50]. Yet, on the two-particle level the situation is inverted: RPA and also QSGW+DMFT [136] yield a conserving density-density response function/polarization. In fully self-consistent GW+DMFT which uses dressed Green functions, on the other hand, gauge invariance for the polarization is not given and its violation can lead to a qualitatively wrong description of, e.g., collective (plasmon) modes [142]. Achieving gauge-invariance respecting one-and two-particle quantities within a GW+DMFT scheme remains a challenge for future works.

Correlations on all time-and length-scales: Ab initio DΓ A 4.1 Method
The essential approximation of DMFT is that the self energy Σ, which is nothing but the one-particle fully irreducible vertex, is local-given by all local skeleton diagrams [24]. We can put this concept on the next level, assuming the locality of the two-particle fully irreducible vertex Λ. 6 This is the dynamical vertex approximation (DΓ A) [60,73]. From the local, fully irreducible vertex Λ, extracted by inverting the local parquet equation of DMFT [143], one can construct-through the self-consistent solution of the parquet equation of the lattice system-the non-local full vertex F , and from that, the non-local DΓ A self energy as well as all physical susceptibilities. This full-fledged parquet DΓ A approach has been employed in Refs. [144,145]. In most calculations however, a restriction to the particle-hole (and transversal particle-hole channel) has been employed. In this so-called ladder DΓ A [60,73], the local vertices irreducible in the two particle-hole channels Γ ph are the starting point and the full vertex F is constructed through the Bethe-Salpeter ladder. This neglects the particleparticle channel which is important, e.g., for superconductivity and weak localization corrections to the conductivity. In both variants the local and non-local self energy is obtained through the Schwinger-Dyson equation of motion, and includes non-local correlation effects such as spin fluctuations and pseudogap physics.
A variety of closely related approaches have been subsequently proposed [62,63,64,65,66]. They all have in common that they include all the local DMFT correlations, and construct additional non-local correlations from the two-particle vertex via Feynman diagrams. The differences are in the details: (i) which two-particle vertex is taken, (ii) whether the real or a dual Green function (subtracting the local Green function) is taken as connection line, (iii) which Feynman diagrams are considered. These diagrammatic extensions of DMFT have been highly successful for studying model systems such as the one-band Hubbard model and we discuss selected results in Section 4.2.
For realistic materials calculations, one might envisage using DΓ A instead of DMFT in a DFT+DΓ A scheme. However, it is more appealing to use the Bethe-Salpeter equation also as a means for calculating the non-local exchange and correlation. This is possible by taking, as the irreducible vertex in the particle-hole channel, the non-local Coulomb interaction V q in addition to the local vertex, see Fig. 8 (b). Besides non-local interactions, such a treatment also allows us to include less strongly correlated orbitals-without the need to calculate the local vertex for them. In the following we discuss this AbinitioDΓ A, while results for SrVO 3 are presented in Section 4.3 The flow diagram of the AbinitioDΓ A algorithm is given in Fig. 8, for a complete presentation and technical details see Ref. [16]: Fig. 8 (a): The first step is to calculate the local generalized susceptibility χ loc via the numerical solution of an Anderson impurity model, 7 and use the local variant of the Bethe-Salpeter equation as well as the bare (bubble) susceptibility χ 0,loc to extract the local irreducible vertex in the particle-hole channel Γ loc (as indicated). This vertex and the local susceptibility depends on three frequencies ν, ν and ω, and four orbitals m, l, m , l . Fig. 8 (b): We supplement this local irreducible vertex Γ loc with the non-local Coulomb interaction V q at momentum q. Together these terms form the AbinitioDΓ A approximation for the irreducible vertex Γ q in the particle-hole channel. Fig. 8 (c): With this Γ q we solve the Bethe-Salpeter equation on the lattice to get the full vertex F q . 8 7 Calculating this vertex by continuous time quantum Monte Carlo simulations [146] is computationally the most demanding step. For getting all components of the vertex, a worm sampling is needed [147]; using an improved estimator [148,149] and vertex asymptotics [150,145,151,149] increase the accuracy and size of the frequency box. 8 As detailed in Ref. [16], besides the displayed particle-hole ladder, also the transversal particle-hole ladder is taken into account, and the double-counted contribution is subtracted. The Bethe-Salpeter equation is formulated in terms of a magnetic and density combination Fig. 8 (d): This F q allows us in turn to calculate the AbinitioDΓ A self energy via the Schwinger-Dyson equation of motion (the second line represents the Hartree and Fock contribution to the self energy).
Self consistency: With a new self energy and local Green function we can, in principle, go back to Fig. 8 (a) and recalculate the local susceptibility and vertex, closing the self-consistency loop.
Before turning to our presentation of selected DΓ A results, let us briefly discuss what kind of physics the AbinitioDΓ A can describe. First of all, we notice that the first, V q term of Fig. 8 (b) yields the RPA screening when inserted into the Bethe-Salpeter equation in the particle hole of Fig. 8 (c). Via Fig. 8 (d) this yields the GW self energy. That is, all GW diagrams are included in AbinitioDΓ A. But on top of GW, there is also the crossing symmetrically related ladder in the transversal particle-hole channel. Second, if we only consider the local Green functions in the Bethe-Salpeter equation Fig. 8 (c), we recover the local F or susceptibility χ of Fig. 8 (a) as well as, via the equation of motion, i.e., the DMFT self energy. In other words, all DMFT diagrammatic contributions are also included. Beyond both, we have more diagrams and physics included, e.g., spin fluctuations. These can be described in weak coupling perturbation theory as the particle-hole and transversal particle-hole ladder with the local bare interaction U as a building block. These diagrams are generated in Fig. 8 (c) when taking the bare U term which is part of Γ loc in Fig. 8 (b). More precisely it is the U contribution to the vertex which is analogous to the second, V k−k term on the right hand side of Fig. 8 (c) which generates the spin fluctuations. Let us emphasize that in DΓ A such spin fluctuations are not restricted to weak coupling, the same kind of diagrams are also generated with the full Γ loc .

Results I: One-band Hubbard model
In the last decade the DΓ A has been intensively applied to study several physical aspects of the single-orbital Hubbard model. On the one hand, this was important to demonstrate the performance of DΓ A-based algorithms to describe intermediateto-strong-coupling parameter regions, hardly accessible to other techniques, in view of subsequent applications to realistic systems. On the other hand, the DΓ A, since its first applications to single-orbital models, has allowed significant progress in the fundamental understanding of important topics in many-body physics. We just mention here, among others, in d = 2 the transformation of the Mott metal-insulator transition into a crossover down to U = 0 and the spin-fluctuation-driven pseudogap [73,72,152,153], and, in d = 3, the critical exponents of the Hubbard model and the breakdown of the paramagnetic Fermi-liquid at low temperatures (T ) because of spin fluctuations [68,153]. Notably, several of these DΓ A findings have been supported [154,69] by complementary results of other powerful diagrammatic-extensions of DMFT, such as the dual-fermion [62] and dual-boson approaches [155], as well as other novel many-body techniques (e.g., the fluctuation diagnostics [156]).
In this section, we will review some of the most recent DΓ A applications to the single-orbital Hubbard model in d = 3, and discuss their possible implications for the development of high-performing multiorbital algorithms. The first DΓ A application we consider is the investigation [71] of the quantum critical properties of the magnetic transition in the three dimensional Hubbard model, as a function of (hole-)doping. As of spins which are not displayed in Fig. 8. Neglecting the second, V k−k term in Fig. 8 (b) simplifies the momentum dependence (Γ qkk → Γ q ) and dramatically reduces the computational effort to solve the Bethe-Salpeter equation. Please note that a corresponding local contribution U is included as part of Γ loc but does not lead to a k, k -dependence. it was also found in DMFT [157], the relatively high antiferromagnetic (AF) ordering temperature (T N ) of the half-filled system is progressively reduced by increasing doping, until at about 20%-doping a quantum critical point (QCP) is found (see Fig. 9). The reduction of T N is associated also to a gradual transformation of the magnetic order from commensurate AF at (π, π, π) to an incommensurate spin-density wave (SDW) at (π, π, Q z < π), see inset of Fig. 9. While the DMFT-description of the related (quantum) critical properties is restricted to mean-field correlations in space, the DΓ A-treatment of both space and temporal correlations on an equal footing yields an improved understanding of the magnetic QCPs in d = 3. In particular, beyond a sizable reduction of T N w.r.t. DMFT throughout the phase-diagram, the finite-T critical exponents γ, ν found in DΓ A for the magnetic susceptibility and correlation length, respectively, are consistent with the 3d-Heisenberg universality class (i.e., γ 1.4, ν 0.7), independently on whether the antiferromagnetism is commensurate or incommensurate. While this already corrects the mere mean-field values of critical exponents found in DMFT (γ = 1, ν = 0.5), the nature of the criticality changes further at the QCP. Here around n ∼ 0.8, the exponents take unexpectedly the values γ 0.7 ÷ 0.8, ν = 1, strongly violating the typical scaling relation γ = 2ν. These values are also incompatible with the standard Hertz-Millis-Moriya theory [158,159] for perturbative QCPs. By means of a complementary semi-analytical analysis [160], the unusual values of the critical exponents have been ascribed to the presence of lines of Kohn's points in the underlying Fermi surface, whose effect is no longer damped by finite-temperature fluctuations at the QCP. The DΓ A results have, thus, identified an additional, important factor controlling the quantum critical properties of correlated systems, hitherto mostly neglected.
The effects of non-local fluctuations are, obviously, not confined to the (quantum) critical properties, as they also affect the spectral properties, in particular at the Fermi-energy. However, while the electronic self energy is significantly corrected w.r.t. the DMFT results, especially at low-T [68,153], a closer inspection [161] reveals that, in d = 3, the intrinsic frequency/momentum structure of the electronic self energy in DΓ A displays specific, important patterns. These, in turn, can be used for devising important simplifications of realistic many-body algorithms for bulk systems, such as GW, GW+DMFT, or the AbinitioDΓ A. Specifically, an inspection of the DΓ A self energy, continued to real frequencies (see Fig. 10) shows that the self energy of the 3d Hubbard model, even in the most correlated low-doping regime (n = 0.9), displays a clear separation in the time/frequency and space/momentum domains: The hallmark of such separation, which extends to a relatively broad frequency interval around the Fermi level, is immediately visible in Fig. 10: in form of the parallel frequency behavior of ReΣ(k, ω) for different k, which reflects a momentum-independent quasi-particle renormalization factor Z k ∼ Z. The DΓ A demonstrates, in fact, that the momentum dependence of Σ(k, ω) is essentially confined to the static sector, which explains the shift among the different parallel self energies. Not surprisingly, the same qualitative behavior, though-quantitatively-less correlated (i.e., with a larger Z) is found in the corresponding GW results, shown for comparison in Fig. 10. It is however noteworthy that the static momentum-dependence is much larger in DΓ A than in GW. This advocates the presence of true non-local correlation effects as opposed to exchange effects that cause a large (static) k-dependence in multi-band GW calculations [129,52]. In fact, lacking spin-fluctuations (both local and non-local), GW-where space-time separation was first evidenced [107,52]-verifies Eq. (14) up to larger frequencies than DΓ A. The validity of Eq. (14) for strongly correlated systems in d = 3 is inspiring for promising algorithmic improvements, potentially applicable to several many-body techniques (e.g., GW, GW+DMFT, AbinitioDΓ A). In all these cases, the assumption of a full time-space separability of Σ(k, ω) would allow to avoid numerically ex- Fig. 11: Momentum-dependence in the k z =0-plane of (a) the real and (b) imaginary part of the electronic self energy at the lowest Matsubara frequency (iν 0 ), computed in AbinitioDΓ A for the V-t 2g orbital 3d xy of SrVO 3 . The corresponding quasi-particle parameters (weight Z k ) and scattering rate γ k , are reported in panels (c) and (d), respectively (reproduced from Ref. [16]).
pensive transfer-momenta/frequency convolutions in the intermediate steps of manybody/diagrammatic calculations, reducing considerably the numerical effort. Further details about such simplifications, and an explicit proposal of a "space-time separated GW" scheme are reported in Ref. [161]. We should also notice, at the end of this section, that complementary simplifications of the momentum structure, were suggested by the recent findings of Ref. [162]: In d = 2 the momentum dependence of Σ(k, ω) can often be approximated by a dependence on the non-interacting dispersion, i.e., Σ(k, ω) → Σ( k , ω).

Results II: SrVO 3
Being a diagrammatic extension of DMFT, the algorithmic implementation of DΓ A is not affected by cluster-size limitations of cluster extensions of DMFT, making possible a systematic generalization of the DΓ A approach to treat realistic multiorbital systems. While the technical aspects of the AbinitioDΓ A [74,16] have been addressed in the previous Section, here we will discuss the physics emerging from the first applications of the AbinitioDΓ A to realistic material calculations. Specifically, we will focus on the very recent AbinitioDΓ A study of Galler et al. [16] performed for the correlated-metal testbed material SrVO 3 . In this compound the 3d−t 2g bands of V are rather well separated from the other bands, which allows for a relatively accurate modelization already in terms of a three-orbital t 2g -only manifold. In fact, this modelization has been exploited in the past for a huge number of many-body calculation, from LDA+DMFT to GW+DMFT (see Section 3), and represents, thus, a sort of drosophila among correlated systems.
The application of the AbinitioDΓ A to SrVO 3 has allowed one of the first nonperturbative analyses of the momentum-dependence of the self energy, and of the spectral function in this compound. The AbinitioDΓ A calculation foots on a DMFT solution of a realistic three-band Hubbard model for the t 2g -orbitals of vanadium. The latter uses a static Hubbard U = 5.0eV and Hund's J = 0.75eV in the rotationallyinvariant Kanamori parametrization. The computation of two-particle quantities resorts to a recent worm algorithm [147]. A sample of the AbinitioDΓ A results, for the d xy -orbital, is shown in Fig. 11, where ReΣ(k, iν) and ImΣ(k, iν) at the lowest Matsubara-frequency (iν = iν 0 = πT ) are reported, as well as the related quasiparticle parameters (Z k ) and (γ k ) extracted from a low-frequency expansion of ImΣ(k, iν n ). These AbinitioDΓ A results show that that a sizable momentum dependence does appear in the electronic self energy of SrVO 3 even if-as in this casenon-local interactions are neglected. These effects thus correct the purely local DMFT results. Interestingly, this k-dependence is mostly confined to ReΣ(k, iν 0 ) (Fig. 11a), where one observes a momentum differentiation larger than 0.2eV, which, however, does not directly mirror the shape of the underlying Fermi-surface. At the same time, the overall k-dependence of Im Σ(k, iν 0 ), and the related quasi-particle coefficient Z k , γ k (Fig. 11b-d) is definitely much weaker (e.g., Z k varies less than 2% over the whole Brillouin zone, with an average value slightly increased w.r.t. DMFT). We note that this behavior matches rather well the conclusions of the space-time separation emerging from previous single band DΓ A calculations [161] discussed above. Moreover, going beyond the single-orbital framework, it is also worth emphasizing that a correlation between the momentum and orbital-dependence is found: the strength of the k-dependence, as computed in Ref. [16], displays the same trends for the orbital dependence of Σ. The latter was found, indeed, to be much more pronounced for ReΣ(k, iν 0 ), than for the other quantities shown in Fig. 11.

Conclusions
One of the main challenges in computational materials science is to predict the properties of materials for which the standard DFT-based methods are not applicable. Present DFT functionals are not reliable if the screening of the electron-electron interaction over different length-and time-scales is insufficient to approximate the electronic exchange and correlation with the common LDA and GGA functionals. This is the case for some correlated semiconductors, most transition metal oxides, as well as heavy fermions. In this paper, we have reviewed the forefront of methodological progress towards a full ab initio treatment of electronic exchange and correlations beyond DFT, and its state-of-the art merger with DMFT. These progresses range from the inclusion of the frequency dependencies in Hedin's GW-scheme for realistic calculations of large systems, to the implementation of a corresponding GW+DMFT algorithm that lifts the quasi-particle approximation in the scheme suggested by van Schilfgaarde and Kotani, and, eventually, to the treatment of correlations beyond the purely local description of DMFT by means of the AbinitioDΓ A. Such advances are of extreme importance, because the new algorithms are conceived to be able to treat all classes of materials, independently of the strength and the range of the screened electronic interaction in the specific compound. While such ambitious goals will certainly require further work in the next decade, the examples we have selected in this work to illustrate the applications of the different method developments already show a promising trend.
Specifically, we have started by discussing the outline of the new GW implementation in the Vienna ab initio simulation Package (VASP), written for massively parallel computations and working-for the first time within the VASP packageon the imaginary time/frequency axis. This allows not only the calculation of full dynamical information at the GW level, but also opens the road for the implementation of a self-consistency at the GW level-beyond the quasi-particle approximation by van Schilfgaarde and Kotani. The applicability of the new implementation has been demonstrated with a calculation of the testbed correlated metal SrVO 3 . The progress in the GW part are also pivotal for allowing a more natural and precise interfacing with DMFT-based algorithms. In particular, after reviewing the generic scheme of the GW+DMFT, where the non-local, but perturbative GW-exchange and correlations are supplemented with the purely local, but non-perturbative ones of DMFT, we have shown self energy results obtained with the G 0 W 0 +DMFT merger of the VASP and the w2dynamics codes, again for the prototype material SrVO 3 . In particular, numerical results obtained by means of different levels of refinement (e.g., retaining/neglecting the dynamical structure of the screened interaction in the DMFT part) have been presented and critically analyzed.
Finally, we have reviewed the most general algorithmic framework in which even the limitations of GW+DMFT can be overcome, i.e., the AbinitioDΓ A approach. This diagrammatic scheme starts with a local irreducible vertex (Γ irr ) obtained by using an impurity-solver (such as w2dynamics) and supplements it with the bare non-local Coulomb interaction. From this starting vertex, ladder diagrams (or for a few orbitals parquet diagrams) are constructed, yielding non-local self energies and correlation functions. While one can easily obtain-within the AbinitioDΓ A formalism-all previously discussed approaches (GW, DMFT, and GW+DMFT) as limiting cases, it also includes non-local correlations beyond these schemes, such as spin fluctuations. For a simple, one-band Hubbard model, we have recapitulated the unexpected properties found at the magnetic QCP in three dimensions, and discussed the space-time separability of the self energy. For realistic multi-orbital materials calculations, we have shown the very first AbinitioDΓ A results for SrVO 3 , and discussed the corrections found w.r.t. DMFT. Indeed, we evidenced a sizable momentum-dependence in the SrVO 3 self-energy even for purely local interactions-an effect well beyond GW approaches [161]. Instead, the momentum dependence in our GW+DMFT results is almost exclusively propelled by exchange contributions originating from non-local interactions-thus far omitted in our AbinitioDΓ A calculations. As a consequence, the results of AbinitioDΓ A and GW+DMFT cannot yet be directly compared. Calculations that include the non-local interaction in AbinitioDΓ A are under way.
As it is typical in physics, and especially true in the case of new algorithmic developments, a significant amount of future work will be inspired by the progress we have reviewed in this paper. In particular, the new method enhancements pave the way towards the implementation of a fully self-consistent, frequency-dependent GW scheme in VASP, while the frequency-dependent treatment of both the GW self energy and the local dynamic interaction of DMFT represents an important step towards the realization of a globally self-consistent GW+DMFT merger of the VASP and the w2dynamic codes. Eventually, the first successful applications of AbinitioDΓ A for treating strong non-local correlations beyond GW+DMFT will encourage further efforts towards a new standard of ab initio materials science calculations for correlated electron systems.