Approximations for many-body Green's functions: insights from the fundamental equations

Several widely used methods for the calculation of band structures and photo emission spectra, such as the GW approximation, rely on Many-Body Perturbation Theory. They can be obtained by iterating a set of functional differential equations relating the one-particle Green's function to its functional derivative with respect to an external perturbing potential. In the present work we apply a linear response expansion in order to obtain insights in various approximations for Green's functions calculations. The expansion leads to an effective screening, while keeping the effects of the interaction to all orders. In order to study various aspects of the resulting equations we discretize them, and retain only one point in space, spin, and time for all variables. Within this one-point model we obtain an explicit solution for the Green's function, which allows us to explore the structure of the general family of solutions, and to determine the specific solution that corresponds to the physical one. Moreover we analyze the performances of established approaches like $GW$ over the whole range of interaction strength, and we explore alternative approximations. Finally we link certain approximations for the exact solution to the corresponding manipulations for the differential equation which produce them. This link is crucial in view of a generalization of our findings to the real (multidimensional functional) case where only the differential equation is known.


I. INTRODUCTION
The one-particle Green's function (GF) 1-3 is a powerful quantity since it contains a wealth of information about a physical system, such as the expectation value of any single-particle operator over the ground state, the ground-state total energy, and the spectral function. In order to access this quantity one can start from its equation of motion: [4][5][6] i ∂ ∂t 1 − h(r 1 ) G(1, 2) where h(r 1 ) is the one-electron part of the many-body Hamiltonian, G 2 (1, 3; 2, 3 + ) is the two-body Green's function, and v(1 + , 3) is the Coulomb potential. The space, spin and time variables are all combined in (1) = (r 1 , σ 1 , t 1 ), and (1 + ) = (r 1 , σ 1 , t + 1 ) with t + 1 = t 1 + δ (δ → 0 + ). Eq. (1) can be manipulated in order to get a more practical expression, by introducing the non-interacting Green's function G 0 with In (3) G 0 determines the appropriate initial condition in time; note that the solutions of (1) and (2) are not unique. Moreover, in order to calculate G, the knowledge of G 2 is required (which in turns requires the knowledge of G 3 , and so on) 4,6 . In order to obtain a closed expression one can generalize G(1, 2) to G(1, 2; [ϕ]), where an external fictitious time-dependent potential ϕ is applied to the system. This allows one to express G 2 as 7 Note that in (4) all Green's functions are generalized to non-equilibrium since they depend on the perturbing potential. The equilibrium G and G 2 in (3) where V H (3) = −i d4 v(3, 4)G(4, 4 + ; [ϕ]) is the Hartree potential. Since the Hartree potential contains the Green's function, this term makes the equations nonlinear. We are interested in the solution of Eq. (5), for ϕ = 0. Its calculation would hence require the solution of a set of coupled, non-linear, first-order differential equations, which is clearly a non trivial task. Moreover one would need a new initial condition to completely define the desired solution of this differential equation, since the derivative δG δϕ has been introduced. Therefore usually another route is taken: one includes the functional derivative in (5) in the definition of a self-energy 4 which, inserted into Eq. (5) for ϕ = 0, gives: This is the Dyson equation for G, where Σ contains all the many-body effects (beyond the Hartree contribution) present in the system. Of course δG δϕ and therefore Σ are still not known and, in practice, Σ has to be approximated. A good starting point is obtained by reformulating the problem in terms of a coupled set of equations containing the one-particle Green's function, the polarizability P , the self-energy Σ, the screened Coulomb interaction W , and the vertex Γ. These equations are most often solved within the so called GW approximation (GW A) 8 , where the vertex Γ is set to unity, resulting in Σ ≈ iGW . Over the last two decades the GW method has become the tool of choice for calculations of quasi-particles (QP) band structures 9,10 of many materials and direct and inverse photo emission spectra (see e.g. Ref. [11][12][13][14] improving substantially over the results provided by static mean-field electronic structure methods. However the GW A suffers from some fundamental shortcomings (see e.g. Refs. [15][16][17][18] ) and with Σ being of first order in W , is not expected to describe strong correlation. Higher orders in W could be added by iterating the equations, but this is technically difficult, and there is no guarantee that results will quickly improve. It is therefore necessary to find guidelines. In the present work we go back to Eq. (5). Our aim is first, to obtain new insights about standard approximations by relating them more directly to the original equations. Second, we want to use Eq. (5) to explore alternative approximations. Finally, it might be interesting to concentrate directly on the set of coupled, non linear, first order functional differential equations for G, Eq (5), although it has been acknowledged that no "practical technique for solving such functional differential equation exactly" 4 is available. However, one may still hope that with new algorithms and the increase in computer power, numerical solutions might become accessible. The present work is hence also meant to explore strategies for, and possible problems of, such a route. In the following we resort to two approximations. First we linearize the set of equations by expanding V H in terms of ϕ. Second, we discretize Eq. (5) and consider in a first instance only one point for each space, spin, and time variable: we will call this latter approximation the "1-point model", as opposed to the full functional problem. The strategy underlying this procedure is the following: for the 1-point model, we can derive the exact explicit solution of the now algebraic differential equation, and solve the initial value problem. One can hence explore approximations to the full solution, which yields valuable insights in the performance of current approaches and suggestions for alternative ones. By determining which manipulations of the differential equation (DE) produce such approximate solutions, one obtains suggestions for analogous manipulations on the differential equation for the full functional problem, which opens the way to translate our model findings to realistic calculations.
The manuscript is structured as follows. In Section II we present the linearized differential equation which can be solved exactly within the 1-point framework. We discuss in particular the initial value problem, and how it can be overcome. In Section III we examine, in the 1-point framework, various common approximations to the solution of the DE: the iteration of Eq. (5) and approximations based on a Dyson equation, in particular different GW flavours. In Section IV we explore other routes to manipulate the initial DE and obtain approximate solutions. We finally give our conclusions and perspectives on future work in Section V.

II. THE SCREENED EQUATION IN A 1-POINT FRAMEWORK
Our first goal is to simplify the equations such that the main physics is retained, but manipulations become more straightforward. To this end we linearize the differential equation with an expansion of the Hartree potential to first order in the external potential ϕ, where G 0 H is a Hartree Green's function containing the Hartree potential at vanishing ϕ,φ = ǫ −1 ϕ is the renormalized external potential, and W = ǫ −1 v is the screened Coulomb potential with ǫ the dielectric function at ϕ = 0 19 . Concerning Eq. (9) three important remarks should be made: first, through the linearization the screened interaction W becomes the central quantity of the equation. This is justified by the physics of extended systems, where screening and plasmons are key concepts. Second, in principle W is the exact screened interaction, which of course is not known. One can however adopt two strategies: either W is considered to be an externally given quantity, obtained within a good approximation, e.g. from a time dependent density functional theory (TDDFT) calculation 20 ; or one could also recalculate W from G[φ] (see in the next section). In this work we will adopt the first strategy, which is illustrated in Fig. 1. Such a philosophy is rigorously justified. In particular, in the framework of the theory of functionals it is possible to pass from the Luttinger-Ward functional (given as functional of G, though indeed one should add the bare Coulomb interaction v as argument) to the socalled Ψ-functional, where v is replaced by W 21 . This explains for example why self-consistency in G only (and not in W) is sufficient to have a conserving GW approximation. Moreover, in practice this is the most current way of proceeding, corresponding e.g. in a GW calculation to the "best G, best W" approach: while the non-interacting G is taken e.g. to be the Kohn-Sham Green's function, W is calculated as accurately as possible, e.g. in the adiabatic local density approximation to TDDFT. Third by approximating the functional derivative δG δφ = −G δG −1 δφ G ≈ GG (which supposes the selfenergy to be independent ofφ) one obtains the Dyson equation for the one-particle Green's function in the GW approximation to the self-energy. The proof is given in App. A. This result shows that, even though the linearization procedure is an approximation, Eq. (9) is still a promising starting point to analyze the different flavours of the GW approximations and to go beyond. After linearizing, the next step consists in discretizing Eq. (9) and then in considering only one value for the space, spin, and time variables respectively (or equivalently concerning space and spin, in considering all Green's functions to be diagonal in a given basis): this is the 1-point model employed throughout the whole manuscript. The 1-point framework has already been used by other authors: in Refs 22, 23 Hedin's equations are combined in one single algebraic differential equation which is solved as a series expansion. This allows the authors to enumerate the diagrams for a certain order of expansion. Several expansion parameters are examined, such as vG 2 H , with v the bare Coulomb potential and G H the Hartree Green's function, vG 2 , with G the exact Green's function, W G 2 , with W the screened Coulomb potential, etc., which shows how at various orders of expansion the number of diagrams decreases by increasing the degree of renormalization. This is also the spirit behind the linearized equation (9), in which the natural ex-pansion parameter would be W G 2 H , where W is treated as an externally given interaction. The advantage of using the 1-point framework is that the equations become algebraic and thus the enumeration of diagrams is facilitated. In Ref. 24 a similar strategy as in Refs 22,23 is used to enumerate diagrams, focusing in particular on the asymptotic behavior of the counting numbers. Moreover Hedin's equations are transformed into a single first oder differential equation for the GF as a function of an interaction parameter, and an implicit solution is obtained. In order to fix the particular solution of this differential equation the initial condition G (v=0) = G 0 is used. Instead here we concentrate on (5), or better its linearized form (9), which is another differential equation for G, as a functional of an external potential. This choice allows us to (i) emphasize the essential physics contained in the screened Coulomb interaction W , (ii) discuss various aspects of the many-body problem in a clear and simple way, (iii) obtain an exact solution of the approximate equation that can be used as a benchmark. Moreover we believe that the 1-point version of Eq. (9) can be a natural starting point for a generalization to the full functional problem. While the equations are easier to manipulate, physical information is of course lost in the 1-point framework. In particular no poles (addition/removal energies) of the GF appear. However, the various aspects that will be explored in the following are intrinsically related to the structure of the equations, and hence exportable also to the full functional problem, in the same spirit as in Refs 22-24 .
A. The 1-point differential equation In the 1-point model Eq. (5) reduces to an algebraic, non-linear, first order differential equation where ϕ → z, G(1, 2; [ϕ]) → y(z), and G 0 (1, 2) → y 0 . Moreover iv(3 + , 4) → −v: this change of prefactor compensates for the time-or frequency integrations that have been dropped in the 1-point model and corresponds to a standard procedure in this context 22,24 . We can now linearize Eq. (10) in the same way as we did starting with Eq. (5) and obtaining Eq. (9). This yields Hence with respect to Eq. (9),φ → x, G 0 H (1, 2) → y H 0 , and iW (3 + , 5) → −u; the subscript u in y u highlights its u dependence. In the following, for simplicity of notation, we denote y 0 H by y 0 unless stated differently. In Appendix B we sketch the main steps to solve Eq. (11), based on the general ansatz y u (x) = A(x) · I(x). With the choice one obtains the equation and the general solution y u (x) reads where C(y 0 , u) is to be set by an initial condition. In the limit x → 0, which is the equilibrium solution we are looking for, Eq. (14) becomes Note that a similar ansatz can also be used for the full functional problem, namely G(1, 2) = d3 A(1, 3) · I(3, 2), in order to get a set of differential equations that are less complicated to manipulate than the original one.

