{'fraction_non_alphanumeric': 0.03431344561804168, 'fraction_numerical': 0.009734513274336283, 'mean_word_length': 4.416048237476809, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Many-body perturbation theory and non-perturbative approaches: the screened interaction as key ingredient

Many-body perturbation theory is often formulated in terms of an expansion in the dressed instead of the bare Green's function, and in the screened instead of the bare Coulomb interaction. However, screening can be calculated on different levels of approximation, and it is important to define what is the most appropriate choice. We explore this question by studying a zero-dimensional model (so called 'one-point model') that retains the structure of the full equations. We study both linear and non-linear response approximations to the screening. We find that an expansion in terms of the screening in the random phase approximation is the most promising way for an application in real systems. Moreover, by making use of the nonperturbative features of the Kadanoff-Baym equation for the one-body Green's function, we obtain an approximate solution in our model that is very promising, although its applicability to real systems has still to be explored.


I. INTRODUCTION
In an interacting electron system, complete knowledge of the many-body wavefunction allows one to calculate all expectation values and therefore to determine all observables. However, in most cases it is neither possible nor desirable to calculate the full wavefunction, as in the calculation of the observables much details of the wavefunction are integrated out. In particular, expectation values of one-body operators and the total energy can be obtained in a straightforward way from the one-body Green's function G(rσt, r σ t ), which is itself an expectation value of an electron creation and an annihilation field operator. This Green's function can be expressed as a series whose terms only include its non-interacting counterpart G 0 (rσt, r σ t ) and the Coulomb potential v c (r, r ).
Such a series, however, can be badly behaved, as, for instance, in the case of extended systems in the thermodynamic limit. In these situations, one can still extract sensible estimates of the Green's function by summing up an infinite subset of terms of the series. Cornerstone of this approach is the Dyson equation, which in a simplified notation reads G = G 0 +G 0 ΣG; here the focus is shifted on the self-energy Σ, whose perturbative expansion is still expressed in terms of G 0 and v c only. Any finite truncation of the series for Σ leads to a resummation of an infinite number of terms in the series of G, and only a few terms of the series of Σ are usually sufficient to greatly improve on truncations of the series of G itself. Still, in order to obtain satisfactory spectra, at least in an extended system it is not sufficient to keep for the self-energy only some diagrams of low order in G 0 and v c . In particular, it has rapidly been realized that one should resum diagrams that dress the bare Coulomb interaction, since the effective interaction between charges is usually strongly screened.
In extended systems at zero temperature, the state of the art approach for the calculation of G, and hence of electron addition and removal spectra, is the so-called GW approximation. In a seminal paper, Hedin formalized this idea by introducing an appropriate screened interaction W and expressing the exchange-correlation contribution Σ xc to the self-energy as a functional of G and W , rather than G 0 and v c . He derived the expression using Schwinger's functional derivative idea, where all the information about exchange and correlation effects on the one-body and higher order Green's functions are contained in a functional differential equation that involves the response of the system to a fictitious external potential 1,2 . From the differential equation, a set of five equations are derived that are known as 'Hedin's equations' 3 . They contain the one-body Green's function G, the self-energy Σ, the irreducible polarizability P which gives rise to the screened Coulomb interaction W , and a vertex function Γ which stems from variations of the self-energy. The GW approximation arises as the first step of an iterative procedure that formally solves Hedin's equations. This first approximation corresponds to setting the vertex Γ = 1. The GW approximation has encountered large success, in particular concerning the calculation of band gaps. However, in cases of failure, such as the description of Mott insulators which come out metallic in GW 4 , or the absence of satellites related to excitations that are not due to the formation of electron-hole pairs 5 , a straightforward further iteration of Hedin's equations has not yet proved to be successful, and, for many applications, it is not even feasible.
The power of the GW approximation is generally attributed to the belief that the screened Coulomb interaction W is weaker than the bare v c , and that therefore a perturbation expansion in W should be more powerful than an expansion in v c . This, however, is not necessarily true. To see the point, it is sufficient to consider the screening due to a single electron in some potential. In this case the density-density response function χ is simply the non-interacting response function χ 0 . The inverse dielectric function −1 = 1 + v c χ , which screens the Coulomb interaction via the relation W = −1 v c , becomes −1 = 1 + v c χ 0 . Since χ 0 is negative, W is indeed often smaller than v c . However, for large enough (negative) χ 0 compared to 1/v c the inverse dielectric function and W become negative, and their absolute value can be arbitrarily large. In this case the argument in favor of a perturbation expansion in terms of W breaks down. Instead, this scenario never happens when W is calculated in the random phase approximation (RPA), where only variations of the Hartree potential screen the interaction. In the RPA the dielectric function is = 1 − v c χ 0 , which is always positive and larger than 1, such that the resulting screened interaction W 0 < v c . This suggests that an expansion in terms of the RPA W 0 might be more powerful than an expansion in W . This indicates a first route to obtain improved expressions, which we follow in the present work.
Instead of searching for a perturbation expansion for the self-energy, one may also go back to the starting point, a functional differential equation. Here we use the Kadanoff-Baym equation (KBE) 6 , which can be reformulated such that screening appears explicitly. The GW approximation is obtained from this equation as a linear-response approximation in conjunction with an approximation on variations of the Green's function with respect to the total classical potential. Also the widely used second-order cumulant approximation can be derived from this equation, again using the linear response approximation, but introducing a decoupling approximation between different states instead of the approximation on the variations of G. In order to go beyond both approximations, it is therefore in particular interesting to investigate contributions beyond linear response.
In the present work, both ideas are explored using a simple model, the "One-Point model" (OPM). Such a mathematical toy model, which has already been used in a similar context 7-9 , represents the 0-dimensional version of a generalized Kadanoff-Baym equation. The model captures many features of the original, multi-dimensional case; it incorporates, for instance, the failure of the skeleton series that was only recently found for some multi-dimensional cases [10][11][12] . The model is exactly solvable, and it reflects many features of findings for real systems, such as the quality of various approximations 13 . Therefore, it constitutes an established first step to test new approximations.
The paper is organized as follows. Section II will be devoted to recall the reader the theoretical framework; our two approaches will be motivated and presented in Section III and IV, respectively; results of tests on the OPM are presented in Section V; feasibility and computational cost of the approaches for real systems will be discussed in Section VI where conclusions will also be drawn.

II. THE EQUATION OF MOTION AND THE FUNCTIONAL APPROACH
In a non-interacting system, one can determine the one-body Green's function G 0 from its equation of motion (EOM), which is a differential equation containing the first-order time derivative of the Green's function and the single-particle Hamiltonian. Together with a boundary condition in time, this fully determines G 0 . In an interacting system, the EOM for the one-body Green's function G contains additional contributions due to the interaction: besides the Hartree potential, a two-body Green's function appears, which contains the information about exchange and correlation between the interacting particles. In turn, the EOM of the two-body Green's function contains the three-body Green's function, and one finds an infinite chain of equations 4 .
The equations can be cast into a more compact form by using the fact that higher-order Green's functions can be expressed as variations of the one-body G with respect to an external potential ϕ. This allows one to write the EOM for the one-body Green's function G as where v c is the Coulomb potential, G 0 (1, 2) = G(1, 2; ϕ = 0) Here, we use (1) = (r 1 , σ 1 , t 1 ) as a short-hand notation to combine the space, spin, and time variables. Moreover, (1 + ) = (r 1 , σ 1 , t + 1 ), where t + 1 = t 1 + η with η → 0 + . In (1) the generalized one-body Green's function G(1, 2; [ϕ]) is a functional of a time-dependent external potential ϕ(1). The equilibrium one-body Green's function of interest is retrieved in the limit of vanishing ϕ(1). Without the last term on the right-hand side of Eq. (1), this equation would be the Dyson equation in the Hartree approximation. All the many-body effects (exchange and correlation (xc)) beyond the Hartree potential are contained in the derivative term.
As pointed out by Baym and Kadanoff in 6 , there is no known way to solve this equation. In particular, such a non-linear multi-dimensional functional integro-differential equation can have many solutions, and it may be difficult to chose which one is physical. Usually, the problem is circumvented by using some kind of iterative approach, starting from the non-interacting or the Hartree solution. A straightforward iteration of Eq. (1) leads to series of terms of increasing order in the interaction and corresponds to a perturbation expansion of G. However, such a series has often bad convergence properties (see also Sec. V for an illustration). Therefore a common way to go is to transform Eq. (1) into a Dyson equation with the exchange-correlation self-energy Here and in the following, we denote integrals by f (n)g(n) ≡ dn f (n)g(n). Once the expression for the self-energy is established, it is enough to solve Eq. 3 for ϕ → 0 in order to obtain the equilibrium Green's function. Of course, since G is unknown so is Σ xc : with the transformation of the EOM into a Dyson equation, the goal becomes to find accurate approximations to the self-energy. The advantage of the Dyson equation is the fact that even a low-order approximation to the self-energy creates in the Green's function terms to infinite order in the interaction, so one may hope that approximating Σ xc is easier than approximating G itself.
One strategy to find approximations for Σ xc has been formalized by Hedin 14 , based on the insight that, at least for extended systems, screening of the interaction plays a key role. The screening is due to the self-consistent reaction of the electron system to the variation of the external potential. This self-consistency is automatically taken into account when the variation of the total, i.e. external plus system-internal, potential is considered. Hedin's formulation uses the classical part of the total potential, ϕ cl = ϕ + V H , to express the exact exchange-correlation self-energy (4) as From the Dyson equation (3), the fifth equation in Hedin's set is derived for the irreducible vertex functionΓ The screened Coulomb interaction W is linked to the density-density response function , and it can also be expressed in terms of the inverse dielectric function By settingΓ ≈ 1, the equations reduce to a closed set of algebraic equations for G, Σ, W and P , leading to the so called GW approximation (GWA) for the self-energy and the Random Phase Approximation (RPA) for P . The equations can now be solved self-consistently. Often, however, some mean-field Green's function is used instead of the self-consistent G to build the self-energy. The mean-field Green's function is usually referred to as G 0 , although it contains some interaction through, e.g., the KS potential or the Fock exchange. This approach is usually called G 0 W 0 . For later reference we highlight the fact that W is calculated in the RPA by using the notation where the RPA functional is Such a functional generalizes the functional form of W obtained from the Random Phase Approximation for P , whereΓ = 1 implies P ≈ −iGG and hence W ≈ The G 0 W 0 scheme has led to good results in many cases, but it suffers from the fact that these results strongly depend on the choice of the mean field in G 0 . For this reason, partially self-consistent schemes have been developed and led to important progress. Technically, even fully self-consistent GW calculations are nowadays feasible 15 . However, it is not clear whether full self-consistency on the GW level is desirable; in particular, it has been pointed out that the resulting χ does not fulfill the f -sum rule 16 . Clearly, vertex corrections, i.e. corrections to the self-energy of higher order in W , are needed to improve this situation.

III. PERTURBATION THEORY IN W0
The success of the GW approximation may be explained by the intuitive picture according to which the screened potential W represents a better parameter than the bare Coulomb interaction v c for a perturbative expansion of Σ xc . This presumes that |W | < v c . Although often true, this is simplistic, and as explained in the introduction |W | can in principle become arbitrarily large. On the other hand, P RP A (ω = 0; [G 0 ]) built with any mean-field G 0 is always negative, which means that W RPA [G 0 ] is always positive and smaller than v c . Moreover, the corresponding density-density response function χ RP A [G 0 ] is well behaved; in particular, it fulfills the f -sum rule.
Based on these considerations, we propose to consider the perturbative expansion of Σ xc in G 0 and W RPA [G 0 ]. This is formally constructed in the following way. First, we use the definition of We then use the well established perturbation expansion of Σ xc in terms of v c and insert the above expression for v c to rewrite the series in terms of W RPA [G 0 ] and, finally, we group terms of the same order in W RPA [G 0 ]. By defining W 0 ≡ W RPA [G 0 ] we can schematically represent the procedure as The first order yields G 0 W 0 , while the second order modifies the Second Born approximation. More in detail, the first order contribution to the self-energy is and the second order contributions read where ) are all second order diagrams of the expansion in v c evaluated with W 0 instead of v c . Note the subtraction of the last diagram, which is necessary in order to avoid double counting of the bubble diagrams.

