Diagrammatic expansion for positive spectral functions beyond GW: Application to vertex corrections in the electron gas

We present a diagrammatic approach to construct self-energy approximations within many-body perturbation theory with positive spectral properties. The method cures the problem of negative spectral functions which arises from a straightforward inclusion of vertex diagrams beyond the GW approximation. Our approach consists of a two-steps procedure: we first express the approximate many-body self-energy as a product of half-diagrams and then identify the minimal number of half-diagrams to add in order to form a perfect square. The resulting self-energy is an unconventional sum of self-energy diagrams in which the internal lines of half a diagram are time-ordered Green's functions whereas those of the other half are anti-time-ordered Green's functions, and the lines joining the two halves are either lesser or greater Green's functions. The theory is developed using noninteracting Green's functions and subsequently extended to self-consistent Green's functions. Issues related to the conserving properties of diagrammatic approximations with positive spectral functions are also addressed. As a major application of the formalism we derive the minimal set of additional diagrams to make positive the spectral function of the GW approximation with lowest-order vertex corrections and screened interactions. The method is then applied to vertex corrections in the three-dimensional homogeneous electron gas by using a combination of analytical frequency integrations and numerical Monte-Carlo momentum integrations to evaluate the diagrams.


I. INTRODUCTION
Many-body perturbation theory (MBPT) has provided a systematic way to study electron-electron (electronphonon) interactions in various systems ranging from molecules to solids. 1,2 Within MBPT the interaction effects are included via a self-energy term which is treated perturbatively. One of the widely used self-energy approximations is the GW approximation 3 which consists of replacing the bare interaction with the screened interaction in the first-order exchange diagram. Diagrammatically the GW approximation can be viewed as an infinite summation of polarization diagrams. It is well known that for solids the GW approximation (usually not implemented self-consistently) tends to give band gap values close to the experimental values, thus improving over the density functional calculations (which instead underestimate the values for the band gaps). 4 In spite of some improvements over complementary theories, the self-consistent GW approximation is known to have a number of deficiencies like the washing out of plasmon features and broadened bandwidths in the electron-gaslike metals. 5 For many decades the common argument has then been that the inclusion of vertex corrections would act as a balancing force for the self-consistency, [6][7][8][9] thus, e.g., hampering the washing out of plasmon satellites. Several people have worked on this issue on various levels, [10][11][12][13][14][15][16] but the most interesting result from our point of view is that the straightforward inclusion of vertex corrections beyond the GW level yields negative spectra in some frequency regions, as first noticed by Minnhagen for the electron gas. 10 This deficiency not only prohibits the usual probability interpretation of the spectral function but also generates Green's functions with the wrong analytic properties. In particular the latter feature prevents an iterative self-consistent solution of the Dyson equation since the analytic properties deteriorate with every self-consistency cycle. This unpleasant situation is not limited to the electron gas as it has also been observed in a study of vertex corrections in finite systems. 17,18 As we will explain in depth in this work the problem lies in the structure of the vertex correction and therefore the solution must be sought in the way we use MBPT. There has been very little work on how to generate positive spectral functions from MBPT. The only work on the issue of positivity that we are aware of has been done by Almbladh but in the context of photoemission. Almbladh showed that the positivity of the photocurrent was guaranteed by expressing the photo-emission triangle diagrams as a square of half-diagrams. 19 This, however, was done only for a certain selection of low-order diagrams and it was not indicated how the idea could be applied to the spectral function.
In this paper we put forward a diagrammatic approach to generate self-energy approximations beyond GW yielding positive spectral functions. We start from an expression of the self-energy derived by Danielewicz 20 and use the Keldysh formalism 21 to extract the lesser/greater components (these components are needed to construct the spectral function). Every lesser/greater diagram is partitioned into different contributions, each with internal times integrated over either the forward or the backward branch of the Keldysh contour. The full lesser/greater diagram corresponds to the sum of all possible partitions. We then factorize each partition into half-diagrams by using the Lehmann representation of the Green's function, 22,23 where the one half of the partition consists of time-ordered quantities and the other half consists of anti-time-ordered quantities. The partitioning can be seen as cutting the diagram in half along the lesser/greater Green's function lines. A similar cutting procedure is used in high-energy physics to calculate the imaginary part of diagrams contributing to the scattering amplitudes. [24][25][26][27][28][29] Our cutting rules agree with those of high-energy physics but our derivation is based on the Keldysh formalism, which allows us to advance the theory further. In fact, the positivity of the spectral function entails that the sum of the products of the half-diagrams is the sum of perfect squares. In some situations the MBPT approximation is already a sum of perfect squares, like for the GW approximation. For other self-energy approximations we instead need to add missing half-diagrams to complete the square, like for the GW GGW self-energy, i.e., the lowest-order vertex correction. We acknowledge here that Danielewicz 30 also studied a cutting procedure for the lesser/greater self-energy diagrams and derived a manifestly positive exact formula for the spectral function in terms of retarded/advanced n-point functions, but the issue of how to cure negative spectral functions of approximate self-energies was not discussed in his work. The focus of our work is to study approximate MBPT self-energies and give simple drawing rules to decide whether or not the approximation generate a positive spectral function. If not we provide extra, but still simple, drawing rules to extend to a minimal set of diagrams the MBPT approximation and turn positive the spectral function. This paper is organized as follows. In section II we derive the Lehmann form of the lesser/greater self-energy and relate it to a diagrammatic representation in terms of half-diagrams. Then we describe how to construct a selfenergy approximation with a positive spectral function from a given MBPT approximation by a minimal selection of additional self-energy diagrams. The theory is developed using noninteracting Green's functions and subsequently extended to self-consistent Green's functions. In section III we illustrate the formalism with text-book examples. In section IV we derive the simplest self-energy with vertex corrections and screened interaction yielding a positive spectral function. We then apply the theory to the three-dimensional homogeneous electron gas in section V. We evaluate the self-energy diagrams by using a combination of analytical frequency integrations and numerical Monte-Carlo momentum integrations, and show how the minimal selection of additional diagrams cures the problem of negative spectra. We finally present our conclusions and outlooks in section VI.