B. The initial value problem
In general in order to set C(y 0 , u), y u (x) has to be known for a given potential x β (i.e. y u (x β ) = y β u ). However it is far from obvious to formulate such a condition in the realistic full functional case; this would indeed require the knowledge of the full interacting G for some given potential ϕ. Therefore the question is whether one can reformulate the condition in a simpler way in order to set C.
To answer this question we expand the exact solution for small values of u, obtaining: When u → 0 the one-body Green's function G has to reduce to the non interacting G 0 , in our framework this translates into: y u u→0 ≡ y 0 . Imposing this condition on which is satisfied if This result for C holds also for u = 0. Indeed it guarantees a non-divergent result for any non-vanishing potential x in (14). Moreover it reproduces the perturbative result which is obtained by iterating Eq. (11); for example the sixth iteration yields This is precisely the same series as the one appearing in Eq. (16) when C(u, y 0 ) is set to −1.

III. ANALYSIS OF COMMON METHODS TO CALCULATE THE ONE-BODY G
In the following we will analyze various established approximations for the calculation of the one particle G, using the knowledge of the exact solution.
A. Iteration of the DE Let us first iterate Eq. (11) starting from y (0) u (x) = y 0 , according to: For x = 0 the first two orders in u read and Eq. (19) for the third order. Results as a function of u are depicted in Fig. 2 together with the exact solution. Two observations can be made: i) very few terms are needed to obtain a good approximation to the exact solution in the small u regime; ii) for a given u = u n , the expansion diverges starting from an order n. The larger is u n , the smaller is n, which limits the precision that can be obtained. As previously mentioned, the iteration coincides with the expansion for small u of the exact solution.
Since the small u expansion is de facto the asymptotic expansion of the error function times an exponential (as can be seen in (16)) the divergent behaviour of the iteration in (20) is not surprising. Divergences of higher orders have been found in perturbation expansions for realistic systems, e.g. for orders higher than 3 in the Møeller-Plesset scheme 25,26 .

