European Theoretical Spectroscopy Facility (ETSF) European Theoretical
Spectroscopy Facility (ETSF) European Theoretical Spectroscopy Facility (ETSF)

# Total energy beyond $GW$: exact results and guidelines for approximations

Abdallah El-Sahili  LSI, CNRS, CEA/DRF/IRAMIS, École Polytechnique, Institut
Polytechnique de Paris, F-91120 Palaiseau, France  abdallah.el-
<EMAIL_ADDRESS>Francesco Sottile  LSI, CNRS, CEA/DRF/IRAMIS, École
Polytechnique, Institut Polytechnique de Paris, F-91120 Palaiseau, France
Lucia Reining  LSI, CNRS, CEA/DRF/IRAMIS, École Polytechnique, Institut
Polytechnique de Paris, F-91120 Palaiseau, France

(August 28, 2024)

###### Abstract

The total energy and electron addition and removal spectra can in principle be
obtained exactly from the one-body Green’s function. In practice, the Green’s
function is obtained from an approximate self-energy. In the framework of
many-body perturbation theory, we derive different expressions that are based
on an approximate self-energy, but that yield nevertheless in principle the
exact exchange-correlation contribution to the total energy for any
interaction strength. Response functions play a crucial role, which explains
why, for example, ingredients from time-dependent density functional theory
can be used to build these approximate self-energies. We show that the key
requirement for obtaining exact results is the consistent combination of
ingredients. Also when further approximations are made, as it is necessary in
practice, this consistency remains the key to obtain good results. All
findings are illustrated using the exactly solvable symmetric Hubbard dimer.

The Dyson equation with an approximate self-energy leads to an approximate
Green’s function and to approximate spectra. However, the same approximate
self-energy yields the exact exchange-correlation energy, provided the latter
is evaluated following a consistent prescription.

##  1 Introduction

Many important properties of materials are linked to observables that can be
expressed in principle as expectation values in the many-body ground state or
in thermal equilibrium. In practice, the use of many-body wavefunctions is
often avoided by rather describing the observables as functionals of more
compact quantities, i.e., quantities that depend on less arguments, such as
the density 1, one-body reduced density matrix 2, 3, 4, or one- or two-body
Green’s functions 5. This represents a trade-off: often, one does not know the
exact functional for an observable in terms of these quantities, and
approximations have to be designed. One important example is the total energy:
it can be straightforwardly formulated in terms of the one-body Green’s
function (GF) 6, whereas no exact explicit expression in terms of the density
is known. The same holds for electron addition and removal spectral functions.
Excitation spectra involving neutral excitations in linear response, instead,
are easily expressed in terms of a two-body Green’s function, but not in terms
of the ground state density or the one-body Green’s function 5. Even when the
expressions are known, one faces another problem: while the use of the compact
quantities carries the promise of reduced computational load, they are
themselves only known explicitly as expectation values involving many-body
wavefunctions. Therefore, nothing is gained, unless one finds ways to
calculate them in a different way, which may be in principle exact, and in
practice, require approximations. Typically, the density is obtained from the
Kohn-Sham equations 7 with an approximate exchange-correlation (xc) potential,
and the GF, from a Dyson equation with an approximate xc self-energy
$\Sigma_{\rm xc}$ 8. It is therefore not always obvious which framework
(Density Functional Theory (DFT), Green’s Function Functional Theory, etc.) is
the best choice to access a given observable.

This holds in particular for the total ground state energy $\text{E}_{0}$.
While the Galitskii-Migdal formula 6 or functional expressions such as the
Luttinger-Ward 9 or Klein 10 functionals yield an in principal exact and, in
the latter two cases, even variational, form in terms of the Green’s function
$G$, the need for approximations to the GF itself strongly impacts the quality
of the results. Together with a computational load far heavier than that of
the most widely used density functionals, this explains why the vast majority
of total energy calculations is performed using DFT, not Green’s functions.
Still, research concerning total energy calculations using Green’s functions
is active and important 11, 12, 13, 14. Besides the - important - fact that in
principle expressions for $\text{E}_{0}$ as functional of $G$ and/or
$\Sigma_{\rm xc}$ are known, the Green’s functions framework benefits from the
existence of powerful approximations. In particular, Many-Body Perturbation
Theory (MBPT) 8 suggests a way to expand the self-energy in diagrams that
carry physical meaning and that are therefore helpful to describe phenomena
such as the van der Waals dispersion interaction 15. For situations that show
only weak to moderate interaction effects, MBPT is often considered to be a
systematic way to proceed, although in practice renormalizations, such as
screening of the Coulomb interaction, are needed. In particular, even the
lowest order of an expansion of the self-energy in terms of the screened
Coulomb interaction $W$, which is the widely used $GW$ approximation 16, has
been very successful for the calculation of the quasi-particle (QP) part of
electron addition and removal spectra in finite and extended systems 17, 18,
19, 20, 21, 22, 23, 24, 5.

However, there are many cases where the $GW$ approximation is not sufficient.
On one hand, the description of QP energies is not always good and certainly
worsens in more strongly correlated systems 25, 5. On the other hand, other
quantities, such as satellite features in the electron addition and removal
spectra, are often less well described, even in absence of strong correlation
26, 27. Most importantly, $GW$ does not necessarily yield total energies of
better quality than currently used density functionals 14, 28, 13. Research on
total energies in terms of GFs goes therefore hand in hand with research on
approximations to the self-energy beyond $GW$. The most straightforward way to
go would be to explore higher orders in $W$, and important research in this
direction is ongoing 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39. In many cases
it cannot, however, bring a practical solution, since the resulting
expressions become quickly very cumbersome and costly, and since perturbation
theory will diverge when the interaction is too strong. Therefore, it would be
desirable to find an efficient way to terminate the perturbation series.

In the various possible ways to express the xc energy contribution to the
total energy such as using the adiabatic connection fluctuation dissipation
theorem 40, 41, the polarizability plays a key role. This suggests to explore
links to other frameworks that are used to access the polarizability, in
particular, Time-Dependent DFT (TDDFT) 42. Indeed, there is a long, and
sometimes very successful, history of attempts to use TDDFT in order to go
beyond $GW$ in terms of vertex corrections based on the xc kernel $f_{\rm xc}$
43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, the functional derivative with
respect to the density of the xc potential $v_{\rm xc}$ of TDDFT 55, or
related linear response kernels that may be closer to the many-body Green’s
functions framework 56, 57, 58, 59. This kind of combination leads to the so-
called $G\tilde{W}$ self-energy, where the Coulomb interaction is screened by
a test charge-test electron (TCTE) dielectric function instead of the test
charge-test charge (TCTC) one used in the $GW$ approximation 60, 61.
Independently of the specific recipe that is used in the various $G\tilde{W}$
expressions, these approaches replace the complicated exact vertex function
$\Gamma$ that depends on three space, spin and time arguments by a two-
arguments function $(1-f_{\rm xc}\chi_{0})^{-1}$, where $\chi_{0}$ is an
independent-particle polarizability. Therefore, the resulting self-energy is
always approximate 49. Nevertheless, using a $G\tilde{W}$ self-energy instead
of $GW$ often improves the QP energies 43, 52, 54. At the same time, the idea
is much less explored when it comes to total energies 52, 62. Moreover, to the
best of our knowledge a systematic study for both total energy and spectra
that would discern the effect of the replacement of the full $\Gamma$ by a
two-arguments vertex on one side, from the effect of approximations to the
$f_{\rm xc}$ itself on the other side, is still missing.

The present work has a focus on the total xc energy, while making a link to
other aspects of the GF when interesting. It addresses the following
questions: Could a self-energy with a two-arguments vertex correction, and in
particular, a TDDFT-derived one, yield in principle exact results? If yes, how
do we have to build the corresponding expressions for the total xc energy? How
do widely used approximations impact the results? And what happens to the
kinetic energy and to spectra?  As we will show, there are indeed different
possibilities to obtain in principle exact expressions for the total xc
energy, which are moreover quite robust when widely used approximations are
made. Consistent combination of ingredients is a key requirement for this to
be true. With these self-energies, the kinetic energy is not exact in
principle, but we examine the possibility to use the virial theorem in order
to overcome this issue. This allows us moreover to make an interesting
comparison to the widely used adiabatic connection approach, which also makes
use of the polarizability, but without involving a self-energy. Spectra are
also approximate in principle when a two-arguments vertex correction is used,
but we find that the $G\tilde{W}$ results still exhibit improvements over
$GW$.

Our investigation and discussion is general, and it is accompanied by an
illustration using the exactly solvable symmetric Hubbard dimer at half-
filling. The paper is organized as follows: the theoretical background is
given in Sec. 2. Theoretical developments made on this basis are presented and
discussed in Sec. 3. The results for the Hubbard dimer are contained in Sec.
4. Conclusions are summarized in Sec. 5.

##  2 Theoretical background

###  2.1 Total energy and spectral function in terms of the GF

The ground state total energy $\text{E}_{0}$ can be expressed in terms of the
time ordered GF 6, 13

$\mathrm{E}_{0}=\underbrace{-i\lim_{t_{2}\rightarrow t^{+}_{1}}\int
dx_{1}\,\bigg{[}-\frac{\nabla^{2}_{r_{1}}}{2}+v_{\mathrm{ext}}(x_{1})\bigg{]}G(x_{1},x_{1};t_{1}-t_{2})}_{\mathrm{E}_{\mathrm{k}}+\mathrm{E}_{\mathrm{ext}}}\\\
\underbrace{-\frac{i}{2}\lim_{t_{2}\rightarrow
t^{+}_{1}}\int{dx_{1}\,v_{\mathrm{H}}(x_{1})G(x_{1},x_{1};t_{1}-t_{2})}}_{\mathrm{E}_{\mathrm{H}}}\\\
\underbrace{-\frac{i}{2}\lim_{t_{2}\rightarrow t^{++}_{1}}\int
dx_{1}dx_{3}dt_{3}\,\Sigma_{\mathrm{xc}}(x_{1},x_{3};t_{1}-t_{3})G(x_{3},x_{1};t_{3}-t_{2})}_{\mathrm{E}_{\mathrm{xc}}}\,,$
(1)

where $x=(r,\sigma)$ stands for position and spin, and $\displaystyle
t^{+}\equiv\lim_{\eta\to 0^{+}}\,(t+\eta)$. Here, we have highlighted the
different contributions to the total energy, namely the kinetic energy
$\mathrm{E}_{\mathrm{k}}$, the contribution $\mathrm{E}_{\mathrm{ext}}$ coming
from the external potential $v_{\mathrm{ext}}$, the Hartree energy
$\mathrm{E}_{\mathrm{H}}$ given in terms of the Hartree potential
$v_{\text{H}}$, and the exchange-correlation energy $\mathrm{E}_{\mathrm{xc}}$
expressed in terms of the exchange-correlation self-energy
$\Sigma_{\mathrm{xc}}$. The last two terms compose the interaction energy
$\mathrm{E}^{\mathrm{inter}}\equiv\mathrm{E}_{\mathrm{H}}+\mathrm{E}_{\mathrm{xc}}$.
Note that here, in the context of MBPT, $\mathrm{E}_{\mathrm{xc}}$ refers
specifically to the exchange-correlation energy of the Coulomb interaction, in
contrast to the DFT framework where the xc energy also includes the
correlation contribution from the kinetic energy. The specific form Eq. (1) of
the Galitskii-Migdal equation is convenient in order to discuss separately the
different contributions to the total energy, and to find specific improvements
for each part. While such a strategy does not benefit from error canceling and
therefore does not necessarily lead to globally improved results, it helps to
obtain deeper insight, and eventually to arrive to the good result for the
good reason.
The main quantity of interest is $\Sigma_{\text{xc}}$, which can be expressed
exactly as

$\Sigma_{\mathrm{xc}}(1,2)=i\int{d(34)G(1,4)W(3,1^{+})\Gamma(4,2,3)}\,,$ (2)

where $1=(x_{1},t_{1})=(r_{1},\sigma_{1},t_{1})$ stands for position, spin and
time. The screened Coulomb interaction $W$ is given by

$W(1,2)=v_{c}(1,2)+\int d(34)\,v_{c}(1,3)\chi(3,4)v_{c}(4,2)\,,$ (3)

with $v_{c}$ the bare Coulomb interaction
$v_{c}(1,2)=\delta(t_{1}-t_{2})\frac{1}{|r_{1}-r_{2}|}$ and $\chi$ the
reducible polarizability. The main complication stems from the vertex function
$\Gamma$, defined as

$\Gamma(4,2,3)=\delta(4,3)\delta(2,3)+\frac{\delta\Sigma_{\text{xc}}(4,2)}{\delta
v_{\text{cl}}(3^{++},3^{+})}\,,$ (4)

where $v_{\rm cl}=v_{\text{H}}+v_{\mathrm{ext}}$ is the total classical
potential. Since $\Gamma$ is in turn determined by the self-energy, in most
cases it cannot be expressed in a closed form nor calculated exactly. To
lowest order in the Coulomb interaction
$\Gamma(4,2,3)\approx\delta(4,3)\delta(2,3)$. Corrections to this result are
called vertex corrections. Neglecting vertex corrections one obtains the $GW$
approximation, where $\Sigma_{\mathrm{xc}}(1,2)=iG(1,2)W(2,1^{+})$16.

Once the self-energy is determined in the chosen approximation, the GF is
obtained by solving the Dyson equation

$G(1,2)=G_{0}(1,2)+\int{d(34)}\,G_{0}(1,3)\bigg{(}v_{\mathrm{H}}(3,4)+\Sigma_{\mathrm{xc}}(3,4)\bigg{)}G(4,2)\,,$
(5)

where $G_{0}$ is the non-interacting GF. Finally, the resulting GF can be used
to calculate $\text{E}_{0}$ from Eq. (1) or to evaluate the spectral function
from the frequency Fourier transform of $G$,

$A(x,x,\omega)=\frac{1}{\pi}|\mathrm{Im}\big{(}G(x,x;\omega)\big{)}|\,.$ (6)

###  2.2 Interaction energy in terms of the polarizability

Our focus is to find accurate expressions for the interaction energy
$\mathrm{E}^{\mathrm{inter}}$. For this purpose, it is useful to express it in
terms of the reducible polarizability $\chi$ 62.

For a system with $N$ electrons in its ground state, the interaction energy is
given by the expectation value of the Coulomb interaction operator $\hat{V}$
in the many-body ground state $\ket{N_{0}}$,

$\mathrm{E}^{\mathrm{inter}}=\braket{N_{0}}{\hat{V}}{N_{0}}=\frac{1}{2}\int{dx_{1}dx_{2}}\,v_{\mathrm{c}}(x_{1},x_{2})\bra{N_{0}}\hat{\psi}^{\dagger}(x_{2})\hat{\psi}^{\dagger}(x_{1})\hat{\psi}(x_{1})\hat{\psi}(x_{2})\ket{N_{0}}\,,$
(7)

