Many-body perturbation theory and non-perturbative approaches: screened interaction as the key ingredient

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({\bf r}\sigma t, {\bf r}'\sigma' 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({\bf r}\sigma t, {\bf r}'\sigma' t')$ and the Coulomb potential $v_{\rm c}({\bf r}, {\bf 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. The cornerstone of this approach is the Dyson equation, which in a simplified notation reads $G=G_0+G_0 \Sigma G$; here the focus is shifted on the self-energy Σ, whose perturbative expansion is still expressed in terms of G₀ 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₀ 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 (GWA). In a seminal paper, Hedin formalized this idea by introducing an appropriate screened interaction W and expressing the exchange-correlation contribution $\Sigma_{\rm xc}$ to the self-energy as a functional of G and W, rather than G₀ 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 GWA arises as the first step of an iterative procedure that formally solves Hedin's equations. This first approximation corresponds to setting the vertex $\Gamma=1$. The GWA 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 GWA 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 $\chi_0$. The inverse dielectric function $\newcommand{\e}{{\rm e}} \epsilon^{-1} = 1+v_{\rm c}\chi$, which screens the Coulomb interaction via the relation $\newcommand{\e}{{\rm e}} W=\epsilon^{-1}v_{\rm c}$, becomes $\newcommand{\e}{{\rm e}} \epsilon^{-1}=1+v_{\rm c}\chi_0$. Since $\chi_0$ is negative, W is indeed often smaller than v_c. However, for large enough (negative) $\chi_0$ compared to $1/v_{\rm 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 $\newcommand{\e}{{\rm e}} \epsilon=1-v_{\rm c}\chi_0$, which is always positive and larger than 1, such that the resulting screened interaction $W_0<v_{\rm c}$. This suggests that an expansion in terms of the RPA W₀ 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 GWA 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 zero-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–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 2 will be devoted to recall the reader the theoretical framework; our two approaches will be motivated and presented in sections 3 and 4, respectively; results of tests on the OPM are presented in section 5; feasibility and computational cost of the approaches for real systems will be discussed in section 6 where conclusions will also be drawn.

In a non-interacting system, one can determine the one-body Green's function G₀ 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₀. 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 function s 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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{KBE} G(1,2;[\varphi]) =&\, G_0(1,2) + \int {\rm d}3 G_0(1,3) \varphi(3) G(3,2;[\varphi])\nonumber \\ & + \int {\rm d}3 G_0(1,3) V_{\rm H}(3;[\varphi]) G(3,2;[\varphi]) \nonumber \\ &+{\rm i}\int {\rm d}3{\rm d}4 \,G_0(1,3) v_{\rm c}(3^+,4) \frac{\delta G(3,2;[\varphi])}{\delta \varphi(4)} , \nonumber \end{align} \tag{ 1 }$$

where v_c is the Coulomb potential, $G_0(1, 2)=$$G(1, 2;\varphi=0)\Big\vert _{v_{\rm c}= 0}$ is the Green's function of the noninteracting system, and the Hartree potential $V_{{\rm H}}(1;[\varphi])$ is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{h.1} V_{{\rm H}}(1;[\varphi]) = -{\rm i}\int {\rm d}2 v_{\rm c}(1,2)G(2,2^+;[\varphi]) . \nonumber \end{align} \tag{ 2 }$$

Here, we use $(1)=({\bf r}_1, \sigma_1, t_1)$ as a short-hand notation to combine the space, spin, and time variables. Moreover, $(1^+) = ({\bf r}_1, \sigma_1, t^+_1)$, where $\newcommand{\e}{{\rm e}} t^+_1 = t_1 + \eta$ with $\newcommand{\e}{{\rm e}} \eta\to 0^+$.

In equation (1) the generalized one-body Green's function $G(1, 2; [\varphi])$ is a functional of a time-dependent external potential $\varphi(1)$. The equilibrium one-body Green's function of interest is retrieved in the limit of vanishing $\varphi(1)$. Without the last term on the right-hand side of equation (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 equation (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 section 5 for an illustration). Therefore a common way to go is to transform equation (1) into a Dyson equation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:dyson} G(1,2;[\varphi]) = G_{\rm H}(1,2;[\varphi]) +{\rm i} \,G_{\rm H}(1,\bar 3;[\varphi]) \Sigma_{\rm xc}(\bar 3,\bar 6;[\varphi]) G(\bar 6,2;[\varphi]) , \nonumber \end{align} \tag{ 3 }$$

with the exchange-correlation self-energy

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Sigma_{\rm xc}(3,6;[\varphi]) \equiv -{\rm i} v_{\rm c}(3^+,\bar 4) G(3,\bar 5;[\varphi])\frac{\delta G^{-1}(\bar 5, 6;[\varphi])}{\delta \varphi(\bar 4)}. \label{eq:sigma} \nonumber \end{align} \tag{ 4 }$$

Here and in the following, we denote integrals by $\newcommand{\e}{{\rm e}} f(\bar n)g(\bar n)\equiv \int {\rm d}n\, f(n)g(n)$.

Once the expression for the self-energy is established, it is enough to solve equation (3) for $\varphi \to 0$ in order to obtain the equilibrium Green's function. Of course, since G is unknown so is $\Sigma_{\rm 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 $\Sigma_{\rm xc}$ is easier than approximating G itself.

One strategy to find approximations for $\Sigma_{\rm xc}$ has been formalized by Hedin [3], 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, $\varphi_{\rm cl}=\varphi + V_{\rm H}$, to express the exact exchange-correlation self-energy (4) as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Sigma_{\rm xc}(3,6;[\varphi])& \equiv -{\rm i} v_{\rm c}(3^+,\bar 4) G(3,\bar 5;[\varphi])\frac{\delta G^{-1}(\bar 5, 6;[\varphi_{\rm cl}])}{\delta \varphi_{\rm cl}(\bar 7)}\frac{\delta \varphi_{\rm cl}(\bar 7)}{\delta \varphi(\bar 4)}\nonumber \\ &\equiv {\rm i}G(3\bar 4)W(3^+,\bar 5)\tilde\Gamma(\bar 4,6;\bar 5) \label{eq:sigma-2} \nonumber \end{align} \tag{ 5 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:WDOP} W(1,2) = v_{\rm c}(1,2) + v_{\rm c}(1,\bar 3)P(\bar 3,\bar 4)W(\bar 4,2) \nonumber \end{align} \tag{ 6 }$$

and the irreducible polarizability

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:PG} P(1,2) = -{\rm i} G(1,\bar 3)G(\bar 4,1)\tilde\Gamma(\bar 3,\bar 4;2). \nonumber \end{align} \tag{ 7 }$$

From the Dyson equation (3), the fifth equation in Hedin's set is derived for the irreducible vertex function $\newcommand{\e}{{\rm e}} \tilde \Gamma \equiv -\delta G^{-1}/\delta \varphi_{\rm cl}$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\sxc}{\Sigma_{\rm xc}} \displaystyle \label{eq:GAMMA} \tilde\Gamma(1,2;3) = \delta(1,2)\delta(1,3)+ \frac{\delta \sxc(1,2)}{\delta G(\bar 4,\bar 5)} G(\bar 4,\bar 6)G(\bar 7,\bar 5)\tilde\Gamma(\bar 6,\bar 7;3). \nonumber \end{align} \tag{ 8 }$$

The screened Coulomb interaction W is linked to the density–density response function

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi(1,2) \equiv \frac{\delta n(1)}{\delta \varphi(2)} \nonumber \end{align} \tag{ 9 }$$

via $W(1, 2) = v_{\rm c}(1, 2) + v_{\rm c}(1, \bar 3) \chi(\bar 3, \bar 4) v_{\rm c}(\bar 4, 2)$, and it can also be expressed in terms of the inverse dielectric function $\newcommand{\e}{{\rm e}} \epsilon^{-1}(1, 2)\equiv \delta(1, 2)+v_{\rm c}(1, \bar 3)\chi(\bar 3, 2)$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle W(1,2) = \epsilon^{-1}(1,\bar 3)v_{\rm c}(\bar 3,2). \nonumber \end{align} \tag{ 10 }$$

By setting $\tilde \Gamma\approx 1$, the equations reduce to a closed set of algebraic equations for G, Σ, W and P, leading to the so called GWA for the self-energy and the 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₀, although it contains some interaction through, e.g. the KS potential or the Fock exchange. This approach is usually called $G_0W_0$. For later reference we highlight the fact that W is calculated in the RPA by using the notation

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \displaystyle \label{oneshotGW} \Sigma_{\rm xc}\approx {\rm i}G_0\Wrpa{G_0} \,\,\,\,\,\,\,\,\,\,\, W\approx \Wrpa{G_0}\ \,\,\,\,\,\,\,\,\,\,\, P \approx -{\rm i}G_0G_0 \,\,\,\,\,\,\,\,\,\,\, \tilde \Gamma \approx 1 , \nonumber \end{align} \tag{ 11 }$$

where the RPA functional is $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \newcommand{\e}{{\rm e}} \Wrpa{F}\equiv (1+{\rm i}FFv_{\rm c}){\hspace{0pt}}^{-1}v_{\rm c}$. Such a functional generalizes the functional form of W obtained from the RPA for P, where $\tilde\Gamma =1$ implies $P\approx-{\rm i}GG$ and hence $W\approx (1+{\rm i}v_{\rm c} GG){\hspace{0pt}}^{-1}v_{\rm c}$.

The $G_0W_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₀. 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 [14]. 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 fulfil the f-sum rule [15]. Clearly, vertex corrections, i.e. corrections to the self-energy of higher order in W, are needed to improve this situation.

The success of the GWA 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 $\Sigma_{\rm xc}$. This presumes that $\vert W\vert < v_{\rm c}$. Although often true, this is simplistic, and as explained in the introduction $\vert W\vert$ can in principle become arbitrarily large. On the other hand, $P^{\rm RPA}(\omega=0;[G_0])$ built with any mean-field G₀ is always negative, which means that $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \Wrpa{G_0}$ is always positive and smaller than v_c. Moreover, the corresponding density-density response function $\chi^{\rm RPA}[G_0]$ is well behaved; in particular, it fulfils the f-sum rule.

Based on these considerations, we propose to consider the perturbative expansion of $\Sigma_{\rm xc}$ in G₀ and $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \Wrpa{G_0}$. This is formally constructed in the following way. First, we use the definition of $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \Wrpa{G_0}$ to invert the map $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} v_{\rm c} \to \Wrpa{G_0}$, which leads to $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} v_{\rm c}=(\Wrpa{G_0}^{-1} + P_0){\hspace{0pt}}^{-1}=$$\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \Wrpa{G_0}(1-{\rm i}\Wrpa{G_0}G_0G_0){\hspace{0pt}}^{-1}$ with $\newcommand{\e}{{\rm e}} P_0 \equiv -{\rm i}G_0G_0$. We then use the well established perturbation expansion of $\Sigma_{\rm xc}$ in terms of v_c and insert the above expression for v_c to rewrite the series in terms of $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \Wrpa{G_0}$ and, finally, we group terms of the same order in $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \Wrpa{G_0}$. By defining $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \newcommand{\e}{{\rm e}} W_0\equiv \Wrpa{G_0}$ we can schematically represent the procedure as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Sigma_{\rm xc}&={\rm i}v_{\rm c} G_0+v_{\rm c} G_0 v_{\rm c} G_0G_0+...\nonumber \\ & = {\rm i} (W_0^{-1} + P_0)^{-1}G_0 +(W_0^{-1} + P_0)^{-1}G_0(W_0^{-1} + P_0)^{-1} G_0 G_0+...\nonumber \\ &={\rm i}W_0G_0+ W_0(G_0 W_0 G_0G_0 -G_0G_0 W_0G_0)+... \nonumber \end{align} \tag{ 12 }$$