B. Self-energy based approximations
In this section the introduction of a self-energy Σ will be discussed along with its most common approximations. The Dyson-like form for Eq. (11), which is the equivalent of Eq. (7), reads: where a self-energy kernel has been defined as With for the vertex function, the selfenergy reads which is the equivalent of Σ = iGW Γ 8 . The Bethe-Salpeter equation for the vertex function Γ is then derived from (23) where y u = y u (x → 0). For x = 0 Eqs. (23), (25), and (27) correspond to a subset of the so-called Hedin's equations 8 , obtained by fixing W . A pictorial representation of this subset for a given W is given in Fig. 1.
In the following we will approximate the equations and the results will be compared to the exact solution of the differential equation, in order to obtain greater insight about these self-energy based techniques. From now on all quantities will hence be understood to be taken at x = 0.

G0W0 and self-consistency
Let us first look at different flavours of the GW approximation 8 . Setting Γ u (y u ) to unity, it follows that Σ u (y u ) = −uy u . Within the initial guess y (0) u = y 0 , one obtains a so-called G 0 W 0 self-energy Σ u = −uy 0 . 27 This is then employed in the Dyson equation (23) in order to get an improved y (1) u . To go beyond this first approximation one can iterate further within the GW approximation, i.e. keeping Γ u = 1. This corresponds to an iteration towards a GW 0 result, since G is iterated towards self-consistency but u, which represents the screened interaction, is kept fixed. We report here the expressions obtained for G 0 W 0 , i.e. the first solution of the Dyson equation, and for three successive loops We call this procedure iterative self-consistent scheme, in contrast with the direct self-consistent scheme where one solves directly the Dyson equation (23), for x = 0, with Σ u = −uy u . In this latter case one gets a second-order equation with two solutions Note that for the full functional problem one would find even more solutions, since a second-order equation has to be solved for each matrix element of G.
In order to choose the physical solution we Taylor expand the square root around u = 0, which leads to Since for u = 0 one has to obtain y u = y 0 , the physical solution is y u = 1 + 4uy 2 0 − 1 2uy 0 . In Fig. 3 we can appreciate how well these GW -based methods are performing against the exact solution in a wide u range.
Interestingly odd iterations quickly converge to the physical solution of the direct sc-GW 0 , while even iterations do also converge but at a slower pace: it can be shown that for u → ∞ their limit forms the sequence of rational numbers · · · which ultimately converges to 0. All the odd iterations have instead the exact large u limit (namely y u = 0 when u → ∞). One might use this property to improve the convergence of the series. An important question is now: does the result of the selfconsistent procedure depend on the starting point of the iteration? Here we have naturally chosen y (0) u = y 0 , but one might fear that this choice is simply lucky. Let us therefore look at the general iterative scheme which is obtained by solving the Dyson equation (23) for x = 0 By starting the iteration with a guess for y u on the right side one obtains For y This contains nothing else but the continued fraction representation for the square root corresponding to the physical solution y u = √ 1 + z − 1 2uy 0 where z = 4uy 2 0 . It converges for all values of the terminator y s . Therefore, this iteration will always converge to the physical solution. Does this mean that there is no risk of running into the unphysical solution? The answer is that it depends on the iterative scheme that is used, and not on the starting point. Look at the following way to re-write the Dyson equation (23): . If we iterate this equation by starting with some y (0) u = y s on the right-hand side we get which, with Eq. (37), is just the continued fraction representation for the unphysical solution y u = (− 1 + 4uy 2 0 − 1)/2uy 0 . In a way, this is good news: usually the iterative scheme adopted in the context of GW calculations is rather the first, safe one. Indeed it has been found empirically that such a scheme leads to self-consistent results independent of the starting point and in reasonably good agreement with experiments (see e.g. [28][29][30] ). However, when one goes beyond GW, higher order equations appear as we will see in the following. There are hence more and more solutions, and more and more ways to iterate the equations. In other words, there will be an increased danger to run into a wrong solution. One should keep this in mind when trying to add vertex corrections beyond GW .