II. THEORETICAL FRAMEWORK
Within the realm of Green's function theory the most common techniques to study equilibrium problems are either the zero-temperature formalism or the Matsubara formalism. These are two special cases of the more general Keldysh formalism which is usually applied in the context of non-equilibrium physics beyond linear response. In this work we show that the Keldysh formalism is also the natural tool to develop a diagrammatic theory for positive-definite spectral functions of systems in equilibrium. We consider a system of interacting fermions with Hamiltonian where the field operatorψ (ψ † ) with argument x = rσ annihilates (creates) a fermion in position r with spin σ. In the Keldysh formalism the field operators evolve on the time-loop contour C shown in Fig. 1. Operators on the minus-branch are ordered chronologically while operators on the plus-branch are ordered antichronologically. Letting z 1 and z 2 be two contour-times, the Green's function G(x 1 z 1 , x 2 z 2 ) can be divided into different components G αβ (x 1 t 1 , x 2 t 2 ) depending on the branch α, β = +/− to which z 1 and z 2 belong. For α = β = − we have the time-ordered Green's function In this expression the average . . . is done with some density matrixρ and T is the time-ordering operator. The subscript "H" attached to a general operatorÔ signifies that that operator is in the Heisenberg picturê whereÛ(t 1 , t 2 ) is the time-evolution operator and t 0 is an arbitrary initial time. Reversing the time arrow the G −− is converted into the anti-time-ordered Green's function Integration contour on complex ⇠-plane.
Integration contour on complex ⇠-plane.
FIG. 1: The closed time-loop contour C. The forward branch is denoted with a "−" label while the backward branch is denoted by a "+" label.
whereT orders the operators anti-chronologically. Finally, choosing z 1 and z 2 on different branches we have . (5b) These two last components are equivalently written as G −+ = G < (lesser Green's function) and G +− = G > (greater Green's function), and describe the propagation of an added hole (G < ) or particle (G > ) in the medium.
A. Positive Semidefiniteness of the Exact Self-Energy In equilibrium the Green's function G αβ depends on the time-difference only. Omitting the dependence on the position and spin coordinates x 1 and x 2 the spectral function is defined according to where here and in the following G αβ (ω) denotes the Fourier transform of the Green's function with respect to the time-difference. From the Lehmann representation it is easy to show that iG > (ω) and −iG < (ω), as matrices in the x-space, are positive semidefinite (PSD). From the Dyson equation on the Keldysh contour one can also show that 2 where are the retarded/advanced Green's function and Σ c is the correlation self-energy. Since G A (ω) = [G R (ω)] † the PSD of ∓iG ≶ implies that ∓iΣ ≶ c is PSD and vice versa. Even though one can prove the PSD property of the self-energy using the corresponding property of the Green's function and the Dyson equation, a direct proof starting from a Lehmann representation is not possible since Σ c is not the average of a correlator. In this section we use the Keldysh formalism to provide an alternative proof of the PSD property of Σ c . For the proof we derive a Lehmann-like representation of Σ c and highlight the connection with the diagrammatic expansion. This connection will be extremely useful to generate approximate PSD self-energies from diagrammatic theory. The starting point is the following expression 2,20 and the subscript "irr" signifies that only one-particle irreducible diagrams should be retained. That is, the expansion of the self-energy (9) contains all the diagrams of the two-particle-one-hole correlation function 31 in which the entrance and exit channels cannot be separated by cutting one Green's function line. 32,33 Unlike the definitions (5) of the Green's functions G ≶ , the equations (9) are not averages of a correlator due to the exclusion of reducible diagrams. Nevertheless a Lehmann-like representation for the self-energy can be derived using diagrammatic methods. Let us study, e.g., Σ < c as the same reasoning applies to Σ > c . For simplicity we restrict the discussion to systems at zero temperature and assume a nondegenerate ground state Ψ 0 . Using Eq. (3) and introducing the short-hand notation 1 = x 1 t 1 , 2 = x 2 t 2 , etc. we can rewrite Eq. (9b) as Next we assume that Ψ 0 can be obtained by evolving backward the noninteracting ground state Φ 0 from a distant future time τ (with τ → ∞) to the arbitrary initial time t 0 using an interaction which is switched-on adiabatically, i.e., |Ψ 0 =Û(t 0 , τ )|Φ 0 where the evolution operator is calculated with the time-dependent interaction e η|t−t0| v (η being an infinitesimal energy). This is the standard assumption of the zero-temperature Green's function formalism. Then Eq. (11) becomes (the limit τ → ∞ is implied) where we inserted a completeness relation i |χ i χ i | = 1 and used the group propertyÛ(t 1 , t 0 )Û(t 0 , τ ) =Û(t 1 , τ ) andÛ † (t 0 , τ )Û(t 0 , t 2 ) =Û(τ, t 2 ). Letĉ k ,ĉ † k denote the annihilation and creation operators of a fermion in the kth eigenstate of the noninteracting problem. In Eq. (12) only states |χ i of the form contribute since the operatorγ (γ † ) annihilates (creates) a fermion. The indices p = (p 1 , . . . , p N +1 ) and pq specify the quantum numbers of theĉ andĉ † operators respectively. We conclude that Eq. (12) does not change under the replacement Here the inner sum denotes integrations or summations over sets of p and q quantum numbers with the restriction that p integration runs over the occupied and q integration runs over the unoccupied states, respectively. The prefactor stems from the inner product of the intermediate states, i.e., where P and Q run over all possible permutations of N +1 and N indices with parities (−) P and (−) Q , respectively. We further denoted the permutation group of N elements by π N . Defining the amplitudes the lesser self-energy takes the following compact form To proceed further we need to analyze the amplitudes S and their complex conjugate S * . Under the adiabatic assumption the evolution of the noninteracting ground state Φ 0 from −τ to τ yields Φ 0 up to a phase factor, i.e.,Û with e iα = Φ 0 |Û(τ, −τ )|Φ 0 . Therefore we can write where the time-argument in the operators specifies their position along the interval (−τ, τ ). The time τ + is infinitesimally greater than τ , which assures the correct ordering of the operators. The amplitude S * can now be expanded in powers of the inter-particle interaction v by means of Wick's theorem, and the generic term of the expansion is a connected diagram of noninteracting time ordered Green's functions g −− with external vertices 1 = x 1 t 1 and p, q at time τ , see left diagram in Fig. 2. Following the same steps it is easy to show that which can also be expanded using Wick's theorem, and the generic term of the expansion is a connected diagram of noninteracting anti-time ordered Green's functions g ++ with external vertices 2 = x 2 t 2 and p, q at time τ , see right diagram in Fig. 2. Let us investigate the result of multiplying a diagram of S * N,pq (1) by a diagram of S N,pq (2) and then sum over p and q. In a diagram of S * N,pq the outgoing Green's functions with q-labels and the ingoing Green's functions with p-labels are calculated at the latest possible time τ . Therefore where (x, t x ) is an internal space-spin-time vertex and we introduced the short-hand notation g ij for the matrix elements of g between two spin-orbital states i and j, hence g αβ (x 1 t 1 , x 2 t 2 ) = g αβ x1x2 (t 1 , t 2 ). Similarly in a diagram of S N,pq the ingoing Green's functions with qlabels and the outgoing Green's functions with p-labels are calculated at the latest possible time τ and therefore g ++ yq (t y , τ ) = g > yq (t y , τ ), g ++ py (τ, t y ) = g < py (τ, t y ). Thus the multiplication of a S diagram by a S * diagram and the subsequent sum over p and q involves the sum over q of g > qx g > yq and the sum over p of g < xp g < py . The noninteracting lesser/greater Green's functions can be expanded in the basis of the noninteracting one-particle eigenstates according to 2 (20) where p is the energy of the one-particle eigenstate |p , f is the zero-temperature Fermi function andf = 1 − f .
When taking the product of the left and the right halfdiagram in Fig. 2 we can use relation (21) to replace each of the N products of g > functions by a single g > connecting two internal times in each half-diagram. Similarly we can use relation (22) to replace each of the N + 1 products of a g < by a single g < . The result of this gluing procedure is a diagram with external vertices 1 and 2.
The structure of the diagram is such that all internal vertices to the left of the glued lines have time labels on the minus branch and all internal vertices to the right of the glued lines have time labels on the plus branch. The gluing procedure can be reversed by cutting all Green's function lines between a vertex labeled − and a vertex labeled +. We will refer to this procedure as the cutting rule for a diagram. At this point we observe that the self-energy in Eq. (15) is not the sum of all possible S-S * diagrams due to the subscript "irr". This means that from the diagrams obtained by the gluing procedure we still have the remove the reducible diagrams, i.e., the diagrams which fall apart in two disjoint pieces by cutting a single Green's function line. An obvious case of a reducible diagram is when there is only a single line to glue in Fig. 2. This happens when N = 0 and we only have the single label p 1 . This case is easily taken care of by letting the sum in Eq. (15) start at N = 1 instead. For N > 1 the gluing procedure leads to reducible diagrams whenever the S-diagram can be disjoint into a piece which contains only vertex 1 and a piece which contains the pq vertices by cutting a single Green's function line. We call these S-diagrams reducible and defineS as the sum of irreducible S-diagrams. Note that theS-diagrams can be disjoint into two pieces by cutting a single Green's function line, but then one of the pieces would contain 1 and some of the pq vertices. For instance half-diagrams with a self-energy insertion on the lines we glue together belong toS. An example is shown in Fig. 3 where a tadpole in inserted in one of the p-lines. This half-diagram can be disjoint by a single cut When glued with some irreducible S * diagram it yields an irreducible self-energy diagram. Notice that this diagram is reducible in the "correlation function" sense because cutting a single Green's function line produces disconnected pieces.
but then the vertex 1 would not be isolated from all the pq. Consequently this half-diagram belongs toS and produces irreducible diagrams for the self-energy. From this analysis we conclude that the self-energy can be written as Equation (23) is the Lehmann-like representation of the self-energy and the main result of this section. The irreducible part of products in Eq. (15) has been transformed into products of irreducible parts in Eq. (23). The product of aS-diagram and aS * -diagram yields an irreducible self-energy diagram in which the internal times are either integrated over the minus-branch or over the plus-branch. We call this product a partition of the self-energy diagram. It is now easy to show that the Fourier transform of −iΣ < c has the PSD property. Fourier transformingS andS * and omitting the dependence on the position-spin variables we find In this equation the Fourier transform ofS(t 2 ) (and similarly ofS * (t 1 )) is performed over the time argument t 2 and not over the time difference t 2 −τ , which is ill-defined for τ → ∞. For the right hand side to depend only on t 1 − t 2 the following property has to be fulfilled. From this property we see that F(ω) ≥ 0 and inserting it back into equation (24) we see that F(−ω) is the Fourier transform of the function −iΣ < c (1, 2) with respect to the time difference t 1 − t 2 . So far we have restricted ourselves to the exact selfenergy. In the following section we will explain how a given approximate diagrammatic expression for the selfenergy can be extended, if necessary, to a similar form as in Eq. (23) by an appropriate selection of additional halfdiagrams, thereby ensuring the positivity of its spectral function.