Being more specific, if the starting G₀ is, for instance, the Hartree Green's function G_H we have that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &\Sigma_{\rm xc}(1,2)=G_{\rm H}(1,2) W_{\rm H}(1,2) +\left(\delta(1,2)W_{\rm H}\right.\nonumber \\ &\left.(1,3)G_{\rm H}(3,4)G_{\rm H}(4,5)G_{\rm H}(5,3)W_{\rm H}(4,5)+\right.\nonumber \\ &\left.+G_{\rm H}(1,3)G_{\rm H}(3,4)G_{\rm H}(4,2)W_{\rm H}(1,4)W_{\rm H}(3,2)\right.\nonumber \\ &\left.+ G_{\rm H}(1,3)G_{\rm H}(3,4)G_{\rm H}(4,2)W_{\rm H}(1,2)W_{\rm H}(3,4)\right)+... \nonumber \end{align} \tag{ 13 }$$

where $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \newcommand{\e}{{\rm e}} W_{\rm H}\equiv \Wrpa{G_{\rm H}}$. The first term represent the GWA, while the others, between brackets, the second order contribution in W_H. Diagrammatically, this is represented by

The previous section designs a way to use a linear response ingredient, namely W₀, 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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{sigma.ini} \Sigma[\varphi] =\varphi + V_{\rm H}[\varphi]+{\rm i}v_{\rm c}\frac{\delta G[\varphi]}{\delta \varphi}G^{-1}[\varphi]. \nonumber \end{align} \tag{ 14 }$$

We assume that $V_{\rm H}[\varphi]$ (defined in (2)) is Taylor expandable:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{VTaylor} V_{\rm H}[\varphi]=-{\rm i}v_{\rm c} n-{\rm i} \varphi v_{\rm c} \frac{\delta G[\varphi]}{\delta \varphi}\Big\vert _{\varphi\to0} -\frac{\rm i}{2}\varphi\varphi v_{\rm c} \frac{\delta^2 G[\varphi]}{\delta \varphi\delta \varphi}\Big\vert _{\varphi\to0}+... \nonumber \end{align} \tag{ 15 }$$