Vertex corrections -First order Γ
We will now analyze the effects of a first order vertex correction which is obtained employing Σ u = −uy u in Eq. (27) 8,31 . Solving for Γ u gives Employing this vertex the self-energy (25) becomes Now two routes can be taken and either a G 0 W 0 Γ (1) (y 0 ) or a self-consistent GW 0 Γ (1) (y u ) calculation can be carried out. The first of the two is once more based on the initial guess for the Green's function y As it can be noticed from the result a cubic equation for the unknown y u had to be solved within this more sophisticated approach. Again the limit of vanishing interaction has been used to pick the physical solution.
In Fig. 4 we can directly compare the two types of vertex corrections. For small u values their performance is similar, however, in a wider u range (see inset), the G 0 W 0 Γ (1) scheme diverges from the exact solution and has the wrong asymptotic limit u → ∞: it hence behaves as the first iteration of the sc-GW 0 approach, which also exhibits the wrong large u limit. Fig. 4 also shows how the GW 0 Γ (1) scheme, for small u values, slightly improves over the sc-GW 0 . However, given the augmented complexity already at this first order of the correction (one could very well iterate further the equations for Γ and Σ and get higher order corrections), the benefits of employing vertex corrections are not obvious. Also note that interestingly, on the scale from u = 0 to u → ∞, the closest curve to the exact one is the sc-GW 0 one.