B. Diagrammatic theory of positive spectral functions
The Lehmann-like representation of Eq. (23) brings to light a general and simple rule to calculate the lesser component of a self-energy diagram. A diagram for Σ < c (1, 2) has two external vertices, one with time t 1 on the minusbranch and the other with time t 2 on the plus-branch, and a certain number of internal vertices with times to be integrated over the Keldysh contour. If we assign to each internal vertex a − or a + sign to signify that the corresponding time is integrated over the minus or plus branch then we obtain a division of the original selfenergy diagram, and the full self-energy diagram is the sum of all of them. Since the two-particle interaction is it only connects two vertices with times on the same branch of the Keldysh contour. Since the external vertices are fixed there are 2 n−2 divisions (n > 1) for a diagram with n interaction lines. However, not all such divisions contribute. As shown in Fig. 2 and in Eq. (23) the only divisions appearing in the expansion of the selfenergy are those in which one side of the diagram only contains "−" vertices and the other side of the diagram only contains "+" vertices. All other divisions must therefore be discarded. These are the divisions for which we get a piece disconnected from the external vertices 1 and 2 upon cutting the +/− g-lines (hence these divisions cannot be written as the product of aS-S * diagram). The number of contributing divisions clearly depends on the topological structure of the diagram. In Section II A we called such divisions partitions. As an example consider the Σ < c diagram shown in Fig. 4. The third division vanishes since we get a piece disconnected from 1 and 2 upon cutting the +/− g-lines. This simple diagrammatic rule to extract the lesser self-energy can be viewed as a generalization of the Langreth rules 34 to diagrams which are neither a product nor a convolution of Green's functions. 13 The same diagrammatic rule can alternatively be derived by working in frequency space and by taking into account the conservation of energy at each vertex. 25,26 There remains one issue to address before introducing our diagrammatic theory of PSD spectral functions: How manyS-S * diagrams do lead to the same partition of a self-energy diagram? To answer this question we need to investigate the expansion ofS in terms of Feynman diagrams. Due to the anti-commutation rules for the creation and annihilation operators it follows from the definition of S N,pq in Eq. (18) that a permutation P of � � the labels p and a permutation Q of the labels q simply changes the sign of S N,pq , and hence alsoS N,pq , by a factor (−1) P +Q . Therefore if we let {D (j) pq } with j ∈ I N be the set of all topologically inequivalent diagrams for S N,pq that differ by more than a permutation of the p or q labels, then we can writẽ where P and Q run over all permutations π N +1 and π N of N + 1 and N indices respectively. Inserting Eq. (25) back into Eq. (23) we find This expression can be simplified further by noticing that the composite permutations P •P i and Q•Q i with i = 1, 2 yield the same contribution as the permutations P i and Q i , i.e., (27) since the effect of P and Q is equivalent to a relabeling of the p and q. There are N !(N + 1)! such relabelings and they all give the same contribution. We can therefore simplify Eq. (26) to Every term of the form pq D (j2) corresponds to a unique partition of a Σ < c diagram. Vice versa, every partition of a Σ < c diagram can be written as the product of a unique D-D * diagram for otherwise there should exist more than one way to cut the self-energy along the +/− g-lines.
Equation (28) is an exact rewriting of Σ < c in terms of D diagrams. MBPT approximations to the self-energy consist of the sum of a (finite or infinite) subset of diagrams. An exotic approximation could be, e.g., the one of Fig. 5. Each diagram contains four interaction lines (n = 4) and, therefore, it is divided in 2 n−2 = 4 different ways. It is a simple exercise to see that the left diagram is the sum of three partitions with three +/− g-lines and one partition with five +/− g-lines whereas the right diagram is the sum of two partitions with three +/− g-lines and two partitions with five +/− g-lines. In Fig. 6 we display, e.g., all partitions with five +/− g-lines (N = 2 in Eq. (28)) and how to write them as D-D * diagrams. There are two different D diagrams, say D (a) and D (b) , which are glued as The products of the half-diagrams is represented in the same order as they appear in this mathematical expression in Fig. 6 after the equality sign. The diagrams D (b) in the first line differ only in a permutation of the labels q 1 and q 2 as shown in the left half-diagrams of Fig. 6. Furthermore, the right half-diagram of the first term (D (a) ) and the left half-diagram of the last term (D (a) * ) have the same topological structure but differ in the labeling of q 1 and q 2 . This example is instructive since it highlights the general structure of a MBPT approximation to the selfenergy. The partitioning leads to an expression of the form of Eq. (28) where the domains of the summation indices and the set of permutations are restricted. The couple (j 1 , j 2 ) runs over a subset I N ⊂ I N × I N of the product set of the topologically inequivalent half-diagrams. In the example of Fig. 6 we have I 2 = {(a, b), (b, a)}. Given the couple (j 1 , j 2 ) ∈ I N the permutations P 1 and Q 1 run over a subset π (j1j2) N +1,p ⊂ π N +1 and π (j1j2) N,q ⊂ π N of the permutation groups π N +1 and π N . In the example  Fig. 6 we considered in detail the D-D * diagrams of the selfenergy of Fig. 5 belonging to the set I 2 . The set I 1 in which we cut three +/− g-lines can be analyzed similarly. All other sets I N with N > 2 are empty in this case. In general, however, we have for an approximate MBPT self-energy where the sets I N may contain zero, a finite or an infinite number of elements. An important remark to be made regarding Eq. (29) is that it does not, in general, fulfill the PSD property. Our analysis shows, however, how to formulate the diagrammatic rules to transform a MBPT self-energy into a PSD self-energy by adding the minimal number of partitions. For the PSD property to be fulfilled the self-energy should have the structure of Eq. (23) with some approx-imateS N,pq . LetĨ N ⊂ I N be the smallest subset such thatĨ andπ N +1,p andπ N,q be the smallest subgroups of the permutation groups π N +1 and π N with the propertỹ Mathemeticallyπ N +1,p andπ N,q are given by the intersection of all subgroups containing the subsets π (j1j2) N +1,p and π (j1j2) N,q . The self-energy contains all partitions of Eq. (29) plus other partitions, and each partition is counted only once. Furthermore, taking into account thatπ N +1,p andπ N,q are two subgroups, Eq. (27) is valid for all P ∈π N +1,p and Q ∈π N,q . Hence we can rewrite Eq. (33) to bring out its product structure. We have where d N +1,p and d N,q are the dimensions of the subgroupsπ N +1,p andπ N,q respectively. The self-energy in Eq. (34) is clearly PSD. It is also clear that any reduction of the setsĨ N ,π N +1,p andπ N,q would either not include the original MBPT diagrams or would not fulfill the PSD property. Thus Eq. (34) contains the minimal number of partitions of self-energy diagrams to correct the MBPT self-energy. This concludes the diagrammatic theory to generate PSD spectral functions. In the next sections we work out explicitly some text-book examples and derive the leading PSD self-energy diagrams with vertex corrections and screened interaction, thus going beyond the GW approximation.

