Solution to the many-body problem in one point

In this work we determine the one-body Greenʼs function as solution of a set of functional integro-differential equations, which relate the one-particle Greenʼs function to its functional derivative with respect to an external potential. In the same spirit as Lani et al (2012 New J. Phys. 14 013056), we do this in a one-point model, where the equations become ordinary differential equations (DEs) and, hence, solvable with standard techniques. This allows us to analyze several aspects of these DEs as well as of standard methods for determining the one-body Greenʼs function that are important for real systems. In particular: (i) we present a strategy to determine the physical solution among the many mathematical solutions; (ii) we assess the accuracy of an approximate DE related to the GW ?> +cumulant method by comparing it to the exact physical solution and to standard approximations such as GW ?> ; (iii) we show that the solution of the approximate DE can be improved by combining it with a screened interaction in the random-phase approximation. (iv) We demonstrate that by iterating the GW ?> Dyson equation one does not always converge to a GW ?> solution and we discuss which iterative scheme is the most suitable to avoid such errors.


Introduction
The one-body Greenʼs function σ σ ′ ′ ′ G t t r r ( , ) contains a wealth of information. Important quantities such as the density, the current density and, in particular, excitation energies such as electron removal and addition energies can be deduced from the one-body Greenʼs function. However, to solve the equation of motion (EOM) of the one-body Greenʼs function one requires knowledge of the two-body Greenʼs function, the EOM of which requires knowledge of the threebody Greenʼs function, etc. The standard procedure to cut short this chain of equations of increasing complexity is to introduce an effective potential, the self-energy, which takes into account all the many-body effects. The goal then becomes to find accurate approximations to the self-energy. Hedin proposed an iterative procedure which systematically improves the self-energy [1]. For example, one iteration of this procedure leads to the GW approximation. However, a second iteration already leads to a self-energy that is too complicated to be calculated for real materials.
This was one of the main motivations of Lani et al [2] to investigate the possibility to calculate the Greenʼs function directly without resorting to the self-energy. The starting point for their work is the Kadanoff-Baym equation (KBE). It is a functional differential equation (DE) for the generalized one-body Greenʼs function φ G (1, 2; [ ]) which is a functional of an external potential φ (1). Here we used σ = t r (1) ( , , ) 1 1 1 as a short-hand notation to combine the space, spin, and time variables. The equilibrium one-body Greenʼs function is retrieved in the limit of vanishing φ (1). The KBE is given by where v c is the Coulomb potential, G 0 is the Greenʼs function of the noninteracting system, and H c is the Hartree potential. Furthermore, we used the notation σ = + + t r (1 ) ( , , ) 1 1 1 , where η = + + t t 1 1 with η → + 0 . All the many-body effects (exchange and correlation (xc)) beyond the Hartree potential are contained in the last term on the right-hand side of equation ( Below we summarize the main questions that we will answer in this work.
• What are the solutions of the full KBE in one point?
• How can we determine the physical solution and is it unique? Solving the KBE leads to a huge number of mathematical solutions. It is therefore crucial that we have a general strategy to determine the physical solution among them. • How does the result obtained from the linearized KBE compare to the exact solution and to well-established approximations such as GW ? The linearized KBE stands at the basis of the GW +cumulant approach which has become increasingly important for the calculation of photoemission spectra [3][4][5][6][7]. It is therefore important to assess its accuracy. • What is the best screened Coulomb interaction to include in the linearized KBE? The linearized KBE can be combined with several approximate screened Coulomb interactions. It is important to verify how the results depend on this choice. • What is the best way to iterate the GW Dyson equation? The GW Dyson equation can be iterated in more than one way and it is not guaranteed that an iteration scheme will converge to the physical solution. Therefore, we will investigate which iteration schemes converge to the physical solution and which do not.
The paper is organized as follows. In section 2 we obtain the general solution of the KBE in one point and determine its physical solution. In section 3 we compare it to the solution of the linearized KBE. We compare the physical solution of the full and linearized KBE to standard approximations such as GW in section 4. In the same section we also investigate several schemes to iterate the GW equations. Finally, we summarize our conclusions in section 5.

KBE in the one-point model
In this section we solve the KBE and demonstrate how to find its physical solution. To transform the KBE in equation (1) to an equivalent and pertinent expression for the one-point model we use the following substitutions: i . As pointed out in [2] the change of prefactors compensates for the time integrals that are dropped in the one-point model. Moreover, this choice guarantees that the following constraints are satisfied when passing from the full functional problem to the one-point model: (i) since the diagonal of G is the density, i.e., ρ − = , the physical solution of the KBE in one point at vanishing potential is real and nonnegative; (ii) for real, nonnegative values of v the inverse dielectric constant in the random-phase approximation (RPA) is real and has values between 0 and 1.
The KBE in one point thus becomes Before we discuss the general solution of equation (3) let us first determine the solution of the KBE in one point for a noninteracting system under the influence of an external potential z. The solution to this problem can be obtained by setting v = 0 in equation (3). We obtain is also a solution of equation (3) for any finite v as can be verified by substitution. This is due to the cancelation of the two terms in equation (3) that contain v which represent Hartree and xc contributions in one point. This cancelation is similar to what occurs for the hole part of the KBE of a one-electron system. However, we are interested in describing using a simple model some aspects of a many-electron system, for which Hartree and xc terms in equation (1) will only partially cancel. To achieve this one could, for example, study a two-point model. In this case, however, the analytical solution of the KBE is, to the best of our knowledge, out of reach. Therefore, here we propose to generalize equation (3) by introducing a parameter λ according to We will now solve equation (5) for two values of λ (λ = 1 and λ = 1 2 ) and find for each its physical solution. We will show that for λ = 1 the physical solution is indeed equal to y z ( ) 0 .
2.1. Solution of the KBE in one point for λ ¼ 1 We start by studying the case λ = 1 as it represents the KBE of equation (1)   However, obtaining the mathematical solutions of the KBE is meaningless if we would not be able to determine the function C v y ( , ) 0 0 that corresponds to the physical solution. It is therefore important to be able to find a strategy to obtain C v y ( , ) 0 0 corresponding to the physical solution that will also be applicable to the general KBE in equation (1). In the next subsection we will present such a strategy and demonstrate that it is essential that C v y ( , ) 0 0 is independent of v.

Physical solution.
It is a major challenge to find constraints that a physical solution should satisfy. We follow the approach of [2] and study the behavior of equation (6) in the limit of vanishing interaction, i.e., → v 0. In this limit a physical solution should reduce to y z ( ) 0 , the solution of the KBE for noninteracting systems, for any z. We will, as a first step, determine . This can be achieved by choosing, for example, = + z i y (1 ) 0 0 . In order to study we make use of the asymptotic expansion of the error function,  The assumption that C v y ( , ) 0 0 has a Taylor expansion around v = 0 is motivated by the fact that the physical result thus obtained reproduces the perturbative result, i.e, substitution of the perturbative series ∑ = ∞ a z v ( ) n n n 0 into equation (3) and solving for the coefficients a n (z) gives back equation (9). The advantage of directly determining C v y ( , ) 0 0 corresponding to the physical solution, is that we immediately obtain the resummed result. This point will become clearer when we will discuss the case λ = 1 2 in the next subsection. We conclude that in the one-point model the Greenʼs function is independent of the interaction (for λ = 1). As mentioned before, we know that in the full functional problem given by equation (1) there is partial cancelation between Hartree and xc terms. For this reason we will now study a case for which λ ≠ 1 to mimic this partial cancelation in one point.

Solution of the KBE in one point for
In order to have partial cancelation we now study the case λ = 1 2 . The rationale behind this choice is related to the fact that in the full functional problem given by equation (1) an important contribution to the last term is the Fock exchange. Since the prefactor of the exchange term in the Fock operator of Hartree-Fock theory is half the prefactor of the term involving the Hartree potential we chose λ =   The derivation of this result can be found in appendix A.2.

Physical solution.
To obtain the physical solution we once again investigate y(z) in the limit → v 0. In this limit a physical solution should reduce to y z ( ) 0 for any z. We will first determine = C v y ( 0, ) 0 0 by studying the limit → + v 0 for potentials z such that < y z ( ) 0 0 2 . Using the asymptotic expansion of the error function (7) we can rewrite the general solution for this case as In that case we obtain It can be verified that the Taylor series of the right-hand side of equation (13) around v = 0 coincides with the perturbative solution of equation (10): v n n 0 0 0 2 As mentioned in the previous subsection the advantage of directly determining C v y ( , ) 0 0 of the physical solution, is that we immediately obtain the resummed result (13).
We see that the physical solution (13) of the KBE in equation (10) depends on the interaction v since there is only partial cancelation between the terms involving v.
In conclusion, we have solved the KBE in one point for λ = 1 and λ = 1 2 and we have determined the physical solution among the many mathematical solutions by demanding that it reduces to the noninteracting solution at vanishing interaction and that it does not contain nonanalytic contributions. This strategy can also be applied to the solution of the full functional problem in equation (1).

Quality of the linear approximation
The GW method is the state-of-the-art approach for calculations of quasiparticle energies and photoemission spectra [8,9]. Recently, the GW +cumulant approach [3][4][5][6][7] has received much attention, since the photoemission spectra obtained with this method give a much better description, both qualitatively and quantitavely, of satellite peaks with respect to those obtained within the GW approach [4]. The GW +cumulant approach is based on the KBE of equation (1) in which the Hartree potential is linearized, i.e., the Hartree potential is replaced by its firstorder Taylor expansion around φ = 0 [4]. It is therefore important to assess how accurate this linearized KBE is with respect to the full KBE of equation (1). Having obtained the exact physical solutions of the KBE in the one-point model we can now compare with the approximate physical solutions of the linearized KBE in one point which have been obtained in [2].
In order to keep this article self-contained as well as to introduce some quantities that will prove useful in the next section we will briefly demonstrate how one obtains the linearized KBE for the one-point model. Let us first rewrite equation (5) The latter two expressions represent the inverse dielectric function, and the (reducible) polarizability in one point, respectively. Introducing the change of variable, leads to the following linear DE, We note that equation (23) is not the same as the DE used in [2]. There, y H 0 and = u z ( 0) were assumed to be given, and independent of the interaction, whereas here we calculate them selfconsistently. This also explains the different choice for the sign of the last term in our Equation (5) with respect to the corresponding equation in [2].
Using a similar strategy as before we obtain the physical solution of equation (23). It is given by ⎛ . For the details on how to solve equation (23) and how to find its physical solution given in equation (25) we refer the reader to [2] in which an equivalent equation to equation (23) is discussed in detail.
To compare y 0 lin to the exact solutions given in equations (9) and (13) (19) and (24). This leads to an explicit expression for y 0 lin in terms of v and y 0 0 , provided that the iterations lead to a converged result. Both approaches lead to the same result. For the calculation of = u z (˜0) we need ′ = y z ( 0). We obtain this quantity from equation (15). In order to numerically compare Greenʼs functions obtained within different approximations it is convenient to plot the scaled Greenʼs function  solution remains much closer to the exact solution than for λ = 1. This is due to the fact that an approximation based on a screened interaction is more suitable to describe an interacting system than a noninteracting system.
In figures 1 and 2 we also report y 0 lin obtained using the exact u and ϵ = − u v RPA RPA 1 . We see that improving the description of the screening does not necessarily improve the solutions since using the exact screening worsens the results. While for λ = 1 the results obtained with an RPA screening deteriorate with respect to y 0 lin , those for λ = 1 2 show an overall improvement when u RPA is employed. We note that tendencies obtained for λ = 1 2 should be given more weight than those obtained for λ = 1 since the former corresponds to a model of an interacting system.
In figures 3 and 4 we compare u 0 lin , i.e., the screened interaction u that corresponds to y 0 lin given by   2 . This shows that there is a considerable cancelation of error when calculating y 0 lin using u RPA . In conclusion, to obtain the best results for interacting systems the linearized KBE should be combined with an RPA screened interaction. We note that an RPA screened interaction is often used in practical calculations. For example, an RPA-type screened interaction was employed when the GW +cumulant method was applied to the description of the photo-emission spectrum of silicon leading to excellent results [4]. Also, the linearization of the KBE could be regarded as an approximation to the vertex in the self-energy expression of Hedinʼs equations. It is known that errors thus introduced in, for example, total energies can be partially canceled by doing a similar approximation of the vertex that appears in the expression for W [10,11].
In order to have a better idea of the relative accuracy of y 0 lin we will now compare it to some approximations obtained from Hedinʼs equations such as the GW approximation.

Solutions of Hedinʼs equations: GW and beyond
As mentioned in the introduction, solving the full functional problem of the KBE in equation (1) is generally considered to be a too complicated task. Therefore, a perturbative approach to obtain the one-body Greenʼs function is usually adopted. In such an approach one defines a selfenergy which contains all the many-body effects of the system, i.e., all the terms involving the interaction. Here we will compare to a perturbative approach due to Hedin [1] in which the perturbation expansion is done in terms of a screened interaction instead of the bare interaction.
In one point Hedinʼs equations are given by = + + y z y z y z zy z y z s y z y z H H H xc λ = s y z y z u z g y z ( ( )) ( ) ( )˜( ( )), where s xc is the exchange-correlation part of the self-energy in one point and p z ( ) and g y z ( ( )) represent the irreducible polarizability and vertex function in one point, respectively. The derivation of the expressions given in equations (27)-(31) can be found in appendix B. Hedinʼs equations form a closed set and are equivalent to the KBE of equation (5). Hedinʼs equations are usually expressed for a vanishing external potential φ. Equivalent expressions can be obtained in one point by setting z = 0 in the above equations.
Since Hedinʼs equations are equivalent to the KBE, their general solution for y(z) is given by equation (6) and its physical solution by equation (9) in the case λ = 1. For λ = 1 2 the general and physical solution are given by equations (11) and (13), respectively. We can now use this   where the functions a(z) and b(z) are defined by We see that for λ = 1 the exchange-correlation part of the self-energy s xc (z) is equal to minus the Hartree potential. This is a consequence of the fact that in one point the physical solution is equal to the noninteracting solution. In this case the interactions included in the Hartree potential that are implicitly contained in y H (z) are spurious. Therefore, the only purpose of s xc (z) is to cancel the effects of the Hartree potential. As a result the Greenʼs function in one point y(z) becomes equal to the noninteracting Greenʼs function as it should. We note that Molinari et al [12,13] have studied a set of Hedinʼs equations (at z = 0) within the one-point model that are slightly modified with respect to those presented in this section. The modification is in the irreducible polarization given in equation (30) to which they added a minus sign, i.e., = − p y g2 . Starting from these modified Hedinʼs equations Pavlyukh and Hübner [14] found an implicit solution for the Greenʼs function. In appendix C we demonstrate how one can obtain an explicit solution.
In the following we will evaluate some approximations to Hedinʼs equations. Since we are mainly interested in the Greenʼs function at vanishing external potential φ, we will only consider the case z = 0 in the remainder of this section.

GW
The GW approximation in one point is obtained by setting = s 0 xc and performing one iteration of Hedinʼs equations starting with the vertex function g. Therefore, in the GW approximation, Hedinʼs equations in one point are written as = g 1, where in the last step we used equation (19).  The other three solutions diverge when → v 0. We note that this is consistent with the first theorem derived in [15]. We obtain similar results for λ = 1 2 .
In figures 1 and 2 we compare y GW to the exact physical solution of the KBE and the approximate physical solution of the linearized KBE. For λ = 1 at small interaction y 0 lin is somewhat better than y GW . Otherwise, we observe that the GW approximation leads to results that are generally very close to those obtained from the linearized KBE. This finding can be surprising since in the full functional problem the linearized KBE contains many more diagrams than GW . However, in the one-point model many features, such as the spectral function, for which these diagrams are important [3][4][5][6][7], are absent. For this reason the difference between GW and the linearized KBE cannot be detected in one point.
In figures 3 and 4 we compare the GW screened interaction u GW , given by GW GW 2 to the exact screening and to u 0 lin . For λ = 1 the GW screening is slightly worse than u 0 lin while in the case of λ = 1 2 the GW screening is somewhat better than u 0 lin . The GW Dyson equation can be rewritten in several ways and therefore many iterative schemes are possible. Let us first study the most common iteration scheme. It is given by  where the Hartree Greenʼs function G n H, and the screened interaction W n are given by To iteratively solve Hedinʼs equations we start from the noninteracting solution, i.e., = y y 0 0 0 , followed by substitution of this expression into the right-hand side of equation (64). This leads to the first-order approximation to the GW solution, This result represents the G W 0 0 solution in one point. By substituting the first-order result into equation (64) one obtains the second-order approximation to y GW and continuation of this iteration scheme will lead to higher-order approximations to y GW .
In figures 6 and 7 we compare several of these approximations to the exact GW result in one point for λ = 1 and λ = 1 2 , respectively. We see that y G W 0 0 is close to y GW for small values of the interaction but largely underestimates it when the interaction becomes stronger. Instead, y G W 1 1 largely overestimates y GW for large v. However, the most interesting result is that beyond a certain interaction strength, i.e., λ > y v [ ] 4 0 0 2 2 , the iterative result does not converge to y GW . In fact, the iterative result oscillates between two values, neither of which is equal to y GW . This conclusion also holds when the starting point is chosen close to y GW .
To explain the discrepancy between the exact and iterative GW result we first note that equation (64)  2 . We note that this is consistent with the second theorem derived in [15].
In [2] it was shown that the solution obtained by solving the GW equations iteratively depends strongly on the iteration scheme used. In their case they solved Hedinʼs equations iteratively for fixed values of y 0 H and u. This leads to two solutions of which only one is a physical solution. They showed that depending on the chosen iteration scheme the solution could be either the physical or the nonphysical solution. Here we show that, moreover, there are iteration schemes that lead to results that are not solutions at all. We also note that we have obtained this result using an iterative scheme that is equivalent to that used in most practical applications [16][17][18].
We obtain a different iteration scheme when, instead of writing the Dyson equation as in equation (61)  One can verify that y GW is an attracting fixed point of equation (67) for all v. Therefore, iteration of equation (67) is guaranteed to converge to the exact GW solution (provided that the initial value lies in the vicinity of y GW ). However, one prefers to iterate the Dyson equation of equation (48) since this form has been implemented in many computer codes. Therefore, we propose an alternative scheme that also converges to the physical solution for all v and entails only a small modification of these codes. The small modification consists of iterating equation (61) for fixed W according to and when convergence is reached for G update W according to where ∞ G represents the converged G for fixed W. Then reiterate equation (68) with W obtained from equation (69), etc, until convergence for both G and W is achieved. In the one-point model this scheme will also lead to the exact GW solution for all v. Finally, we note that in the absence of analytical results there does not yet exist a general strategy which permits one to know whether the converged result of a given iteration scheme corresponds to a physical or a nonphysical solution. Therefore, the study of models for which analytical results exist give important insights on how to solve this problem. where we used equation (19). We therefore have to solve the following nonlinear equation With this self-energy we will obtain the direct solutions of equation (72) as well as its iterative solution using equation (71) in the next two subsections. where, for notational convenience, we dropped the superscript (2). The above equation for y (2) cannot be solved analytically. We can, however, solve it numerically for any given value of v.
We obtain eight solutions of which only one tends to y 0 0 in the limit → v 0. This is the physical solution   to the exact screening, to u 0 lin and to u GW . We see that for λ = 1 2 the Γ GW screening provides the best approximation to the exact screening for small interactions.

Iterative solution.
The Γ GW Dyson equation can be rewritten in even more ways than the GW Dyson equation and therefore a large number of iteration schemes are possible. Here we will report the results of the iteration scheme that is similar to the most common GW iteration scheme given in equation (61), i.e., In figure 8 we compare the results obtained with this iteration scheme to the physical Γ GW solution for λ = 1. Once more we observe that the iterative result does not converge to the physical solution beyond a certain interaction strength. As for GW it oscillates between two values. Moreover, compared to the iterative GW result, the range over which the iteration scheme converges to the physical solution is reduced by about 30%. We have obtained similar results for λ = 1 2 . This seems to indicate that using a more complex self-energy reduces the convergence range. In other words, beyond GW the choice of the iteration scheme becomes more crucial.

Conclusions
We solved the KBE in one point. We obtained a family of solutions and we showed that only one is a physical solution. Moreover, we presented a strategy to obtain the physical solution that is, in principle, not limited to the one-point model.
We showed that in the one-point model the physical solution is equal to the noninteracting solution. This is due to the fact that the terms in the KBE that contain the interaction cancel. Therefore, we proposed a slightly modified KBE in which this cancelation is partial, just as in the full functional problem. We also solved this DE and found its unique physical solution.
We assessed the accuracy of an approximate linearized KBE that lies at the heart of the GW +cumulant method. We compared it to the exact physical solution of the full KBE in one point as well as to standard approximations such as GW . The solutions of the linearized KBE are nearly exact for small interactions but deviate from the exact results for larger interactions. We showed that the solution of the linearized KBE can be improved by employing an RPA screened interaction. Adding low order vertex corrections does not seem to lead to significant improvement.
Much of the results are governed by the structure of the equations, which is the same in the one-point model and in real systems. One can hence extrapolate the findings, as long as one is careful not to over-interpret. Of course, important questions such as the position of poles corresponding to excitation energies can by definition not be treated in the one-point model.
Finally, we demonstrated that by iterating the GW Dyson equation in the usual way one does not always converge to a GW solution. We proposed a practical iterative scheme that leads to the physical GW solution for all interaction strengths.
We will now address the two cases, λ = 1 and λ = The general solution of a The solution of the homogeneous equation corresponding to the above equation is with K a constant. We can then vary the constant to obtain the solution of the inhomogeneous equation. It is given by ⎛ We can thus rewrite equation (A.6) according to

. A particular solution of this second-order DE is
which can be verified by substitution into equation (A.8). A second independent solution is then given by [19] The general solution for u(x) can therefore be written as ⎡ where for notational convenience we have suppressed the dependence of C 1 and C 2 on v and y 0 0 . The derivative of u(x) is given by ⎡ Substituting these results into equation (A.7) then leads to the general solution for y(x) ( ) In the limit → + v 0 the first term on the right-hand will tend to ∓ ∞ i while the second term on the right-hand side will tend either to zero (for ). This means that y(z) will tend to ∓ ∞ i for any