IV. NONPERTURBATIVE LINEAR RESPONSE
The previous section designs a way to use a linear response ingredient, namely W 0 , in order to build more powerful perturbation series. In the present section we will discuss a way to make non-linear response appear. Since the expressions become quickly very clumsy, we suppress in the presentation all arguments except for the functional dependence on the fictitious external potential, which is crucial for the derivations.
Starting from (1), we focus on the self-energy We assume that V H [ϕ] (defined in (2)) is Taylor expandable: where n(1) ≡ G(1, 1 + ; [0]) is the charge density. As first approximation, we truncate the above series to first order which leads to This rewriting suggests to renormalize (or screen) the external potential as Since (19) implies If one now approximates the functional derivative as one obtainsΣ Since the derivative part has been approximated, this equation is no longer a functional-differential equation but a Dyson equation, which can be solved directly at ϕ =φ = 0. The equilibrium self-energy is then obtained as Here the RPA screened interaction appears, since These are, respectively, the self-consistent GW approximation with W ≈ W RPA [G], and the so-called G 0 W 0 approximation, where the expressions are evaluated with non-interacting Green's functions; in practice, most often some mean-field solution such as a Kohn-Sham Green's function is used for G 0 . The above derivation of the G 0 W 0 approximation relies on two ordered approximations: (17) and (25), that are tied to each other via the definition of the screened perturbing potentialφ. The first one makes the problem linear, while the second one removes the differential character and leads to an algebraic equation. If one would apply (25) directly on (15), without going through (17) and the definition ofφ, one would get to which is simply the Hartree-Fock approximation.
The derivation seems to suggest that a better approximation may come from improving on (17), namely taking more orders of the Taylor expansion of V H [ϕ]. However, solving the Dyson equation at equilibrium kills all possible improvements. Such a mechanism is clear already when keeping the second order in (16). Let us first write Now we introduce the screened perturbing potential (in the even shorter notation for which functional derivatives for which ]. This implies that which, in combination with givesΣ which at equilibrium is justΣ It follows that the only way to include non-linear screening effects is to go beyond the second approximation (25,35). This suggests that in order to go beyond the GWA one should keep the functional derivative, and attempt to solve the linearized functional-differential equation itself, as it has been suggested in 9 and pursued in 17,18 in the context of cumulant approximations. According to the order of the truncation of (16), we get different approximate KBEs, which we shall collectively refer to as 'hierarchy of linearized KBEs', since these differential equations are linear in the variable G[ϕ]. The coefficients of the Taylor expansion can in principle be estimated within Time-dependent Density Functional Theory, or evaluated self-consistently once the differential equation is solved. 19 The hope is that one recovers the original KBE by climbing the hierarchy of linearized KBEs, namely by considering more and more terms of (16), and that the resulting hierarchy of physical solutions to these equations smoothly converges to the exact Green's function.