C. Self-consistency
Before we discuss some examples in detail we would first like to address the issue of self-consistency which plays an important role in so-called conserving approximations. 2,35 So far we used the noninteracting Green's functions g to evaluate the diagrams. Suppose that a specific selection of partitions guarantees the positivity of the spectral function for Σ ≶ c,PSD [g], where we indicate explicitly the functional dependence of the self-energy on g. Then we can use this self-energy in the Dyson equation to evaluate a new Green's function, which we may call G to distinguish it from the noninteracting g. In the next step we can evaluate our approximate diagrammatic expression for the self-energy in terms of G, i.e., we evaluate Σ ≶ c,PSD [G]. The natural question to ask then is whether Σ ≶ c,PSD [G] still has the PSD property. We now demonstrate that this is indeed the case by a modification of the derivation in Sections II A and II B.
The largest modification involves relations (21) and (22) since they are not valid anymore for general dressed Green's functions. This is, for example, easily demonstrated for the exact interacting lesser Green's function of an N -particle system with ground state energy E 0 , which has the Lehmann representation Here α labels the many-body eigenstates Ψ α,N −1 with energy E α,N −1 of the system with N − 1 particles, Ω α = E α,N −1 − E 0,N are the removal energies and f α (x) = Ψ α,N −1 |ψ(x)|Ψ 0 the so-called Dyson orbitals. Although the Lehmann representation (35) looks formally identical to the expansion in Eq. (19) for g, it does not allow us to derive Eq. (21) anymore since the nonvanishing f α (x) form an overcomplete and nonorthonormal one-particle basis-set (in a noninteracting system most f α (x) are zero and the nonvanishing ones form an orthonormal basis-set). Our strategy, therefore, is to replace Eqs. (19) and (20) with a different relation that still allows us to formulate a cutting rule. Crucial in our reasoning is that a PSD self-energy generates a PSD spectral function for G, as was explained just below Eq. (8). This implies that G < has the form where the removal-part of the spectral function A < xx (ω) ≡ f (ω)A xx (ω) is a self-adjoint and PSD matrix in the one-particle indices for every ω. Denoting by x|a i (ω) the eigenstates of A < xx (ω) with eigenvalue a i (ω), the spectral representation of this matrix can be written as Without loss of generality we shall assume that the eigenstates |a i (ω) form an orthonormal basis-set for every ω.
Due to the PSD property of A < the eigenvalues a i ≥ 0, and therefore we can define the square root of the spectral function according to where for notational convenience we suppressed the ωdependence of the eigenvalues and eigenvectors. Correspondingly we can define the square root of the lesser Green's function according to This function has the property that as follows from a quick calculation using the definitions above. Similarly G > is the integral over a positive spectral function, and we can therefore define the square root √ G > with the property that The relations (36) and (37) provide a new cutting rule for a self-energy diagram with a dressed Green's function. Whenever we cut a self-energy diagram we obtain half-diagrams with outgoing lines √ G < and √ G > . Using these modified half-diagrams we obtain an equation that is identical in structure to Eq. (33). The only thing that changes in Eq. (33) is that the sums over pq are replaced by integrals over y's andt's. However this does not change the quadratic structure of the equation. We therefore conclude that Σ < c,PSD [G] also is PSD. From this new self-energy we can use the Dyson equation to calculate a new Green's function which has again a PSD spectral function and can be decomposed as in Eqs. (36) and (37), thereby yielding yet another PSD self-energy. By repeating the procedure we obtain a series of PSD Green's functions. If this series converges to a limiting Green's function then we have solved the Dyson equation self-consistently for our approximate Σ < c,PSD in which both the Green's function and the self-energy are PSD.
The conclusion of this analysis is that our diagrammatic approach to PSD spectral functions also applies to self-consistent perturbation theory. Of course, in order to avoid double countings, the partitions of which Σ < c,PSD is made of should be skeletonic, i.e., should not contain selfenergy insertions. 2 An important approximation that has this structure is the GW approximation which we study in more detail below. In general, however, it should be emphasized that a PSD self-energy made exclusively of all the partitions of conserving diagrams is rare. In fact, approximations which are simultaneously conserving and PSD are exceptional.