where $\hat{\psi}$ and $\hat{\psi}^{\dagger}$ are the annihilation and
creation field operators, respectively. On the other hand, the reducible
polarizability $\chi$ is defined as

$\chi(x_{1},t_{1};x_{2},t_{2})=-iG(x_{1},t_{1};x_{1},t_{1}^{+})G(x_{2},t_{2};x_{2},t_{2}^{+})\\\
-i\bra{N_{0}}\hat{T}\big{[}\hat{\psi}^{\dagger}(x_{1},t^{+}_{1})\hat{\psi}(x_{1},t_{1})\hat{\psi}^{\dagger}(x_{2},t_{2}^{+})\hat{\psi}(x_{2},t_{2})\big{]}\ket{N_{0}}\,.$
(8)

In the limit $t_{2}=t_{1}^{++}$,

$\chi(x_{1},t_{1};x_{2},t_{1}^{++})=-iG(x_{1},t_{1};x_{1},t_{1}^{+})G(x_{2},t_{1}^{++};x_{2},t_{1}^{+++})\\\
-i\braket{N_{0}}{\hat{\psi}^{\dagger}(x_{2})\hat{\psi}(x_{1})}{N_{0}}\delta(x_{2}-x_{1})-i\braket{N_{0}}{\hat{\psi}^{\dagger}(x_{2})\hat{\psi}^{\dagger}(x_{1})\hat{\psi}(x_{1})\hat{\psi}(x_{2})}{N_{0}}\\\
=i\,n(x_{1})n(x_{2})-i\,\rho(x_{1},x_{2})\delta(x_{2}-x_{1})-i\braket{N_{0}}{\hat{\psi}^{\dagger}(x_{2})\hat{\psi}^{\dagger}(x_{1})\hat{\psi}(x_{1})\hat{\psi}(x_{2})}{N_{0}}\,,$
(9)

where we used the anticommutation relation
$\\{\psi(x_{2}),\psi^{\dagger}(x_{1})\\}=\delta(x_{2}-x_{1})$, and where we
introduced the one-body reduced density-matrix
$\rho(x_{1},x_{2})=\bra{N_{0}}\hat{\psi}^{\dagger}(x_{2})\hat{\psi}(x_{1})\ket{N_{0}}=-iG(x_{1},t,x_{2},t^{+})$,
with the electron density $n(x_{1})=\rho(x_{1},x_{1})$. The last term in Eq.
(9) enters the definition of the interaction energy in Eq. (7). The
interaction energy can therefore be expressed in terms of the polarizability
$\chi$ as

$\mathrm{E}^{\mathrm{inter}}=\frac{1}{2}\int{dx_{1}dx_{2}}v_{\mathrm{c}}(x_{1},x_{2})n(x_{1})n(x_{2})\\\
+\frac{i}{2}\int{dx_{1}dx_{2}v_{\mathrm{c}}(x_{1},x_{2})\chi(x_{1},t_{1};x_{2},t_{1}^{++})}-\frac{1}{2}\int{dx_{1}dx_{2}v_{\mathrm{c}}(x_{1},x_{2})\rho(x_{1},x_{2})\delta(x_{2}-x_{1})}\,,$
(10)

where the first term is the Hartree energy and the last two terms are the
exchange-correlation energy. This formulation of $\text{E}_{\rm xc}$ is not
directly suitable for practical purposes, since it consists of terms
containing a divergence that cancels in the sum. It is, however, a good
starting point for the developments in the next section.

##  3 Theoretical developments

###  3.1 A freedom of choice

In order to eliminate the problematic last term in Eq. (10), we introduce a
generalized independent-particle polarizability defined as
$\chi_{0}(1,2)\equiv-i\bar{G}(1,2^{+})\bar{G}(\bar{2},1^{+})$. Its time
diagonal is

$\displaystyle\chi_{0}(x_{1},t,x_{2},t^{++})$
$\displaystyle=-i\bar{G}(x_{1},t;x_{2},t^{+++})\bar{G}(x_{2},t^{++};x_{1},t^{+})$
(11)
$\displaystyle=\bar{\rho}(x_{1},x_{2})\bigg{(}-i\bra{\bar{N}_{0}}\hat{\psi}(x_{2})\hat{\psi}^{\dagger}(x_{1})\ket{\bar{N}_{0}}\bigg{)}$
(12)
$\displaystyle=-i\bar{\rho}(x_{1},x_{2})\bra{\bar{N}_{0}}\delta(x_{2}-x_{1})-\hat{\psi}^{\dagger}(x_{1})\hat{\psi}(x_{2})\ket{\bar{N}_{0}}$
(13)
$\displaystyle=-i\bar{\rho}(x_{1},x_{2})\delta(x_{2}-x_{1})+i\bar{\rho}(x_{1},x_{2})\bar{\rho}(x_{2},x_{1})\,,$
(14)

where $\ket{\bar{N}_{0}}$ is the many-body ground state corresponding to a
system that could be the true interacting system or an auxiliary interacting
or non-interacting system. $\bar{G}$ and $\bar{\rho}$ are the corresponding GF
and the corresponding density matrix, respectively.

In the last term of Eq. (10), only the diagonal of the density matrix is
needed. In order to replace this term, we can therefore consider all systems
that yield the exact density $\bar{\rho}(x,x)=n(x)$, such as the true
interacting system, or the Kohn-Sham auxiliary system. This leaves
considerable freedom, which we can use to derive different exact expressions
for $\mathrm{E}^{\mathrm{inter}}$ and to design efficient approximations.
Indeed, when $\bar{\rho}(x,x)=n(x)$ we have, from Eq. (14)

$n(x_{1})\delta(x_{2}-x_{1})=i\chi_{0}(x_{1},t;x_{2},t^{++})+\bar{\rho}(x_{1},x_{2})\bar{\rho}(x_{2},x_{1})\,,$
(15)

which, replaced in Eq. (10), yields

$\displaystyle\text{E}_{\text{xc}}$
$\displaystyle=\begin{multlined}-\frac{1}{2}\int{dx_{1}dx_{2}v_{\mathrm{c}}(x_{1},x_{2})\bar{\rho}(x_{1},x_{2})\bar{\rho}(x_{2},x_{1})}\\\
+\frac{i}{2}\int{dx_{1}dx_{2}}v_{\mathrm{c}}(x_{1},x_{2})\bigg{(}\chi(x_{1},t_{1};x_{2},t_{1}^{++})-\chi_{0}(x_{1},t_{1};x_{2},t_{1}^{++})\bigg{)}\end{multlined}-\frac{1}{2}\int{dx_{1}dx_{2}v_{\mathrm{c}}(x_{1},x_{2})\bar{\rho}(x_{1},x_{2})\bar{\rho}(x_{2},x_{1})}\\\
+\frac{i}{2}\int{dx_{1}dx_{2}}v_{\mathrm{c}}(x_{1},x_{2})\bigg{(}\chi(x_{1},t_{1};x_{2},t_{1}^{++})-\chi_{0}(x_{1},t_{1};x_{2},t_{1}^{++})\bigg{)}$
(18)
$\displaystyle=\bar{\mathrm{E}}_{\mathrm{x}}+\bar{\mathrm{E}}_{\mathrm{c}}=\bar{\mathrm{E}}_{\mathrm{x}}+\mathrm{E}_{\mathrm{c}}+\left(\mathrm{E}_{\mathrm{x}}-\bar{\mathrm{E}}_{\mathrm{x}}\right)\,.$
(19)

The first term in Eq. (18) is $\bar{\mathrm{E}}_{\mathrm{x}}$, the exchange
energy corresponding to $\ket{\bar{N}_{0}}$. Since the derivation shows that
the sum of all terms is the exact exchange-correlation energy, the second term
in Eq. (18) $\bar{\mathrm{E}}_{\mathrm{c}}$ contains the exact correlation
energy plus a correction that compensates the error of
$\bar{\mathrm{E}}_{\mathrm{x}}$ with respect to the exact exchange energy
$\mathrm{E}_{\mathrm{x}}$. _It is crucial to note that one can use any system
defined by a ground state $\ket{\bar{N}_{0}}$, as long as this system yields
the exact density: this will yield the exact Coulomb interaction energy,
although $\bar{\rho}$ is not the density matrix of the true interacting
system._ Two most obvious choices are either the true many-body (MB) system
with $\bar{G}=G$, which leads to $\chi_{0}^{\mathrm{MB}}\to-iGG$ and
$\bar{\rho}(x_{1},x_{2})=\rho(x_{1},x_{2})$ the true density matrix, or the
Kohn-Sham (KS) system with $\chi_{0}\to\chi_{0}^{\mathrm{KS}}\equiv-
iG^{\mathrm{KS}}G^{\mathrm{KS}}$ the independent-particle polarizability built
with the Kohn-Sham Green’s function, and $\bar{\rho}\to\rho^{\text{KS}}$ the
KS density matrix.

###  3.2 Exact exchange-correlation energy from approximate self-energies

Our next goal is to make a self-energy appear in the expression of
$\mathrm{E}_{\mathrm{xc}}$. To this aim, we rewrite Eq. (18) as

$\displaystyle\mathrm{E}_{\mathrm{xc}}$ $\displaystyle=\frac{1}{2}\int
dx_{1}d2\,\bar{G}(1,2^{+})\bar{G}(2,1^{+})v_{\mathrm{c}}(2,1)+\frac{i}{2}\int
dx_{1}d3\,v_{\mathrm{c}}(3,1)\left[\chi(1,3^{++})-\chi_{0}(1,3^{++})\right]$
(20) $\displaystyle=\begin{multlined}\frac{1}{2}\int
dx_{1}d2\,\bar{G}(1,2^{+})\bar{G}(2,1^{+})v_{\mathrm{c}}(2,1)\\\
+\frac{i}{2}\int
dx_{1}(234)\,\chi_{0}(1,2)\left[v_{\mathrm{c}}(2,4)+\bar{f}_{\mathrm{xc}}(2,4)\right]\chi(4,3^{++})v_{\mathrm{c}}(3,1)\,,\end{multlined}\frac{1}{2}\int
dx_{1}d2\,\bar{G}(1,2^{+})\bar{G}(2,1^{+})v_{\mathrm{c}}(2,1)\\\
+\frac{i}{2}\int
dx_{1}(234)\,\chi_{0}(1,2)\left[v_{\mathrm{c}}(2,4)+\bar{f}_{\mathrm{xc}}(2,4)\right]\chi(4,3^{++})v_{\mathrm{c}}(3,1)\,,$
(23)

where we have introduced the generalized exchange-correlation kernel
$\bar{f}_{\mathrm{xc}}$ that, once a choice for $\chi_{0}$ is made, is defined
from the Dyson-like equation

$\chi(1,2)=\chi_{0}(1,2)+\int
d(34)\chi_{0}(1,3)\bigg{(}v_{\mathrm{c}}(3,4)+\bar{f}_{\mathrm{xc}}(3,4)\bigg{)}\chi(4,2)\,,$
(24)

keeping in mind that $\chi$ is always the exact reducible polarizability. When
$\chi_{0}$ is chosen to be the KS independent particle polarizability,
$\bar{f}_{\rm xc}=f_{\rm xc}$, the xc kernel of linear response TDDFT55, but,
as pointed out above, other choices are possible.111Note that here we have
given the equations in terms of time-ordered quantities, whereas TDDFT is
usually causal. One has to pay attention to be consistent when combining the
GFFT and TDDFT frameworks in practice. By using the definition of $\chi_{0}$,
given in the beginning of 3.1, Eq. (23) can be written as

$\displaystyle\mathrm{E}_{\mathrm{xc}}$
$\displaystyle=\int{dx_{1}d2\,\bar{G}(1,2^{+})\bar{G}(2,1^{+})\bigg{(}v_{\mathrm{c}}(2,1)+\int{d(34)\,\big{(}v_{\mathrm{c}}(2,4)+\bar{f}_{\mathrm{xc}}(2,4)\big{)}\chi(4,3^{++})v_{\mathrm{c}}(3,1)}}\bigg{)}$
(25)
$\displaystyle=\int{dx_{1}d2\,\bar{G}(1,2)\bar{G}(2,1^{++})\bigg{(}v_{\mathrm{c}}(2,1^{+})+\int{d(34)\,\big{(}v_{\mathrm{c}}(2,4)+\bar{f}_{\mathrm{xc}}(2,4)\big{)}\chi(4,3^{++})v_{\mathrm{c}}(3,1^{+})}}\bigg{)}$
(26) $\displaystyle=\frac{1}{2}\int
dx_{1}d2\,\bar{G}(1,2)\bar{\tilde{W}}(2,1^{+})\bar{G}(2,1^{++})\,,$ (27)

where we have defined the generalized TCTE screened Coulomb interaction222The
double infinitesimals in $\chi(4,3^{++})$ do not change the spectrum of
$\bar{\tilde{W}}$, but we keep them here explicitly since they give a
straightforward prescription for the contour integral in frequency space
yielding $E_{\rm xc}$.

$\bar{\tilde{W}}(2,1)=v_{\mathrm{c}}(2,1)+\int{d(34)}\bigg{(}v_{\mathrm{c}}(2,4)+\bar{f}_{\mathrm{xc}}(1,4)\bigg{)}\chi(4,3^{++})v_{\mathrm{c}}(3,1)\,.$
(28)

In this way, the exact exchange-correlation energy takes a form analogous to
the last term in the Galitskii-Migdal expression Eq. (1):

$\mathrm{E}_{\mathrm{xc}}=-\frac{i}{2}\int
dx_{1}d2\,\bar{\Sigma}_{\mathrm{xc}}(1,2)\bar{G}(2,1^{++})\,,$ (29)

with an exchange-correlation self-energy

$\bar{\Sigma}_{\mathrm{xc}}(1,2)\equiv
i\bar{G}(1,2)\bar{\tilde{W}}(2,1^{+})\,.$ (30)

The important point to stress here is the fact that the exact
$\mathrm{E}_{\mathrm{xc}}$ is obtained with an approximate self-energy Eq.
(30). This approximation is often called $G\tilde{W}$. It is usually derived
43 by replacing $\Sigma_{\mathrm{xc}}(4,2)$ in the functional derivative in
Eq. (4) with a local $\delta(4,2)\bar{v}_{\rm xc}(4)$. Most often,
$\bar{v}_{\rm xc}\equiv v_{\rm xc}$, the KS xc potential of TDDFT, is chosen
and the resulting $f_{\rm xc}$ is approximated, for example, in the adiabatic
local density approximation. In our derivation, $\bar{f}_{\mathrm{xc}}$ does
not have to be a functional derivative, since it is defined by Eq. (24), which
generalizes the definition of $\bar{\tilde{W}}$. This gives a rigorous
foundation to attempts to use $\bar{f}_{\mathrm{xc}}$ other than approximate
TDDFT ones in order to approximate vertex corrections to the self-energy, in
particular, the so-called nanoquanta kernel and approximations to it:56, 57,
58, 59, 44, 63, 64, 65, 66, 67 the only requirement is that $\chi_{0}$
corresponds to the correct density. It should, however, be noted that a
$\bar{f}_{\mathrm{xc}}$ fulfilling Eq. (24) does not necessarily exist for
every $\chi_{0}$. We will give an illustration below in the Hubbard dimer.