V. TESTS ON THE OPM
Solving the functional-differential equation (1), even with the simplest approximations for V H [ϕ], is a tough problem. First attempts in this direction were carried out by Lani and collaborators 9 , who found it fruitful to study the equation in its zero-dimensional version, in which case the equation becomes a simple differential equation for a scalar function of one variable. Later on, Berger and collaborators 13 generalized such a differential equation in a way to describe more physics by keeping the same simple framework. The resulting equation, usually referred to as the 'one-point model' (OPM), encodes highly nontrivial features of the functionals discussed in Section II, which otherwise could be studied only in their perturbative regime and only for a very limited order of expansion, due to the factorially increasing number of terms produced by perturbation theory. Despite its simplicity, the OPM has proved to be a surprisingly reliable benchmark for qualitative analysis of approximation schemes, as proved by the analyses performed by Berger et al. and, later, Stan et al., who also showed that the model can be used to study pathologies of Perturbation Theory such as the misconvergence of the skeleton series 10,12 . This makes the model particularly appropriate for the study of the nonperturbative features of the equation of motion and its approximations, discussed in the previous section. Moreover, one can also use it to test perturbative schemes like the one based on an expansion in the RPA screened interaction, with the possibility of exploring at glance a large number of terms of a series and getting qualitative insights concerning its convergence.
In the next subsections, we shall recall the relevant ingredients of the OPM (V A); show the performance of the OPM equivalent of GWA (V B), which will serve as benchmark for our suggested approximations; derive and study the series in the equivalent of W RPA [G 0 ] and compare it with the series in the exact W (V C); solve the hierarchy of the OPM equivalent of the linearized KBE equations (V D); and, finally, compare all approaches (V E).

