Partial self-consistency and analyticity in many-body perturbation theory: particle number conservation and a generalized sum rule

We consider a general class of approximations which guarantees the conservation of particle number in many-body perturbation theory. To do this we extend the concept of $\Phi$-derivability for the self-energy $\Sigma$ to a larger class of diagrammatic terms in which only some of the Green's function lines contain the fully dressed Green's function $G$. We call the corresponding approximations for $\Sigma$ partially $\Phi$-derivable. A special subclass of such approximations, which are gauge-invariant, is obtained by dressing loops in the diagrammatic expansion of $\Phi$ consistently with $G$. These approximations are number conserving but do not have to fulfill other conservation laws, such as the conservation of energy and momentum. From our formalism we can easily deduce if commonly used approximations will fulfill the continuity equation, which implies particle number conservation. We further show how the concept of partial $\Phi$-derivability plays an important role in the derivation of a generalized sum rule for the particle number, which reduces to the Luttinger-Ward theorem in the case of a homogeneous electron gas, and the Friedel sum rule in the case of the Anderson model. To do this we need to ensure that the Green's function has certain complex analytic properties, which can be guaranteed if the spectral function is positive semi-definite.The latter property can be ensured for a subset of partially $\Phi$-derivable approximations for the self-energy, namely those that can be constructed from squares of so-called half-diagrams. In case the analytic requirements are not fulfilled we highlight a number of subtle issues related to branch cuts, pole structure and multi-valuedness. We also show that various schemes of computing the particle number are consistent for particle number conserving approximations.


I. INTRODUCTION
Many-body perturbation theory using Green's functions is a successful and powerful method for studying interacting systems in and out of equilibrium.Its most general non-equilibrium version 1 is routinely used for very diverse situations such as quantum transport (e.g. through single molecules 2 ), cold atoms in optical lattices 3 and transient dynamics (e.g.transient photoabsorption 4 ).In the limiting case of equilibrium systems it can also be used to study standard photo-absorption and photo-emission spectroscopies 5,6 .Moreover the nonequilibrium formalism even in this limiting case has found useful applications to derive approximations that guarantee positive spectral distributions as it allows for a direct expansion for those quantities [7][8][9] .
One of the reasons of the successes of many-body perturbation theory is the use of diagrammatic techniques, where one can sum up diagrams to infinite order in a practical way.However, if partial resummations are used, when one sums over a subset over all diagrams, important conservation laws can be violated.The most important laws are the fulfillment of the continuity equation, momentum and angular momentum conservation, and energy conservation.Approximations which fulfill all of these conservation laws are referred to as conserving 10 .
Baym 11 considered approximate generating functionals for the single-particle Green's function G, the Φfunctional, from which the self-energy Σ could be constructed via Σ = δΦ/δG.If G is obtained from the selfconsistent solution of Dyson's equation for such a Σ, the approximation is called Φ-derivable and is automatically conserving.Commonly used Φ-derivable many-body approximations (fully self-consistent) are the Hartree-Fock approximation, the Second Born Approximation (2BA), the GW approximation and the T-matrix approximation.
While fully self-consistent schemes have many conceptual advantages, they usually carry a high computation cost, as well as other features that are not desirable, both in and out of equilibrium.In fully self-consistent calculations, it has been shown 12,13 that strong timedependent fields can yield artificial steady states in finite systems.Partial self-consistency lessened this effect.Another method related to partial self-consistency is the Generalized Kadanoff-Baym Ansatz (GKBA), where one approximates certain time-nondiagonal elements of the Green's function and can be made conserving, which is also free of artificial damping [14][15][16] .
0][21][22][23][24][25] for examples of fully selfconsistent GW calculations).These types of approximations tend to yield accurate band gaps and spectral properties.Furthermore, for GW for the electron gas, full self-consistency worsens the spectral features compared to non-selfconsistent calculations 20,26 , even if the scheme is energy conserving.While these drawbacks will be cured by vertex corrections to GW , there has been a long discussion in the literature if one should focus on non-selfconsistent vertex corrections, or full selfconsistency (see, for example, [27][28][29] ).Moreover, partially self-consistent approximations are expected to yield reasonable energies, especially if one uses the Luttinger-Ward functional for calculating the total energy 30 .
Thus, while the concept of a Φ-derivable approximation is extremely useful, it is quite restrictive.As soon as a scheme is not fully self-consistent, the concept of Φ-derivability is lost, together with a convenient way of seeing if the chosen approximation is particle number conserving.The class of approximations that do conserve particle number is larger than those generated by Φ-functionals, however.For example, the partially self-consistent GW 0 approximation (where the screened interaction W is kept fixed to be W 0 during the selfconsistency cycle) and the non-selfconsistent G 0 W 0 approximation are not Φ-derivable.However, GW 0 is particle number conserving (first shown for the equilibrium electron gas 26,31 , and later more generally 32 ) but not fully conserving 33 , while G 0 W 0 is not particle number conserving 32,34 .
Another property that can be violated by summing over subclasses of diagrams is the positivity of the spectral functions of the Green's function and the polarizability.Positivity is an important ingredient for two reasons.First of all, it is related to a probability interpretation of photo-emission and photo-absorption processes, which naturally requires positive probabilities.Second, it guarantees correct analyticity and causality properties of the various propagators, which is a necessary condition for performing self-consistent calculations, as the analytic properties in general will deteriorate every iteration cycle.
The conserving and analytic properties are related but not equivalent.A situation in which both of these properties play a role is in the derivation of sum rules for the particle number, such as the Luttinger-Ward 35 for the electron gas and the Friedel sum rule 36,37 for the Anderson model.Another occasion is the derivation of the Ward identity and the frequency sum rule for the density response function.Conserving approximation in connection with the Luttinger-Ward sum rule were also studied in the context of the non-perturbative approximation of dynamical mean field theory in Ref. 38 .
In this work we generalize the concept of Φ-derivability to partially self-consistent schemes.While the resulting approximations are not fully conserving, a large class of many-body approximations used in the literature is shown to be included in this larger class of approximations.The formalism provides an easy way to see if an approximation conserves particle number or not.Since we are using NEGF, the concept of partially self-consistent Φ-functionals is valid both in and out of equilibrium.Using the new formalism, we study particle number sum rules in equilibrium in their most general settings, and highlight the importance of particle conservation.We state sufficient conditions for the sum rules to be valid in a given diagrammatic approximation.In particular, we stress that the sum rules are valid provided that G has the correct analytical properties, and we study the consequences of incorrect analytical properties, including issues related to pole structure and multivaluedness of the logarithm.We exemplify the formalism and sum rules in model systems of quantum transport.
The paper is structured as follows: In the first part of this work, Sec.II, we will introduce the concept of partial Φ-derivability.Subsequently in Sec.III we will discuss the complex analytic properties for approximate Green's function and self-energies and the way these properties affect the calculation of the particle number.Then finally, in Sec.IV, we discuss the derivation of a generalized sum rule for the particle number and illustrate the concepts by explicit numerical calculations.We end with our conclusions in Sec.V.