III. EXAMPLES
In this section we apply the formalism developed in Section II to some illustrative examples. a. Single bubble diagram Let us first consider the first bubble diagram shown in Fig. 7(A). The lesser component of this self-energy reads which can be partitioned in only one way. Upon cutting along the +/− g-lines we find where the D diagram reads Equation (39) is already of the form in Eq. (33) and therefore the first bubble diagram produces a PSD spectrum.
b. Second-order exchange The exchange diagram of the first bubble diagram is shown in Fig. 7(B) and the application of the formalism is more interesting. The lesser component of the self-energy now reads . This diagram too can be partitioned in only one way. However, upon cutting the +/− g-lines we find a product of D diagrams differing by a permutation P (p 1 , p 2 ) = (p 2 , p 1 ) p1p2q1 is the same as in Eq. (40). The smallest subgroup of π 2 which contains P is π 2 itself and the domain j 1 = a and j 2 = a is already of the form I 1 ×Ĩ 1 . Therefore we can form a PSD self-energy by tak-ingĨ 1 = {a},π 2,p = π 2 andπ 1,q = {1} = π 1 . In this way we end up with the diagrams shown in Fig. 7(C). This is the second-Born (2B) approximation, which we have now shown to give a PSD spectrum. The secondorder exchange diagram is particularly instructive since it shows that the PSD outcome of the sum of MBPT diagrams is not the sum of the PSD outcome of the MBPT diagrams taken separately. Indeed the PSD outcome of the first-bubble diagram is the first bubble-diagram itself whereas the PSD outcome of the second-order exchange diagram is the 2B approximation. This implies that if we had summed the PSD outcomes of these two separate diagrams we would have counted the first-bubble diagram twice. No double counting occurs if we apply the rules in Eqs. (30)(31)(32) to the 2B approximation; the PSD outcome would be the 2B approximation itself.
c. Two bubbles As a third example let us consider the sum of the first and second bubble diagrams as shown in Fig. 8(A). After distributing the pluses and minuses over the internal vertices we find that the lesser selfenergy diagram can be written in terms of two types of D-diagrams The three D-D * diagrams in this expression are represented in the bottom line of Fig. 8(A). We observe that the only permutation appearing is the identity permutation. Therefore, according to the rules in Eqs. (30)(31)(32), we only need to find the smallestĨ 1 such that , a), (a, b), (b, a)}. This is simplyĨ 1 = {a, b}. Thus, the PSD outcome of the self-energy is (1) and the corresponding diagrams are shown in Fig. 8(B). In accordance with our notation a diagram with no plus/minus on the internal vertices represents the full diagram, i.e., the sum of all possible partitions. d. GW and T-matrix approximations With similar arguments it is easy to show that the GW self-energy has the PSD property. We take the RPA W = v +vP W with P x1x2 (z 1 , z 2 ) = −ig x1x2 (z 1 , z 2 )g x2x1 (z 2 , z 1 ). Then where we took into account that if P ++ appears to the left of P < and/or P −− appears to the right of P < then the division is not a partition since a cut along the +/− g-lines generates a disconnected + and/or − island, see Section II B. It is a matter of simple algebra to prove that (43) where the diagrams D (j) are defined as in Fig. 9(A). Equation (43) is clearly of the form in Eq. (33). In a similar fashion it can be shown that the symmetrized version of the T-matrix approximation also has the PSD property since it can be written as in the second row of Eq. (43) with D diagrams given in Fig. 9(B).