A. The One-Point Model
We here summarize the relevant features of the OPM. The reader interested to more details is referred to the work of Berger et al. 13 . At the basis of the model we find the differential equation: with y 0 (x) ≡ y 0 0 /(1 − y 0 0 x), x ∈ R, λ > 0, v ≥ 0, y 0 0 > 0. For λ = 1, the above equation corresponds to the zero-dimensional version of the equation of motion (1), provided that For this specific value of λ, the exact problem is characterized by the trivial physical solution y(x) = y 0 (x). A much reacher result is obtained by considering λ = 1. Here we take λ = 1/2, namely The choice λ = 1 reflects the fact that in the real system Hartee and exchange contributions do not cancel perfectly, whereas this would be the case in the model for λ = 1. In this sense, λ = 1/2 is linked to the prefactor 1/2 that is often used in real systems for exchange terms when spin is not explicitly considered.
In analogy with the original problem, one can use the equation to get the perturbative expansion of y(x) in v without actually solving the differential problem. In fact, the resulting series: can be summed up to the function which is the analytic continuation of the series (41). Since (42) is also solution to (40), we identify this as the physical solution of the problem, all other solutions to (40) being labeled as unphysical ones. Such a function is a model of the functional that in the standard many-body problem realizes the map ϕ → G by analytic continuation of the series produced by Perturbation Theory. The Green's function at equilibrium, namely G[ϕ = 0], is then represented by which is plotted in fig. 1. The function y 0 (x), to which y p (x) reduces for v → 0, is then identified as the out-of-equilibrium non-interacting Green's function, the corresponding version at equilibrium being simply y 0 (0) = y 0 0 . Provided with these functions, we can define the physical self-energy Σ, the Hartree self-energy Σ H , the exchange-correlation self-energy Σ xc , the screened potential w and the RPA function(al) w rpa (f ) as In the rest of the section we shall consider Σ H as given once for all. This is different from the work of Berger et al. where Σ H was calculated self-consistently. We make this choice since it will allow us to better represent standard many-body techniques, and G 0 W 0 in particular, in which Σ H is usually a given input coming from a DFT calculation. Moreover it allows us to discuss the two proposed schemes in a consistent way.