II. PARTIAL Φ-DERIVABILITY
Our main object of study is the single-particle Green's function G (1, 2), where 1 = (r 1 , σ 1 , z 1 ) is a collective space-spin-time variable.The Green's function is defined on a time contour γ that has a forward, backward, and Matsubara (imaginary-time) branch, 1 Tr T e −i γ dz Ĥ(z) , where ψ ( ψ † ) are fermionic field destruction (creation) operators, and T is the time-ordering operator on the contour γ.The time-dependent Hamiltonian Ĥ(z) = Ĥ0 (z) + Ĥint (z) is given by where h(x, z) is the single-particle part and v(x, x , z) is the two-body interaction.The non-equilibrium formalism allows for both time-dependent single-particle potentials, such as electric fields, as well as time-dependent two-body interactions, such as an adiabatic switch-on.
To keep the discussion general, the time-dependence will not be specified further at this point.We also define, for convenience, x = (r, σ) and to highlight the temporal dependence without being dependent on the spatial basis.
The reason for considering the contour-ordered G is that general time-dependent systems at finite temperature can be considered, in and out of equilibrium.When the contour times z are on the Matsubara branch, we obtain the Matsubara Green's function as ĜM (τ 1 , τ 2 ) = Ĝ(t 0 − iτ 1 , t 0 − iτ 2 ).For general time-dependent systems, we obtain the lesser Ĝ< and greater Ĝ> Green's functions from Ĝ< (t 1 , t 2 ) = Ĝ(t − 1 , t + 2 ), and Ĝ> (t , where z = t − /t + is on the forward/backward branch.All single-particle quantities can be obtained by the contourordered Green's function.For example, if we take t 1 = t 2 , we obtain the time-dependent density n(x, t) and the time-dependent current density j(x, t) as n(x, t) = −iG < (x, t; x, t) (5) The equations of motion for G(1, 2) is given by the Kadanoff-Baym equations 1,39 , with the Kubo-Martin-Schwinger (KMS) boundary con- is the self-energy, which has to be approximated.An equivalent way of writing the equations of motion is given by the Dyson equations, where G0 (1, 2) is the corresponding Green's function with Ĥint = 0 that satisfies the KMS boundary conditions.In a given approximation to Σ, diagrammatic or nor, there is no guarantee that the resulting scheme will be conserving.Since we in this work will focus on the fulfillment of the continuity equation, we here briefly describe it.We subtract Eq. ( 8) from Eq. ( 7), and put 2 = 1 + , where 1 + denotes (x 1 , z + 1 ) and z + 1 denotes a time infinitesimally later than z 1 on the contour.We then obtain Eq.( 10) is a continuity equation on the contour, with a source/drain term caused by the interactions.Thus, approximations that guarantee (local) particle conservation fulfill for all times, which is then a condition on the approximate Σ (1, 2).If this condition is fulfilled, the total particle number is conserved in time.
Particle number conservation is also relevant in equilibrium situations, as can be seen as follows.Let us consider a finite system in which the particles are non-interacting at t 0 with a switch-on of the interaction for times t > t 0 .If the continuity equation is fulfilled for all times, the switching of the interaction cannot change the number of particles.In particular, this is true for an adiabatic switch on of the interaction.Thus, particle conservation is an issue also in equilibrium calculations, a point also stressed in Refs. 32,33.For example, within a conserving approximation the particle number can not depend on the bond length of a molecule 33 .
Baym 11 introduced a convenient way to generate approximations that are automatically conserving, with the use of a diagrammatically defined functional Φ[G].If Φ[G] is invariant under gauge transformations, translations, rotations and time, the self-energy given via will yield a conserving approximation when the Kadanoff-Baym equations, Eq. ( 7) and Eq. ( 8), are solved selfconsistently with this Σ.In particular, the source term in the continuity equation, Eq.( 10), will vanish.The scheme is then said to be Φ-derivable.We now extend the idea of Φ-derivable approximations to include partially self-consistent schemes, and find sufficient criteria for the fulfillment of the continuity equation.We define a partially dressed functional Φ[G, G 0 ], in which the Green's function lines can be dressed with either G or additional, fixed, Green's functions G 0 .The fixed Green's functions could be outputs from a Hartree-Fock or density-functional theory calculation, but this is in no way necessary.In fact, sets of different fixed Green's functions can be used.We define the corresponding self- Solving the Kadanoff-Baym equations, Eq.( 7) and Eq.( 8), with this self-energy defines a partially self-consistent scheme.We will refer to such schemes as partially Φ- In FIG. 1 we show the example of a ring Φ-diagram dressed in various ways.By dressing only one Green's function line with the full G, we obtain non-selfconsistent (one-shot, or single-shot) approximations.Combining FIG. 1 a) and FIG.2d), we obtain one-shot 2nd Born.Moreover, FIG. 1 a) is part of the G 0 W 0 approximation.The G 0 W 0 approximation can be obtained by considering all higher order ring diagrams dressed with only one G.If all lines except one are dressed with G, as in FIG. 1  c), we obtain several types of diagrams.In one diagram the baseline in Σ is fully dressed, while in the other diagrams the baseline is not.Any one of the Σ-diagram in FIG. 1 c) is by itself not partially Φ-derivable, since closing it with a G-line and subsequently differentiating yields also the two other Σ-diagrams.The first of the Σ diagrams in FIG. 1 c) is part of the G 0 W approximation.It is then readily seen that the G 0 W approximation for Σ is not partially Φ-derivable.The same can be seen in FIG.3b), where we get two non-equivalent classes of Σdiagrams from one Φ-diagram.Each separate Σ-diagram is not partially Φ-derivable, but their sum is.
An important class of diagrams is obtained when each loop in Φ is dressed with either all G or all G 0 separately.We will refer to this as a consistent dressing of the loops.Cutting a G-line from consistently dressed Φ-diagrams yields Σ[G, G 0 ]-diagrams with fully dressed baselines.An example of a consistently dressed Φ-diagram is shown in FIG. 1 b).This diagram is part of the GW 0 approximation.The remaining terms in the GW 0 approximation can be obtained by considering all higher order ring diagrams with only one loop consistently dressed with G, see also FIG. 3a).In the next section, we will show that all approximations coming from consistently dressed Φfunctionals are number conserving.
In a given Φ-derivable approximation, the diagrammatic structure of Φ[G, G 0 ] is similar to the one for the fully dressed Φ where c nk is a symmetry prefactor, n is the number of interaction lines and k labels the Σ diagrams.For a collection of Σ k diagrams obtained from the same Φ-diagram, the number c nk is generally given by one divided by the number of G-lines in Φ.The proof of this statement is analogous to the proof that c nk = 1/2n in case that all 2n G-lines are dressed 1 .As an example, in FIG.3a) the GW 0 diagram has a prefactor of 1/2, while the T -matrix diagram has 1/4.In fact, for GW 0 each Φ-diagram has the same prefactor 1/2, independent on the order of the diagram.Note that the Feynman rules for Φ include a symmetry prefactor as well, related to the dimension of the symmetry group of each diagram 1 .

−→
Examples of partially dressed ring Φ-functionals (left) and the corresponding Σ = δΦ δG (right).a) A single partially dressed loop, thus Φ is not gauge invariant and the approximation is not number conserving.This diagram is part of one-shot 2nd Born, as well as G0W0.b) A single, fully dressed loop, which gives a gauge-invariant, and thus particle-conserving approximation.This diagram is part of the GW0 approximation.c) Partially dressed loops, not gauge invariant.Several types of self-energy diagrams are obtained, where the first is part of the G0W approximation.d) Fully dressed, thus conserving.Part of 2nd Born, and GW .