The important message of this section is that _the exact exchange-correlation
energy can be obtained with an approximate self-energy $\bar{\Sigma}_{\rm xc}$
and with an approximate Green’s function $\bar{G}$ which is not the solution
of the Dyson equation using $\bar{\Sigma}_{\mathrm{xc}}$, but which has been
chosen from the beginning. The two important requirements are consistency of
the ingredients used in Eq. (29), and the fact that they stem from the real or
from an auxiliary system yielding the exact density._ In the following, we
will call this a consistent scheme, as opposed to a non-consistent scheme
where different GFs are used in
$\mathrm{E}_{\mathrm{xc}},\bar{\Sigma}_{\mathrm{xc}},\bar{\tilde{W}}$. Here,
we have shown that there is more than one possible consistent choice, which
may help to design efficient approximations.

###  3.3 The kinetic energy

The TCTE screened self-energy $\bar{\Sigma}_{\rm xc}$ does in general not
correspond to the exact self-energy, and therefore one does not have access to
the exact GF nor to the exact density matrix. As a consequence, the kinetic
energy $\mathrm{E}_{\mathrm{k}}$ cannot be computed exactly. However, with the
exact Coulomb interaction energy $\mathrm{E}^{\mathrm{inter}}$ at hand, this
problem can in principle be overcome by using the virial theorem for the
electron system 68, 69,

$2\mathrm{E}_{\mathrm{k}}+\mathrm{E}^{\mathrm{inter}}=\int{d^{3}{\bf
r}\,n({\bf r}){\bf r}\cdot\nabla v_{\mathrm{ext}}(r)}\equiv S_{VT}\,.$ (31)

Using the virial theorem requires in principle to know the exact density. This
is not an additional requirement here, since it was already assumed throughout
the above derivations. Moreover, research in the framework of DFT shows that
errors induced by approximate functionals are often predominantly due to the
form of the functional, whereas in many cases errors due to an approximate
density are small 70. Therefore, the use of the virial theorem is a promising
route to take when, as it is the case here, one can expect to access the
interaction energy with good accuracy.

###  3.4 Comparison to the adiabatic connection

Finally, it is interesting to compare our equations to the adiabatic
connection (AC) approach 40, 41, 71. In principle, this approach yields the
exact full correlation energy, which encompasses correlations arising from
both kinetic and Coulomb interaction energies, as well as the difference
between the exchange energy calculated with the true and the KS density
matrix, respectively.

Also in this approach, the correlation energy is expressed in terms of $\chi$
and $\chi_{0}$, but with an integration over a coupling constant $\lambda$
that scales the Coulomb interaction and modifies $v_{\mathrm{ext}}$ such that
the density remains constant,

$\text{E}^{\text{full}}_{\text{c}}=\frac{i}{2}\int_{0}^{1}{d\lambda}\int{dx_{1}dx_{2}}\,v_{\mathrm{c}}(x_{1},x_{2})\bigg{(}\chi^{\lambda}(x_{1},t_{1};x_{2},t_{1}^{++})-\chi_{0}^{\mathrm{KS}}(x_{1},t_{1};x_{2},t_{1}^{++})\bigg{)}\,.$
(32)

Since the structure of the expression is the same as that of
$\bar{\mathrm{E}}_{\mathrm{c}}$ in Eq. (19), one can express also the AC
result in terms of an effective self-energy,

$\text{E}^{\text{full}}_{\rm
c}=-\frac{i}{2}\int{dx_{1}d3\,\Sigma^{\text{eff}}_{\text{c}}(1,3)G^{\mathrm{KS}}(3,1^{++})},$
(33)

where

$\Sigma^{\text{eff}}_{\rm
c}(1,3)=iG^{\mathrm{KS}}(1,3)\int_{0}^{1}{d\lambda\,\tilde{W}_{\rm
pol}^{\lambda}(3,1^{+})}\,,$ (34)

with $\tilde{W}_{\rm pol}^{\lambda}$ the polarization contribution to the
$\lambda$-dependent KS TCTE screened interaction. $\Sigma^{\text{eff}}_{c}$ is
an effective correlation self-energy that contains kinetic and interaction
contributions. Here, we have worked with the KS scheme, since the AC
expression is often (though not always, see, e.g., 38) used in the framework
of KS-DFT. Analogous expressions are obtained for other allowed choices of
$\chi_{0}$, e.g., stemming from a generalized KS scheme. Of course, this self-
energy yields the exact full correlation energy, while it is not meant to be
used in a Dyson equation to yield the GF.

It is interesting to compare the errors to be expected in practice from the AC
approach on one side, and, with the errors of the approach discussed here,
i.e., the combination of the calculation of $\mathrm{E}^{\mathrm{inter}}$ plus
the use of the virial theorem. For this estimate, we suppose the virial term
$S_{\rm VT}$ in Eq. (31) to be known with an error that is negligible with
respect to the error $\Delta\mathrm{E}^{\mathrm{inter}}$ stemming from
approximations to $\chi$. This is consistent with the fact that we suppose the
density to be known with good accuracy. Using the virial theorem
$2\mathrm{E}_{\mathrm{k}}+\mathrm{E}^{\mathrm{inter}}=S_{\rm VT}$, the error
in the kinetic energy will be
$\Delta\mathrm{E}_{\mathrm{k}}=-\frac{\Delta\mathrm{E}^{\mathrm{inter}}}{2}$,
leading to a total error of
$\Delta\text{E}=+\frac{\Delta\mathrm{E}^{\mathrm{inter}}}{2}$. In the case of
the AC, the error is determined entirely by the integral over response
functions Eq. (32). Since the non-interacting $\chi_{0}$ is subtracted, it is
reasonable to suppose that the dominant contribution is linear in $\lambda$.
Evidence that this is true can be found for small systems in Ref. 72. Assuming
linearity in $\lambda$, one obtains the _same_ error
$\Delta\text{E}=+\frac{\Delta\mathrm{E}^{\mathrm{inter}}}{2}$ as in our
alternative scheme. Whether higher orders in $\lambda$ will rather reduce or
increase this result depends on whether $\chi_{\lambda}$ is convex or concave.
In any case, this discussion suggests that similar errors are to be expected,
while the $\lambda$-integration is avoided in the approach using the virial
theorem.

###  3.5 Shortcomings of the TCTE self-energy

While different flavors of the TCTE screened $G\tilde{W}$ self-energies yield
the exact xc energy, they will in general not yield the correct spectral
function calculated from the solution of the Dyson equation. One may expect
some improvement with respect to the $GW$ approximation for the quasiparticle
(QP) energies, since the use of $\bar{f}_{\rm xc}$, which is negative, reduces
the polarization contribution and therefore approximates one important effect
of the full vertex corrections, which is to reduce self-polarization 25.
However, one may expect that it will not be sufficient to bring significant
correction to the satellites, which are in general poorly described by the
$GW$ approximation. One reason lies in the fact that the $GW$ self-energy is
of first order in $W$, but used in the solution of the Dyson equation plus
infinite order. Another reason is the following: the poles of the exact
Green’s function are the total energy differences $\pm(E_{N\pm 1,s}-E_{N})$,
where $N$ is the particle number and $s$ labels a ground ($s=0$) or excited
state $s$. This can be written as $\pm(E_{N\pm 1,s}-E_{N\pm 1,0})\pm(E_{N\pm
1,0}-E_{N})$, i.e., the excitation energy of the $N\pm 1$-electron system plus
the chemical potential for electrons or holes. This means that satellites of
the QP, that lies at the respective chemical potential, must be found at a
distance equal to the excitation energies of the $N\pm 1$-electron system, and
not, as it would be the case in the $GW$ approximation for small systems with
a discrete spectrum, at a distance close to the excitation energies of the
$N$-electron system (plus differences in input and output QP energies when $G$
is not calculated self-consistently). This shortcoming cannot be overcome by a
$f_{\mathrm{xc}}(\omega)$ that depends on a single frequency and multiplies
$\chi(\omega)$ in frequency space: such a structure cannot shift the poles of
$\chi(\omega)$. This could only be achieved by a frequency integration, as it
is the case when the true three-times vertex correction is used. One should
therefore at best expect corrections of the intensities of the satellites when
moving from $GW$ to $G\tilde{W}$.

##  4 Illustrations

In order to illustrate our main findings and suggestions, we will use a simple
exactly solvable model, the symmetric Hubbard dimer25, 73, 74, 75, 76. Its
hamiltonian reads77, 78

$\hat{H}=\sum_{i,\sigma}\epsilon_{0}\hat{n}_{i\sigma}-\sum_{<i,j>,i\neq
j,\sigma}t\hat{c}^{\dagger}_{i\sigma}\hat{c}_{j\sigma}+U\sum_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}.$
(35)

where $i,j$ denote the sites $1,2$, the spin $\sigma=\uparrow,\downarrow$, the
external on-site potential is $\epsilon_{0}$, and $t$ is the hopping that is
linked to the kinetic energy. $U$ is the onsite Coulomb repulsion, and
$\hat{n}_{i\sigma}=\hat{c}^{\dagger}_{i\sigma}\hat{c}_{i\sigma}$ is the
particle number operator, where $\hat{c}$ and $\hat{c}^{\dagger}$ annihilate
and create a fermion, respectively. Using this simple model allows us to
explore the full range of correlation, which can be quantified by the ratio
$U/t$, and to have an unambiguous benchmark. We will use it at half filling,
i.e., with two electrons, which yields non-trivial electron removal and
addition features, and we will set $U=$4 eV throughout the illustrations. One
limitation of the model is the fact that the density is trivial and always
exact in all methods that conserve symmetry and particle number. Since in the
present work we suppose to know the exact KS ingredients, this is not a main
drawback. Moreover, asymmetry in the potential removes degeneracy and
therefore has a tendency to decrease correlation effects. The symmetric dimer
is therefore the most critical test case. Exploring density-driven errors
would be interesting, but beyond the scope of this work.

Figure 1: Symmetric Hubbard dimer at half filling and $U=4$ eV: error of the
exchange-correlation energy as a function of the hopping $t$.
$\mathrm{E}_{\mathrm{xc}}$ is obtained from
$\mathrm{E}_{\mathrm{xc}}=-\frac{i}{2}\int
G^{\mathrm{out}}\Sigma_{\mathrm{xc}}[G^{\mathrm{in}}]$ for different
$G^{\mathrm{in}},G^{\mathrm{out}}$, which are, respectively, the input GF used
to build the self-energy, and the GF that is usually the output of the Dyson
equation, but for which we have more options here. $\Sigma_{\mathrm{xc}}$ is a
$G\tilde{W}$ self-energy, built with $G^{\mathrm{in}}$ and using the
consistently chosen xc kernel. For a compact notation, we denote this by
$G^{\mathrm{out}}|G^{\mathrm{in}}|\bar{f}_{\mathrm{xc}}(\omega)$. The black
and sky blue solid lines with dot markers result from a self-energy built with
KS ingredients and integrated, respectively, consistently with
$G^{\mathrm{out}}=G^{\mathrm{in}}=G^{\mathrm{KS}}$ or, inconsistently, with
the $G^{\mathrm{out}}=G^{\mathrm{D}}$ resulting from the Dyson equation. The
red and green dashed lines with triangle markers result from a self-energy
built with QP ingredients and integrated, respectively, consistently with the
$G^{\mathrm{out}}=G^{\mathrm{in}}=G^{\mathrm{QP}}$ or, inconsistently, with
the $G^{\mathrm{out}}=G^{\mathrm{D}}$ resulting from the Dyson equation. Left
panel: results using the exact consistent $\bar{f}_{\mathrm{xc}}(\omega)$.
Middle panel: results using the adiabatic approximation $\omega=0$ for
$\bar{f}_{\mathrm{xc}}$. Right panel: results obtained by neglecting
$\bar{f}_{\mathrm{xc}}$ completely, which corresponds to a $GW_{0}$
approximation, where $W_{0}$ is calculated in the RPA and $G=G^{\mathrm{in}}$.

The exact analytical expressions for the time-ordered Green’s function and
self-energy are given in Supporting Information. For the approximate Green’s
functions we have solved the Dyson equation numerically. Our code uses
retarded quantities 79, 80 which is more reliable, since the numerical results
of the time-ordered calculations suffer from instabilities for some
approximations in the small range of $t$ ($t\rightarrow 0$). Computational
details are given in Supporting Information.

As worked out in Subsec. 3.1, different choices for $\chi_{0}$ are possible.
The simplest choice is to build $\chi_{0}$ with KS Green’s functions. In this
case, the corresponding xc kernel $f_{\mathrm{xc}}^{\mathrm{KS}}$ is the one
defined in TDDFT. In the symmetric Hubbard dimer the KS xc potential is a
number that we set by constraining the highest occupied level (HOMO) energy of
the KS system to yield the exact ionization potential. In this way we obtain
the KS Green’s function and $\chi_{0}^{\mathrm{KS}}$, and subsequently
$f_{\mathrm{xc}}^{\mathrm{KS}}$ by inversion of Eq. (24). This inversion is
not unique in the symmetric Hubbard dimer, because both the exact $\chi$ and
$\chi_{0}^{\mathrm{KS}}$ have only one non-zero element, which is the
antibonding/antiboding one (see Supporting Information). Therefore, as already
pointed out in 75, only the antibonding/antibonding matrix element of the
resulting $f_{\mathrm{xc}}$ is defined. The other elements are arbitrary, but
their choice has no impact on the results, since $f_{\mathrm{xc}}$ appears
only in the combination $\chi_{0}f_{\mathrm{xc}}\chi$.