B. The GWA
To assess the quality of the proposed approximations, we shall compare them with the OPM equivalent of the GWA. A thorough analysis of such an approximation was worked out by Berger et al. 13 , who considered a self-consistently calculated Σ H . As pointed out above, here we make an analogous analysis, but with a Σ H which is supposed to be given. The self-consistent GW approximation corresponds to solving the set of equations This set of equations admits 3 solutions, of which only the following reduces to y 0 0 in the noninteracting limit: For real systems a direct solution of the GW equations is, in general, not feasible. Therefore, in practice, they are usually solved iteratively. Since the GW equations can be rewritten in several ways, various iterative schemes are possible 13 . Interestingly, the above solution can be obtained by the following iterative procedure: starting from y 1 = y 0 0 , which is the standard iterative scheme used in GW calculations for real systems. We note that the above iterative scheme converges to the physical solution for all v thanks to the fact that we consider Σ H fixed. Otherwise the above iterative scheme does not converge for large v 13 .
Since we can also calculate the exact screened potential W , it is worthwhile to explore also a GW self-energy built with the exact W , which could in principle be obtained from TDDFT. In this case, we solve the Dyson equation Eqn.(45c) with the exact screened potential W = w of (44f). Of the two solutions, the only one that goes to zero in the limit v → 0 is In figure 2, the OPM equivalent of the self-consistent GW approximation and self-consistent GW with W exact are shown in comparison with the equivalent of G 0 W 0 . In analogy with what is found in realistic calculations, selfconsistent GW does not improve systematically on G 0 W 0 and is clearly worse for larger values of the interaction. The function y p (0), which represent the equilibrium interacting Green's function, is written in terms of two parameters: y 0 0 , which represents the equilibrium non-interacting Green's function, and v, which parametrizes the interaction. Following Section III, we can rewrite it in terms of y 0 0 and w rpa ≡ W rpa (y 0 0 ) by simply using the definition of w rpa to express v, which leads to In figure 3, we compare approximations to y p (0) resulting from the straight perturbative expansion (51) and from the solution to the Dyson equation with self-energy approximated by the expansion in (52). All hierarchies of approximations smoothly converge to the correct result. It is worth noticing that, even though the leading order expansion of Σ xc gives a better approximation for large values of v compared with the other leading orders, the straight expansion of y p (0) converges much faster than the others. Finally, we would like to reply to the question: is w rpa a better perturbative parameter than the exact W ? To do that, one has to invert the expression (44g) to get the function corresponding to the map W → v. There are three solutions to the equation for v, of which only one goes to 0 when W goes to 0, namely v = Truncations of the above series are plotted in figure 5. Few orders of the series are representative of the overall bad behaviour of the series. This suggests that the exact W does not represent a good parameter for a perturbation expansion. Indeed, a glance at the comparison of the RPA and the exact inverse dielectric function in fig. 4 shows that while the RPA value is steadily decreasing with increasing interaction, the tendency is inverted starting from a critical value of the interaction in the exact case, making the screened interaction too large to be a useful parameter for perturbation theory.