A. Number conservation
As was shown by Baym 11 for fully dressed Φfunctionals, particle number conservation is guaranteed by the inherent gauge invariance in Φ [G].Here, we repeat the same steps for partially dressed Φ[G, G 0 ]-functionals.
For an arbitrary variation in G(1, 2), by definition the corresponding Φ-functional will change according to For an infinitesimal gauge, the variation in Φ becomes For gauge-invariant Φ-functionals, δΦ = 0 for any gauge Λ, and thus we have This is the source term which appeared in the continuity equation, Eq. (10).Thus, a gauge-invariant Φ leads to an approximation that conserves the particle number.If the approximate Φ used is not gauge invariant, the source term will in general be non-zero, leading to a violation of the particle number.Since the resulting Σ from a gauge-invariant Φ has a fully dressed baseline and consistently dressed loops, we can also infer that Σ (2) .G Λ and Σ[G Λ , G 0 ] fulfills the Kadanoff-Baym equations, Eq. ( 7) and Eq. ( 8), with a Hamiltonian shifted with the gauge.One can then use the same considerations as in Ref. 1 to show that the Ward identity is satisfied, which then also leads to the fulfillment of the frequency sum rule for the density response function.As in the fully conserving case, this is however only valid under additional assumptions on the correct analytical structure of the response functions 1,40 .
An approximate, gauge invariant, Φ-functional is guaranteed to yield an approximation that fulfill the continuity equation for all times, independent on the particular shape of the time dependence in the Hamiltonian.For example, the relation is valid for a time-dependent singleparticle potential, such as a time-dependent electric bias, relevant for quantum transport.The importance of using particle conserving approximations is exemplified by nonselfconsistent 2nd Born and G 0 W 0 , which can severely violate the continuity equation in quantum transport. 2,41or gauge-dependent approximations, particle conservation can be violated already in ground-state calculations.A clear example can be found in Ref. 33 , where the H 2 molecule was considered, starting with 2 particles in the non-interacting ground state.After an adiabatic switching, it was found that the particle number depended on the molecular bond length for G 0 W 0 , while it did not for GW 0 .
We also stress that any linear combination of consistently dressed Φ-functionals generates particle number conserving approximations.One can then generate particle number conserving approximations when combining many-body methods, such as GW + T-matrix 42 and the FLEX approximation 43,44 , provided that each Φ-diagram is consistently dressed.
Finally, we stress that gauge-invariant partial Φderivability guarantees number conservation, but not other types of conservation laws ensured by a fully conserving scheme, such as momentum, angular momentum, and energy conservation 11 .For example, the conservation of energy depends on the time-invariance of all four Green's function lines joining an interaction line, and as such it is not enough for each loop to be dressed consistently separately.Furthermore, no mention has been done about how large the violations of the conservation laws are, since this also depend on how strongly correlated the system is 34 .In equilibrium, the violation of particle number seem to be small in some systems 32,33 , but larger for other systems, 34 see also Sec.IV C. In biased systems in quantum transport, the violation in the particle current can be as large as the current itself 2,41 .

A. Analytical continuation and causality
In the previous section we discussed particle number conserving approximations.In our discussion the particle number was defined by Eq.( 12) using the Green's function in real time.In the following we will discuss only the equilibrium situation.In that case there are at least two others ways to calculate the particle number.The first one is from an integral over the spectral function, as we will discuss in much more detail below.Another way is from the derivative of the grand canonical potential Ω with respect to the chemical potential µ as Baym showed that for a conserving approximation the latter equation yields the same result as Eq.( 12).In his derivation he, however, implicitly assumed some analytic requirements of the Green's function, which are not guaranteed to be fulfilled.In the reminder of the paper we will show that all three ways of calculating the particle number yield the same results for partially Φ-derivable approximations provided that they preserve the correct complex analytic properties.We will further discuss a generalized sum rule for the particle number, the validity of which depends crucially on the analyticity of the Green's function.For these reasons we will in this section first review the analytic properties of the Green's function and how they relate to the positive semi-definiteness of the spectral function.
For systems in equilibrium at inverse temperature β, the integrals over the γ contour reduces to integrals over the imaginary axis.We can then make use of the Matsubara representation ĜM with Matsubara frequencies ω m = 2m+1 −iβ .The exact Green's function can be rewritten using the Lehmann representation 45 as where Â(ω) is the spectral function.
We are interested in the analytical continuation of this function to the complex frequency plane.It is not obvious that such a continuation is unique, since ĜM (ω m ) is only defined on isolated points.Nevertheless, if we restrict the continuation to be analytic and bounded at infinity, the continuation is unique, and is given by We use ζ as a complex frequency, while real frequencies are denoted by ω.We note that Eq.( 23) immediately implies that ĜM is analytic away from the real axis.Eq.( 23) can be used to define the commonly used retarded and advanced Green's functions in their respective half-planes as The real-frequency retarded and advanced Green's functions are then obtained as ĜM (ω ± iη) = ĜR/A (ω + µ), where ( ĜA ) † (ω) = ĜR (ω) and η a positive infinitesimal.ĜM (ζ) is discontinuous when crossing the real axis, with the discontinuity given by ĜR as can be deduced from Eq.( 22).The functions ĜR/A (ζ) are analytic in their respective half-planes, which is what we will refer to as the correct analytical properties.The discussion so far has been independent of the selfenergy.However, in practice, we often start from an equation of the self-energy.We therefore want to relate the analytical properties of G to those of Σ.This connection is provided by the Dyson equation, Eq.( 9), which in terms of Matsubara notation becomes ĜM (ω m ) = 1 where for operators Â we use the notation Â−1 = 1/ Â. 46 The same considerations for ĜM (ζ) also applies to the Matsubara self-energy.Its analytic continuation is where Γ(ω) is the rate operator.ΣR/A (ζ) are defined as in Eq.( 24).Furthermore, ΣM (ω ±iη) = ΣR/A (ω +µ), and The discussion so far considered the exact Green's function and self-energy.However, in practice we typically approximate the self-energy and solve the Dyson equation.In a given approximation for Σ the Green's function does not necessarily have the correct analytical properties.For instance, the denominator in Eq.( 28) may have zeros away from the real axis in which case ĜM (ζ) has poles away from the real axis and consequently a representation as in Eq.( 22) does not exist.However, such non-analyticities are not possible if the approximate selfenergy is of the form of Eq.( 27) in which the rate operator Γ(ω) is Positive Semi-Definite (PSD).With this we mean that ϕ| Γ(ω)|ϕ ≥ 0 for any one-particle basis function ϕ.In this case one can show that ĜM (ζ) is analytic away from the real axis and has a representation as in Eq.( 22) with Â(ω) being PSD as well 789 .For this reason a diagrammatic perturbation method based on so-called half-diagrams has been devised in order to guarantee that the rate function is always PSD 789 .
In case an approximate form for ΣM does give rise to poles for ĜM a generalized spectral representation, as in Eq.( 22), is possible.In Appendix C we show that if ĜM is analytic except for simple poles away from the real axis we can write ĜM where αl is the residue matrix, ξ R l is the location of a pole in the upper half-plane, and the spectral function Â(ω) is also in this case given by Eq.( 25).If ĜM is obtained from Eq.(28) then from taking the limit |ζ| → ∞ (assuming that the approximate Σ is bounded at infinity) we see from where 1 is the unit operator.Thus, the presence of poles can change the norm of the spectral function.
While discussing the analytic properties, we would like to clear up a possible point of confusion regarding the analytic continuation of the retarded and advanced functions G R/A (ζ) from the upper/lower half planes in which they are analytic, across the real axis to the other half planes.Such continuations are useful for derivations with the help of Cauchy's theorem and have also been useful for interpreting quasi-particle life-times in terms of analytic properties of the analytically continued G R/A (ζ).
Due to the discontinuity in G M (ζ), the analytical continuation from the upper half-plane, where G M (ζ) = G R (ζ), into the lower half-plane yields a completely different function than G A (ζ).The analytically continued G R (ζ) is used in several textbooks (see, for example, Refs. 45,47) when the quasi-particle peak of the spectral function has a Lorentzian shape.A typical example is the homogeneous electron gas in which case we can regard the spectral function as a scalar function rather than a matrix, as it is diagonal in momentum space.We consider a Lorentzian with width Γ as We obtain G R (ζ) = 1 ζ+iΓ , Im ζ > 0 from Eq.( 23).The analytical continuation is given by the same expression, but defined for all ζ and has a pole in the lower halfplane.A useful application of the analytic continuation is the calculation of the retarded Green's function in timespace, by the following procedure.For t < 0, G R (ζ)e −iζt decays exponentially in the upper half-plane, and we can integrate around a closed half-circle in the upper half-plane.Since G R (ζ) has no pole there, G R (t) = 0 for t < 0, befitting of a retarded function.This also shows that the existence of poles in the upper half-plane will in general violate the causality property of G R (t).For t > 0, we integrate around a half-circle in the lower half-plane and find Thus, the location of the pole in the lower half-plane yields the lifetime of an added particle.This interpretation was entirely based on the Lorentzian form of the quasi-particle peak of the spectral function.However, it is known for the homogeneous electron gas that the moment dω ω 2 A(ω) of the spectral function is finite, which precludes a Lorentzian form 48 .A study of the short-time properties of the spectral function 48,49 shows that for this system a more realistic form of the retarded Green's function in real time is: with γ and τ parameters.This Green's function has a Gaussian short-time behavior and an exponential longtime behavior.We would like to demonstrate that such a short-time behavior can be described by a simple spec- , where a > 0 determines the width of the single peak, and concentrate on the analytical properties.From Eq.( 23), we obtain where Erfi is the imaginary error function, analytic in the whole complex plane.As it should, which is analytic for all ζ in the complex plane.Thus, for a Gaussian spectral function, G R (ζ) does not have any poles in the lower half-plane, as opposed to the case with a Lorentzian spectral function.
for |ζ| → ∞.This affects the calculation for G R (t) using contour integration.We can close a contour in the upper half-plane for t < 0 yielding G R (t) = 0, but we cannot close the contour in the lower half-plane for t > 0, since the integral on the half-circle diverges, see FIG. 4. Nevertheless, G R (t) is well-defined for all t, and can be obtained by Fourier transforming G R (ω) directly, yielding We see that, similarly to the Lorentzian case, the width of the spectral function determines the lifetime of an added particle excitation.However, in this case the lifetime is unrelated to the properties of G R (ζ) analytically continued to the lower half-plane.The main message that we want to give in this section is therefore that the analytic continuations of the retarded and advanced functions beyond their original analytic domains in general can give rise to a richer analytical structure than just simple poles.This not only includes unbounded analytic functions, but, for instance, also functions with algebraic or logarithmic branch cuts.