Another natural choice would be to use the exact Green’s function $G$ to build
$\chi_{0}$, since it also yields the exact density, as required. However,
interestingly there is no solution to the inversion of Eq. (24) in this case.
The reason is that also the bonding/bonding element of this $\chi_{0}=-iGG$ is
non-vanishing. Further analysis shows that this stems from the satellite
contributions to $G$, which are not canceled by proper vertex corrections.
This is a nice illustration for one of the problems of this ill-behaved
polarizability which also, for example, does not fulfill the $f$-sum rule 81.
We will instead use $\chi_{0}^{\mathrm{QP}}\equiv-
iG^{\mathrm{QP}}G^{\mathrm{QP}}$. It is built with the QP approximation
$G^{\mathrm{QP}}$ to the exact $G$, where satellites are neglected and the
remaining intensities normalized to 1. This can be seen as a realization of a
generalized KS Green’s function, stemming from a potential that is non-local
in space but instantaneous in time. Such a potential can lead to accurate QP
energies 82, but not to satellites. In a real material, the widely used hybrid
functionals 83 fall into this class. Also many scalar long-range kernels are
designed to be used on top of a $\chi_{0}^{\mathrm{QP}}$. It should again be
stressed that both kernels, whether the one of the KS or the one of the QP
scheme, can be called “exact”, as long as they are used consistently in
conjunction with $\chi_{0}$ built with the corresponding Green’s functions.

###  4.1  Results using exact xc kernels

In the following we will focus on the results obtained with the two kernels
$f_{\mathrm{xc}}^{\mathrm{KS}}(\omega)$ and
$f_{\mathrm{xc}}^{\mathrm{QP}}(\omega)$, without approximating them further.
This will allow us to illustrate the effect of using an $f_{\mathrm{xc}}$ to
simulate the full three-argument vertex of many-body perturbation theory,
without further approximations.

####  4.1.1 Exchange-correlation energy

First, let us examine the xc contribution to the total energy, given by Eq.
(29)
$\mathrm{E}_{\mathrm{xc}}=-\frac{i}{2}\int\bar{G}\bar{\Sigma}_{\mathrm{xc}}$.
As pointed out above, here $\bar{G}$ should not be the exact Green’s function
nor the one resulting from the Dyson equation with
$\bar{\Sigma}_{\mathrm{xc}}$, which we will call $G^{\mathrm{D}}$ in the
following, but $\bar{G}$, which is the one used to build the $G\tilde{W}$
self-energy $\bar{\Sigma}_{\mathrm{xc}}$. This point is important since in
practical applications, using $G^{\mathrm{D}}$ would often seem to be a
natural choice, being the best available Green’s function, i.e. the one
closest to the exact $G$. We will therefore compare these choices in the
following, by evaluating $\mathrm{E}_{\mathrm{xc}}=\int
G^{\text{out}}\bar{\Sigma}_{\mathrm{xc}}[G^{\text{in}}]$. Here,
$G^{\mathrm{in}}$ is the input GF used to build the $G\bar{\tilde{W}}$ self-
energy $\bar{\Sigma}_{\mathrm{xc}}$, and the $G^{\mathrm{out}}$ is either the
output of the corresponding Dyson equation $G^{\mathrm{D}}$, or equal to
$G^{\mathrm{in}}$. In all cases, $\bar{\Sigma}_{\mathrm{xc}}$ is built with
the xc kernel $\bar{f}_{\mathrm{xc}}$ that is consistent with
$G^{\mathrm{in}}$.

For a compact notation, we use
$G^{\mathrm{out}}|G^{\mathrm{in}}|\bar{f}_{\mathrm{xc}}(\omega)$. For example,
$G^{\mathrm{D}}|G^{\mathrm{KS}}|f_{\mathrm{xc}}^{\mathrm{KS}}(\omega)$ stands
for $\mathrm{E}_{\mathrm{xc}}=-\frac{i}{2}\int
G^{\mathrm{D}}\Sigma_{\mathrm{xc}}[G^{\mathrm{KS}}]$, where the
$G\bar{\tilde{W}}$ self-energy is built using the KS Green’s function and KS
xc kernel. The Dyson equation is then solved using this self-energy, and the
resulting Green’s function $G^{\mathrm{D}}$ is used in the integral. Note that
while $G^{\mathrm{D}}$ is not the same in the KS and QP frameworks, we do not
highlight this difference in the notation, since it is clear from the context.
Comparison of the various flavors allows us to illustrate the importance of
the consistency requirement advocated in Sec. 3.2. For subsequent
investigation, we also indicate by $|\bar{f}_{\mathrm{xc}}(\omega)$ whether
the exact consistent $\bar{f}_{\mathrm{xc}}(\omega)$ is used or further
approximations are made, e.g., $|\bar{f}_{\mathrm{xc}}(\omega=0)$. Fig. 1
shows the difference to the exact xc energy $\mathrm{E}_{\mathrm{xc}}$. The
results in the left panel were obtained using the exact consistent
$\bar{f}_{\mathrm{xc}}(\omega)$. As predicted by Eq. (29), the two consistent
calculations
$G^{\mathrm{KS}}|G^{\mathrm{KS}}|f_{\mathrm{xc}}^{\mathrm{KS}}(\omega)$ and
$G^{\mathrm{QP}}|G^{\mathrm{QP}}|f_{\mathrm{xc}}^{\mathrm{QP}}(\omega)$ both
yield the exact result. Instead, when solution of the Dyson equation is used
for $G^{\text{out}}$ we obtain
$G^{\mathrm{D}}|G^{\mathrm{KS}}|f_{\mathrm{xc}}^{\mathrm{KS}}(\omega)$ and
$G^{\mathrm{D}}|G^{\mathrm{QP}}|f_{\mathrm{xc}}^{\mathrm{QP}}(\omega)$, which
are both inconsistent and therefore not exact. The error of the former is
larger than that of the latter. This can be understood, since the difference
between $G^{\mathrm{D}}$ and the input Green’s function is larger in the case
of KS than in the case of the QP input. In all cases, errors are vanishing for
large $t$, whereas they increase in the inconsistent calculations with
decreasing $t$. Even closer to the atomic limit, all errors tend to zero.
Nevertheless, the importance of consistency is nicely illustrated by this
result.

Figure 2: Kinetic energy errors as a function of the hopping parameter $t$.
Left panel: $\mathrm{E}_{\mathrm{k}}$ is calculated with the Green’s function
resulting $G^{\mathrm{D}}$ from the Dyson equation with a $G\tilde{W}$ self-
energy (red with triangles and sky blue with dots) or with a $GW^{\rm TCTC}$
self-energy, where $W^{\rm TCTC}$ is the exact TCTC screened Coulomb
interaction (violet with dots). The self-energy is built with KS ingredients
(sky blue and violet) or QP ingredients (red). Right panel: $f_{\mathrm{xc}}$
is approximated adiabatically (orange with dots for KS ingredients, black with
triangles for QP ingredients) or completely neglected (green with dots for KS
ingredients, dark blue with triangles for QP ingredients).

####  4.1.2 Kinetic energy

While an approximate self-energy used in the consistent scheme yields exact
results for $\mathrm{E}_{\mathrm{xc}}$, no such scheme exists for the kinetic
energy. Instead, by definition the result of the Dyson equation
$G^{\mathrm{D}}$ is used to determine the density matrix and hence, the
kinetic energy. We will therefore examine the error introduced by various
flavors of the self-energy, starting with those that can yield the exact
$\mathrm{E}_{\mathrm{xc}}$. The left panel of Fig. 2 shows the results for
$G^{\mathrm{D}}|G^{\mathrm{KS}}|f_{\mathrm{xc}}^{\mathrm{KS}}(\omega)$ and
$G^{\mathrm{D}}|G^{\mathrm{QP}}|f_{\mathrm{xc}}^{\mathrm{QP}}(\omega)$. Both
show errors that only vanish at large $t$ and for $t\to 0$. The KS flavor
converges more quickly to the exact result with increasing $t$ than the QP
version. This favors the use of the $G\tilde{W}$ self-energy built with KS,
rather than QP, ingredients. Still, the error is significant. However, as
noted in Subsec. 3.3, with an exact interaction energy one can, in principle,
also obtain the exact kinetic energy by using the virial theorem. This allows
one to overcome the problem of not knowing the exact density matrix.

Figure 3: Error of the QP energies as a function of the hopping $t$. The LUMO
and HOMO energy errors are shown in the upper and lower panels, respectively.
Left column: The result of the Dyson equation with a $G\tilde{W}$ self-energy
with KS ingredients (blue with dots) or QP ingredients (red with triangles),
or with a $GW$ self-energy using the exact TCTC screened Coulomb interaction
(violet with dots) is shown. Right column: The result of the Dyson equation
with a $G\tilde{W}$ self-energy with KS ingredients using a static
$f_{\mathrm{xc}}(\omega=0)$ (orange with dots) or neglecting $f_{\mathrm{xc}}$
(green with dots), or with QP ingredients using a static
$f_{\mathrm{xc}}(\omega=0)$ (black with triangles) or neglecting
$f_{\mathrm{xc}}$ (blue with triangles) is shown.

####  4.1.3 Spectra

The situation is different for spectral properties: here, the shortcomings of
an approximate Green’s function cannot be overcome easily. As for the kinetic
energy, the result of the Dyson equation is used to calculate the spectra. We
will explore which of the flavors of the self-energy that gives an in
principle exact total energy will yield the best spectral properties, and what
are the remaining problems.

Let us first look at the QP peaks of the spectral function Fig. 3 shows the
error of the position of the HOMO (lower panel) and of the lowest unoccupied
state (LUMO) (upper panel) as a function of $t$. The two panels in the left
column contain the HOMO and LUMO energy errors obtained with the exact KS or
QP ingredients to build the $G\tilde{W}$ self-energy. While the KS and QP
flavor perform very similarly for the HOMO, with small errors at larger $t$
and a significant deviation from the exact result for small $t$ that goes up
to U/2 for $t\to 0$, the LUMO is relatively well described for all $t$, and
the error vanishes for $t\to 0$. Results for the LUMO are particularly
satisfying when KS ingredients are used, in which case the error does not
exceed 0.1 eV for any $t$. For larger $t$, above 2.5 eV, the errors become
small for both HOMO and LUMO, especially in the KS flavor, where they remain
well below 0.1 eV and quickly move into the meV range (see insets).

Beyond the QP features, Fig. 4 shows the entire spectral functions for $t=0.5$
eV. We will concentrate on the satellites. They are due to the peaks in the
imaginary part of the self-energy, which are in turn determined by the peaks
of $\tilde{W}$: the poles of ${\rm Im}\,\bar{\Sigma}_{\mathrm{xc}}$ are
situated at energies $\bar{\varepsilon}_{i}\pm\omega_{j}$, where
$\bar{\varepsilon}_{i}$ is a removal/addition pole of the
$G^{\text{in}}=\bar{G}$ used to build the self-energy, and $\omega_{j}$ is a
pole of $\chi$. Not all poles are visible in all matrix elements: in the
symmetric Hubbard dimer, the bonding (antibonding) matrix element of the self-
energy is dominated by the addition (removal) part of $G^{\mathrm{in}}$. The
bonding (antibonding) matrix element satellites are therefore found at
energies higher (lower) than the LUMO (HOMO). In many real materials, all
parts of the Green’s function contribute to all matrix elements of the self-
energy, and satellites are found on both sides of a QP. In this sense, the
Hubbard dimer is an extreme case, where a given matrix element selects just
one particular excitation, that may moreover not be the intuitively expected
one. This does not influence our conclusions, but it is interesting to note.

The most obvious feature in Fig. 4 is the fact that satellites are not well
described in general when the exact $\bar{f}_{\mathrm{xc}}(\omega)$ is used.
Their position at $\bar{\varepsilon}_{i}\pm\omega_{j}$ combines two errors:
the fact that the excitation energy $\omega_{j}$ of the $N$ electron system is
used (see Subsec. 3.5), and the fact that $\bar{\varepsilon}_{i}$ can be
different from the true QP energy. Since in our case the antibonding matrix
element of the self-energy is dominated by the HOMO $\bar{\varepsilon}$, the
exact QP energy is used in all cases studied here and the error is entirely
due to the difference between the (too high) excitation energy of the $N$
electron system with respect to the $N-1$ electron one. For the bonding matrix
element, instead, the LUMO $\bar{\varepsilon}$ is used, which is exact when QP
ingredients are used, but which is too low in the KS case. This adds to the
error of the (too high) excitation energy of the $N$ electron system with
respect to the $N+1$ electron one. Since the two errors are of opposite sign,
the KS ingredients yield the best result for the bonding matrix element. The
difference between the $N$ and $N\pm 1$ excitation energies should be of
particular importance in finite systems, but an analogous error might also
impact results in infinite systems with localized electrons.84 Note, that the
problem discussed here is different from another issue in extended systems,
where the satellite position can be spoiled by the appearance of a plasmaron,
a spurious solution of the QP condition that is found at some distance from
the peak in the imaginary part of the self-energy 85, 86, 87. In a discrete
system such as the Hubbard dimer, instead, the satellites are always found
close to the position of peaks of the imaginary part of the self-energy, and
the point here is that this position is calculated with the wrong number of
electrons.

The biggest effect of $f_{\mathrm{xc}}$ is to decrease screening, which
remedies the self-screening problem for the QPs25, 88, 89, but which also
decreases the satellite intensity because, as can be seen in Fig. 5,
$f_{\mathrm{xc}}$ is always negative. Indeed, the $G\tilde{W}$ satellites in
Fig. 4 are of much too weak intensity. KS ingredients do a bit better than QP
ones in this respect, since in this case a weaker $\bar{f}_{\mathrm{xc}}$ is
used (see Fig. 5), which leads to a smaller decrease of the satellite
intensities, but the result is still unsatisfactory. This dilemma cannot be
solved with such a simple vertex correction that is multipicative in
frequency. In other words and as expected, $G\tilde{W}$, even with exact KS or
QP ingredients, cannot yield reliable satellites.

Figure 4: Bonding-bonding (left panel) and antibonding-antibonding (right
panel) matrix elements of the spectral function for $U=4$ eV and $t=0.5$ eV.
The continuous black curves are the exact result. The result of the Dyson
equation using a $G\tilde{W}$ self-energy with KS ingredients and the exact
$f_{\mathrm{xc}}^{\mathrm{KS}}$ is in sky blue. While the result for
$f_{\mathrm{xc}}^{\mathrm{KS}}(\omega)=0$ is in green. The result of the Dyson
equation using QP ingredients and the exact $f_{\mathrm{xc}}^{\mathrm{QP}}$ is
in dashed red, while the result for $f_{\mathrm{xc}}^{\mathrm{QP}}(\omega)=0$
is in dashed blue. The red arrow indicates the position of the very weak
satellite obtained when QP ingredients are used. Moreover, the result of the
Dyson equation using a $GW$ self-energy with the exact TCTC screening is shown
in violet. The exact HOMO is situated at 0. Note that the corresponding
satellites are found at higher energies. The exact QP of the LUMO is situated
at 3.5 eV, with satellites in the electron removal energy range.

###  4.2 Impact of approximating $f_{\mathrm{xc}}$