IV. PSD SELF-ENERGY BEYOND GW
The GW self-energy is the leading order approximation in the screened interaction W . Despite the numerous successful applications of the GW approximation in reproducing experimental spectra there is also a commensurable number of examples for which the GW approximation fails, pointing out the importance of including vertex corrections. The next to leading order approximation in W is represented by the diagram Σ vert of Fig. 10 and contains a vertex correction. Unfortunately the selfenergy Σ GW+v ≡ Σ GW + Σ vert is not PSD and the resulting spectral function has the undesired feature of being negative in some region of the frequency space. In this section we use the rules of Eqs. (30)(31)(32) to identify the minimal number of additional partitions for constructing the leading order beyond-GW self-energy with the PSD property.
On the right hand side of Fig. 10 the self-energy Σ vert is written as the sum of four partitions. Let us determine the underlying D diagrams. We observe that the screened interactions W ++ or W −− can be partitioned in only one way with all internal polarizations P ++ or P −− . Indeed a partition of W ++ in which appears a P > does necessarily contain also a P < and hence the cut along  the +/− g-lines would generate a disconnected island. A similar argument applies to W −− . Accordingly the first diagram on the right hand side of Fig. 10 is the sum of the partitions shown in Fig. 11(A), and can be written as the sum of products between the D (j) diagrams of the GW approximation, see Fig. 9(A), and theD (j) diagrams displayed in Fig. 11(B). The same is true for the second diagram of Fig. 10 provided we exchange D (j) with D (j) . The partitions of the third diagram of Fig. 10 are very simple due to the presence of only W ++ and W −− , see Fig. 12. The result is the sum of products between D (j) diagrams and permuted D (j) diagrams. Finally the fourth diagram of Fig. 10 is partitioned as illustrated in Fig. 13(A), and can be written as the sum of products between two D (ij) diagrams, see Fig. 13(B), where i refers to the number of top bubbles and j to the number of bottom bubbles. In conclusion (omitting the dependence on 1 and 2) Thus Σ GW+v is the sum of a contribution Σ 3 containing partitions with three +/− g-lines and a contribution Σ 5 (the last term in Eq. (44)) containing partitions with five +/− g-lines. From the rules in Eqs. (30)(31)(32) we see that Σ 5 has already a PSD structure. Instead Σ 3 should be corrected according to Σ 3 → Σ 3,PSD with This self-energy contains the original four partitions of Σ 3 plus four more partitions arising from the permutation of the two p dangling lines. The latter are explicitly worked out in Fig. 14. Collecting all results together we conclude that the leading order self-energy diagrams with screened interactions and vertex corrections yielding a PSD spectral function are those in Fig. 15. Here we recall that a diagram with no +/− on the internal vertices is the full diagram (hence the sum over all possible partitions). Noteworthy the minimal completion of the square requires a fourth order diagram in W .