B. PSD and analyticity
One way of ensuring the correct analytical properties of G M (ζ) is to ensure that the chosen diagrammatic approximation retain the PSD property.If the diagrammatic structure has this property, then the correct analytical properties are guaranteed 50 .Recently, it has been shown 7,9 that by choosing diagrams in a specific way, by only keeping diagrams that can be built up from so-called half-diagrams, the PSD property and thus the correct analytical properties are guaranteed at zero temperature.Moreover, the PSD properties depend only on the diagrammatic structure, but not on the dressing of the loops provided that we dress with a G coming from a PSD approximation.It is thus important to note that conservingness and PSD are two completely different properties, a point which we feel is not stressed enough in the literature.For example, G 0 W 0 is not conserving, but PSD, while the fully dressed second-order exchange approximation (see FIG. 2b)) is conserving, but not PSD. 7Thus, we are guaranteed to have the correct analytical properties for the former approximation, but not the latter.
Non-PSD spectral functions has been observed 7,51 in studies of the homogeneous electron gas, when one goes beyond GW and takes vertex corrections into account.The lowest order vertex correction is the second-order exchange diagram, where the interaction lines have been dressed with the screened interaction.Moreover, it has also been observed 51 that in this approximation a pole can be created in the lower half-plane for G A (ζ), violating the analytical properties.Also in finite systems, restricting to certain classes of diagrams can violate positivity for the photo-absorption spectra. 40,52o elucidate the above considerations, and to illustrate Eq.( 29), we consider a very simple model of an interacting self-energy.We assume that the analytically continued Σ R (ζ) has a single pole on the real axis, with residue Σ R (ζ) is analytic in the upper half-plane, and α > 0(< 0) corresponds to a positive (negative) definite rate operator Γ which has two poles, at is analytic in the upper and lower half-planes, but otherwise we obtain two poles in the separate half-planes.The imaginary part is If G R (ζ) is non-analytic in the upper half-plane, it is impossible to write G R (ζ) in the form of Eq. ( 23).Nevertheless, we can define the spectral function as A(ω) = − 1 π Im G R (ω) due to Eq. ( 29), which will be the sum of two delta-functions in case the poles are on the real axis.Now, however, the integral over the spectral function will depend on the position on the pole.The integral can be obtained by integrating G R (ζ) over a semi-circle in the upper half-plane.In this model, G R (ζ) → 1/ζ, ζ → ∞, and thus the half-circle integration yields πi.If 4α < −b 2 , then there is a pole ω + in the upper half-plane, with residue This gives the integral as The spectral function is plotted in FIG. 5 for different values of α.As can be seen, if α > 0, A(ω) is positive definite.If −b 2 < 4α < 0, A(ω) has one negative and one positive peak.In this simple model, we can thus identify three regions.
• α > 0. The rate operator is positive definite, and thus the spectral function is positive definite.This leads to that G R (ζ) is analytic in the upper halfplane, and dωA(ω) = 1.
• 4α < −b 2 .The rate operator is negative definite, and A(ω) = 0.There is a pole in the upper halfplane for G R (ζ) and is of the form of Eq.( 29) in the upper half-plane; where Keeping in mind that Σ R (ζ) is a very simple model of an interacting self-energy, the case α > 0 corresponds to when we can write the self-energy as a sum of squares of half-diagrams 9,50 .The case α < 0 corresponds to when the self-energy is not a sum of squares of half-diagrams.However, we see that it is not guaranteed that the correct analytic properties are broken.Taking α as a measure of the correlation effects, for α < 0 we find a critical point when the analytical properties are broken.While being a simple model, the same trend can be observed in the Anderson model, discussed in Sec.IV C.