D. Solutions to the Hierarchy of linearized KBE
We will now move on to the second route beyond the standard approximations, namely the expansion of the Hartree potential introduced in Section IV. To this end, we define the the 'perturbed' Hartree potential as V H (x) ≡ −vy(x), which leads to rewrite (40) as and consider the hierarchy of approximations generated by truncating its Taylor expansion in the perturbing potential The leading term V 0 simply corresponds to the Σ H previously defined and assumed to be known, while the other terms involve derivatives of y(x) evaluated at zero: V 1 = −vy (0), V 2 = −vy (0)/2, etc. which are not known a priori. Just like for V 0 , we could consider them as given inputs, since for real systems they can in principle be obtained within the framework of TDDFT. In practice, however, we can only get approximations to their value, on which we have no control. We therefore prefer to calculate them self-consistently within our scheme. To simplify this task, we express them in terms of y 0 0 , v and the unknown y(0), which is achieved by setting x to 0 in the equation (58) and its derivatives. This leads to: ...
In previous works 9, 13 , the zeroth and first order approximations to V H , namely V H ≈ V 0 and V H (x) ≈ V 0 + xV 1 , were considered, although, as pointed out earlier, in slightly different versions from what we are now studying. 20 We will now establish a general procedure that solves the equation for an arbitrarily high order truncation of (59). We start by noticing that, since V H (x) is a given function, the equation (58) is in the linear form For α(x) = 0, the family of solutions of the above equation is This suggests to use the ansatz y(x) =ȳ(x)e β(x)dx when α(x) = 0. When inserted in (61), we obtain the following equation forȳ(x):ȳ which is solved by the family of functions where the lower bound is a fixed constant, while C parametrizes the family, leading to The value of a must be chosen carefully once the integrand β(x) is specified, in order to avoid possible divergences of the integral. When applied to (58), we have the exact result Since a has the dimension of a potential, we chose a = y 0 0 −1 , which will turn out to create no problems of convergence of the integral. We finally arrive at which is the family of solutions to (58), parametrized by the constant C.
Now that the family of solutions is established, the next task is to identify the physical one. It is possible to prove that the choice C = 0 is sufficient to pick a solution that reproduces the correct noninteracting function y 0 (x) in the limit v → 0, therefore we shall refer to that as the approximate physical solution. It should be noticed that C = 0, however, is not the only possible choice. This is because C is a constant with respect to x, the variable in the differential problem (58), but it can still depend on v and y 0 0 via the dimensionless quantity vy 0 0 2 ; it follows that any choice of such a dependence such that e −f (x) C → 0 for v → 0 would be allowed. This is a fundamental problem of the framework set by the Kadanoff-Baym equation, which will not be addressed here, but is left to future work. For the moment it should suffice to say that an arbitrary choice of such a dependency identifies a hierarchy of approximate physical solutions that do not, in general, converge to the exact one; our choice C = 0, on the other hand, seems to create a hierarchy with the correct convergence. Although we do not have an analytic proof for this fact, we do have strong numerical evidence, as we shall soon see.
Since each truncation of (59) corresponds to a different approximation of the original equation (58), we label the corresponding physical solution with an index (n), where n is the order of the truncation: which at equilibrium simplifies to Exact expressions for the integral in (70) can be found for n = 0 and n = 1, while for higher orders the integral can be evaluated numerically. As a matter of fact, equation (70) only gives an implicit definition of y p (0) to VH ≈ V0 + xV1 and so on.