where $\newcommand{\e}{{\rm e}} n(1)\equiv G(1, 1^+;[0])$ is the charge density. As first approximation, we truncate the above series to first order

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{approx1} V_{\rm H}[\varphi]\approx -{\rm i}v_{\rm c} n-{\rm i} \varphi v_{\rm c} \frac{\delta G[\varphi]}{\delta \varphi}\Big\vert _{\varphi\to0} \nonumber \end{align} \tag{ 16 }$$

which leads to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{sigma.1app} \Sigma[\varphi] \approx -{\rm i}v_{\rm c} n+\varphi \left(\delta-{\rm i} v_{\rm c} \frac{\delta G[\varphi]}{\delta \varphi}\Big\vert _{\varphi\to0}\right) +{\rm i}v_{\rm c}\frac{\delta G[\varphi]}{\delta \varphi}G^{-1}[\varphi]. \nonumber \end{align} \tag{ 17 }$$

This rewriting suggests to renormalize (or screen) the external potential as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{barphi} \varphi\to \bar{\varphi}\equiv \varphi \left(\delta-{\rm i} v_{\rm c} \frac{\delta G[\varphi]}{\delta \varphi}\vert _{\varphi\to0}\right). \nonumber \end{align} \tag{ 18 }$$

Since (18) implies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\delta G[\varphi]}{\delta \varphi} =\frac{\delta G[\varphi[\bar{\varphi}]]}{\delta \bar{\varphi}} \frac{\delta \bar{\varphi}[\varphi]}{\delta \varphi} = \frac{\delta \bar{G}[\bar{\varphi}]}{\delta \bar{\varphi}} \left(\delta-{\rm i}v_{\rm c}\frac{\delta G[\varphi]}{\delta \varphi}\Big\vert _{\varphi\to 0} \right) \nonumber \end{align} \tag{ 19 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Rightarrow \frac{\delta G[\varphi]}{\delta \varphi}\Big\vert _{\varphi\to 0} = \frac{\delta \bar{G}[\bar{\varphi}]}{\delta \bar{\varphi}}\Big\vert _{\bar{\varphi}\to 0 \Leftrightarrow \varphi\to 0} \left(\delta-{\rm i}v_{\rm c}\frac{\delta G[\varphi]}{\delta \varphi}\Big\vert _{\varphi\to 0} \right) \nonumber \end{align} \tag{ 20 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Rightarrow \frac{\delta G[\varphi]}{\delta \varphi}\Big\vert _{\varphi\to 0} = \left(\delta+{\rm i}v_{\rm c}\frac{\delta \bar{G}[\bar{\varphi}]}{\delta \bar{\varphi}}\Big\vert _{\bar{\varphi}\to 0} \right)^{-1} \frac{\delta \bar{G}[\bar{\varphi}]}{\delta \bar{\varphi}}\Big\vert _{\bar{\varphi}\to 0} \nonumber \end{align} \tag{ 21 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Rightarrow \frac{\delta G[\varphi]}{\delta \varphi} = \frac{\delta \bar{G}[\bar{\varphi}]}{\delta \bar{\varphi}} \left(\delta-{\rm i}v_{\rm c} \left(\delta+{\rm i} v_{\rm c} \frac{\delta \bar{G}[\bar{\varphi}]}{\delta \bar{\varphi}}\Big\vert _{\bar{\varphi}\to 0} \right)^{-1} \frac{\delta \bar{G}[\bar{\varphi}]}{\delta \bar{\varphi}}\Big\vert _{\bar{\varphi}\to 0} \right) \nonumber \end{align} \tag{ 22 }$$

with $\newcommand{\e}{{\rm e}} \bar{G}[\bar{\varphi}]\equiv G[\varphi[\bar{\varphi}]]$, we can write (17) as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{intermedio}& \bar{\Sigma}[\bar{\varphi}]\equiv \Sigma[\varphi[\bar{\varphi}]] \approx -{\rm i}v_{\rm c} n+\bar{\varphi}+{\rm i}v_{\rm c} \frac{\delta \bar{G}[\bar{\varphi}]}{\delta \bar{\varphi}}\nonumber \\ &\left(\delta-{\rm i}v_{\rm c} \left(\delta+{\rm i} v_{\rm c} \frac{\delta \bar{G}[\bar{\varphi}]}{\delta \bar{\varphi}}\Big\vert _{\bar{\varphi}\to 0} \right)^{-1} \frac{\delta \bar{G}[\bar{\varphi}]}{\delta \bar{\varphi}}\Big\vert _{\bar{\varphi}\to 0} \right)G^{-1}[\varphi]. \nonumber \end{align} \tag{ 23 }$$

The functional derivative can be approximated with the leading order of its perturbative expansion in $\bar{G}[\bar{\varphi}]$, as⁵

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{approx2} \frac{\delta \bar{G}[\bar{\varphi}]}{\delta \bar{\varphi}}\approx \bar{G}[\bar{\varphi}]\bar{G}[\bar{\varphi}]. \nonumber \end{align} \tag{ 24 }$$

Combining (23) and (24) one finds

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &\bar{\Sigma}[\bar{\varphi}] \approx -{\rm i}v_{\rm c} n+\bar{\varphi}+{\rm i}v_{\rm c} \bar{G}[\bar{\varphi}]\bar{G}[\bar{\varphi}] \nonumber \\ &\left(\delta-{\rm i}v_{\rm c} \left(\delta+{\rm i} v_{\rm c} \bar{G}[0]\bar{G}[0] \right)^{-1} \bar{G}[0]\bar{G}[0] \right)G^{-1}[\varphi]. \nonumber \end{align} \tag{ 25 }$$

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 $\varphi=\bar{\varphi}=0$. The equilibrium self-energy is then obtained as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{sigma.almostrpa} \bar{\Sigma}[0]=\Sigma[0]\approx V_{\rm H}[0] +{\rm i}v_{\rm c} \left(1 -{\rm i} G[0] G[0]v_{\rm c} \left(1+{\rm i} v_{\rm c} G[0]G[0] \right)^{-1} \right) G[0] \nonumber \end{align} \tag{ 26 }$$

since $\bar{G}[0]=G[0]$. Here the RPA screened interaction appears, since

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{W}_{\rm RPA}[\,f]&=v_{\rm c} \epsilon^{-1}_{\rm RPA}[\,f]=v_{\rm c}(1+v_{\rm c} \chi_{\rm RPA}[\,f])\nonumber \\ &\equiv v_{\rm c}\left(\delta + v_{\rm c} \left(-{\rm i}ff(\delta+{\rm i}ff v_{\rm c})^{-1}\right)\right) \nonumber \end{align} \tag{ 27 }$$

where f denotes a generic function on which the operator $\mathcal{W}_{\rm RPA}$ acts. Since $\bar{G}[0]=G[0]$ the self-energy (26) becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{sigma.AD} \Sigma[0] \approx V_{\rm H}[0]+{\rm i}\mathcal{W}_{\rm RPA}[G[0]] G[0] \approx V_{\rm H} +{\rm i}\mathcal{W}_{\rm RPA}[G_0] G_0. \nonumber \end{align} \tag{ 28 }$$

These are, respectively, the self-consistent GWA with $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} W\approx\Wrpa{G}$, and the so-called $G_0W_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₀.

The above derivation of the $G_0W_0$ approximation relies on two ordered approximations: (16) and (24), that are tied to each other via the definition of the screened perturbing potential $\bar{\varphi}$. 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 (24) directly on (14), without going through (16) and the definition of $\bar{\varphi}$, one would get to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Sigma[\varphi]\approx \varphi + V_{\rm H}[\varphi]+{\rm i}v_{\rm c} G[\varphi] \nonumber \end{align} \tag{ 29 }$$

which is simply the Hartree–Fock approximation.

The derivation seems to suggest that a better approximation may come from improving on (16), namely taking more orders of the Taylor expansion of $V_{\rm H}[\varphi]$. However, solving the Dyson equation at equilibrium kills all possible improvements. Such a mechanism is clear already when keeping the second order in (15). Let us first write

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Sigma[\varphi]\approx &-{\rm i}v_{\rm c} n+\varphi \left(\delta-{\rm i} v_{\rm c} \frac{\delta G[\varphi]}{\delta \varphi}\Big\vert _{\varphi\to0} -\frac{\rm i}{2}\varphi v_{\rm c} \frac{\delta^2 G[\varphi]}{\delta \varphi\delta \varphi}\Big\vert _{\varphi\to0}\right) \nonumber \\ &+{\rm i}v_{\rm c}\frac{\delta G[\varphi]}{\delta \varphi}G^{-1}[\varphi]. \nonumber \end{align} \tag{ 30 }$$

Now we introduce the screened perturbing potential (in the even shorter notation for which functional derivatives $\delta F[\varphi]/\delta \varphi$ are denoted by $F'[\varphi]$)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{\varphi}\equiv \varphi(1-{\rm i}v_{\rm c}G'[0]-\frac{\rm i}{2}v_{\rm c}\varphi G''[0]) \nonumber \end{align} \tag{ 31 }$$

for which

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G'[\varphi]=\tilde{\varphi}'[\varphi]\tilde{G}'[\tilde{\varphi}]=(1-{\rm i}v_{\rm c}G'[0]-\frac{\rm i}{2}\varphi v_{\rm c} G''[0])\tilde{G}'[\tilde{\varphi}] \nonumber \end{align} \tag{ 32 }$$

where $\newcommand{\e}{{\rm e}} \tilde{G}[\tilde{\varphi}]\equiv G[\varphi[\tilde{\varphi}]]$. This implies that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G'[0]=\tilde{G}[0](1+{\rm i}v_{\rm c}\tilde{G}'[0])^{-1} \nonumber \end{align} \tag{ 33 }$$