IV. GAUGE INVARIANCE AND SUM RULES
Gauge invariant Φ-functionals leads to approximations that obey well-known sum rules in equilibrium.Such sum rules have been discussed before (see, e.g. 35,53,54).Here, we repeat the derivations for approximations that are partially Φ-derivable, but we keep a general structure and focus on the importance of the analytical properties.
We consider the number of particles in equilibrium, given by the Matsubara Green's function as where η → 0 + after the summation.According to our previous discussion, we can analytically continue ĜM (ω m ) to ĜM (ζ) when Im ζ = 0, for any finite temperature.Here, we wish to point out an issue regarding the analytical properties.Matsubara sums can be rewritten according to the following rule (see, for example, Ref. 45 ), which follows from considering two half-circles, in the upper and lower halfplane, when Q where f (ω) = 1/(e βω + 1) is the Fermi function.By taking Q = ĜM and using Eq.( 46), we can write N as would have simple poles, we must add their residues to Eq.( 48) as these poles give an extra contribution in the derivation of Eq.47.We conclude, that if ĜM (ζ) is non-analytic, then −iTr ĜM (τ, τ + ) and This shows yet another inconsistency if one works with approximations that do not guarantee the correct analytical properties.Note, however, that no assumption was made on the existence of a generating Φ-functional.
Yet another way of obtaining the particle number is via the grand canonical potential, as was considered by Baym 11 for conserving schemes, via Eq.(20).For a partially Φ-derivable approximation with the correct analytical properties, this give the same particle number as the two other definitions discussed.For details, see Appendix D.
We now focus on rewriting Eq.( 46).Since the continuation of ĜM (ω m ) to ĜM (ζ) is analytic at the Matsubara frequencies, we can perform derivatives around ω m .We can thus write Eq.( 46) as where we have defined I 1 and I 2 , and taken the branch cut of the complex logarithm to be the negative real axis.
The relation used for the matrix logarithm, is valid for all Q(ω) that are smooth and invertible but not necessarily diagonalizable, see Appendix B. The reason for introducing I 1 and I 2 , is that at zero temperature I 2 is related to variations of the Φ-functional, and is zero for number-conserving approximations.Furthermore, I 1 can be integrated analytically.Eq.( 49) simplifies at zero temperature.We will first discuss I 2 .For small temperature T , the distance between the Matsubara frequencies, dω = 2πi/β, becomes small.In the limit T → 0, a Matsubara sum can be written as an integral over the upper and lower half-plane, lim where the principal value of the integral is Eq. ( 51) is valid if Q(ω m ) is not too singular around 0, see Appendix A, and also the discussion in Ref. 55 .With this relation, we can write I 2 as T → 0 as In Eq.( 52) we have made the assumption of integrability and replaced the principal value of the integral by the integral itself.We will here discuss these assumptions.ĜM (ζ) is analytic and thus integrable for Im ζ = 0, and has a discontinuity around ζ = 0 given by the spectral function at µ, Eq.( 25).Thus, assuming that Â(µ) is finite, ĜM (ζ) is integrable.By the same argument, ΣM (ζ) and ∂ ΣM ∂ζ are integrable for Im ζ = 0. From Eq.( 27), from which we see that if Γ(ω) decays quickly enough at µ, ∂ Σ ∂ζ is continuous also at ζ = 0.For Fermi liquids, Γ(µ) ∼ (ω − µ) 2 , for ω close to µ, 56 and as such Eq.( 53) is integrable.However, for other systems, such as the one-dimensional electron gas 57 , or in Mott insulators 58 , the rate operator can decay slower or even diverge at µ. Excluding such cases, the integral in Eq.( 52) is welldefined, and equal to its principal value.
I 2 , also referred to as the Luttinger integral, can be related to the Φ-functional under special types of variations in G, namely those that shift the Matsubara frequencies.As an example, consider the consistently dressed ring di- where spatial integration is suppressed in the second line, to focus on the ω-dependence.The Kronecker delta represents frequency conservation at every vertex, due to the time-locality of the interaction.Φ evaluated at a G which has shifted frequencies, G M (ω m +δω), (defined via its analytical continuation, Eq.( 23), as G M (ζ)| ζ=ωm+δω ), would in general be different from the unshifted case.
For the special case of a shift corresponding to one (or more) Matsubara frequency, i.e. δω = 2πi β which gives ω m + δω = ω m+1 , we can relabel the sum for m and n in Eq.( 54).Since we perform the same relabeling for both m and n, the Kronecker delta does not change, which leaves Φ invariant.If the loop was not consistently dressed, Φ would change since the Kronecker delta would change.The argument works for a general diagram, and thus we see that if the loops are consistently dressed, δΦ = 0 when ω m → ω m+1 in G M .Otherwise, δΦ = 0 in general.
The reason why consistently dressed Φ-diagrams are invariant under frequency shifts is due to the interaction being time-local.The loops in Φ can then be considered separately.This type of argument can not, however, be used for other types of conservation laws.For example, momentum conservation comes from the invariance of Φ under spatial translations in G 11 , for interactions that depend only on the interparticle distance.For such a shift in G in Eq.( 54), by a variable transform we only shift one index in the interaction, and as such Φ is not invariant.For a Φ-diagram to be conserving, we need, in general, full Φ-derivability.
Consider again the change in Φ when we make frequency shift in G.This is given by Eq.( 16), which can be written using Matsubara Green's functions as Let us fix the shift to δω = 2πi β , and take the zerotemperature limit of Eq.( 55).Here, we cannot use Eq.( 51) directly, since the integrand is too singular around 0 (see Appendix A).The limit is δΦ → + lim By making use of partial integration, we obtain Thus, at zero temperature, we can relate I 2 to special types of variations of Φ.For a partially Φ-derivable approximation, I 2 = 0 if Φ is gauge invariant.Otherwise, I 2 = 0 in general.We now derive expressions for I 1 and I 2 involving integrals over real frequencies at zero temperature.Using Eq.( 47) we obtain I 2 as Eq.( 59) can also be derived from Eq.( 58) by considering the contour integral around two quarter-circles in the upper and lower half-plane.
If one or more of the eigenvalues traverse the branch cut at (−∞, 0) as function of ω, a unique logarithm cannot be found, and the argument principle (see, e.g., Ref. 59 ) has to be used instead of Eq.( 61).Crossing the branch cut changes the logarithm with 2πi if it is crossed from above, and −2πi if it is crossed from below.However, if the rate operator is PSD no eigenvalue can cross the real axis, as we will now show.
Eq.( 60) is a sum rule valid for those approximations that give analytic ĜM (ζ) and ΣM (ζ).Furthermore, if the approximation comes from a gauge-invariant Φ, I 2 = 0.For example, the PSD approximations GW and GW 0 will fulfill Eq.( 60) with I 2 = 0, while the PSD approximation G 0 W 0 fulfills Eq.( 60) with I 2 = 0 in general.We now give four different examples and limits of the sum rule.