IV. EXPLORING OTHER APPROXIMATIONS FOR G
In this section we will explore alternative approximations to the exact solution of the 1-point DE and the corresponding manipulations of the initial differential equation producing them. Here we will report in particular approximations that might be eventually transposed to the full functional framework.

A. Continued fraction approximation
A well known approximation for the error function is its continued fraction representation 32 . The exact expression for y u (Eq. (15)) transforms into We will now show how one can obtain Eq. (45) starting simply from the initial DE in Eq. (11), equivalent to (9), without any information about its exact solution. Beginning with Eq. (11) and taking successively higher order derivatives of the equation, one obtains : and so on. Neglecting derivatives e.g. from the 4 th order on and then setting x = 0, this truncation allows us to solve all the above equations, beginning with Eq.
which is precisely the result obtained by approximating the exact solution with a continued fraction expression for the error function ( Eq. (45) ). We will name this manipulation limited order differential equation. In Fig. 5 we compare the different orders of this approximation to the exact expression for y u . The approximation gets rapidly closer and closer to the exact solution by including higher derivatives. However, also for this continued fraction, odd and even orders, converge towards the exact result with a different speed. In analogy with the continued fraction of Eq. (37), even iterations have the correct large u limit, while the odd ones don't, although they do eventually approach it for a very large number of steps. We notice that the above continued fraction converges slower than the one arising from the sc-GW 0 ; however, the former will eventually converge towards the exact solution, whereas the latter only to the self-consistent GW 0 result. It is therefore interesting to note that such a procedure can in principle be used also in the full functional framework (see related manipulations e.g. in 6,33 ), where the functional differential equation can be differentiated to an arbitrary order and the corresponding approximated G obtained. For example, differentiating Eq. (9) with respect to the external potentialφ one gets Truncating the highest order derivative δ 2 G δφ 2 and solving for ϕ = 0 (which means alsoφ = 0) gives: which reinserted in Eq. (9) yields: Like in the 1-point model, this first step simply provides the one-particle GF in the G 0 H W 0 approximation to the self-energy. One can go further: differentiating Eq. (50) with respect toφ and neglecting the third order derivative δ 3 G δφ 3 yields: which is a four-point quantity of a similar complexity as the Bethe-Salpeter equation 5 . Indeed in the full functional problem the equations become quite involved since terms like uy 2 0 correspond to large matrices. However the approach doesn't require self-consistency. This might turn out to be a significant advantage, compared to vertex corrections to Σ, as we have discussed in the previous subsection concerning self-consistency. More details about the derivation are given in App. C.