Understanding the impact of replacing the full vertex corrected self-energy
with a $G\tilde{W}$ one is of fundamental interest. For practical
applications, one also has to face the problem that the exact
$f_{\mathrm{xc}}(\omega)$ is in general not known. Therefore, we also briefly
examine the impact of two widely used approximations: either a complete
neglect of $f_{\mathrm{xc}}$, which brings us back to the $GW$ approximation
with an RPA $W=W_{0}$, or at least an adiabatic approximation where only
$f_{\mathrm{xc}}(\omega=0)$ is used, since the frequency dependence of
$f_{\mathrm{xc}}(\omega)$ is notoriously difficult to approximate. As we will
see, these approximations do not have the same impact according to the flavor
(KS or QP) that is chosen, and according to the combination of ingredients.

####  4.2.1 Exchange-correlation energy: impact of approximations

Let us first look at the quantity that is obtained exactly when $G\tilde{W}$
is used consistently, namely, the xc contribution $\mathrm{E}_{\mathrm{xc}}$
to the total energy. The middle panel of Fig. 1 compares results using the
adiabatic approximation $\bar{f}_{\mathrm{xc}}(\omega=0)$ and combining the
ingredients in a consistent or inconsistent way, respectively. Similarly,
results in the right panel were obtained by completely neglecting
$f_{\mathrm{xc}}$. In all cases, the consistent results now show an error, but
it is smaller than that of the corresponding inconsistent results, which
demonstrates that a consistent choice of ingredients remains essential to
obtain good total energies. The impact of neglecting $f_{\mathrm{xc}}$ is
smaller when KS ingredients are used. The best results are obtained using the
consistent KS flavor. When the adiabatic approximation is used, the fact that
the performance of KS remains good can be explained by the fact that the
quadratic frequency dependence of the kernel, which is a universal feature of
$f_{\mathrm{xc}}$ 90, is milder in the KS than in the QP case, as shown in
Fig. 5. Although approximate, the benefit of using $f_{\mathrm{xc}}$ remains
very important, as can be seen by comparing the middle panel and the right
panel, where results on the $GW$ level with an RPA $W_{0}$ are given. The
$GW_{0}$ results tend to the exact result very slowly with increasing $t$, and
a part from the consistent KS flavor, they deviate significantly from the
exact result in the atomic limit. The
$G^{\mathrm{D}}|G^{\mathrm{KS}}|f_{\text{xc}}(\omega)=0$ flavor tends to
$U/2$,while both consistent and non-consistent QP cases tend to $U/4$.

Figure 5: Antibonding matrix element of $f_{\text{xc}}(\omega)$ as a function
of frequency for $U=4$ eV and $t=3$ eV. The light blue curve shows the xc
kernel corresponding to the KS system, whereas the dashed orange result is the
xc kernel that is consistent with QP ingredients.

####  4.2.2 Kinetic energy: impact of approximations

The kinetic energy is never exact in $G\tilde{W}$, as explained above and as
illustrated in Fig. 2. The right panel of Fig. 2 also shows the impact on the
kinetic energy of approximations to $f_{\mathrm{xc}}$. The adiabatic
approximation $f_{\mathrm{xc}}(\omega=0)$ has a very moderate effect, with a
tendency that is rather towards improving the results. The reason for this is
the fact that the kinetic energy suffers from the underestimate of the
satellite intensity discussed above, which is improved when $f_{\mathrm{xc}}$
is weaker. Neglecting the quadratic frequency-dependence of the kernel shown
in Fig. 5 is therefore rather beneficial for the kinetic energy. A complete
neglect of $f_{\mathrm{xc}}$, instead, spoils results in the moderate to large
$t$-range, while further slightly improving the small-$t$ regime, where the
satellites are important. Overall, KS flavors perform slightly better than QP
ones. Finally, we also show in the left panel the result of a
$GW^{\text{TCTC}}$ calculation, where the self-energy is of $GW$ form using KS
ingredients and the exact $\chi$ and therefore the exact $W$ is used, but
where the vertex $\Gamma=1$ in the self-energy, i.e., the functional
derivative in Eq. (2) is set to 1. This means that the exact test-charge test-
charge (TCTC) screening is used instead of the TCTE one that is used in the
$G\tilde{W}$ approximation. Indeed, it would be tempting to think that a very
good $W$ used in $GW$ could improve results. However, with respect to a
standard $GW_{0}$ calculation using an RPA $W=W_{0}$, where
$f_{\mathrm{xc}}=0$ also in $\chi$, the results are worse, especially in the
moderate to large-$t$ regime. It has been pointed out that vertex corrections
in the polarizability and in the self-energy tend to cancel partially91, 47:
the present result is a good illustration.

####  4.2.3 Use of the adiabatic connection versus virial theorem

Finally, we can examine the quality of the result that can be obtained by
using the virial theorem, instead of approximating the kinetic energy
directly, and compare to the results obtained using the AC fluctuation-
dissipation theorem discussed in Sec. 3.4. Both approaches are in principle
exact, but might react differently to approximations.

Fig. 6 gives the errors of the full correlation energy including interaction
and kinetic contributions, obtained using an adiabatic kernel,
$f_{\mathrm{xc}}(\omega=0)$ and KS ingredients. In order to use the virial
theorem, one has to determine the term $S_{\rm VT}$ in Eq. (31). We bypass the
difficulty to adapt this equation to the Hubbard dimer by using the fact that
here we work with the exact density in all cases, which allows us to use the
exact $S_{\rm VT}$, which we obtain from the exact solution as $S_{\rm
VT}\equiv 2\mathrm{E}_{\mathrm{k}}+\mathrm{E}^{\mathrm{inter}}$ for all values
of the hopping $t$. The resulting $S_{\rm VT}$ is then used in place of the
right hand side of Eq. (31) in order to obtain
$\mathrm{E}_{\mathrm{k}}=(S_{\rm VT}-\mathrm{E}^{\mathrm{inter}})/2$ for a
given approximation to $\mathrm{E}^{\mathrm{inter}}$. This procedure gives the
light blue curve (VT) in Fig. 6. As predicted in Sec. 3.4, the error is
similar to the one of the AC approach using the same approximation
$f_{\text{xc}}(\omega=0)$. This is indeed due to the fact that the integrand
of the full correlation energy depends approximately linearly on the coupling
constant $\lambda$, as one can see in the inset of Fig. 6. The difference of
the correlation energy $\text{E}^{\text{full}}_{\text{c}}$ is very small
around $t=3$ eV where the behaviour is almost exactly linear, while the
deviation is larger at the smaller $t=0.5$ eV, where a quadratic
$\lambda$-dependence is clearly visible. In this small-$t$ regime, where the
function is convex, the approach using the virial theorem performs better,
while also avoiding the need for the $\lambda$-integration.

Figure 6: Error of the full correlation energy (kinetic and interaction
contributions) as a function of the hopping $t$, when the adiabatic
approximation $f_{\mathrm{xc}}(\omega=0)$ is made and KS ingredients are used:
comparison of the adiabatic connection result (AC, in magenta) with the result
obtained using the virial theorem (VT, in cyan). Insets: $\lambda$-resolved
exact full correlation energy
$\text{E}_{\text{c}}^{\text{full}}=\int_{0}^{1}{d\lambda\text{E}_{\text{c}}^{\lambda}}$,
evaluated within the KS scheme, as a function of $\lambda$, for $t=0.5$ eV
(upper inset) and $t=3$ eV (lower inset).

####  4.2.4 Spectra

The fact that $G\tilde{W}$ does not yield the correct spectral properties
cannot be overcome, but it is still interesting to examine the effect of
approximations made in practice. This is done in the right panels of Fig. 3
for the QP energies, and in Fig. 4 for the satellites. For the LUMO position,
both the adiabatic approximation and neglecting $f_{\mathrm{xc}}$ completely
lead to significant worsening of the result in the small-$t$ regime, the worst
results being obtained with KS and QP $GW_{0}$, i.e, $f_{\mathrm{xc}}=0$,
which also slightly deteriorates results at larger $t$. It is interesting to
note that keeping the exact $f_{\mathrm{xc}}(\omega)$ in $W$ alone, i.e.,
using the exact $W^{\rm TCTC}$ instead of the RPA $W=W_{0}$, does not fix any
of these problems, as one can see in the left panel of Fig. 3 for the LUMO.
The same is true also in the case of the HOMO. These findings are in line with
observations on real systems 92. Concerning the other approximations for the
HOMO, shown in the lower right panel of Fig. 3, the observation concerning the
$GW$ approximation is similar to the LUMO for moderate to large $t$, whereas
the adiabatic approximation is rather beneficial, especially for smaller $t$.
Also a complete neglect of $f_{\mathrm{xc}}$, i.e., the $GW$ solution with RPA
$W=W_{0}$, decreases the error for small $t$, and when KS ingredients are
used, the $GW_{0}$ results even reaches the correct $t\to 0$ limit. However,
in this case the improvement is limited to a very small range of $t$ close to
the atomic limit. The observed trends highlight the fact that the effect of
including $f_{\mathrm{xc}}$ is beneficial for the LUMO at all $t$ and for the
HOMO at large $t$, but too strong for the HOMO at small $t$. Since, as
discussed in 4.1.3, matrix elements of the self-energy are quite particular in
the Hubbard dimer, this observation should not be generalised and further
analysis will be needed to eventually turn these findings into a systematic
correction, which is beyond the scope of the present work.

Finally, Fig. 4 illustrates that including $f_{\mathrm{xc}}$ in $W$ alone,
i.e., performing a $GW^{\rm TCTC}$ calculation, rather worsens the satellites
as compared to a $GW$ result obtained with RPA $W_{0}$, which illustrates
again the error canceling. Therefore, in Fig. 4 the best satellite results are
obtained using the $GW_{0}$ with KS ingredients and a complete neglect of
$\bar{f}_{\mathrm{xc}}$. Note, however, that this is not a general finding for
all values of $U/t$.

##  5 Conclusions

In conclusion, the _exact_ exchange-correlation contribution to the total
interaction energy can be calculated using an _approximate_ self-energy of the
form $G\tilde{W}$. Here, $\tilde{W}$ is a test-charge test-electron screened
Coulomb interaction, which replaces the RPA or the TCTC screened interaction
that are commonly used in the $GW$ approximation. Different choices for
$\tilde{W}$ are possible, one of them being the traditionally used KS scheme,
which adds an xc kernel $f_{\rm xc}(\omega)$ from linear response TDDFT to the
bare Coulomb interaction in the dielectric function. For all choices the
condition is that the GF and xc kernel used to build the self-energy are
consistent and yield the correct density. On top of the KS choice, we have
examined the case where the GF is built with the exact QP energies. For all
possible choices, it holds that the exact xc energy is obtained by integrating
the approximate self-energy with the very same GF that was used to build it.
Instead, when the approximate self-energy is used in a Dyson equation and
integrated with the resulting GF, the results carry an error. The importance
of consistency between the GF used to build the self-energy and the GF used
for the integration may explain the success of self-consistent $GW$ total
energy calculations, which indeed fulfill the requirement that the self-energy
is integrated with the GF that is used to build it. Here, we show that one can
obtain good quality results by being consistent without carrying out self-
consistent calculations.

The exact correlation contribution to the kinetic energy cannot be accessed in
the same way. Instead, we propose to use the virial theorem. We have studied
the impact of widely used approximations to this approach, and compared with
the use of the adiabatic connection fluctuation dissipation theorem. Our
general derivation predicts that the final errors are similar, without the
need of a coupling constant integration in the present approach.

Using the approximate self-energies in the Dyson equation leads to approximate
GFs and therefore, to approximate spectral functions. Still, $G\tilde{W}$
yields overall better QPs than $GW$, and since the computational difficulty is
similar, it should be preferred. The satellite problem, instead, cannot be
fixed in this way.

All statements have been illustrated for the symmetric half-filled Hubbard
dimer, confirming our conjectures and highlighting the fact that results
obtained using KS ingredients are overall superior and less impacted by
additional approximations with respect to results obtained using QP
ingredients. While the Hubbard dimer is a simple model, our findings relie on
derivations that are valid for the general case, and they should open the way
for interesting applications to more realistic systems.

The authors acknowledge the fruitful discussions with Steffen Backes, Fabien
Bruneval, Kieron Burke and Steven Crisostomo.
This project has received funding from the European Union’s Horizon 2020
research and innovation programme under grant agreement No 800945 — NUMERICS —
H2020-MSCA-COFUND-2017. The exact and $GW$ solutions for the symmetric Hubbard
dimer model at half-filling (two electrons) are given in literature 25. In
this section, we provide the solutions for the model within the $GW$ and
$G\tilde{W}$ approximations using both the KS and QP flavors.

##  1 $GW$ solutions

The two ingredients needed to calculate $\Sigma_{\mathrm{xc}}$ at the $GW$
level are the GF and the screened Coulomb interaction.

  1. 1.

The exact Kohn-Sham GF and the exact Quasi-Particle (QP) GF in the dimer sites
basis read respectively

$G^{\mathrm{KS}}_{ij\sigma}(\omega)=\frac{1}{2}\bigg{(}\frac{1}{\omega-(\epsilon_{0}+t-(c-U)/2)-i\eta}+\frac{(-1)^{i-j}}{\omega-(\epsilon_{0}+3t-(c-U)/2)+i\eta}\bigg{)},$
(36)
$G^{\mathrm{QP}}_{ij\sigma}(\omega)=\frac{1}{2}\bigg{(}\frac{1}{\omega-(\epsilon_{0}+t-(c-U)/2)-i\eta}+\frac{(-1)^{i-j}}{\omega-(\epsilon_{0}-t+(c+U)/2)+i\eta}\bigg{)},$
(37)

where $c=\sqrt{16t^{2}+U^{2}}$. The KS GF ($G^{\mathrm{KS}}$) is obtained by
introducing an energy shift to the poles of the non-interacting GF such that
HOMO energy becomes exact93. The QP GF $(G^{\mathrm{QP}}$) equals the exact GF
without the satellite contributions and with the quasiparticle intensities set
to $1$.

  2. 2.

We use both the exact screened Coulomb interaction $W$ and approximations
denoted $W_{0}^{\mathrm{KS}}$ or $W_{0}^{\mathrm{QP}}$, depending on the
choice of the GF used to compute the polarizability. The $W_{0}^{\mathrm{KS}}$
and $W_{0}^{\mathrm{QP}}$ are calculated within the Random Phase Approximation
(RPA), using the following irreducible polarizabilities

$P_{\mathrm{RPA}}^{\mathrm{KS}}(1,2)=-iG^{\mathrm{KS}}(1,2^{+})G^{\mathrm{KS}}(2,1^{+}),$
(38)

and

$P_{\mathrm{RPA}}^{\mathrm{QP}}(1,2)=-iG^{\mathrm{QP}}(1,2^{+})G^{\mathrm{QP}}(2,1^{+}),$
(39)

respectively for the KS and QP cases. They have the following analytical
expressions