A. Luttinger-Ward theorem
In the case of the homogeneous electron gas, Eq.( 60) leads immediately to what is known as the Luttinger-Ward theorem. 35See Ref. 57 for a more detailed discussion of this important sum rule in the electron gas.For homogeneous systems, the Hamiltonian and self-energy are diagonal in the momentum-basis, and Σ R (p, µ) = Σ A (p, µ).Assuming a number-conserving approximation, I 2 = 0 and using where θ is the step function, we obtain the Luttinger-Ward theorem where p are the single-particle energies, V is the volume and the factor of 2 accounts for spin.Thus, the Luttinger-Ward theorem is valid for those numberconserving approximations which retain the correct analytical properties.

B. The Friedel sum rule
Eq.( 60) also allows for a useful relation in a completely different context, namely the Friedel sum rule in quantum transport.We consider a small region in space (e.g. a molecular region), where it is important to consider interaction effects carefully, coupled to macroscopic reservoirs (e.g.metallic leads) where a mean-field description is satisfactory.This implies that the frequency-dependent part of the self-energy is non-zero only in the central region.
The sum rule of Eq.( 60) still applies for the total number of particles, but usually we are more interested in the properties of the central region.The trace Tr C over the central region gives the number of particles N C in that region.By the use of embedding techniques 60,61 the Green's function in the central region can be written as ĜR,A (ω) = 1 where the embedding self-energy Σemb takes the environment into account in an exact way, and ĥC is the Hamiltonian for the disconnected central region.The trace over the basis of the finite region gives the same type of sum rule as in Eq.( 60), provided we replace Σ by Σ + Σemb : ΣR emb (ω) enters I 2 , which now has the form We can thus split the contributions as I 2 = I 2,M B + I 2,emb .Now, I 2 = 0 in general, even for non-interacting systems.For a number conserving scheme, I 2,M B = 0.I 1 can be written in terms of the eigenvalues λ k (ω) of − ĜR (ω).Since we will make use of this, we write out the form explicitly.
To recover the usual formulation of the Friedel sum rule, we take our central system to consist of a single interacting site with energy , and use the wide-band limit approximation, Σ R/A emb (ω) = ∓iΓ.Eq.( 68) becomes we obtain By using a gauge independent Φ, I 2,M B = 0, and we obtain the well-known Friedel sum rule 36,37 .Thus, we can regard Eq.( 68) as a generalized Friedel sum rule.

C. Analytical properties in the Anderson model
In this section, we show explicitly, for the singleimpurity Anderson model in the zero-temperature limit, that even conserving approximations can have spectral functions that can become negative.The Anderson model is a simple model in quantum transport, which consists of a single interacting dot in contact to a featureless reservoir with coupling strength Γ.The Hamiltonian of the dot is where nσ is the number operator for spin σ and U is the interaction strength.Furthermore, we consider nonmagnetic situations, and as such where Γ is the coupling to the reservoir.As our approximation to Σ R (ω), we take various 2nd order diagrams; the single-shot 2nd Born approximation FIG.1a) + FIG.2a), the single-shot 2nd order exchange approximation FIG.2a), and the self-consistent 2nd order exchange approximation FIG.2b).The 2nd Born approximation is PSD, while the two 2nd order exchange approximations are not.We write here the explicit expression for Σ(z 1 , z 2 ) for a general system with local interactions, since we will make use of larger systems in Sec.IV D. 1, and c = 1 for 2nd Born.This simple structure is due to the interaction being space-local.Furthermore, since 2nd Born is a PSD approximation, we have immediately that the rate operator of 2nd order exchange is negative semi-definite.
The explicit expression for ΣR kl (ω) in equilibrium is obtained from the Langreth rules 36 , and Fourier transforming (for more details, see Ref. 62 .) where the lesser and greater Green's functions in equilibrium are given by the fluctuation-dissipation theorem 1 For single-shot 2nd order exchange, Σ R (ω) was evaluated using a Green's function from a self-consistent Hartree-Fock calculation.The self-consistent 2nd order exchange calculations did not converge if we would start from the Hartree-Fock Green's functions, however (see below).Instead, as in Ref. 63 , we performed selfconsistency with smaller values of U , and then starting new calculations with this initial guess.
In FIG. 6, we show A(ω) for the three approximations.As expected, the 2nd order exchange approximation yields spectral functions which can become negative.In the one-shot case, one can by increasing U make the spectral function negative in a large region.
dωA(ω) = −0.2 for one-shot 2nd order exchange, indicating a pole in the upper half-plane of G R (ζ).For the other approximations, dωA(ω) = 1.Similar to the case of α in Sec.III B, we find a crossover (not shown) with respect to U above which a pole moves into the upper half-plane.Note also that since the particle number is , it can become negative.However, the particle number defined via the spectral function can be unphysical due to the pole structure, see the discussion of the next section.
We now discuss the Friedel sum rule, Eq.( 68), for the Anderson model.The single-shot 2nd Born is PSD and thus fulfills the sum rule, but since the approximation , and ending at µ = 0, for the same situation as in FIG. 6.The 2nd order exchange approximations are not PSD, and thus −G R (ω) can cross the real axis.Single-shot 2OE breaks analyticity as well as the validity of the sum rules.Self-consistent 2OE does not.The single-shot 2nd Born approximation is PSD, and thus never crosses the real axis.For illustrative purposes, we have divided G R from 2OE0 with 3 and multiplied with 2 for 2BA0.
is not particle conserving, I 2 = 0.The 2nd order exchange approximations are not PSD, and are thus not guaranteed to fulfill the sum rule.For our parameters, we found that self-consistent 2OE fulfilled the sum rule with I 2 = 0, while for single-shot 2OE, we find that To understand the term −2, we show in FIG.7 −G R (ω) as parametrized by ω from ω = −∞ to ω = µ = 0.As can be seen, −G R (ω) from 2nd Born never crosses the real axis, which is guaranteed by its PSD property.However, −G R (ω) obtained from single-shot 2OE crosses the negative real axis, and thus breaks the assumptions in the derivations of the sum rules.Since we change Riemann sheet, but still use the principal branch of the logarithm, we miss an additional term of 2πi.This trans-lates into a term −2 in I 1 , as can be seen from Eq.( 70).The self-consistent 2OE also cross the real axis, but not at the branch cut.In fact, if we push the parameters further (increase U ), analyticity breaks during the selfconsistent cycle, and we are unable to reach convergence.This again underscores the importance of PSD approximations in self-consistent approximations.
Finally, we note that using only the second-order exchange diagram can seem like a constructed example.However, the same considerations apply if the diagram is dressed with the screened interaction W .The second order diagram would then correspond to vertex corrections from GW .This approximation can yield a non-PSD spectral function 7 , and here we show that the same considerations apply in a model in quantum transport.