V. VERTEX CORRECTIONS IN THE HOMOGENEOUS ELECTRON GAS
The numerical implementation of the full self-energy of Fig. 15 is rather demanding. We observe, however, that the exclusion of the last partition of Fig. 10 leads to a much simpler PSD self-energy since no permuted D diagrams appear. It is straightforward to show that in this case the rules of Eqs. (30-32) lead to the following PSD self-energy where the corresponding D diagrams are defined in Fig. 9(A) and Fig. 11(B). The diagrammatic representation of Eq. (46) is shown in Fig. 16. This self-energy too goes beyond the GW approximation, but the vertex correction is only partial. The three-dimensional homogeneous electron gas (3d HEG) is one of the most studied correlated many-body system. 36 We still lack detailed knowledge of one directly observable quantity -the spectral function A(k, ω)when departing from the on-shell energy, i.e., when ω ≈ k . Discrepancies with experimental measurements contributed to the debates on the position of satellites, 37,38 bandwidth of simple metals, 14,39 cancellation of vertex function and self-consistency effects, 5,40,41 and the spectral function shape. 42,43 Despite numerous efforts there are just a few results that go beyond the GW approximation in the study of quasiparticle properties. Analytically, these are the results of Onsager et al. 44 on the second-order exchange energy and of Ziesche 45 and of Glasser and Lamb 46 on the second-order exchange self-energy. These works only contain analytic results for on-shell self-energy and only a contribution of the bare Coulomb interaction is included. The screened Coulomb interaction is possible to treat numerically. The non self-consistent calculations were performed by Hedin 3 and Lundqvist. [47][48][49] It took three decades to implement the same approach partially or fully self-consistently. This was achieved in works by von Barth and Holm. 5,50 The non-positivity of the spectral function first observed by Minnhagen 10,40 hindered systematic diagrammatic explorations and stimulated development of synthetic approaches: analyzing real time Kadanoff-Baym dynamics, 51 neglecting the incoherent part of the electron spectral function, 12 employing the Ward identities and a model form of the exchangecorrelation kernel, 14,15,52,53 or the self-consistent cumulant expansion. 37 Using the presented formalism it is now possible to pursue the diagrammatic route. In this section we present the results of a non self-consistent calculation. Thus, we evaluate the diagrams in Fig. 16 for 3d HEG using the analytical frequency and numerical Monte-Carlo momentum integrations. The method was developed in a prior publication, 43 however, especially the analytical frequency integration part had to be substantially extended. For HEG the bare time-ordered Green's function can be written as a function of frequency ω and momentum k as For self-consistent calculations this equation can be extended to include multiple poles (e.g., to describe quasiparticle satellites) for each momentum value. In this work we perform a one-shot calculation and therefore set B(k) = n k and A(k) = 1 − n k with 0 ≤ n k ≤ 1 denoting the occupation number and k the energy of the state with momentum k. For the energy dispersion we used the prescription by Hedin 2,3,54 to put the pole of the dressed G = G 0 + G 0 Σ c [G 0 ]G with Fermi momentum k = k F correctly at the Fermi surface, i.e., k = k 2 /2 + µ − F where µ is the chemical potential and F = k 2 F /2 is the Fermi energy. Analogous definitions have been used for the anti-time-ordered, lesser and greater Green's functions.
Generally, each interaction line in Fig. 16 can designate either the bare Coulomb interaction or include scattering with generation of a plasmon or a particle-hole pair. On the RPA level the plasmon oscillator strength t(q) vanishes when its dispersion curve Ω(q) enters the particlehole continuum at the critical wave number q c . Because plasmonic excitations exist on the bounded momentum interval they are especially convenient for numerics: momentum integrations need to be performed on the finite interval. Physically, plasmons also dominate the particlehole response in electron liquids with typical metallic densities. Therefore, in present work we only treat the plasmonic contributions. The screened Coulomb inter-action is treated in the plasmon pole approximation, i.e. W (k, ω) ≈ W 0 (k, ω) (cf. Eq. 25.11 of Ref. 54): where we denote w(q) = t(q)Ω 2 (0)/Ω(q) and 0 ≤ t(q) ≤ 1 is the plasmonic spectral weight with t(0) = 1 and t(q c ) = 0. Analogous defininitions have been used for the other Keldysh components of W 0 . Our numerical approach also allows us to include contributions from the particle-hole continuum to exhaust the f -sum rule for the dielectric function ε(q, ω): where Ω(0) = 4 αrs 3π F , with α = [4/(9π)] 1/3 ≈ 0.521 and r s = 1/(αk F ) the Wigner-Seitz radius. Monte-Carlo momentum integration of these terms is, however, more involved as it requires an extra integration for each interaction line.
The frequency integrations can be done completely analytically (facilitated by the mathematica computer algebra system) whereas for the remaining momentum integrations one has to rely on numerics. The frequency integration is implemented for a general case of timeordered G 0 (k, ω) and W 0 (k, ω). Since each correlator in the self-energy expression comprises two terms the total number of terms grows geometrically with the diagrams order. In order to optimize the calculation we implemented an approach where for each ω-integration the program closes the integration contour in such a way that the least number of terms is generated. One might argue that the Hedin set of equations gives a prescription on which half of the complex plane the integration should be closed (see for example Eqs. (A30) of Ref. 3). However, it is easy to verify that this is only relevant for the first order term. Higher order terms decay faster than 1/ω at infinity (ω here denotes the frequency to be integrated over) and, therefore, the choice of the half-plane is not important. Finally, we notice that it is sufficient to consider only the lesser self-energy. The greater selfenergy component required to describe properties of the states above the Fermi level is obtained from symmetry consideration.
Our analytic approach covers the case of Fig. 16(A) diagram where all partitions are included. Other diagrams on this figure contain only a subset of the partitions of the MBPT diagrams from which they originate. The corresponding analytic expressions can be easily extracted from the general result by analyzing the frequency dependence. The overall ω-dependence can readily be obtained on paper by integrating δ-functions in lesser and greater correlators. Despite the omission of several partitions the final result is bulky and requires an additional simplification using A(k) + B(k) = 1.  15: (Color online) Leading order beyond-GW self-energy with the PSD property. Thick wavy lines denote the screened Coulomb interaction in the random phase approximation.  (47) and (48) and notations below Eqs. (49). For instance g5(z −y3), connecting vertices with "−" labels as in C) and D), corresponds to the time-ordered Green's function (see also Eq. (2)): where we define In these equations we adopt the following notations: The quantities i and Ω i are energies labeled by the momenta, as shown in Fig. 16 This leads to the following density-dependent prefactors: where i = 1, 2, 3 denotes the diagram's order. From the particle-hole symmetry we can obtain the greater self-energy (Σ > ). For this it is sufficient to replace B i with A i and vice versa, and to change the sign of each Ω i . Upon a close inspection it is evident that the four terms in Eq. (49) can be combined together as follows: The first order integrand (cf. Eq. (49a)) is multiplied by a full square. Since we know that the GW self-energy fulfills the PSD property the same holds for the sum of the four terms in Fig. 16. This analytic conclusion is numerically confirmed. In order to make Eqs. (49) suitable for Monte-Carlo integration the following has to be done: (i) integrate in spherical coordinates with zenith direction along the z vector, (ii) one angular integration is trivially done (2π) because the system is isotropic, (iii) map integration variables to the interval [0, 1]. The speed and quality of a pseudo-random number generator is very important for the present calculations. We use the Mersenne twister 19937 generator as implemented in gsl library combined with a highly efficient jump ahead method 57 due to Haramoto et al. 55 for parallelization. The method is roughly 3 times faster than the standard fortran implementation. We used roughly 10 12 Monte-Carlo real-izations to get the full frequency dependence for each momentum value. Real self-energy parts were computed using the Hilbert transform: where Σ x (k) is the frequency independent exchange self-energy and Γ(k, is the rate operator. Each term of Eq. (49) was separately computed as shown in Figs. 17,18 for two values of k. The first order self-energy is well understood. One of its marked features is the existence of logarithmic singularities at ω = k ± Ω(0). 54 These singularities have never been observed in an experiment and are believed to be washed out by higher order contributions. Although our calculations do not preclude this conclusion we found that our   Fig. 17 for a different momentum k = 1.2kF (the energy ω is measured with respect to µ). The inset magnifies the region of the logarithmic singularity. Notice the almost complete cancellation of Σ > c (k, ω) for ω > k + Ω(0).
selection of second-order terms makes the singularities even more pronounced. The real part of the retarded self-energy of Eq. (51) is displayed in Fig. 19. We observe that the real part of the first order and the complete self-energy cross the y-axis at almost the same point which is equal to µ/ F − 1. This implies that the higher order corrections do not change the chemical potential appreciably. As expected, the rate operator is everywhere positive despite a large negative contribution of the second order term. We notice an almost complete cancellations of different order  2 Γ(k, ω) = −Im Σ R c (k, ω) of the homogeneous electron gas at the density of rs = 4 and k = kF (the energy ω is measured with respect to µ). Shaded dashed curves denote the first order calculation. Full lines denote contribution from diagrams shown in Fig. 16 plus the first order contribution from the particle-hole excitations. The horizontal line bounding the shaded area in the graph of Re Σ R crosses the y-axis at the value of the exchange self-energy Σx(kF )/ F = − 2α π rs ≈ −1.327.
terms beyond the singularities, i.e. ω > k + Ω(0) for particle (k > k F ) and ω < k − Ω(0) for hole (k < k F ) states. High accuracy of the Monte-Carlo integration was required to get the cancellations properly. This is especially important at metallic densities where different orders have comparable contributions. Due to the density scaling (see Eq. (50)) the first order self-energy becomes dominant at large densities (r s → 0), while the third order is largest in the correlated low density regime (r s → ∞).
The selection of diagrams of Fig. 16 and computed in this work in the plasmon pole approximation describes scattering processes accompanied by the emission or absorption of one plasmon. Correspondingly, the scattering operator has pronounced (more narrow in comparison with the first order result) peaks at ω = k ± Ω(0). But how important are the remaining contributions? The simplest first order term, absent in the plasmon pole approximation, but included in the results of Fig. 19, involves generation of a single particle-hole pair. Due to less restrictions on the available phase space for scattering it is important for ω → ∞ and also determines the life-time of quasiparticles in the vicinity of the Fermi energy. It also gives rise to secondary peaks in Fig. 19. These, however, are not related to scattering mechanisms involving generation of two plasmons. Inclusion of such processes is important for the interpretation of multiple satellites in the spectral function. To lowest order they result from the partition of the second-order self-energy diagram included in Fig. 15, but omitted in Fig. 16. Respective calculations are on the way and will be the subject of a forthcoming publication.

VI. CONCLUSIONS
Approximations of MBPT to the self-energy can lead to unphysical density of states, with a negative spectral weight in some frequency region. This undesired feature entails unphysical results on the system properties and makes self-consistent calculations impossible due to a progressive deterioration of the analytic properties of the Green's function. In 1985 Almbladh proposed a diagrammatic perturbation theory of the photoemission current, 56 and in a subsequent paper 19 he elaborated on the theory and gave a prescription how to combine diagrams of different order to get a physically sensible result: the positive-definite photoemission current. These ideas are precursors for our theory. Using the Green's function of the Keldysh formalism we developed a method to construct manifestly PSD spectral functions. The method becomes particularly lucid when expressed in diagrammatic language as it amounts to apply a few simple drawing rules.
We derive a Lehmann-like representation of the exact self-energy and show that it is given by the sum of squares of irreducible correlators. We then elucidate the connection between the diagrammatic expansion of the irreducible correlators and MBPT. Any lesser/greater self-energy diagram can be partitioned into two halves with internal time-vertices on opposite branches of the Keldysh contour. Thus, by simply drawing diagrams and assigning a sign to the internal vertices we are able to extend to a minimal set of diagrams any MBPT approximation and to generate PSD spectral functions. Several important MBPT approximations, such as the GW or T-matrix approximations, do not require any corrections. Our theory applies equally well to diagrammatic expansions with noninteracting and with self-consistent Green's functions because PSD self-energies do preserve the correct analytic structure.
In standard MBPT approximations the straightforward inclusion of vertex corrections inevitably ruin the PSD property and, hence, our additional diagrams must be included. Remarkably, these diagrams are of higher order. For instance, the inclusion of the full first-order vertex leads to diagrams of the fourth order in the screened interaction. Required computational power to numerically evaluate them is immense. Fortunately, excluding some partitions allows us to construct an approximation containing diagrams of maximally third order. They are feasible for numerics as our calculations for the 3d HEG demonstrate.
Even though we only presented in detail the formalism for the spectral function, the same ideas apply to the spectrum of the density response function. This extension, however, goes beyond the scope of the present work and will be presented elsewhere.

GS acknowledges funding by MIUR FIRB Grant
No. RBFR12SW0J. YP acknowledges support by DFG-SFB762. AMU would like to thank the Alfred Kordelin Foundation for support. RvL would like to thank the Academy of Finland for support.