which, in combination with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{approx2.1} \frac{\delta \tilde{G}[\tilde{\varphi}]}{\delta \tilde{\varphi}}\approx \tilde{G}[\tilde{\varphi}]\tilde{G}[\tilde{\varphi}], \nonumber \end{align} \tag{ 34 }$$

gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & \bar{\Sigma}[\bar{\varphi}]\equiv \Sigma[\varphi[\bar{\varphi}]]\approx -{\rm i}v_{\rm c} n+\tilde{\varphi}+{\rm i}v_{\rm c}\nonumber \\ &\left(1 -{\rm i}\tilde{G}[0]\tilde{G}[0](1+{\rm i}v_{\rm c} \tilde{G}[0]\tilde{G}[0])^{-1}v_{\rm c}\right) \tilde{G}[\tilde{\varphi}] \nonumber \end{align} \tag{ 35 }$$

which at equilibrium is just

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \bar{\Sigma}[0]=\Sigma[0]\approx -{\rm i} v_{\rm c} \tilde{G}[0]+{\rm i}\mathcal{W}_{\rm RPA}[\tilde{G}[0]] \tilde{G}[0]. \nonumber \end{align} \tag{ 36 }$$

It follows that the only way to include non-linear screening effects is to go beyond the second approximation (24) and (34). 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 [16, 17] in the context of cumulant approximations.

According to the order of the truncation of (15), 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[\varphi]$. The coefficients of the Taylor expansion can in principle be estimated within time-dependent density functional theory (TDDFT), or evaluated self-consistently once the differential equation is solved⁶. The hope is that one recovers the original KBE by climbing the hierarchy of linearized KBEs, namely by considering more and more terms of (15), and that the resulting hierarchy of physical solutions to these equations smoothly converges to the exact Green's function.

Iterating the KBE (1) starting from the (solvable) case $v_{\rm c}=0$ gives rise to the perturbative expansion of the Green's function. Apart from that, there are no standard techniques to solve the equation, preventing us to extract information on the Green's function in its nonperturbative regime. Even worse, the equation admits an infinity of solutions and it is unclear what sort of condition should be added to identify the desired one. In recent years, it has been suggested to reduce the complexity of the problem by looking at the zero-dimensional version of the equation [9]. The equation would retain its differential character, but rather than functionals one would have to deal with simple functions of one variable. In particular, with the substitution $G(1, 2;[\varphi])\to {\rm i}y(x)$, $G_0(1, 2[\varphi])\to y_0(x)$, $\varphi \to -{\rm i}x$, $v_{\rm c}\to {\rm i}v$ the equation reads $y(x) = y_0(x) + y_0(x) (vy(x) + vy(x)y(x){\hspace{0pt}}^{1}) y(x)$. The interacting zero-temperature equilibrium Green's function would then be represented by the value of the function(al) y at vanishing external field x, namely $y(0)$. Such an equation is, however, too simple: its iteration leads to the solution y(0)  =  y₀(0), i.e. the interacting Green's function equals the noninteracting one. A way to add some complexity of the original multidimensional problem is to avoid the perfect cancellation occurring between the Hartree and the exchange-correlation self-energy in zero dimensions. This is achieved by introducing a parameter λ: $y(x) = y_0(x) + y_0(x) vy(x) + \lambda vy(x)y(x){\hspace{0pt}}^{1} y(x)$. We refer to the resulting equation as the one-point model (OPM). Inspired by Hartree–Fock, in which only like spins are affected by exchange interaction, the specific value $\lambda=1/2$ was considered [13]. With such a choice the model can be solved explicitly: Iteration of the equation with $\lambda=1/2$ leads to a perturbative expansion that is non-trivial but still simple enough to be summed up in a closed form, as we shall later see. The resulting function $y(x)$, which represents the OPM version of $G(1, 2;[\varphi])$, is now complex enough to discuss several aspects of the original multidimensional problem. In fact, despite its simplicity, the OPM with $\lambda=1/2$ has proved to be a surprisingly reliable benchmark for qualitative analysis of approximation schemes, as proved by the analyses performed by Berger et al [13] and, later, Stan et al [12], who also showed that the model can be used to study pathologies of perturbation theory such as the misconvergence of the skeleton series [10]. 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 (section 5.1); show the performance of the OPM equivalent of GWA (section 5.2), which will serve as benchmark for our suggested approximations; derive and study the series in the equivalent of $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \Wrpa{G_0}$ and compare it with the series in the exact W (section 5.3); solve the hierarchy of the OPM equivalent of the linearized KBE (section 5.4); and, finally, compare all approaches (section 5.5).