D. Example of the generalized sum rule
In the last example, we consider a 3-site linear chain at zero temperature to elucidate the generalized Friedel sum rule, Eq.( 68), in a finite region.The Hamiltonian is where h ii = , h ij = −1 for nearest neighbors, and zero otherwise.As contacts, we choose two one-dimensional tight-binding leads with band-widths of 4, attached to the first and third site, respectively, with coupling strength −1.The 3 × 3 Green's function matrix for the central region is 1][62] .
We considered four different approximation to ΣR (ω), with various amounts of dressing: The self-consistent ring diagram FIG.1d), the partially self-consistent ring diagram FIG.1b), the single-shot ring diagram FIG.1a), and the single-shot 2nd order exchange diagram FIG. 2. Of these, the first three approximations are PSD, and as such the total number of particles in the central region N C is given by the generalized Friedel sum rule, N C = I 1 + I 2 , where I 1 can be conveniently written in terms of the eigenvalues λ k (ω) to − ĜR (ω), using Eq.( 70), and I 2 = I 2,M B +I 2,emb is given by Eq.( 69).As discussed before, for the fully self-consistent and the consistently dressed partially self-consistent approximations I 2,M B = 0, while for the single-shot cases I 2,M B = 0, in general.
To illustrate the sum rule, we also calculate N C according to N C = 0 −∞ dωTr Â(ω) .To explore a large parameter range, we sweep with the gate voltage , symmetrically from = −U/2.N C , I 1 , I 2,emb and I 2,M B are shown in FIG. 8, for all the different approximations.We note that the correction from I 2,emb is small, but becomes larger as we move away from = −U/2.Note also that the contribution from I 2,emb would increase with more sites connected to leads and lead coupling strength.We also see, as expected, that I 2,M B = 0 for the consistently dressed approximations.For the single-shot case, I 2,M B can be very large.At = −U/2, I 2,M B = 0, but is large away from this point.We plot I 1 , which is the resulting particle number one would obtain if I 2 is not taken into account.The magnitude of I 2,M B makes N C and I 1 very different, and in fact I 1 increases dramatically when raising the gate potential.Not taking I 2,M B into account gave a similar erroneous trend for the Anderson model using single-shot 2nd Born in Ref. 64 , as was also pointed out by Ref. 63 .We also note that in the single-shot ring approximation, N C increases slightly as increases.This non-physical behavior is due to only taking the ring diagram into account.Single-shot 2nd Born (i.e. the ring diagram and the 2nd order exchange diagram) does not exhibit this non-physical behavior (not shown).
Finally, in the non-PSD 2nd order exchange approximation N C behaves non-physically and can even become negative.The norm of the spectral function (not shown) also becomes negative for certain parameter values, similarly to the case of the Anderson model.The change of norm indicates poles in the upper half-plane for ĜR (ζ), according to Eq. (31).We also note that the particle number we calculate is N A = µ −∞ dωTr Â(ω) .According to the discussion we had above, when there are poles in the upper half-plane for ĜR (ζ), with residue a R k = 0 and location ξ R k , N A does not have to equal N = −iTr ĜM (τ, τ + ) , since we have to add the residues in Eq. (48).For simplicity, we discuss the situation in a single-site case.From Eq.( 31), we have in the single-site case that At zero temperature, a pole can only contribute to the particle number if it is located to the left of µ, Assuming that the jump in N A ( ) occurs at the same time as when Re ξ R ( ) equals µ, we found that the jump was canceled by the contribution from the residue, calculated from the norm in Eq.( 81).This suggests that N is a more physical quantity than N A , in general.Numerically, we find N C − I 1 − I 2 = 0 to a high degree of accuracy for all the ring approximations.For the non-PSD 2nd order exchange approximation, we find instead FIG. 8. Number of particles NC in a 3-site system, for different dressings of ring approximations and single-shot 2nd order exchange.We show the self-consistent ring approximation (FIG.1d)), partially dressed ring diagram, (FIG.1b)), single-shot ring diagram (FIG.1a)), and single-shot 2nd order exchange (FIG.2a)).I 2,emb and I2,MB are the embedding and many-body contributions to the Luttinger integral, Eq.( 69), respectively.Note the different scale for I2 in the non-selfconsistent calculations.
that N C − I 1 − I 2 varies, depending on the parameters.This can be understood by studying how the eigenvalues of − ĜR (ω) change as ω goes from −∞ to µ.To illustrate, we track λ k (ω for the PSD self-consistent ring approximation, shown in FIG. 9, and the non-PSD single-shot 2nd order exchange approximation, shown in FIG.10. The eigenvalues λ k (ω) of the PSD ring approximation in FIG. 9 do not cross the real axis when we vary ω, as expected.The eigenvalues of the non-PSD approximation in FIG. 10 do, however.For = −U/2, all three eigenvalues cross the branch cut.This means that, identical to the discussion we had above, for each crossing eigenvalue we obtain an additional factor of 2 from I 1 , which explains the result N C −I 1 −I 2 = −6 in FIG. 8 for 2OE0.Depending on , three, two, one, or no eigenvalues can cross the branch cut, the number of net crossings being (I 1 + I 2 − N C )/2.As a side remark, we noted that eigenvalues could become degenerate for some frequencies (not shown).

V. CONCLUSIONS
We have extended the notion of Φ-derivable schemes to include partially self-consistent and non-selfconsistent approximations.For such partially Φ-derivable schemes.we have shown that by dressing the loops in Φ consistently, a gauge invariant approximation is obtained.This implies a fulfillment of the continuity equation at all times, which in turn implies number conservation.By making use of the non-equilibrium Green's function formalism, the resulting partially Φ-derivable scheme can  be applied in equilibrium as well as out of equilibrium, in finite and infinite systems.In equilibrium, partial Φderivability allows for deriving a generalized sum rule for the particle number, where we have stressed the importance of the correct analytical properties of the Green's functions.If the approximation is not gauge invariant, as is the case for many approximations used in the literature, we obtain another term I 2 in the sum rule.The rule was also applied to systems which allow for a partitioning into smaller subsystems, allowing for a generalized Friedel sum rule.The frequency-dependent embedding self-energies need to be taken into account in I 2 .
We have elucidated the known fact that PSD approximations automatically yield correct analytical properties for G R (ω).The converse is not true, which we have ex-emplified in several systems; a spectral function that can become negative can still come from a G R (ω) that has the correct analytical properties.We have proven that the validity of the Luttinger-Ward theorem and the Friedel sum rule depend crucially on the correct analytical properties of the Green's function.These considerations show that the simplest vertex corrections to GW can lead to a violation of common sum rules.Extra diagrams need to be added to obtain a PSD spectral function, thus guaranteeing the satisfaction of the sum rules.
While a given approximation can be shown to be number conserving without making use of partial Φderivability, we believe that this formalism is convenient and has several advantages.It is easy to see whether an approximation will be particle conserving or not, and also if the corresponding sum rules will be fulfilled.Combined with the knowledge of when the approximation has the correct analytical properties, these considerations could serve as a helpful tool in determining which diagrams to choose in a given situation.][9] Appendix B We made use of several properties of matrix logarithms in the derivations of the sum rules, which will be more explained here.The key point is that matrix functions are not defined according to their Taylor expansions, since this assumes a radius of convergence.However, we cannot use a spectral representation, since Green's functions are not Hermitian, and as such not guaranteed to be diagonalizable.
We assume that the Green's functions can be written as L × L matrices.We also assume that the matrices are invertible, which they have to be in order to define a matrix logarithm.We make use of the fact that a general complex matrix can be written in its Jordan normal form, Q = Ẑ Ĵ Ẑ−1 , where Ĵ is block diagonal with M blocks Ĵi corresponding to non-degenerate eigenvalues.Each m i ×m i block has the eigenvalue λ i of Q on the diagonal, with multiplicity m i .The superdiagonal of each block is equal to 1.
Eq.(B3) immediately yields Eq.( 62), where it is explicit that Q must be invertible for the matrix logarithm to exist, since then all λ k (ω) = 0.If we furthermore assume that Q(ω) and the transformation matrices are smooth functions of ω we obtain Eq.( 50) as These relations can be combined to derive Eq.( 61) as (B6) The integrals on the right hand side are called logarithmic integrals. 59If the λ k (ω):s do not cross the branch cut at (−∞, 0), unique primitives can be found for each eigenvalue, yielding Eq.( 61):