B. Large u expansions
Perturbation theory usually deals with weak interactions, hence the small u limit. However, it is also very interesting to examine the large u limit for several reasons: i) this is the regime of strong correlation, where current approximations exhibit failures; ii) the large u expansion of the exact solution gives a convergent series (being a product of two convergent Taylor expansions, one for the exponential and the other one for the error function) and one can, for instance, obtain a better approximation to the exact solution by adding higher order terms (which instead does not improve the result for the small u expansion of the solution); iii) excellent approximations for the exact solution are Padé approximants 34 , which have to be constructed using both the small and large u limit. In this subsection we will present two possible routes to approach this limit: the first is a straightforward large u expansion of the exact solution for y u , while the second combines the latter with the large u expansion for the Dyson equation.
one obtains for the different orders of the full solution (62) Fig. 6 shows how these different expansions perform versus the exact result. Overall their behaviour is very good for large u and few orders are sufficient to get a good approximation over a wide u range (which is our ultimate goal), however for u = 0 they all diverge.

Large u expansion for yu and for the Dyson equation
When u gets larger, also Σ u increases. This implies that, using the Dyson equation for the one-particle Green's function y u = y −1 0 − Σ u −1 one could expand y u as Hence to lowest order y u ≈ −Σ −1 u or This simple relation allows us to use the large u expansion of the exact solution for y u to approximate Σ u for large u; we can then use this approximate Σ u in the Dyson equation to recalculate y u . For example, using the lowest order of the large u expansion of the exact y u one gets the following self-energy: which reinserted in the Dyson equation y u = (y −1 0 − Σ u ) −1 gives: In Fig. 7 the performance of this approximation for y u is plotted against two orders of the straightforward large u expansion for the Green's function, G 0 W 0 and the exact solution. The "large Σ" approach shows an overall good agreement (generally better than G 0 W 0 ) with the exact solution and it has the desirable property of being exact in the small and large u limits, mending the divergence of all orders of the straightforward expansion for y u . At higher orders of the approximation this property remains true, although undesired poles appear. In conclusion the methodology is promising and worthwhile to be explored further. The main difficulty is that in the framework of a large u expansion, without knowing the exact solution, one would not straightforwardly know how to set the constant C, i.e. how to pick the physical solution: this issue requires further analysis.

C. Self-consistent calculations of the Hartree Green's function and of the screened interaction
In the above discussions we have treated the Hartree Green's function and the screened interaction as externally given quantities. This is justified by the fact that realistic calculations are most often following such a pragmatic ansatz. In principle these quantities should be part of Hedin's self-consistent cycle. A fully self-consistent treatment, in the full functional framework, is today out of reach. In the 1-point model, however, it is possible to go beyond this limitation and indeed, the implicit solution of Hedin's equation that has been achieved in the work of 24 contains all quantities calculated on the same footing. Also in the linearized version, that is employed in the present work, one can obtain the Hartree Green's function and the screened potential consistently from the equations, as we will discuss in the following. Let us first turn to the Hartree Green's function y 0 H . In terms of the truly non-interacting Green's function y 0 it reads in other words, it depends (through the density) on the solution y u at vanishing external potential. In a selfconsistent scheme this y 0 H should then replace y 0 in the solution Eq. (15), which leads to an implicit equation for y u . For a self-consistent treatment of the screened interaction we can use of the fact that the 1-point differential equation can be solved for dy u dx , and insert the result into the expression for the screened interaction u in terms of the bare v, which reads u = v + v dy u dx v. Two routes can be taken. The first one is based on the linearized equation (11) where the interaction is already screened from the very beginning. This leads to a quadratic equation for u, with two solutions The physical solution is the one of the positive square root, since it approaches the bare v in the limit of vanishing interaction, hence vanishing screening. The second route consists in calculating dy u dx from the initial equation (10), where the bare y 0 and the interaction v appear. This yields In both cases, the solution for u should be used into Eq. (15), which again makes the expression for the GF implicit. One may argue about which of the two ways to calculate u self-consistently is more adequate. In a realistic calculation one would probably use the former approach, in an iterative way: after calculating the GF as a functional of the external potential for a given initial interaction in the linearized DE, one would recalculate the W from the functional derivative, and so on. Whatever choice, however, does not influence the main conclusions that can be drawn from the above considerations. Specifically: i) a self-consistent calculation leads to an implicit solution (like in the work of 24 ) which however would not be identical to theirs because of our linearization procedure; ii) the behaviour for the small interaction limit is unchanged by the self-consistent treatment, as one can verify from equations (67),(68) and (69); this means in particular that the constant C is chosen in the same way as before. iii) Finally also the discussion about the limit of large interaction is unchanged: by making the ansatz that to lowest order y u ∝ 1 √ u one finds consistency.
Altogether, this shows that the linearization of the equations does not imply necessarily that one has to treat the Hartree Green's function and the screened interaction as externally given quantities. It also shows that a more refined, self-consistent treatment does not change the overall behaviour of the solution.