$P_{\mathrm{RPA},ij\sigma}^{\mathrm{KS}}(\omega)=\frac{(-1)^{i-j}}{4}\bigg{(}\frac{1}{\omega-2t+i\eta}-\frac{1}{\omega+2t-i\eta}\bigg{)},$
(40)
$P_{\mathrm{RPA},ij\sigma}^{\mathrm{QP}}(\omega)=\frac{(-1)^{i-j}}{4}\bigg{(}\frac{1}{\omega+(2t-c)+i\eta}-\frac{1}{\omega-(2t-c)-i\eta}\bigg{)}.$
(41)

For the exact $W$, we use the exact reducible polarizability $\chi$, which is
related to the 2-particle Green’s Function ($G_{2}$)

$\chi(1,2)=-iG(1,1^{+})G(2,2^{+})+iG_{2}(1,2,1^{+},2^{+}),$ (42)

where $G_{2}$ is the 2-GF. In the Hubbard dimer site ($ij$) basis we have

$\chi_{ij\sigma_{1}\sigma_{2}}(\omega)=\sum_{s\neq
0}\bigg{[}\braket{N_{0}}{\hat{c}^{\dagger}_{i\sigma_{1}}\hat{c}_{i\sigma_{1}}\ket{N_{s}}\bra{N_{s}}\hat{c}^{\dagger}_{j\sigma_{2}}\hat{c}_{j\sigma_{2}}}{N_{0}}\frac{1}{\omega+(\text{E}^{N}_{0}-\text{E}^{N}_{s})+i\eta}\\\
-\braket{N_{0}}{\hat{c}^{\dagger}_{j\sigma_{2}}\hat{c}_{j\sigma_{2}}\ket{N_{s}}\bra{N_{s}}\hat{c}^{\dagger}_{i\sigma_{1}}\hat{c}_{i\sigma_{1}}}{N_{0}}\frac{1}{\omega-(\text{E}^{N}_{0}-\text{E}^{N}_{s})-i\eta}\bigg{]}\,,$
(43)

which leads to the following solutions

$\chi_{ij\uparrow\uparrow}(\omega)=\frac{(-1)^{i-j}}{2a^{2}}\bigg{(}\frac{1}{\omega-(c+U)/2+i\eta}-\frac{1}{\omega+(c+U)/2-i\eta}\bigg{)}\\\
+(-1)^{i-j}\frac{16t^{2}}{2a^{2}(c-U)^{2}}\bigg{(}\frac{1}{\omega-(c-U)/2+i\eta}-\frac{1}{\omega+(c-U)/2-i\eta}\bigg{)},$
(44)
$\chi_{ij\uparrow\downarrow}(\omega)=\frac{(-1)^{i-j}}{2a^{2}}\bigg{(}\frac{1}{\omega-(c+U)/2+i\eta}-\frac{1}{\omega+(c+U)/2-i\eta}\bigg{)}\\\
-(-1)^{i-j}\frac{16t^{2}}{2a^{2}(c-U)^{2}}\bigg{(}\frac{1}{\omega-(c-U)/2+i\eta}-\frac{1}{\omega+(c-U)/2-i\eta}\bigg{)},$
(45)

where $a^{2}=2\bigg{(}\frac{16t^{2}}{(c-U)^{2}}+1\bigg{)}$ and
$\chi_{ij\uparrow\uparrow}=\chi_{ij\downarrow\downarrow}$,
$\chi_{ij\uparrow\downarrow}=\chi_{ij\downarrow\uparrow}$.

The spin-independent $\chi$ matrix in the site basis is a sum over spins, i.e.
$\chi_{ij}(\omega)=\chi_{ij\uparrow\uparrow}+\chi_{ij\uparrow\downarrow}+\chi_{ij\downarrow\uparrow}+\chi_{ij\downarrow\downarrow}$.
In the bonding and anti-bonding (b/ab) basis, $\chi$ reads

$\begin{pmatrix}0&0\\\ 0&\chi_{\text{ab-ab}}(\omega)\end{pmatrix},$ (46)

where $\chi_{\text{ab-ab}}(\omega)=2(\chi_{11}(\omega)+\chi_{22}(\omega))$. We
can now write the screened Coulomb interaction $W$ in the different flavors.
The exact screened Coulomb interaction is

$\displaystyle W(1,2)$
$\displaystyle=v_{\mathrm{c}}(1,2)+\int{d(34)v_{\mathrm{c}}(1,3)P(3,4)W(4,2)},$
$\displaystyle=v_{\mathrm{c}}(1,2)+\int{d(34)v_{\mathrm{c}}(1,3)\chi(3,4)v_{\mathrm{c}}(4,2)}\,,$

or, in the site basis

$\displaystyle W_{ij}(\omega)$
$\displaystyle=U\delta_{ij}+U\sum_{k\sigma}P_{ik\sigma}(\omega)W_{kj}(\omega),$
$\displaystyle=U\delta_{ij}+U^{2}\sum_{\sigma\sigma^{\prime}=\uparrow,\downarrow}\chi_{ij\sigma\sigma^{\prime}}(\omega)\,,$

which leads to

$W_{ij}(\omega)=U\delta_{ij}+(-1)^{i-j}\frac{2U^{2}}{a^{2}}\bigg{(}\frac{1}{\omega-(c+U)/2+i\eta}-\frac{1}{\omega+(c+U)/2-i\eta}\bigg{)}\,.$

By using Eq.s (40)(41), we find

$W_{0,ij}^{\mathrm{KS}}(\omega)=U\delta_{ij}+(-1)^{i-j}\frac{U^{2}t}{h}\bigg{(}\frac{1}{\omega-h+i\eta}-\frac{1}{\omega+h-i\eta}\bigg{)}\,,$
(47)

where $h=\sqrt{4t^{2}+4Ut}$, and

$W_{0,ij}^{\mathrm{QP}}(\omega)=U\delta_{ij}+(-1)^{i-j}\frac{U^{2}(c/2-t)}{h^{\prime}}\bigg{(}\frac{1}{\omega-h^{\prime}+i\eta}-\frac{1}{\omega+h^{\prime}-i\eta}\bigg{)}\,,$
(48)

where $h^{\prime}=\sqrt{(2t-c)^{2}+4U(c/2-t)}$. Now, by using
$G^{\mathrm{KS}}$, $G^{\mathrm{QP}}$, $W_{0}^{\mathrm{KS}}$,
$W_{0}^{\mathrm{QP}}$ and $W$ we calculate the different flavors of
$\Sigma_{\mathrm{xc}}$ by integrating in frequency space. Finally, we convert
$\Sigma_{\mathrm{xc}}$ in Eq. (27) to frequency space. Note that the use of
multiple infinitesimals in Eq. (27) is not important for the self-energy
itself, but for the calculation of $\mathrm{E}_{\mathrm{xc}}$, as the
respective weight of the infinitesimals in the different contributions
indicates the contour that is to be used in the integral. So

$\Sigma_{\mathrm{xc},ij\sigma}(\omega)=\frac{i}{2\pi}\int{d\omega^{\prime}\,G_{0,ij\sigma}(\omega^{\prime}+\omega)W_{0,ji}(\omega^{\prime})e^{3i\omega^{\prime}\eta}}\,,$
(49)

where $G_{0}$ can be $G^{\mathrm{KS}}$ or $G^{\mathrm{QP}}$, and $W_{0}$ can
be $W_{0}^{\mathrm{KS}}$, $W_{0}^{\mathrm{QP}}$ or $W$.

The solutions of the different $GW$ flavors are

$\Sigma^{G^{\mathrm{KS}}W_{0}^{\mathrm{KS}}}_{\text{xc},ij\sigma}(\omega)=-\frac{U}{2}\delta_{ij}+\frac{U^{2}t}{2h}\bigg{(}\frac{1}{\omega-(\epsilon_{0}+3t-(c-U)/2+h)+i\eta}\\\
+\frac{(-1)^{i-j}e^{-3i\omega\eta}}{\omega-(\epsilon_{0}+t-(c-U)/2-h)-i\eta}\bigg{)}\,,$
(50)

$\Sigma^{G^{\mathrm{KS}}W}_{\text{xc},ij\sigma}(\omega)=-\frac{U}{2}\delta_{ij}+\frac{U^{2}}{a^{2}}\bigg{(}\frac{1}{\omega-(\epsilon_{0}+3t+U)+i\eta}+\frac{(-1)^{i-j}e^{-3i\omega\eta}}{\omega-(\epsilon_{0}+t-c)-i\eta}\bigg{)}\,,$
(51)

$\Sigma_{\text{xc},ij\sigma}^{G^{\mathrm{QP}}W_{0}^{\mathrm{QP}}}=-\frac{U}{2}\delta_{ij}+\frac{U^{2}(\frac{c}{2}-t)}{2h^{\prime}}\bigg{(}\frac{1}{\omega-(\epsilon_{0}-t+(c+U)/2+h^{\prime})+i\eta}+\\\
\frac{(-1)^{i-j}e^{-3i\omega\eta}}{\omega-(\epsilon_{0}+t-(c-U)/2-h^{\prime})-i\eta}\bigg{)}\,,$
(52)

and,

$\Sigma_{\text{xc},ij\sigma}^{G^{\mathrm{QP}}W}=-\frac{U}{2}\delta_{ij}+\frac{U^{2}}{a^{2}}\bigg{(}\frac{1}{\omega-(\epsilon_{0}-t+c+U)+i\eta}+\frac{(-1)^{i-j}e^{-3i\omega\eta}}{\omega-(\epsilon_{0}+t-c)-i\eta}\bigg{)}.$
(53)

The non-interacting $\chi_{0}(1,2)=-iG(1,2^{+})G(2,1^{+})$ does not have the
same structure as the exact $\chi$ for the Hubbard dimer, when $G$ is the
exact GF . In fact,

$\chi^{GG}_{0,ij}(\omega)=(-1)^{i-j}\frac{\big{(}1+\frac{4t}{c-U}\big{)}^{4}}{2a^{4}}\times\bigg{(}\frac{1}{\omega+2t-c+2i\eta}-\frac{1}{\omega-2t+c-2i\eta}\bigg{)}+\\\
\frac{\big{(}1+\frac{4t}{c-U}\big{)}^{2}\big{(}1-\frac{4t}{c-U}\big{)}^{2}}{a^{4}}\times\bigg{(}\frac{1}{\omega-c+2i\eta}-\frac{1}{\omega+c-2i\eta}\bigg{)}\\\
+(-1)^{i-j}\frac{\big{(}1-\frac{4t}{c-U}\big{)}^{4}}{2a^{4}}\times\bigg{(}\frac{1}{\omega-2t-c+2i\eta}-\frac{1}{\omega+2t+c-2i\eta}\bigg{)}\,,$
(54)

which yields, in the b/a-b basis

$\chi_{0}^{GG}(\omega)=\begin{pmatrix}2C_{2}&0\\\
0&2C_{1}+2C_{3}\end{pmatrix}\,,$ (55)

where $C_{1}$, $C_{2}$ and $C_{3}$ correspond to the first, second and last
term in Eq. 54, respectively. The fact that the bonding-bonding matrix element
does not vanish, contratry to the exact interacting $\chi$ Eq. 46, explains
why no $f_{\mathrm{xc}}$ can be found that would link $\chi_{0}^{GG}$ and
$\chi$ in a Dyson equation.

##  2 $G\tilde{W}$ solutions

The test-charge test-electron screened interaction is defined as

$\tilde{W}(1,2)=v_{\mathrm{c}}(1,2)+\int{d(34)}\bigg{(}v_{\mathrm{c}}(1,3)\\\
+f_{\mathrm{xc}}(1,3)\bigg{)}\chi(4,2)v_{\mathrm{c}}(4,2)\,.$ (56)

The two $f_{\mathrm{xc}}$ kernels that we used in the main text are given by
the matrix equations below

$f_{\mathrm{xc}}^{\mathrm{KS}}(\omega)=\left[\chi_{0}^{\mathrm{KS}}(\omega)\right]^{-1}-\left[\chi(\omega)\right]^{-1}-v_{\mathrm{c}}\,,$
(57)
$f_{\mathrm{xc}}^{\mathrm{QP}}(\omega)=\left[\chi_{0}^{\mathrm{QP}}(\omega)\right]^{-1}-\left[\chi(\omega)\right]^{-1}-v_{\mathrm{c}}\,.$
(58)

In the bonding-antibonding basis, because of Eq. (46) the KS $\chi_{0}$ cannot
be inverted and $f_{\mathrm{xc}}^{b-b}$ is not determined. Instead,
$f_{\mathrm{xc}}^{\text{ab-
ab}}=\frac{1}{2\chi_{0,11}}-\frac{1}{2\chi_{11}}-U$, where
$\chi^{\mathrm{KS}}_{0,11}(\omega)=\chi^{\mathrm{KS}}_{0,11\uparrow}(\omega)+\chi^{\mathrm{KS}}_{0,11\downarrow}(\omega)$
and
$\chi_{11}(\omega)=\chi_{11\uparrow\uparrow}(\omega)+\chi_{11\uparrow\downarrow}(\omega)+\chi_{11\downarrow\uparrow}(\omega)+\chi_{11\downarrow\downarrow}(\omega)$.
$f_{\text{xc}}^{\text{b-b}}$ and $f_{\text{xc}}^{\text{ab-ab}}$ are the
bonding-bonding and antibonding-antibonding elements of the $f_{\text{xc}}$
matrix. This leads to

$f_{\text{xc}}^{\mathrm{KS,ab-
ab}}(\omega)=\omega^{2}\bigg{(}\frac{1}{4t}-\frac{a^{2}}{4(c+U)}\bigg{)}-t+\frac{a^{2}(c+U)}{16}-U\,,$
(59)

and

$f_{\text{xc}}^{\mathrm{QP,ab-
ab}}(\omega)=\omega^{2}\bigg{(}\frac{1}{2c-4t}-\frac{a^{2}}{4(c+U)}\bigg{)}+2t-c+\frac{a^{2}(c+U)}{16}-U\,,$
(60)

respectively for the KS and QP cases. In the symmetric Hubbard dimer,
$f_{\mathrm{xc}}$ does not have poles. Its frequency dependence is quadratic.
Comparison of Eq. (59) with Eq. (60) shows that, since
$\frac{1}{2c-4t}-\frac{a^{2}}{4(c+U)}>\frac{1}{4t}-\frac{a^{2}}{4(c+U)}$,
$f_{\mathrm{xc}}^{\mathrm{QP}}$ varies more strongly than
$f_{\mathrm{xc}}^{\mathrm{KS}}$ with $\omega$, as it is also shown in Fig. 5.

The self-energies based on these kernels read