5.1. 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 [13]. At the basis of the model we find the differential equation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y(x)=y_0(x)+y_0(x)\left(-v y(x)+\lambda v y'(x)y(x)^{-1}\right)y(x) \nonumber \end{align} \tag{ 37 }$$

with $\newcommand{\e}{{\rm e}} y_0(x)\equiv y_0^0/(1-y_0^0 x)$, $x\in \mathbb{R}$, $\lambda>0$, $v\geqslant 0$, ${y_{0}^{0}}>0$. For $\lambda=1$, the above equation corresponds to the zero-dimensional version of the equation of motion (1), provided that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G(1,2;[\varphi]) \to {\rm i}y(x) , \nonumber \end{align} \tag{ 38a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G_0(1,2) \to {\rm i}y_0^0 , \nonumber \end{align} \tag{ 38b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \varphi(1) \to -{\rm i}x , \nonumber \end{align} \tag{ 38c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v_{\rm c} \to {\rm i}v . \nonumber \end{align} \tag{ 38d }$$

For this specific value of λ, the exact problem is characterized by the trivial physical solution y(x)  =  y₀(x). A much richer result is obtained by considering $\lambda \neq 1$. Here we take $\lambda=1/2$, namely

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{opmkbe} y(x)=y_0(x)+y_0(x)\left(-v y(x)+\frac{1}{2} v y'(x)y(x)^{-1}\right)y(x). \nonumber \end{align} \tag{ 39 }$$

The choice $\lambda \neq 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 $\lambda =1$. In this sense, $\lambda=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:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{opmpt} y(x)=y_0(x)-\frac{1}{2} v y_0(x){\hspace{0pt}}^3+\frac{1}{4} v^2 y_0(x){\hspace{0pt}}^5-\frac{1}{8} v^3 y_0(x){\hspace{0pt}}^7+\frac{1}{16} v^4 y_0(x){\hspace{0pt}}^9+... \nonumber \end{align} \tag{ 40 }$$

can be summed up to the function

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{opmphys} y_p(x)\equiv \frac{y_0(x)}{\frac{1}{2} v y_0(x){\hspace{0pt}}^2+1} \nonumber \end{align} \tag{ 41 }$$

which is the analytic continuation of the series (40). Since (41) is also solution to (39), we identify this as the physical solution of the problem, all other solutions to (39) being labeled as unphysical ones. Such a function is a model of the functional that in the standard many-body problem realizes the map $\varphi \to G$ by analytic continuation of the series produced by perturbation theory. The Green's function at equilibrium, namely $G[\varphi = 0]$, is then represented by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{exacty} y_p(0)=\frac{y_0^0}{\frac{y_0^0{\hspace{0pt}}^2 v}{2}+1}, \nonumber \end{align} \tag{ 42 }$$

which is plotted in figure 1. The function y₀(x), to which y_p(x) reduces for $v\to 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 $\newcommand{\sh}{{\Sigma_{\rm H}}} \sh$, the exchange-correlation self-energy $\newcommand{\sxc}{\Sigma_{\rm xc}} \sxc$, the screened potential w and the RPA function(al) $w_{\rm rpa}(\,f)$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Sigma \equiv y_0(0)^{-1}-y_p(0)^{-1}=-\frac{y_0^0 v}{2} \nonumber \end{align} \tag{ 43a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\sh}{{\Sigma_{\rm H}}} \displaystyle \sh \equiv -vy_p(0)=-\frac{y_0^0 v}{\frac{y_0^0{\hspace{0pt}}^2 v}{2}+1} \nonumber \end{align} \tag{ 43b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\sh}{{\Sigma_{\rm H}}} \newcommand{\sxc}{\Sigma_{\rm xc}} \displaystyle \sxc \equiv \Sigma -\sh =\frac{y_0^0 v \left(2-y_0^0{\hspace{0pt}}^2 v\right)}{2 y_0^0{\hspace{0pt}}^2 v+4} \nonumber \end{align} \tag{ 43c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p \equiv y_p'(0) \nonumber \end{align} \tag{ 43d }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \epsilon^{-1} \equiv 1-v p \nonumber \end{align} \tag{ 43e }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle w \equiv \epsilon^{-1}v=\frac{v \left(3 y_0^0{\hspace{0pt}}^4 v^2+4\right)}{\left(y_0^0{\hspace{0pt}}^2 v+2\right){\hspace{0pt}}^2}\label{Wexact} \nonumber \end{align} \tag{ 43f }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle W_{\rm rpa}(\,f) \equiv \frac{v}{f^2 v+1}.\label{wrpadef} \nonumber \end{align} \tag{ 43g }$$

Zoom In Zoom Out Reset image size

In the rest of the section we shall consider $\Sigma_{\rm H}$ as given once for all. This is different from the work of Berger where $\Sigma_{\rm 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 $\Sigma_{\rm 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.

5.2. The GW approximation (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 [13], who considered a self-consistently calculated $\newcommand{\sh}{{\Sigma_{\rm H}}} \sh$. As pointed out above, here we make an analogous analysis, but with a $\newcommand{\sh}{{\Sigma_{\rm H}}} \sh$ which is supposed to be given. The self-consistent GWA corresponds to solving the set of equations

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P \approx y^2 \nonumber \end{align} \tag{ 44a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle W \approx v-vPW \nonumber \end{align} \tag{ 44b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\sh}{{\Sigma_{\rm H}}} \displaystyle y \approx {y_{0}^{0}}+{y_{0}^{0}}(\sh+\frac{1}{2}y W)y\label{scGWdyson}. \nonumber \end{align} \tag{ 44c }$$

This set of equations admits three solutions, of which only the following reduces to ${y_{0}^{0}}$ in the noninteracting limit:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{ygw} y_{gw}=\frac{\frac{3\ 3^{2/3} v \left(3 y_0^0{\hspace{0pt}}^2 v \left(y_0^0{\hspace{0pt}}^2 v \left(y_0^0{\hspace{0pt}}^2 v-8\right)-12\right)-16\right)}{\sqrt[3]{B}}+3 \sqrt[3]{3} \sqrt[3]{B}+9 y_0^0 v \left(y_0^0{\hspace{0pt}}^2 v+2\right)}{18 v \left(3 y_0^0{\hspace{0pt}}^2 v+2\right)} \nonumber \end{align} \tag{ 45 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle A \equiv \sqrt{3 y_0^0{\hspace{0pt}}^2 v \left(3 y_0^0{\hspace{0pt}}^2 v \left(y_0^0{\hspace{0pt}}^2 v \left(3 y_0^0{\hspace{0pt}}^2 v \left(y_0^0{\hspace{0pt}}^2 v \left(2 y_0^0{\hspace{0pt}}^2 v+13\right)+80\right)+488\right)+416\right)+464\right)+256} \nonumber \end{align} \tag{ 46 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq47} B \equiv &\, 2 \sqrt{3} A v^{3/2} \left(3 y_0^0{\hspace{0pt}}^2 v+2\right)+9 y_0^0{\hspace{0pt}}^9 v^6+216 y_0^0{\hspace{0pt}}^7 v^5\nonumber \\ &+648 y_0^0{\hspace{0pt}}^5 v^4+576 y_0^0{\hspace{0pt}}^3 v^3+144 y_0^0 v^2. \nonumber \end{align} \tag{ 47 }$$

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:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_i = y_i^2 \nonumber \end{align} \tag{ 48a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle W_i = v-vP_iW_i \nonumber \end{align} \tag{ 48b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\sh}{{\Sigma_{\rm H}}} \displaystyle y_{i+1} ={y_{0}^{0}}+{y_{0}^{0}}(\sh+\frac{1}{2}y_i W_i)y_{i+1} \nonumber \end{align} \tag{ 48c }$$

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 $\Sigma_{\rm 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 (44c) with the exact screened potential W  =  w of (43f). Of the two solutions, the only one that goes to zero in the limit $v \to 0$ is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq49} y_{\rm gwexact}=\frac{\left(y_0^0{\hspace{0pt}}^2 v+2\right) \left(3 y_0^0{\hspace{0pt}}^2 v-\sqrt{-6 y_0^0{\hspace{0pt}}^6 v^3+9 y_0^0{\hspace{0pt}}^4 v^2+4 y_0^0{\hspace{0pt}}^2 v+4}+2\right)}{y_0^0 v \left(3 y_0^0{\hspace{0pt}}^4 v^2+4\right)}. \nonumber \end{align} \tag{ 49 }$$

In figure 2, the OPM equivalent of the self-consistent GWA and self-consistent GW with W exact are shown in comparison with the equivalent of $G_0W_0$. We include the case derived here, where the self-consistency is performed keeping v_H fixed where v_H is calculated self-consistently. We notice that this latter, which obeys to

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \displaystyle \label{scGW} \left\{\begin{array}{@{}l@{}} y_{\rm sc}=y_{\rm H}+y_{\rm H}\Sigma_{GW}[y_{\rm sc}]y_{\rm sc} \nonumber \\ \Sigma_{GW}[y]=\frac{1}{2}\Wrpa{y_{\rm sc}}y_{\rm sc} \nonumber \\ y_{\rm H}=y_0^0+y_0^0(-v y_{\rm H})y_{\rm H} \end{array}\right. \;\;\;{\rm and}\;\;\; \lim_{v\to 0}y_{\rm sc}=y_0^0 \nonumber \end{align} \tag{ 50 }$$

and was studied in [13], leads to the best of all GW-like approximations.

Zoom In Zoom Out Reset image size

5.3. Perturbation theory in W₀ and W

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 3, we can rewrite it in terms of ${y_{0}^{0}}$ and $\newcommand{\e}{{\rm e}} w_{\rm rpa}\equiv W_{\rm rpa}({y_{0}^{0}})$ by simply using the definition of w_rpa to express v, which leads to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{ywrpa} y_p(0)=\frac{y_0^0}{1-\frac{y_0^0{\hspace{0pt}}^2 w_{{\rm rpa}}}{2 \left(y_0^0{\hspace{0pt}}^2 w_{{\rm rpa}}-1\right)}}=y_0^0-\frac{y_0^0{\hspace{0pt}}^3}{2}w_{{\rm rpa}}-\frac{1}{4} y_0^0{\hspace{0pt}}^5 w_{{\rm rpa}}^2-\frac{1}{8} y_0^0{\hspace{0pt}}^7 w_{{\rm rpa}}^3+... \nonumber \end{align} \tag{ 51 }$$

Similarly, we can express Σ and $\newcommand{\sxc}{\Sigma_{\rm xc}} \sxc$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{sigmawrpa} \Sigma = \frac{y_0^0 w_{{\rm rpa}}}{2 y_0^0{\hspace{0pt}}^2 w_{{\rm rpa}}-2}=-\frac{y_0^0}{2}w_{{\rm rpa}}-\frac{1}{2} y_0^0{\hspace{0pt}}^3 w_{{\rm rpa}}^2-\frac{1}{2} y_0^0{\hspace{0pt}}^5 w_{{\rm rpa}}^3-\frac{1}{2} y_0^0{\hspace{0pt}}^7 w_{{\rm rpa}}^4+... \nonumber \end{align} \tag{ 52 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\sxc}{\Sigma_{\rm xc}} \displaystyle \sxc = \frac{y_0^0 w_{{\rm rpa}} \left(2-3 y_0^0{\hspace{0pt}}^2 w_{{\rm rpa}}\right)}{2 y_0^0{\hspace{0pt}}^2 w_{{\rm rpa}} \left(y_0^0{\hspace{0pt}}^2 w_{{\rm rpa}}-3\right)+4}=\frac{y_0^0}{2}w_{{\rm rpa}}-\frac{1}{4} y_0^0{\hspace{0pt}}^5 w_{{\rm rpa}}^3-\frac{3}{8} y_0^0{\hspace{0pt}}^7 w_{{\rm rpa}}^4+...\label{sxcg0w0} \nonumber \end{align} \tag{ 53 }$$

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 $\newcommand{\sxc}{\Sigma_{\rm xc}} \sxc$ 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.

Zoom In Zoom Out Reset image size

If the starting point of the GWA is not the completely noninteracting Green's function, but rather the Hartree one, here defined by $y_{\rm H}=y_{\rm H}+y_0^0(-vy_{\rm H})y_{\rm H}$, then the expansion is performed in $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \newcommand{\e}{{\rm e}} w_{\rm H}\equiv \Wrpa{y_{\rm H}}$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y_p(0) = y_{\rm H}+\frac{1}{2} w_{\rm H} y_{\rm H}^3+\frac{1}{4} w_{\rm H}^2 y_{\rm H}^5-\frac{7}{8} w_{\rm H}^3 y_{\rm H}^7-\frac{63}{16} w_{\rm H}^4 y_{\rm H}^9+... \nonumber \end{align} \tag{ 54 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Sigma = -\frac{w_{\rm H} y_{\rm H}}{2}-w_{\rm H}^2 y_{\rm H}^3-2 w_{\rm H}^3 y_{\rm H}^5-4 w_{\rm H}^4 y_{\rm H}^7+\frac{87}{32} w_{\rm H}^5 y_{\rm H}^9+... \nonumber \end{align} \tag{ 55 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Sigma_{\rm xc} = \frac{w_{\rm H} y_{\rm H}}{2}-w_{\rm H}^3 y_{\rm H}^5-3 w_{\rm H}^4 y_{\rm H}^7+\frac{119}{32} w_{\rm H}^5 y_{\rm H}^9+...\label{pertWh} \nonumber \end{align} \tag{ 56 }$$

which also shows a good convergence, as reported in figure 4.

Zoom In Zoom Out Reset image size

Finally, we would like to reply to the question: is the exact W a better parameter than w_rpa? Even though for real systems the exact W is not accessible, an affirmative reply to this question would justify the search for better approximations by incorporating vertex corrections into W. To find an answer within our model, we have to invert the expression (43f) to get the function corresponding to the map w, the exact screening of the OPM, to v. There are three solutions to the equation for v, of which only one goes to 0 when W goes to 0, namely

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v=\frac{y_0^0{\hspace{0pt}}^2 w \left(\sqrt[3]{36 \sqrt{3} A+B}+y_0^0{\hspace{0pt}}^2 w+36\right)+\left(36 \sqrt{3} A+B\right)^{2/3}-36}{9 y_0^0{\hspace{0pt}}^2 \sqrt[3]{36 \sqrt{3} A+B}} \nonumber \end{align} \tag{ 57 }$$

with $\newcommand{\e}{{\rm e}} A\equiv \sqrt{y_0^0{\hspace{0pt}}^2 w \left(y_0^0{\hspace{0pt}}^2 w \left(2 y_0^0{\hspace{0pt}}^2 w+83\right)-36\right)+12}$ and $\newcommand{\e}{{\rm e}} B\equiv w y_0^0{\hspace{0pt}}^2 (432 + w y_0^0{\hspace{0pt}}^2 (54 + w y_0^0{\hspace{0pt}}^2))$. Provided with such a map, we can calculate again the expansion of y_p(0), $\newcommand{\sxc}{\Sigma_{\rm xc}} \sxc$ and Σ to be

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{sigmawexact} y_p(0) = y_0^0-\frac{y_0^0{\hspace{0pt}}^3 w}{2}-\frac{y_0^0{\hspace{0pt}}^5 w^2}{4}+\frac{y_0^0{\hspace{0pt}}^7 w^3}{8}+\frac{13 y_0^0{\hspace{0pt}}^9 w^4}{16}+\frac{43 y_0^0{\hspace{0pt}}^{11} w^5}{32}+... \nonumber \end{align} \tag{ 58 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Sigma = -\frac{y_0^0 w}{2}-\frac{y_0^0{\hspace{0pt}}^3 w^2}{2}-\frac{y_0^0{\hspace{0pt}}^5 w^3}{4}+\frac{5 y_0^0{\hspace{0pt}}^7 w^4}{8}+\frac{33 y_0^0{\hspace{0pt}}^9 w^5}{16}+... \nonumber \end{align} \tag{ 59 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\sxc}{\Sigma_{\rm xc}} \displaystyle \sxc = \frac{y_0^0 w}{2}-\frac{y_0^0{\hspace{0pt}}^5 w^3}{2}-y_0^0{\hspace{0pt}}^7 w^4-\frac{5 y_0^0{\hspace{0pt}}^9 w^5}{8}+... \nonumber \end{align} \tag{ 60 }$$

Truncations of the above series are plotted in figure 6. 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 figure 5 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.4. 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 4. To this end, we define the the 'perturbed' Hartree potential as $\newcommand{\e}{{\rm e}} V_{\rm H}(x)\equiv -vy(x)$, which leads to rewrite (39) as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{opmkbe2} y(x)={y_{0}^{0}} +{y_{0}^{0}} \left(x+V_{\rm H}(x)+\frac{1}{2}vy'(x)y(x)^{-1} \right) y(x), \nonumber \end{align} \tag{ 61 }$$

and consider the hierarchy of approximations generated by truncating its Taylor expansion in the perturbing potential

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{vhseries} V_{\rm H}(x)=V_0+V_1 x+V_2 x^2+.... \nonumber \end{align} \tag{ 62 }$$

The leading term V₀ simply corresponds to the $\newcommand{\sh}{{\Sigma_{\rm H}}} \sh$ previously defined and assumed to be known, while the other terms involve derivatives of $y(x)$ evaluated at zero: V₁  =  −vy'(0), V₂  =  −vy"(0)/2, etc which are not known a priori. Just like for V₀, 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 (61) and its derivatives. This leads to:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{vhy0} V_1 =\frac{2 y_0^0{\hspace{0pt}}^3 v-6 y_0^0{\hspace{0pt}}^2 v y(0)+4 y_0^0-4 y(0)}{y_0^0{\hspace{0pt}}^3 v+2 y_0^0} \nonumber \end{align} \tag{ 63a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{vhy1} V_2 =\frac{y_0^0 \left(v \left(3 y_0^0 y(0) \left(y_0^0{\hspace{0pt}}^2 v \left(y_0^0{\hspace{0pt}}^2 v-2\right)-4\right)-2 y(0)^2 \left(y_0^0{\hspace{0pt}}^2 v+2\right) \left(3 y_0^0{\hspace{0pt}}^2 v+2\right)+2 y_0^0{\hspace{0pt}}^2 \left(3 y_0^0{\hspace{0pt}}^2 v+8\right)\right)+8\right)-8 y(0)}{y_0^0{\hspace{0pt}}^2 v \left(y_0^0{\hspace{0pt}}^2 v+2\right){\hspace{0pt}}^2} \nonumber \end{align} \tag{ 63b }$$

$$\begin{align} &... \end{align} \tag{ 63c }$$

In previous works [9, 13], the zeroth and first order approximations to V_H, namely $V_{\rm H}\approx V_0$ and $V_{\rm H}(x)\approx V_0+xV_1$, were considered, although, as pointed out earlier, in slightly different versions from what we are now studying⁷.

We will now establish a general procedure that solves the equation for an arbitrarily high order truncation of (62). We start by noticing that, since $V_{\rm H}(x)$ is a given function, the equation (61) is in the linear form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{lindiffeqn} y'(x)=\alpha(x)+\beta(x) y(x). \nonumber \end{align} \tag{ 64 }$$

For $\alpha(x)=0$, the family of solutions of the above equation is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y(x)=C{\rm e}^{\int \beta(x) {\rm d}x}. \nonumber \end{align} \tag{ 65 }$$

This suggests to use the ansatz $y(x)=\bar{y}(x){\rm e}^{\int \beta(x) {\rm d}x}$ when $\alpha(x)\neq 0$. When inserted in (64), we obtain the following equation for $\bar{y}(x)$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \bar{y}'(x)=\alpha(x){\rm e}^{-\int \beta(x) {\rm d}x}, \nonumber \end{align} \tag{ 66 }$$

which is solved by the family of functions

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \bar{y}_C(x)=C+\int_a^x \alpha(x'){\rm e}^{-\int \beta(x') {\rm d}x'}{\rm d}x', \nonumber \end{align} \tag{ 67 }$$

where the lower bound is a fixed constant, while C parametrizes the family, leading to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{sollineqdiff} y_C(x)={\rm e}^{\int \beta(x) {\rm d}x}\left(C+\int_a^x \alpha(x'){\rm e}^{-\int \beta(x') {\rm d}x'}{\rm d}x'\right). \nonumber \end{align} \tag{ 68 }$$

The value of a must be chosen carefully once the integrand $\beta(x)$ is specified, in order to avoid possible divergences of the integral. When applied to (61), we have the exact result

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \beta(x)=-\frac{2 (y_0^0 {V_{\rm H}}(x)+y_0^0 x-1)}{y_0^0 v}\;\;\;\;{\rm and}\;\;\;\;\alpha(x)=-\frac{2}{v}. \nonumber \end{align} \tag{ 69 }$$

Since a has the dimension of a potential, we chose $a={y_{0}^{0}}{\hspace{0pt}}^{-1}$, which will turn out to create no problems of convergence of the integral. We finally arrive at

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{opmansatz} y_C(x)={\rm e}^{-f(x)}\left(C-\frac{2}{v}\int_{{y_{0}^{0}}{\hspace{0pt}}^{-1}}^{x}{\rm d}x'{\rm e}^{\,f(x')}\right) \nonumber \end{align} \tag{ 70 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f(x)&=\int {\rm d}x\left(\frac{2 \left(y_0^0 \left(V_0+V_1 x+V_2 x^2+...\right)+y_0^0 x-1\right)}{y_0^0 v}\right)\nonumber \\ & =\frac{2 \left(y_0^0 V_0 x+\frac{1}{2} y_0^0 V_1 x^2+\frac{1}{3} y_0^0 V_2 x^3+...+\frac{y_0^0 x^2}{2}-x\right)}{y_0^0 v} \nonumber \end{align} \tag{ 71 }$$

which is the family of solutions to (61), 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₀(x) in the limit $v\to 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 (61), but it can still depend on v and $y_0^0$ via the dimensionless quantity $v{y_{0}^{0}}{\hspace{0pt}}^2$; it follows that any choice of such a dependence such that ${\rm e}^{-f(x)}C\to 0$ for $v\to 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 (62) corresponds to a different approximation of the original equation (61), we label the corresponding physical solution with an index $(n)$, where n is the order of the truncation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{opmypnx} y^{(n)}_p(x)=&\,-\frac{2}{v} {\rm e}^{-\frac{2 \left(y_0^0 V_0 x+\frac{1}{2} y_0^0 V_1 x^2+\frac{1}{3} y_0^0 V_2 x^3+...+\frac{1}{n+1} y_0^0 V_n x^{n+1}+\frac{y_0^0 x^2}{2}-x\right)}{y_0^0 v}}\nonumber \\ &\times\int_{{y_{0}^{0}}{\hspace{0pt}}^{-1}}^{x}{\rm d}x' {\rm e}^{\frac{2 \left(y_0^0 V_0 x'+\frac{1}{2} y_0^0 V_1 x'^2+\frac{1}{3} y_0^0 V_2 x'^3+...+\frac{1}{n+1} y_0^0 V_n x'^{n+1}+\frac{y_0^0 x'^2}{2}-x'\right)}{y_0^0 v}}, \nonumber \end{align} \tag{ 72 }$$

which at equilibrium simplifies to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{opmypn} y^{(n)}_p(0)={y_{0}^{0}}\int_{0}^{1}\!{\rm d}\xi\; \frac{2}{v{y_{0}^{0}}{\hspace{0pt}}^2}\; {\rm e}^{\frac{2 \left(y_0^0 V_0 \xi+\frac{1}{2} V_1 \xi^2+\frac{1}{3} y_0^0{\hspace{0pt}}^{-1} V_2 \xi^3+...+\frac{1}{n+1} y_0^0{\hspace{0pt}}^{1-n} V_n \xi^{n+1}+\frac{\xi^2}{2}-\xi\right)}{v{y_{0}^{0}}{\hspace{0pt}}^2}}. \nonumber \end{align} \tag{ 73 }$$

Exact expressions for the integral in (73) 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 (73) only gives an implicit definition of $y^{(n)}_p(0)$, since V₁, V₂, ... also depend on $y^{(n)}_p(0)$ via ((63a), (63b),...) where $y(0)$ has to be replaced by $y^{(n)}_p(0)$. Although already for n  =  1 we cannot find a closed, explicit form for $y^{(n)}_p(0)$, these functions can be evaluated numerically⁸. In figure 7 the first six functions $y^{(n)}_p(0)$, namely $n=0, 1, 2, 3, 4, 5$, are shown as function of v in comparison with the exact physical solution y_p(0) (equation (73)).

Zoom In Zoom Out Reset image size

5.5. Comparing all approximations

A comparison between all different approximations is reported in figure 8. For the expansion in $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \Wrpa{G_{\rm H}}$ we consider a $\newcommand{\sxc}{\Sigma_{\rm xc}} \sxc$ approximated with the first nontrivial two terms after the one prescribed by $G_0W_0$ (which are of third and fourth order in $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \Wrpa{y_{\rm H}}$, the second being identically zero). We also report the second order expansion of y_p(0) in $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \Wrpa{y_0^0}$, which in this model for certain regimes is a better estimate than the expansion of $\Sigma_{\rm xc}$ in $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \Wrpa{y_{\rm H}}$ at lower computational cost, since it involves a lower number of diagrams. For the linearized KBE, we consider the simplest approximation $V_{\rm H}(x)\approx V_{\rm H}(0)$. We see that the expansion in $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \Wrpa{G_{\rm H}}$ leads to an improvement over $G_0W_0$ and self-consistent GW. However, in the strongly coupled regime, the best approximation is given by the solution to the linearized KBE $y^{(0)}_p(0)$.

Zoom In Zoom Out Reset image size

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_0W_0$ can be regarded as the leading order expansion of $\newcommand{\sxc}{\Sigma_{\rm xc}} \sxc$ in G₀ and $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \newcommand{\e}{{\rm e}} W_0\equiv \Wrpa{G_0}$. Since there are reasons to believe that W₀ 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_{\rm H}[\varphi]$ in the KBE. We consider the hierarchy of linearized KBEs arising by higher and higher order approximations to $V_{\rm H}[\varphi]$.

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 $\newcommand{\Wrpa}[1]{\mathcal{W}_{\rm RPA}[#1]} \Wrpa{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.

For a real system, the first approach presents a theoretically clear path to follow. The way to build the hierarchy of approximations is well defined and, being an expansion in the bare Green's function, there are no problems of spurious solutions. The only limitation is the computational cost of the approach, which is the same as that of an expansion in G₀ and v_c. The second approach is more sophisticated from a mathematical perspective, but at the same time more clear from a physics point of view: the only approximation is made in the response of the Hartree potential to a perturbing potential. However, the approach is theoretically incomplete, for a lack of mathematical tools for a systematic treatment of the pertinent equations. Our findings in the OPM may give us some guidance, but a complete picture of the procedure for a real system is still to be achieved. Nonetheless, our preliminary tests show such a striking improvement compared with competing perturbative methods that more work in this direction is certainly motivated.

The research leading to these results has received funding from the European Research Council under the European Unions Seventh Framework Programme (FP/20072013) / ERC Grant Agreement no. 320971. BSM acknowledges the Laboratoire des Solides Irradiés (Ecole Polytechnique, Palaiseau, France) for the support and hospitality during a sabbatical year. BSM acknowledges partial support from CONACYT-México Grant 153930.