(B7)
As a side remark, we note that we can write a logarithmic integral using determinants, In this case, however, it is not possible to find a unique primitive in general, since the product of λ i (ω) can cross the branch cut even if the individual λ i (ω):s do not.This can be easily seen in the case of a noninteracting 3-site model.We consider the model in Sec.IV D, and put U = 0, = 0 and use the wideband limit approximation.The eigenvalues of − ĜR (ω) are , neither of which cross the real axis.The multiplication of them , does cross the real axis, at ω = ±1.Thus, if µ > −1, we cannot find a unique primitive logarithm in Eq.(B8).The situation is depicted in FIG.11.

Appendix C
In Sec.III, we made use of the spectral representation for ĜM (ζ), Eq. ( 23), and ΣM (ζ), Eq. ( 27).Using the Lehmann representation, one can show that the exact ĜM (ζ) and ΣM (ζ) can be written in this form.In this Appendix, we will show that any analytic operator with certain boundedness properties at infinity can be written in a spectral form.Furthermore, we will discuss how the spectral representation changes if the functions are non-analytic.We first consider scalar functions, and then generalize to operators.
Suppose that we are given a function f M (ζ), analytic except for single poles away from the real axis and decaying to zero at infinity, i.e. f (ζ) → 0 for |ζ| → ∞.
The poles have residue a k located at ξ k in the upper and lower half-plane.We define functions f R (ζ) and f Consider the integral along a closed half-circle in the upper half-plane.From the residue theorem, for Im ζ > 0, we get ) is bounded at infinity, the integral on the half-circle vanishes.Thus, we can write, for Im ζ > 0, where we have analytically continued f R (ω) = lim η→0 f R (ω + iη) in order to perform the integration on the real axis.In the same way, for ζ < 0, where we also analytically continued f A (ζ) to the real axis.Consider now the function This function satisfies Eq. (C1), which can be seen by considering f M (ζ) − f R/A (ζ), for Im ζ > 0 and Im ζ < 0, respectively.For Im ζ > 0, for example, By integrating around a small circle around a chosen ξ R m , which is in the upper half-plane, we obtain that Since a R m = 0, we must have that a R m = (a A m ) * and ξ R m = (ξ A m ) * .From Eq.(C10) it then immediately follows that A(ω) = A * (ω).We can then write f M (ζ), Eq.(C9), as By considering f R/A (ω) = lim η→0 f M (ω ± iη) in Eq.(C12), we find that ∓ 1 π Im f R/A (ω) = A(ω).While A(ω) is real, it is not necessarily positive and can have any sign.
From Eq. (C12) we can see the equivalence between functions analytic for Im ζ = 0 and functions that can be written as an integral over a spectral function.This gives justification to the use of the spectral representation for ĜM (ζ), Eq.( 23), and ΣM (ζ), Eq.( 27), for approximations that guarantee the correct analytical properties.The full justification follows when we consider operators below.When the spectral function is PSD, − ĜR (ζ) is a function that maps the upper half-plane to itself.This is a so-called Nevanlinna function.Such functions are considered in, e.g., 66 , where it is shown that such functions can be written in a spectral representation.However, the notion of spectral function can be generalized to measures, which means that also non-differentiable spectral functions, such as the ones consisting of delta functions, can be considered.
If f M (ζ) is not analytic away from the real axis, but has simple poles, the spectral function can still be defined via the real-frequency retarded or advanced Green's function.We can thus regard Eq.(C12) as a generalized spectral representation.If we would define the particle number as N A = µ −∞ A(ω), we would miss the contribution from the poles if they exist, and different definitions of particle number would yield different results.
We now generalize to operators F M (ζ), which we take to have finite dimension, and impose the condition [ F M (ζ * )] † = F M (ζ).We apply the same considerations as above to each matrix element F M ij (ζ).Thus, Eq.(C12) holds for each matrix element, yielding where Â † (ω) = Â(ω).Each matrix element can have their own set of poles.We now define Using this rewriting, we obtain the spectral representation in operator form We see, that even in the presence of poles, we can define the spectral function operator as F R (ω)− F A (ω) = −2πi Â(ω).We have thus shown that analytic operators, whose matrix elements decay to zero at infinity, can be written in a spectral representation.This discussion fully justifies the use of the spectral representation for ĜM (ζ), Eq.( 23), and ΣM (ζ), Eq.( 27), under the assumption of analyticity.

Appendix D
The particle number can be obtained from the diagonal of the Green's function, Eq. (12).In this Appendix, we will show that we can also obtain the particle number from the grand canonical potential Ω in equilibrium, and that both definitions agree for a partially Φ-derivable approximation.For the exact case, 35 as well as for any conserving approximation, 11 Ω is given by the Luttinger-Ward functional, which we write in the Klein form as and the contour is taken to be the Matsubara contour, γ d1G(1, 1 + ) = −i β 0 dτ dx G(x, τ ; x, τ + ).Note that the assumption of correct analytical properties is already implicit in Eq.(D1) since the logarithm has to be singlevalued (see also the discussion below Eq.( 62)).
A partially Φ-derivable scheme which is non-conserving will not have unambiguous total energies, since different energy functionals give different values. 13The Luttinger-Ward functional, however, is variational with respect to G, meaning that δΩ δG = 0 when evaluated at a G coming from the Dyson equation from a conserving approximation.][69][70] The Luttinger-Ward functional has recently been evaluated at finite temperature in extended systems 71 .
Following Baym, 11 we vary Eq.(D1) with respect to the chemical potential.For a partially Φ-derivable approximation, we have that Note that the derivation of Eq.(D9)only makes use of partial Φ-derivability, and that the assumption of gauge invariance is not needed.

REFERENCES
derivable.Dressing all Green's function lines in Φ with G gives back the original Φ[G]-functional as defined by Baym.Graphically, we obtain Σ[G, G 0 ] from cutting away a G-line from all Φ[G, G 0 ]-diagrams in every possible way.Examples of partially dressed Φ-functionals and the corresponding self-energies are shown in FIG. 1, FIG. 2 and FIG. 3.

FIG. 3 .
FIG. 3. Examples of partially dressed, gauge invariant, Φfunctionals (left) which fulfill particle number conservation, and their corresponding Σ = δΦ δG (right).a) A 4th order GW0 diagram.b) Dressing two loops yields two non-equivalent classes of diagrams.c) A 4th order diagram in the particleparticle T-matrix approximation.Note that we would get the same diagram if we dress only the inner loop.
where the indices k, l label the interacting sites.For the Anderson model, k = l = 1.The first term is the Hartree-Fock contribution, with n k (z 1 ) = −iG kk (z 1 , z + 1 ).The coefficient c in Eq.(75) is c = −1 for the 2nd order exchange diagram, FIG. 2, c = 2 for the ring diagram, FIG.