V. CONCLUSIONS AND OUTLOOK
In this paper we explore several aspects of the set of first order nonlinear coupled differential equations which are conventionally solved perturbatively in order to calculate the one-particle Green's function. After the linearization of the Hartree potential with respect to the external one, we employ a 1-point model where the set of -now linear-differential equations reduces to a 1 st order algebraic differential equation, that can be solved exactly. This provides insights into the structure of the general family of solutions, and on how to determine which amongst them corresponds to the physical one. Within the model we study the performance of established approaches over the whole range of interaction strength: we find that iterations towards self-consistency in the GW scheme sensibly improve on the one-shot (G 0 W 0 ) calculation and that including first order vertex corrections improves the self-consistent GW 0 results only slightly and only for small u. We also find that in case of selfconsistent GW 0 two solutions are possible, of which only one is physical. We show that the standard iterative scheme will always converge to the physical solution, although other schemes may yield different results. This is an important finding: when going beyond GW both the number of possible solutions for the Green's function and the number of possible ways to iterate the equations increase, creating a danger to run into a wrong solution. We finally explore other approximations to the exact solution that might be transposed to the full functional framework, namely a continued fraction approximation and the expansion for large interaction, and we relate these approximations to the corresponding manipulations of the differential equation which produce them. These links are crucial to prepare a generalization of the approach to the full functional framework. full functional problem. A general ansatz for the structure of y u (x) is: where the only restriction is that A and I are not zero. Substituting the ansatz in the DE (11) gives: The idea is now to solve two separate, simpler with respect to the initial one, DEs for A(x) and I(x). Putting together the left-hand side and the second and third terms of the right-hand side of Eq. (B2) one obtains: We can choose the solution which will then determine I(x). One is now left with the equation for I(x) reading Plugging in the expression for A(x) previously obtained and integrating on both sides one obtains: The integral on the right-hand side is: has been made, and the lower limit of the last integral has been chosen to be zero, which requires to set a constant C(u, y 0 ). Hence The exact solution y u (x) = A(x) · I(x) is given in Eq.  (19)). All the three orders are close to the exact solution for small u values, whereas when a given order of the series starts to diverge, the lower orders of the expansion reproduce the exact results better. For each curve C(u, y0) = −1, and we arbitrarily set y0 = 1.  (15)) and different flavours of the GW approximation. In general the self-energy based approximations perform better than the iteration of the DE shown in Fig. 2. In the main panel the self-consistent GW0 (black stars, Eq. (32)) is the best approximation to the exact result. Iterations starting from G = G0 converge towards the self-consistent result (the 2 nd iteration is represented with light blue triangles, the 3 rd with green circles and the 4 th with grey empty triangles). However, analyzing a larger u range (inset), one observes that odd iterations approach the exact u = ∞ limit, while the even ones don't seem to. It can be shown that they also do, however in a very slow fashion and according to the following sequence y and sc-GW0 (black stars, Eq. (32)) is shown. In this range of u, adding a vertex correction, no matter if within a selfconsistent scheme or not, improves over the simpler selfconsistent GW0. However, analyzing a wide u range (inset, semi-logarithmic plot), gives a different perspective: the first iteration of G0W0Γ (1) clearly exhibits the wrong u → ∞ limit and the sc-GW0 scheme becomes the closest approximation to the exact result.  (57-59)). We also report the G0W0 result (blue dots, Eq. (28)) as an example of a small u expansion. Over a wide u range the large u expansions are very satisfactory.  66)). We observe that the latter approximation performs extremely well over the range of interaction examined, being even exact both in the large and small u limits.