E. Comparing all approximations
A comparison between all different approximations is reported in fig. 7. For the expansion in W RPA [G 0 ] we consider a Σ xc approximated with the first nontrivial term after the one prescribed by G 0 W 0 (which is of third order in W RPA [y 0 0 ], the second being identically zero); we also report the second order expansion of y(0) in W RPA [y 0 0 ], actually a better estimate than the expansion of Σ xc at lower computational cost. For the linearized KBE, we consider the simplest approximation V H (x) ≈ V H (0). We see that the expansion in W RPA [G 0 ] leads to an improvement over G 0 W 0 , as hoped. However, in the strongly coupled regime, the best approximation is given by the solution to the linearized KBE y  (43): the self-consistent GW (46); G0W0, which is obtained by taking only the leading order in (53); the next order expansion of Σxc in W0, obtained by taking the leading and nextto-leading order in (53); the second and third order expansion of G in W0 from (51); and, finally, the physical solution to the simplest approximation to the KBE, namely y (0) p from (70). Top panel: overall behavior; bottom panels: ratio between approximate and exact functions in the weakly (left) and strongly (right) coupling regimes.

VI. CONCLUSIONS
We presented two approaches aimed at improving on the GWA outside of the framework set by Hedin in his original formulation. The first one relies on the fact that G 0 W 0 can be regarded as the leading order expansion of Σ xc in G 0 and W 0 ≡ W RPA [G 0 ]. Since there are reasons to believe that W 0 may be a better parameter for such a perturbative expansion compared to the bare Coulomb potential v c and even the exact screened one W , we then propose to consider higher order terms in such an expansion.
The second approach relies on a Taylor expansion of the Hartree potential in terms of the perturbing potential V H [ϕ] in the KBE. We consider the hierarchy of linearized KBEs arising by higher and higher order approximations to V H [ϕ].
Both schemes have been tested on the one-point model. This is not a 'model system', in the sense of a model for physical systems, but rather as a 'model functional', in the sense of a model for the functionals used in the Many-body Perturbation Theory. The OPM allows to avoid certain theoretical limitations of our current perturbative tools and explore aspects of the theory that would be computationally too demanding to study in real systems. This is a particularly suitable framework for our study, since, as mentioned, there is not a known systematic way to deal with equations such as the KBE. Nonetheless, the procedure we used for the OPM is a necessary proof of principle of the feasibility of the approach and gives leads on how to tackle the full many-body problem.
Results concerning the quality of the approximations in comparison with the OPM equivalent of GWA are encouraging. The series in W RPA [G 0 ] seems to be very well behaved and provide good estimates of Green's functions in the strong coupling regimes already with a few orders. This contrasts with the expansion in W , which presents quite an erratic behaviour. The second approach works better, in fact beyond expectations. The solution to even the simplest linearized KBE leads to an approximate Green's function that outperform any other perturbative approach in the strong coupling regime. This suggests that the nonperturbative character of the KBE makes the equation a particularly promising starting point for methods in the strong coupling regime.