Vertex corrections for positive-definite spectral functions of simple metals

We present a systematic study of vertex corrections in the homogeneous electron gas at metallic densities. The vertex diagrams are built using a recently proposed positive-definite diagrammatic expansion for the spectral function. The vertex function not only provides corrections to the well known plasmon and particle-hole scatterings, but also gives rise to new physical processes such as generation of two plasmon excitations or the decay of the one-particle state into a two-particles-one-hole state. By an efficient Monte Carlo momentum integration we are able to show that the additional scattering channels are responsible for the bandwidth reduction observed in photoemission experiments on bulk sodium, appearance of the secondary plasmon satellite below the Fermi level, and a substantial redistribution of spectral weights. The feasibility of the approach for first-principles band-structure calculations is also discussed.

Despite these surpluses we still have a poor knowledge of the energy-and momentum-resolved spectral function A(k, ω) away from the on-shell manifold, i. e., when ω ε k . In angular resolved photoemission this is the regime where electrons with reduced energy (as compared to the prediction based on band-structure and energy balance) are observed. Selfconsistent (sc) perturbation theory, e.g., sc-GW [17], accurately predicts total energies [18][19][20] and it is fully conserving at the one-particle level, a crucial property in the description of transport phenomena [21]. However, for spectral properties sc schemes do not show the expected improvement over simpler one-shot calculations [22][23][24]. In fact, they suffer from serious drawbacks: the incoherent background in the spectral function gains weight at the expenses of the qp peak [25], the qp energy does not agree with experiments (overestimating the bandwidth of simple metals) [26,27], and the screened interaction does not obey the f -sum rule [28,29]. It was then proposed that self-energy (SE) diagrams with vertex corrections may cancel the spurious sc effects [25,[30][31][32][33]. This fueled a number of notable attempts to include the vertex function in a model fashion: using the plasmon model for the screened interaction [34], neglecting the incoherent part of the electron spectral function [35], employing the Ward identities and a model form of the exchange-correlation kernel [27,[36][37][38], or the sc cumulant expansion [39][40][41]. Although these methods clarified a number of issues, they did not provide an exhaustive picture [42]. The major obstacle for a full-fledged vertex calculation, besides numerical complexity, is the issue of negative spectral densities, first noted by Minnhagen [34,43] and only recently solved by us using a positive-definite diagrammatic expansion [44,45]. Our solution merges many-body perturbation theory (MBPT) and scattering theory, thus returning a positive-semidefinite (PSD) spectral function by construction.
With the PSD tool at our disposal, in this Letter we investigate the influence of vertex corrections on the spectral function of the homogeneous electron gas (HEG), paving the way towards first-principles correlated calculations of bandstructures. We demonstrate that the vertex function leads to a number of novel physical phenomena which cannot be reduced to mere self-consistency cancellations. Stochastic methods, long been used to describe integral properties, are shown to be well suited for the calculation of spectral features too. In fact, our Monte Carlo momentum integration is so efficient that the numerical part of the calculations does not pose any difficulty.
Let us motivate and discuss the PSD approximation used in this work. In terms of dressed electronic propagators and screened interaction W there is a single second order SE diagram Σ (2) = . Its straightforward inclusion, how-  ever, yields negative spectra in some frequency regions. This prohibits the usual probability interpretation and, even worse, it jeopardizes sc calculations since the resulting Green's function (GF) has the wrong analytic structure [46]. The key idea of the PSD scheme [44,45] consists in (1) writing a SE diagram as the sum of its partitions, i.e., diagrams with particle and hole propagators, (2) bisect each partition into two halfdiagrams, (3) add the missing half-diagrams to form a perfect square, and (4) glue the half-diagrams back. For Σ (2) the half-diagrams, see Fig. 1(a-c), contain scatterings with up to three particles and two holes in the final state [44]. The SE partitions stem from the interference between these scatterings and after the PSD treatment one obtains partitions up to the fourth order in W, see Ref. 44 for the full list. Among them there are three which deserve special attention. Σ aa in Fig. 1(d) results from the interference of scattering (a) with itself. As illustrated in Fig. 1(g) Σ aa involves a particlehole (ph) pair (orange area) or a plasmon in the final state (this is the first order effect described by the GW SE). The plus and minus vertices in the SE partitions have the purpose of distinguishing the constituent half-diagrams (resulting from the cut of all propagators with +/− and −/+ vertices). Σ aā in Fig. 1(e) is formed by the interference between the scattering (a), leading to two-particle-one-hole (p f 1 -p f 2q f ) final state, and the same scattering with interchanged particle momenta (indicated withā), see Fig. 1(h). Finally, Σ cc in Fig. 1(f) is formed by the interference between the scattering (c), in which a particle loses its energy by exciting 2 ph pairs, 2 plasmons or a mixture of them, and the same scattering with intechanged particle and hole momenta, see Fig. 1(i). Plasmon generation is a dominant second order scattering process although it has a severely limited phase-space (dark blue line and light-blue area) due to energy and momentum conservation. Higher order terms in W (up to fourth order) arise from other interferences and are needed to assure the overall positivity [44]. In general the PSD procedure leads to a manifestly positive Fermi Golden rule form of the SE, , where the sum runs over all final states of energy E (n) f with (n + 1)-particles and n-holes (r s being Wigner-Seitz radius). The role of high order diagrams is two-fold: they bring new scattering mechanisms into play (hence new rates Γ (n) ) and renormalize them through the perturbative corrections γ (n) i . We already mentioned that one of the motivations for including diagrams beyond GW is the excessive broadening of the spectral features when the level of self-consistency is increased, e.g., G (0) W (0) → GW (0) → GW. As a full sc calculation of the PSD SE of Ref. [44] is out of reach, we partially account for self-consistency by using a G (0) W (0) GF (finite qp-broadening and plasmon satellites) and an RPA screened interaction. Our calculations indicate that higher-order diagrams (aside from bringing in new spectral features) counteract the undesired sc effects, thus suggesting the occurrence of sizable cancellations. We then explore the possibility of producing the PSD results with less diagrams and bare GF's. In the bare GF the chemical potential µ is iteratively adjusted by imposing that the energy of states on the Fermi sphere (where the discontinuity in the momentum distribution n k occurs) is exactly equal to µ [47]. In Fig. 2(a) we compare the rate as obtained from the PSD diagrams of Ref. [44] with G (0) W (0) GF and from the much simpler Σ = Σ aā + Σ cc + Σ aa with qp GF (in both cases we used an RPA W). The left flank and the hight of the peak are in perfect agreement. At energies in the region of plasmon satellites the full PSD rate decays faster but the trend is similar and the impact of this discrepancy on the spectral function is only minor. More calculations at different k (not shown) confirm the agreement between the two SEs. We 10  therefore infer that the relevant scattering mechanisms for a positive-conserving, leading-order vertex correction are those of Fig. 1(g-i). This reduction of diagrams represents an important advance in view of correlated band-structure calculations of solids. In the following we use the vertex correction of Fig. 1(d-f) to calculate qp and plasmon energy dispersions, spectral function, scattering rate, renormalization factor and momentum distribution. a. Results The electron density and the dimensionality completely determine the properties of the HEG; in the 3d case they fix the Fermi momentum and energy as follows: k F = (αr s ) −1 and F = 1/2(αr s ) −2 with α = [4/(9π)] 1/3 . We consider the case of metallic densities r s = 4.0 a B appropriate for, e. g., bulk Na metal. Angle-resolved photoemission experiments have pointed out a substantial narrowing of the occupied band in sodium [48]. In Ref. [26] the origin of this narrowing was ascribed to the effects of electron correlations on the unoccupied bands. Such conclusion is by no means obvious as fully sc-GW results [25] indicate the opposite trend (bandwidth larger by 20% as compared to the noninteracting electron dispersion). As it was pointed out in Ref. [27], vertex corrections rather than screening are crucial to reproduce the experimentally observed dispersion. Our diagrammatic approximation confirms this fact, providing a bandwidth reduction of 27.5% (as compared to the noninteracting case), see Fig. 2(b). Fig. 2(c) also shows the dispersion of plasmon satellites (red and blue curves). In the calculations with vertex corrections (solid) the high-energy plasmon branch smears out and, unlike in the G (0) W (0) approximation (dotted), there is no real solution to the Dyson equation. We further observe an upward renormalization of the G (0) W (0) low-energy plasmon branch, in agreement with experiment [49], as well as the emergence a second branch at lower energy (orange).
Vertex corrections have a sizable impact on the energy and momentum resolved spectral function A(k, ω) = i[G R − G A ](k, ω). In Fig. 3 we display the color plot of A(k, ω) without (top) and with (bottom) vertex corrections. The dotted lines denote the solutions of the real part of the Dyson equation (used to produce the curves in Fig. 2 Appearance of the 2nd plasmon satellite below µ, redistribution of the spectral weight between 1st and 2nd satellite, and further broadening of plasmonic spectral features above µ are the most important findings of this work. They confirm the plasmon-pole model analysis of Ref. [35] that predicted only hole satellites and much broader particle features. Our results call into question the cumulant parameterization of the spectral function in Ref. [39] where no distinction between hole and particle features is made. It is interesting to notice that vertex corrections make the qp-peak sharper. This can be inferred from the explicit SE expression or from the rate Γ of Eq. (1) which we plot in Fig. 4(a-c) for three different values of the momentum. Plasmons do not contribute to the on-shell properties at energies around µ because they carry finite energy at zero momentum (ω pl (q = 0) = 1.881 F for r s = 4). Instead, the life-time of qps in the vicinity of the Fermi sphere is mainly determined by Σ aa (GW SE) involving ph production or by Σ aā (qp → 3qp). The latter, shown as yellow shaded curve in Fig. 4(a-c), contributes with negative sign and leads to the observed reduction of Γ (hence an enhancement of the qp peak) [50]. Such a behavior (alternating series in αr s ) is typical of many perturbation theories. Notice that Σ aa + Σ aā also dominates the asymptotic (ω → ∞) behavior. The scattering with generation of 2 plasmons, contained in Σ cc , plays a crucial role for the off-shell properties as it gives rise to new spectral peaks, see green shaded curve.
In the vicinity of a qp or plasmon (pl) peak the spectral function acquires the form [51]: where α = qp, pl and ω (qp) k = k whereas ω (pl) k = ω pl (k) is the dispersion of plasmon satellites. This expression contains two quantities of physical interest that we computed using our vertex function: the renormalization factor and the broadening of the qp or pl excitations 1/τ (α) The renormalization factor is shown in Fig. 4(d). At the band bottom (k = 0) G (0) W (0) gives only one plasmon satellite Z (pl) = 0.382 whereas our vertex approximation gives two satellites with comparable weight Z (pl) = 0.217 and Z (2pl) = 0.207. Furthermore, the qp weight is reduced from Z (qp) = 0.578 (in G (0) W (0) ) to Z (qp) = 0.550 indicating that the incoherent part of the spectrum gains weight. These two effects cannot be seen in the cumulant expansion scheme [39] which suppresses the Z of higher plasmon satellites according to the Poissonian distribution [49] and, due to the neglect of the coupling between particle and hole seas, yields the same Z (qp) as in G (0) W (0) .
For k = k F vertex corrections reduce only slightly the G (0) W (0) qp renormalization factor. It is known that sc-GW overestimates (Z (qp) = 0.793) the already good G (0) W (0) value Z (qp) = 0.638 (our calculation) or Z (qp) = 0.646 (Hedin [17]). The proposed approximation to the vertex gives Z (qp) = 0.628, which thus remains rather close to QMC results 0.64 to 0.69 (at r s = 3.99 a B ) [6]. At the Fermi momentum Z (qp) can also be deduced from the discontinuity of the momentum distribution function n k [52][53][54]. In Fig. 4(e) we show n k as obtained by a straightforward integration of the smooth part of the spectral function, n k = µ −∞ dω 2π A(k, ω), and by adding the singular contributions analytically. The G (0) W (0) and vertex results are almost indistinguishable.
We finally analyze in Fig. 4(f) the quasiparticle life-time, a measure of electronic correlations [55,56]. In ab initio calculations for realistic systems this quantity is typically estimated using the G 0 W 0 approximation [57][58][59]. However, there have also been attempts to go beyond this level of theory by, e.g., including T -matrix diagrams. In Ref. [60] a reduction of qp life-time (increase of Γ by 50% (70%) in relation to GW for r s = 2.07(4.86)) has been predicted and explained by "the multiple scattering". Our findings show the opposite trend, i.e., an increase of the qp life-time (reduction of Γ by −50% in relation to GW for r s = 4).
b. Conclusions Numerous authors emphasized that the inclusion of the vertex function should remedy the drawbacks of self-consistent calculations [25,27,30,61]. Using our recently proposed diagrammatic analysis we have been able to confirm these expectations and show that this is only a part of the whole picture. Additionally, other second-order processes appear. They can be best described in the language of scattering theory with the link provided by the PSD formalism [44,45].
We fully characterized the spectral function of 3d HEG in the k − ω plane. The main original features that we found are: a second plasmon satellite for holes, redistribution of the spectral weight between hole satellites, reduction of the plasmon spectral weight for particles, bandwidth reduction of the main qp-band. So far these effects have only been partially captured by other, non-diagrammatic methods. Our proposed approach has a universal character and can be extended to firstprinciple calculations of metals. In fact, in going from continuous to discrete translational symmetry the functions simply turn into matrix functions (e. g., Σ(k, ω) → Σ GG (k, ω)), something which does not pose any conceptual difficulties for Monte Carlo momentum integration [62]. Alkali metals for which a vertex function was partially included (typically using a model exchange-correlation kernel [63,64]) is a logical next step for our method.