$\Sigma_{\text{xc},ij\sigma}(\omega)=i\int{\frac{d\omega_{1}}{2\pi}G_{ij\sigma}(\omega_{1}+\omega)W_{ji}(\omega_{1})e^{3i\omega_{1}\eta}}\\\
+iU\sum_{m=1,2}\int{\frac{d\omega_{1}}{2\pi}G_{ij\sigma}(\omega_{1}+\omega)f_{xcjm}(\omega_{1})\chi_{mi}(\omega_{1})e^{3i\omega_{1}\eta}}\,,$
(61)

which yields

$\Sigma^{G^{\mathrm{KS}}\tilde{W}^{\mathrm{KS}}}_{\text{xc},ij\sigma}(\omega)=-\frac{U}{2}\delta_{ij}+(-1)^{i-j}\bigg{(}\frac{U}{4}-\frac{U(c+U)}{4ta^{2}}\bigg{)}e^{-3i\omega\eta}+\\\
\frac{U}{4a^{2}}\times\frac{(c+U)^{2}/4-4t^{2}}{t}\bigg{(}\frac{1}{\omega-(\epsilon_{0}+3t+U)+i\eta}+\frac{(-1)^{i-j}e^{-3i\omega\eta}}{\omega-(\epsilon_{0}+t-c)-i\eta}\bigg{)}\,,$
(62)
$\Sigma^{G^{\mathrm{QP}}\tilde{W}^{\mathrm{QP}}}_{\text{xc},ij\sigma}(\omega)=-\frac{U}{2}\delta_{ij}+(-1)^{i-j}\bigg{(}\frac{U}{4}-\frac{U(c+U)}{a^{2}(2c-4t)}\bigg{)}e^{-3i\omega\eta}\\\
+\frac{U}{a^{2}}\times\frac{(c+U)^{2}/4-(2t-c)^{2}}{(2c-4t)}\bigg{(}\frac{1}{\omega+(\epsilon_{0}-t+c+U)+i\eta}+\frac{(-1)^{i-j}e^{-3i\omega\eta}}{\omega-(\epsilon_{0}+t-c)-i\eta}\bigg{)}\,,$
(63)

in which $\tilde{W}$ is the TCTE screened interaction based on the exact
$\chi$ and consistent $f_{\mathrm{xc}}$ kernel.

When we use the the adiabatic approximation for $f_{\mathrm{xc}}$ to evaluate
$\chi$ and $\Sigma_{\mathrm{xc}}$, the Dyson equation becomes

$\chi_{\mathrm{adiab}}(\omega)=\bigg{(}\chi_{0}^{-1}(\omega)-f_{\text{xc}}(\omega=0)-v_{\mathrm{c}}\bigg{)}^{-1}\,.$
(64)

We calculate two different $\chi_{\mathrm{adiab}}$ depending on the choice of
$\chi_{0}$ and corresponding $f_{\mathrm{xc}}(\omega=0)$. So, we have
$\chi_{\mathrm{adiab}}^{\mathrm{KS}}$ and
$\chi_{\mathrm{adiab}}^{\mathrm{QP}}$ when $\chi_{0}^{\mathrm{KS}}$,
$f_{\mathrm{xc}}^{\mathrm{KS}}(\omega=0)$ and $\chi_{0}^{\mathrm{QP}}$,
$f_{\mathrm{xc}}^{\mathrm{QP}}(\omega=0)$ are used respectively in the
equation above. This yields

$\Sigma^{G^{\mathrm{KS}}\tilde{W}^{\mathrm{KS}}_{\text{adiab}}}_{\text{xc},ij\sigma}(\omega=0)=-\frac{U}{2}\delta_{ij}+\\\
\frac{Ut}{2\omega_{1}}\bigg{(}\frac{1}{\omega-(\epsilon_{0}+3t-(c-U)/2+\omega_{1})+i\eta}+\frac{(-1)^{i-j}e^{-3i\omega\eta}}{\omega-(\epsilon_{0}+t-(c-U/2)-\omega_{1})-i\eta}\bigg{)}\,,$
(65)

where $\omega_{1}=\sqrt{4t^{2}+2tf_{\text{Hxc}}^{\mathrm{KS}}(\omega=0)}$,
where
$f_{\text{Hxc}}^{\mathrm{KS}}(\omega=0)=f_{\text{xc},11}^{\mathrm{KS}}(\omega)-f_{\text{xc},12}^{\mathrm{KS}}(\omega=0)+2U=-2t+\frac{a^{2}(c+U)}{8}$
For the QP ingredients, we have, similarly

$\Sigma^{G^{\mathrm{QP}}\tilde{W}^{\mathrm{QP}}_{\text{adiab}}}_{\text{xc},ij\sigma}(\omega=0)=-\frac{U}{2}\delta_{ij}+\frac{U(2c-4t)}{8\omega_{2}}\\\
\times\bigg{(}\frac{1}{\omega-(\epsilon_{0}-t+(c+U)/2+\omega_{2})+i\eta}+\frac{(-1)^{i-j}e^{-3i\omega\eta}}{\omega-(\epsilon_{0}+t-(c-U/2)-\omega_{2})-i\eta}\bigg{)}\,,$
(66)

where
$\omega_{2}=\sqrt{(2t-c)^{2}+(2c-4t)f_{\text{Hxc}}^{\mathrm{QP}}(\omega=0)}$,
with
$f_{\text{Hxc}}^{\mathrm{QP}}(\omega=0)=\frac{2t-c}{2}+\frac{a^{2}(c+U)}{16}$.
Note that $\tilde{W}_{\text{adiab}}$ is the TCTE screened interaction that
includes $f_{\text{xc}}(\omega=0)$ and $\chi_{\mathrm{adiab}}(\omega)$ within
the two different schemes.

##  3 Total energy contributions for the Hubbard dimer

The xc and kinetic energy contributions to the total energy, given in the
Galitskii-Migdal formula in Eq. (1) are written in the site basis and
frequency space of the Hubbard dimer respectively as follows

$\text{E}_{\text{xc}}=-\frac{i}{2}\sum_{ij\sigma}\int_{-\infty}^{+\infty}{\frac{d\omega}{2\pi}\,\Sigma_{\text{xc},ij\sigma}(\omega)G_{ij\sigma}(\omega)e^{2i\omega\eta}}\,,$
(67)

$\text{E}_{\text{k}}=it\sum_{ij,i\neq
j,\sigma}\int_{-\infty}^{+\infty}{\frac{d\omega}{2\pi}\,G_{ij\sigma}(\omega)e^{i\omega\eta}}\,.$
(68)

##  4 Computational details

The entire framework for this work has been developed using an in-house code,
using the Julia programming language.94 For the purpose of performing the
energy integrals, we use the ‘quadgk’ library, that relies on Gauss-Kronod
quadratures.95 Additionally, to visualize our findings effectively, we rely on
the ‘Plots.jl’ library,96 coupled with the GR backend. While we provide the
analytic solutions and equations for the time-ordered quantities, the
numerical calculations have been performed using the retarded Green’s function
framework,79, 80 which yields numerically stable results for small $t$.
Retarded $G^{R}$ and $\Sigma^{R}$ are obtained from the above equations with
the usual sign changes of the imaginary infinitesimals. The total energy
contributions defined in Eq. (67) and Eq. (70) become

$\text{E}_{\text{xc}}=-\frac{1}{2\pi}\sum_{ij\sigma}\int_{-\infty}^{\mu}{d\omega\,\text{Im}(\Sigma^{R}_{\text{xc},ij\sigma}(\omega)G^{R}_{ij\sigma}(\omega))}\,,$
(69)

$\text{E}_{\text{k}}=\frac{t}{\pi}\sum_{ij,i\neq
j,\sigma}\int_{-\infty}^{\mu}{d\omega\,\text{Im}G^{R}_{ij\sigma}(\omega)}\,,$
(70)

where $\mu$ is the chemical potential.

The code of this project, called “Symmetric Hubbard Dimer”, is available at
the following address: https://gitlab.com/tsg1860938/symmetric-hubbard-dimer

## References

  * Hohenberg and Kohn 1964 Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. _Phys. Rev._ 1964, _136_ , B864–B871
  * Gilbert 1975 Gilbert, T. L. Hohenberg-Kohn theorem for nonlocal external potentials. _Phys. Rev. B_ 1975, _12_ , 2111–2120
  * Donnelly and Parr 2008 Donnelly, R. A.; Parr, R. G. Elementary properties of an energy functional of the first‐order reduced density matrix. _The Journal of Chemical Physics_ 2008, _69_ , 4431–4439
  * Levy 1979 Levy, M. Universal variational functionals of electron densities, first-order density matrices, and natural spin-orbitals and solution of the ¡i¿v¡/i¿-representability problem. _Proceedings of the National Academy of Sciences_ 1979, _76_ , 6062–6065
  * Martin et al. 2016 Martin, R.; Reining, L.; Ceperley, D. _Interacting Electrons: Theory and Computational Approaches_ ; Cambridge University Press, 2016
  * V.M. Galitskii 1950 V.M. Galitskii, A. M. Application of Quantum Field Theory Methods to the Many Body Problem. _JETP_ 1950, _7_ , 96
  * Kohn and Sham 1965 Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. _Phys. Rev._ 1965, _140_ , A1133–A1138
  * Mahan 1990 Mahan, G. _Many-particle physics_ ; Plenum Press: New York, 1990
  * Luttinger and Ward 1960 Luttinger, J. M.; Ward, J. C. Ground-State Energy of a Many-Fermion System. II. _Phys. Rev._ 1960, _118_ , 1417–1427
  * Klein 1961 Klein, A. Perturbation theory for an infinite medium of fermions. II. _Phys. Rev._ 1961, _121_ , 950–956
  * Holm 1999 Holm, B. Total Energies from $\mathit{GW}$ Calculations. _Phys. Rev. Lett._ 1999, _83_ , 788–791
  * García-González and Godby 2001 García-González, P.; Godby, R. W. Self-consistent calculation of total energies of the electron gas using many-body perturbation theory. _Phys. Rev. B_ 2001, _63_ , 075112
  * Caruso et al. 2012 Caruso, F.; Rinke, P.; Ren, X.; Scheffler, M.; Rubio, A. Unified description of ground and excited states of finite systems: The self-consistent $GW$ approach. _Phys. Rev. B_ 2012, _86_ , 081102
  * Bruneval et al. 2021 Bruneval, F.; Rodriguez-Mayorga, M.; Rinke, P.; Dvorak, M. Improved One-Shot Total Energies from the Linearized GW Density Matrix. _Journal of Chemical Theory and Computation_ 2021, _17_ , 2126–2136, PMID: 33705127
  * García-González and Godby 2002 García-González, P.; Godby, R. W. Many-Body $\mathit{GW}$ Calculations of Ground-State Properties: Quasi-2D Electron Systems and van der Waals Forces. _Phys. Rev. Lett._ 2002, _88_ , 056406
  * Hedin 1965 Hedin, L. New Method for Calculating the One-Particle Green’s Function with Application to the Electron-Gas Problem. _Phys. Rev._ 1965, _139_ , A796–A823
  * Hybertsen and Louie 1985 Hybertsen, M. S.; Louie, S. G. First-Principles Theory of Quasiparticles: Calculation of Band Gaps in Semiconductors and Insulators. _Phys. Rev. Lett._ 1985, _55_ , 1418–1421
  * Godby et al. 1987 Godby, R. W.; Schlüter, M.; Sham, L. J. Trends in self-energy operators and their corresponding exchange-correlation potentials. _Phys. Rev. B_ 1987, _36_ , 6497–6500
  * Godby et al. 1988 Godby, R. W.; Schlüter, M.; Sham, L. J. Self-energy operators and exchange-correlation potentials in semiconductors. _Phys. Rev. B_ 1988, _37_ , 10159–10175
  * Blase et al. 1995 Blase, X.; Rubio, A.; Louie, S. G.; Cohen, M. L. Quasiparticle band structure of bulk hexagonal boron nitride and related systems. _Phys. Rev. B_ 1995, _51_ , 6868–6875
  * van Schilfgaarde et al. 2006 van Schilfgaarde, M.; Kotani, T.; Faleev, S. V. Adequacy of approximations in $\mathit{GW}$ theory. _Phys. Rev. B_ 2006, _74_ , 245125
  * Kotani et al. 2007 Kotani, T.; van Schilfgaarde, M.; Faleev, S. V. Quasiparticle self-consistent $GW$ method: A basis for the independent-particle approximation. _Phys. Rev. B_ 2007, _76_ , 165106
  * Reining 2018 Reining, L. The GW approximation: content, successes and limitations. _WIREs Computational Molecular Science_ 2018, _8_ , e1344
  * Bruneval et al. 2021 Bruneval, F.; Dattani, N.; van Setten, M. J. The GW Miracle in Many-Body Perturbation Theory for the Ionization Potential of Molecules. _Frontiers in Chemistry_ 2021, _9_
  * Romaniello et al. 2009 Romaniello, P.; Guyot, S.; Reining, L. The self-energy beyond GW: Local and nonlocal vertex corrections. _The Journal of Chemical Physics_ 2009, _131_ , 154111
  * Aryasetiawan et al. 1996 Aryasetiawan, F.; Hedin, L.; Karlsson, K. Multiple Plasmon Satellites in Na and Al Spectral Functions from Ab Initio Cumulant Expansion. _Phys. Rev. Lett._ 1996, _77_ , 2268–2271
  * Guzzo et al. 2011 Guzzo, M.; Lani, G.; Sottile, F.; Romaniello, P.; Gatti, M.; Kas, J. J.; Rehr, J. J.; Silly, M. G.; Sirotti, F.; Reining, L. Valence Electron Photoemission Spectrum of Semiconductors: Ab Initio Description of Multiple Satellites. _Phys. Rev. Lett._ 2011, _107_ , 166401
  * Stan et al. 2009 Stan, A.; Dahlen, N. E.; van Leeuwen, R. Levels of self-consistency in the GW approximation. _The Journal of Chemical Physics_ 2009, _130_ , 114105
  * Minnhagen 1975 Minnhagen, P. Aspects on diagrammatic expansion for models related to a homogeneous electron gas. _J. Phys. C: Solid State Phys._ 1975, _8_ , 1535
  * Bobbert and van Haeringen 1994 Bobbert, P. A.; van Haeringen, W. Lowest-order vertex-correction contribution to the direct gap of silicon. _Phys. Rev. B_ 1994, _49_ , 10326–10331
  * Shirley 1996 Shirley, E. L. Self-consistent GW and higher-order calculations of electron states in metals. _Phys. Rev. B_ 1996, _54_ , 7758–7764
  * Grüneis et al. 2014 Grüneis, A.; Kresse, G.; Hinuma, Y.; Oba, F. Ionization Potentials of Solids: The Importance of Vertex Corrections. _Phys. Rev. Lett._ 2014, _112_ , 096401
  * Hinuma et al. 2014 Hinuma, Y.; Grüneis, A.; Kresse, G.; Oba, F. Band alignment of semiconductors from density-functional theory and many-body perturbation theory. _Phys. Rev. B_ 2014, _90_ , 155405
  * Ren et al. 2015 Ren, X.; Marom, N.; Caruso, F.; Scheffler, M.; Rinke, P. Beyond the $GW$ approximation: A second-order screened exchange correction. _Phys. Rev. B_ 2015, _92_ , 081104
  * Kutepov 2016 Kutepov, A. L. Electronic structure of Na, K, Si, and LiF from self-consistent solution of Hedin’s equations including vertex corrections. _Phys. Rev. B_ 2016, _94_ , 155101
  * Kutepov 2017 Kutepov, A. L. Self-consistent solution of Hedin’s equations: Semiconductors and insulators. _Phys. Rev. B_ 2017, _95_ , 195120
  * Pavlyukh et al. 2016 Pavlyukh, Y.; Uimonen, A.-M.; Stefanucci, G.; van Leeuwen, R. Vertex Corrections for Positive-Definite Spectral Functions of Simple Metals. _Phys. Rev. Lett._ 2016, _117_ , 206402
  * Maggio and Kresse 2017 Maggio, E.; Kresse, G. GW Vertex Corrected Calculations for Molecular Systems. _Journal of Chemical Theory and Computation_ 2017, _13_ , 4765–4778, PMID: 28873298
  * Wang and Ren 2022 Wang, Y.; Ren, X. Vertex effects in describing the ionization energies of the first-row transition-metal monoxide molecules. _The Journal of Chemical Physics_ 2022, _157_ , 214115
  * Langreth and Perdew 1975 Langreth, D.; Perdew, J. The exchange-correlation energy of a metallic surface. _Solid State Communications_ 1975, _17_ , 1425–1429
  * Langreth and Perdew 1977 Langreth, D. C.; Perdew, J. P. Exchange-correlation energy of a metallic surface: Wave-vector analysis. _Phys. Rev. B_ 1977, _15_ , 2884–2901
  * Runge and Gross 1984 Runge, E.; Gross, E. K. U. Density-Functional Theory for Time-Dependent Systems. _Phys. Rev. Lett._ 1984, _52_ , 997–1000
  * Del Sole et al. 1994 Del Sole, R.; Reining, L.; Godby, R. W. GW$\Gamma$ approximation for electron self-energies in semiconductors and insulators. _Phys. Rev. B_ 1994, _49_ , 8024–8028
  * Reining et al. 2002 Reining, L.; Olevano, V.; Rubio, A.; Onida, G. Excitonic Effects in Solids Described by Time-Dependent Density-Functional Theory. _Phys. Rev. Lett._ 2002, _88_ , 066404
  * Overhauser 1971 Overhauser, A. W. Simplified Theory of Electron Correlations in Metals. _Phys. Rev. B_ 1971, _3_ , 1888–1898
  * Petrillo and Sacchetti 1988 Petrillo, C.; Sacchetti, F. Electron-gas self-energy at metallic density. _Phys. Rev. B_ 1988, _38_ , 3834–3840
  * Mahan and Sernelius 1989 Mahan, G. D.; Sernelius, B. E. Electron-electron interactions and the bandwidth of metals. _Phys. Rev. Lett._ 1989, _62_ , 2718–2720
  * Hybertsen and Louie 1986 Hybertsen, M. S.; Louie, S. G. Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies. _Phys. Rev. B_ 1986, _34_ , 5390–5413
  * Hindgren and Almbladh 1997 Hindgren, M.; Almbladh, C.-O. Improved local-field corrections to the ${G}_{0}W$ approximation in jellium: Importance of consistency relations. _Phys. Rev. B_ 1997, _56_ , 12832–12839
  * Schmidt et al. 2017 Schmidt, P. S.; Patrick, C. E.; Thygesen, K. S. Simple vertex correction improves $GW$ band energies of bulk and two-dimensional crystals. _Phys. Rev. B_ 2017, _96_ , 205206
  * Hung et al. 2016 Hung, L.; da Jornada, F. H.; Souto-Casares, J.; Chelikowsky, J. R.; Louie, S. G.; Öğüt, S. Excitation spectra of aromatic molecules within a real-space $GW$-BSE formalism: Role of self-consistency and vertex corrections. _Phys. Rev. B_ 2016, _94_ , 085125
  * Olsen et al. 2019 Olsen, T.; Patrick, C. E.; Bates, J. E.; Ruzsinszky, A.; Thygesen, K. S. Beyond the RPA and GW methods with adiabatic xc-kernels for accurate ground state and quasiparticle energies. _npj Computational Materials_ 2019, _5_ , 106
  * Chen and Pasquarello 2015 Chen, W.; Pasquarello, A. Accurate band gaps of extended systems via efficient vertex corrections in $GW$. _Phys. Rev. B_ 2015, _92_ , 041115
  * Tal et al. 2021 Tal, A.; Chen, W.; Pasquarello, A. Vertex function compliant with the Ward identity for quasiparticle self-consistent calculations beyond $GW$. _Phys. Rev. B_ 2021, _103_ , L161104
  * Gross and Kohn 1985 Gross, E. K. U.; Kohn, W. Local density-functional theory of frequency-dependent linear response. _Phys. Rev. Lett._ 1985, _55_ , 2850–2852
  * Tokatly and Pankratov 2001 Tokatly, I. V.; Pankratov, O. Many-Body Diagrammatic Expansion in a Kohn-Sham Basis: Implications for Time-Dependent Density Functional Theory of Excited States. _Phys. Rev. Lett._ 2001, _86_ , 2078–2081
  * Bruneval et al. 2005 Bruneval, F.; Sottile, F.; Olevano, V.; Del Sole, R.; Reining, L. Many-Body Perturbation Theory Using the Density-Functional Concept: Beyond the $GW$ Approximation. _Phys. Rev. Lett._ 2005, _94_ , 186402
  * Gatti et al. 2007 Gatti, M.; Olevano, V.; Reining, L.; Tokatly, I. V. Transforming Nonlocality into a Frequency Dependence: A Shortcut to Spectroscopy. _Phys. Rev. Lett._ 2007, _99_ , 057401
  * Botti et al. 2007 Botti, S.; Schindlmayr, A.; Sole, R. D.; Reining, L. Time-dependent density-functional theory for extended systems. _Reports on Progress in Physics_ 2007, _70_ , 357
  * Singhal and Callaway 1976 Singhal, S. P.; Callaway, J. Exchange correction to the dielectric function in the local exchange approximation. _Phys. Rev. B_ 1976, _14_ , 2347–2351
  * Hybertsen and Louie 1987 Hybertsen, M. S.; Louie, S. G. Ab initio static dielectric matrices from the density-functional approach. I. Formulation and application to semiconductors and insulators. _Phys. Rev. B_ 1987, _35_ , 5585–5601
  * Hellgren and Baguet 2023 Hellgren, M.; Baguet, L. Strengths and limitations of the adiabatic exact-exchange kernel for total energy calculations. _The Journal of Chemical Physics_ 2023, _158_ , 184107
  * Sottile et al. 2003 Sottile, F.; Olevano, V.; Reining, L. _Phys. Rev. Lett._ 2003, _91_ , 056402
  * Sharma et al. 2011 Sharma, S.; Dewhurst, J. K.; Sanna, A.; Gross, E. K. U. Bootstrap Approximation for the Exchange-Correlation Kernel of Time-Dependent Density-Functional Theory. _Phys. Rev. Lett._ 2011, _107_
  * Adragna et al. 2003 Adragna, G.; Del Sole, R.; Marini, A. _Phys. Rev. B_ 2003, _68_ , 165108
  * Marini et al. 2003 Marini, A.; Del Sole, R.; Rubio, A. _Phys. Rev. Lett._ 2003, _91_ , 256402
  * Rigamonti et al. 2015 Rigamonti, S.; Botti, S.; Veniard, V.; Draxl, C.; Reining, L.; Sottile, F. Estimating Excitonic Effects in the Absorption Spectra of Solids: Problems and Insight from a Guided Iteration Scheme. _Phys. Rev. Lett._ 2015, _114_ , 146402
  * Levy and Perdew 1985 Levy, M.; Perdew, J. P. Hellmann-Feynman, virial, and scaling requisites for the exact universal density functionals. Shape of the correlation potential and diamagnetic susceptibility for atoms. _Phys. Rev. A_ 1985, _32_ , 2010–2021
  * Jiang et al. 2020 Jiang, K.; Mosquera, M. A.; Oueis, Y.; Wasserman, A. Virial relations in density embedding. _International Journal of Quantum Chemistry_ 2020, _120_ , e26204
  * Kim et al. 2013 Kim, M.-C.; Sim, E.; Burke, K. Understanding and Reducing Errors in Density Functional Calculations. _Phys. Rev. Lett._ 2013, _111_ , 073003
  * Ren et al. 2012 Ren, X.; Rinke, P.; Joas, C.; Scheffler, M. Random-phase approximation and its applications in computational chemistry and materials science. _Journal of Materials Science_ 2012, _47_ , 7447
  * Savin et al. 2001 Savin, A.; Colonna, F.; Allavena, M. Analysis of the linear response function along the adiabatic connection from the Kohn–Sham to the correlated system. _The Journal of Chemical Physics_ 2001, _115_ , 6827–6833
  * Romaniello et al. 2012 Romaniello, P.; Bechstedt, F.; Reining, L. Beyond the $GW$ approximation: Combining correlation channels. _Phys. Rev. B_ 2012, _85_ , 155131
  * Carrascal et al. 2015 Carrascal, D. J.; Ferrer, J.; Smith, J. C.; Burke, K. The Hubbard dimer: a density functional case study of a many-body problem. _Journal of Physics: Condensed Matter_ 2015, _27_ , 393001
  * Aryasetiawan and Gunnarsson 2002 Aryasetiawan, F.; Gunnarsson, O. Exchange-correlation kernel in time-dependent density functional theory. _Phys. Rev. B_ 2002, _66_ , 165119
  * Coveney and Tew 2023 Coveney, C. J. N.; Tew, D. P. A Regularized Second-Order Correlation Method from Green’s Function Theory. _Journal of Chemical Theory and Computation_ 2023, _19_ , 3915–3928
  * J.Hubbard 1963 J.Hubbard, Electron correlations in narrow energy bands. _Proc. R. Soc. Lond. A_ 1963, _276_ , 238
  * J.Hubbard 1964 J.Hubbard, Electron correlations in narrow energy bands. _Proc. R. Soc. Lond. A_ 1964, _277_ , 237
  * Spataru et al. 2004 Spataru, C. D.; Benedict, L. X.; Louie, S. G. Ab initio calculation of band-gap renormalization in highly excited GaAs. _Phys. Rev. B_ 2004, _69_ , 205204
  * Honet et al. 2022 Honet, A.; Henrard, L.; Meunier, V. Exact and many-body perturbation solutions of the Hubbard model applied to linear chains. _AIP Advances_ 2022, _12_ , 035238
  * Holm and von Barth 1998 Holm, B.; von Barth, U. Fully self-consistent $\mathrm{GW}$ self-energy of the electron gas. _Phys. Rev. B_ 1998, _57_ , 2108–2117
  * Perdew et al. 2017 Perdew, J. P.; Yang, W.; Burke, K.; Yang, Z.; Gross, E. K. U.; Scheffler, M.; Scuseria, G. E.; Henderson, T. M.; Zhang, I. Y.; Ruzsinszky, A.; Peng, H.; Sun, J.; Trushin, E.; Görling, A. Understanding band gaps of solids in generalized Kohn–Sham theory. _Proceedings of the National Academy of Sciences_ 2017, _114_ , 2801–2806
  * Heyd et al. 2003 Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. _The Journal of Chemical Physics_ 2003, _118_ , 8207–8215
  * Aryasetiawan et al. 2012 Aryasetiawan, F.; Sakuma, R.; Karlsson, K. $GW$ approximation with self-screening correction. _Phys. Rev. B_ 2012, _85_ , 035106
  * Hedin et al. 1967 Hedin, L.; Lundqvist, B.; Lundqvist, S. New structure in the single-particle spectrum of an electron gas. _Solid State Communications_ 1967, _5_ , 237 – 239
  * Bergerse.B et al. 1973 Bergerse.B,; Kus, F. W.; Blomberg, C. SINGLE-PARTICLE GREENS FUNCTION IN ELECTRON-PLASMON APPROXIMATION. _Canadian J. Phys._ 1973, _51_ , 102–110
  * Guzzo et al. 2014 Guzzo, M.; Kas, J. J.; Sponza, L.; Giorgetti, C.; Sottile, F.; Pierucci, D.; Silly, M. G.; Sirotti, F.; Rehr, J. J.; Reining, L. Multiple satellites in materials with complex plasmon spectra: From graphite to graphene. _Phys. Rev. B_ 2014, _89_ , 085425
  * Nelson et al. 2007 Nelson, W.; Bokes, P.; Rinke, P.; Godby, R. W. Self-interaction in Green’s-function theory of the hydrogen atom. _Phys. Rev. A_ 2007, _75_ , 032505
  * Fernandez 2009 Fernandez, J. J. $GW$ calculations in an exactly solvable model system at different dilution regimes: The problem of the self-interaction in the correlation part. _Phys. Rev. A_ 2009, _79_ , 052513
  * Botti et al. 2005 Botti, S.; Fourreau, A.; Nguyen, F. m. c.; Renault, Y.-O.; Sottile, F.; Reining, L. Energy dependence of the exchange-correlation kernel of time-dependent density functional theory: A simple model for solids. _Phys. Rev. B_ 2005, _72_ , 125203
  * Minnhagen 1974 Minnhagen, P. Vertex correction calculations for an electron gas. _Journal of Physics C: Solid State Physics_ 1974, _7_ , 3013
  * Lewis and Berkelbach 2019 Lewis, A. M.; Berkelbach, T. C. Vertex Corrections to the Polarizability Do Not Improve the GW Approximation for the Ionization Potential of Molecules. _Journal of Chemical Theory and Computation_ 2019, _15_ , 2925–2932
  * Perdew and Norman 1982 Perdew, J. P.; Norman, M. R. Electron removal energies in Kohn-Sham density-functional theory. _Phys. Rev. B_ 1982, _26_ , 5445–5450
  * Bezanson et al. 2017 Bezanson, J.; Edelman, A.; Karpinski, S.; Shah, V. B. Julia: A fresh approach to numerical computing. _SIAM Review_ 2017, _59_ , 65–98
  * Laurie 1997 Laurie, D. P. Calculation of Gauss-Kronrod Quadrature Rules. _Mathematics of Computation_ 1997, _66_ , 1133–1145
  * Christ et al. 2023 Christ, S.; Schwabeneder, D.; Rackauckas, C.; Borregaard, M. K.; Breloff, T. Plots.jl – a user extendable plotting API for the julia programming language. 2023, _11_ , 5