# Correlation energy for the homogeneous electron gas: exact Bethe-Salpeter
solution and new approximate evaluation

Emanuele Maggio  Georg Kresse<EMAIL_ADDRESS>University of Vienna,
Faculty of Physics and Center for Computational Materials Science, Sensengasse
8/12, A-1090 Vienna, Austria

###### Abstract

The correlation energy of the homogeneous electron gas is evaluated by solving
the Bethe-Salpeter equation (BSE) beyond the Tamm-Dancoff approximation for
the electronic polarisation propagator. The BSE is expected to improve upon
the random phase approximation, owing to the inclusion of exchange diagrams.
For instance, since the BSE reduces in second order to Møller-Plesset
perturbation theory, it is self-interaction free in second order. Results for
the correlation energy are compared with Quantum Monte Carlo benchmarks and
excellent agreement is observed. For low densities, however, we find imaginary
eigenmodes in the polarisation propagator. To avoid the occurrence of
imaginary eigenmodes, an approximation to the BSE kernel is proposed, which
allows to completely remove this issue in the low electron density region. We
refer to this approximation as the random phase approximation with screened
exchange (RPAsX). We show that this approximation even slightly improves upon
the standard BSE kernel.

yyy

###### pacs:

xxx

##  I Introduction

The study of the homogeneous electron gas (HEG), as a model in condensed
matter theory, has a long tradition. In fact, it allows to focus on the
properties of the many electron system without complications arising from the
discretised lattice symmetry. Instead, the neutralising positive charges are
uniformly spread out and non-responsive. Even if interactions among electrons
are not present in the model, its pair correlation function is not constant.
This is a consequence of the Pauli principle affecting same-spin electrons.
The exchange energy introduced this way promotes a spin polarisation within
the electron gas and, at the same time, increases the electrons’ kinetic
energy Giuliani and Vignale (2005). As the Coulomb repulsion among electrons
is turned on, correlations between particles (prevalently of opposite spin)
Endo  _et al._ (1999) set in. This mechanism produces even stronger deviations
from the free-electron case. One of the reasons for a continued interest in
this model is that it provides density functional theory with a natural
starting point for the unknown exchange and correlation potential.

An attempt to supersede the mean-field description of the HEG dates back to
the seminal work by Hubbard Hubbard (1958). There, a local field factor was
introduced to compute the system’s response functions. Albeit progress has
been made to evaluate this quantity, whose asymptotic behaviour is known
exactly Gori-Giorgi  _et al._ (2000), its complete characterisation is still
lacking Sturm and Gusarov (2000). On the other hand, the field factor can be
related to the irreducible electron-hole scattering amplitude. This appears in
the integral (and recursive) equation for the polarisation propagator, known
as the Bethe-Salpeter equation (BSE). Therefore by solving the BSE one can
gain access to exchange and correlation properties. The BSE kernel $I$ itself
is related to the irreducible self-energy through the relation
$I=-\frac{1}{i}\frac{\delta\Sigma}{\delta G}$. Hence, the degree of
sophistication in the BSE approach can be systematically improved by including
more self-energy diagrams at the many-body perturbation theory (MBPT) level.
In this study, we implement a computational scheme to solve the BSE when its
kernel is derived from the $GW_{0}$ approximation for the self-energy. We then
evaluate the correlation energy of the HEG and compare it against Quantum
Monte Carlo benchmarks.

The $GW_{0}$ and the Random Phase Approximation (RPA) include topologically
equivalent energy diagrams Caruso  _et al._ (2013). In both cases, a sub-class
of self-energy diagrams, the so-called ”bubble” diagrams are summed up to
infinite order Stefanucci and van Leeuwen (2013). While the RPA approach has
witnessed alternate fortune over the years, there has been a strong resurgence
of interest in the field recently Furche (2001a); Miyake  _et al._ (2002);
Niquet  _et al._ (2003); Jiang and Engel (2007); Furche (2008); Scuseria  _et
al._ (2008); Harl and Kresse (2008); Harl  _et al._ (2010); Toulouse  _et al._
(2009, 2010); Jansen  _et al._ (2010). Application of the RPA method to
molecules Furche (2001a); Olsen and Thygesen (2013); Heßelmann and Görling
(2010); Macher  _et al._ (2014); Teale  _et al._ (2010); Bleiziffer  _et al._
(2012) and solid state systems Miyake  _et al._ (2002); Harl  _et al._ (2010);
Kaltak  _et al._ (2014); Grüneis  _et al._ (2009); Schimka  _et al._ (2010);
García-González  _et al._ (2007) has highlighted its capabilities but also its
limitations. The main advantage of RPA related schemes stems from the seamless
inclusion of long-range dispersion (typically not included in standard DFT
potentials) Harl and Kresse (2008). On the other hand, systematic
underestimation of binding energies Furche (2001a); Harl and Kresse (2009);
Ren  _et al._ (2011) and non-physical features for dissociation curves in
diatomic systems Mori-Sánchez  _et al._ (2012); Henderson and Scuseria (2010)
are some shortcomings of the RPA. To improve on this, the BSE has been
exploited to describe hydrogen dissociation Olsen and Thygesen (2014), to gain
access to optical properties of semiconductors and insulators Albrecht  _et
al._ (1998); Rohlfing and Louie (2000); Yan  _et al._ (2012) and to study
excitonic effects in extended Gatti and Sottile (2013) and molecular systems
Del Puerto  _et al._ (2006); Duchemin and Blase (2013). To the best of our
knowledge, the vast majority of previous studies not concerned with exciton
characterisation have been carried out on DFT reference states; this is
rectified herein.

To assess the correlation energy of the HEG the adiabatic connection (AC)
formalism will be employed in this study. From its inception, the AC formalism
has been applied to obtain the correlation energy starting from reference
states evaluated either using Hartree-Fock Linderberg (1967) or Kohn-Sham
Gunnarsson and Lundqvist (1976); Langreth and Perdew (1977); Harris and Jones
(1974) frameworks. In either case, the fluctuation-dissipation theorem is
exploited to relate the ground state correlation energy to the system’s linear
response functions integrated over the AC path.

The study of the HEG is challenging for many reasons. Among others,
correlations associated with large momentum transfer interactions Singwi  _et
al._ (1968) and between spin-parallel electrons Endo  _et al._ (1999) are
overestimated by RPA. The computational scheme employed here addresses these
issues by including screening effects and exchange contributions to infinite
order in perturbation theory.

The manuscript is structured as follows: we first derive an expression for the
correlation energy in terms of the polarisation propagator. Then in subsection
II.2 the implementation of the BSE as a non-Hermitian generalised eigenvalue
problem is described. Properties peculiar to the Bloch representation are
discussed at the end of the theory section. Computational details are reported
in section III; finally in section IV the computational scheme proposed and
related approximations are put to fruition to assess the correlation energy of
the HEG.

##  II Theory

###  II.1 Exchange and Correlation Energy evaluated along the AC path

The adiabatic connection formalism was originally introduced in the DFT
framework Gunnarsson and Lundqvist (1976); Langreth and Perdew (1977); Harris
and Jones (1974) in order to provide a compact expression for the exchange and
correlation energy. In DFT this quantity is given by the sum of two
contributions: the electron-electron interaction and the difference in the
kinetic energy between the ‘physical’ system and the Kohn-Sham system. These
terms are related to the two body density matrix and to the one body density
matrix respectively. Similarly, we can also decompose the exchange-correlation
energy into one- and two-body contributions. The two-body part of the
correlation energy can be estimated solving the BSE, as described in
subsection II.2. The impact of one-body contributions has been the subject of
recent studies Klimeš  _et al._ (2015); Ren  _et al._ (2011). In this section
we compactly derive an expression for the correlation energy and discuss the
main assumptions required, the details are reported in Appendix A.

The system’s Hamiltonian in second quantised form along the AC path reads:

$\hat{\mathcal{H}}_{\lambda}=\int
d\mathbf{x}d\mathbf{x}^{\prime}\hat{\psi}^{\dagger}(\mathbf{x})\langle\mathbf{x}|\hat{h}|\mathbf{x^{\prime}}\rangle\hat{\psi}(\mathbf{x}^{\prime})+\frac{\lambda}{2}\int
d\mathbf{x}d\mathbf{x}^{\prime}v(\mathbf{x},\mathbf{x}^{\prime})\hat{\psi}^{\dagger}(\mathbf{x})\hat{\psi}^{\dagger}(\mathbf{x}^{\prime})\hat{\psi}(\mathbf{x}^{\prime})\hat{\psi}(\mathbf{x})+\int
d\mathbf{x}d\mathbf{x}^{\prime}\hat{\psi}^{\dagger}(\mathbf{x})\langle\mathbf{x}|\hat{V}^{\rm{xc}}_{\lambda}|\mathbf{x^{\prime}}\rangle\hat{\psi}(\mathbf{x}^{\prime}).$
(1)

The field operators $\hat{\psi}(\mathbf{x})$ and their adjoint are given in
the ordinary Schrödinger picture; the Coulomb potential is defined in the
spin-space basis as
$v(\mathbf{x,x^{\prime}})=\frac{\delta_{\sigma,\sigma^{\prime}}}{|\mathbf{r-r^{\prime}}|}$
with the usual space-time notation $(\mathbf{x},t)=(\mathbf{r},\sigma,t)$ for
the space $(\mathbf{r})$, spin $(\sigma)$ and time $(t)$ variables. The one-
particle operator $\hat{V}^{\rm{xc}}_{\lambda}$ can be any approximate
exchange correlation potential, even one that is not diagonal in real space,
for instance, an Hermitian energy independent approximation for the self-
energy $\Sigma_{\rm{xc}}(\omega)$.

Along the AC path the many body interactions are progressively switched on,
inducing changes in the correlations between electrons previously captured
only by the exchange-correlation potential $\hat{V}^{\rm{xc}}_{0}$. The latter
also varies with varying coupling constant and different switching can be
realised for $\hat{V}^{\rm{xc}}_{\lambda}$. These are either designed to keep
the electron density fixed along the AC path Langreth and Perdew (1977), or
assume a linear dependence on $\lambda$ Harris (1984); van Aggelen  _et al._
(2014). Here we use the linear switching
$\hat{V}^{\rm{xc}}_{\lambda}=(1-\lambda)\hat{V}^{\rm{xc}}_{0}$.

The correlation energy can be defined in full generality as the difference
between the expectation values of the interacting Hamiltonian acting on its
ground state and on the ground state of the non-interacting system:

$E_{c}=\langle\Psi_{1}|\hat{\mathcal{H}}_{1}|\Psi_{1}\rangle-\langle\Psi_{0}|\hat{\mathcal{H}}_{1}|\Psi_{0}\rangle.$
(2)

If the symmetry of the ground state does not change along the AC path, then
$|\Psi_{0}\rangle$ evolves into $|\Psi_{1}\rangle$ for non-degenerate ground
states thanks to the Gell-Mann and Low theorem. It is easy to see that with
$E_{\lambda}\equiv\langle\Psi_{\lambda}|\mathcal{H}_{\lambda}|\Psi_{\lambda}\rangle$,
$\langle\Psi_{0}|\mathcal{H}_{1}|\Psi_{0}\rangle=E_{0}+\langle\Psi_{0}|\frac{1}{2}v-\hat{V}^{\rm{xc}}_{0}|\Psi_{0}\rangle$.
Then the correlation energy can be expressed as (see also Klimeš  _et al._
(2015); Toulouse  _et al._ (2010)):

$\displaystyle E_{c}$
$\displaystyle=E_{1}-E_{0}-\bigg{\langle}\Psi_{0}\left|\frac{\partial{\mathcal{H}_{\lambda}}}{\partial{\lambda}}\right|\Psi_{0}\bigg{\rangle}$
$\displaystyle=\int_{0}^{1}d\lambda\frac{dE_{\lambda}}{d\lambda}-\left.\frac{dE_{\lambda}}{d\lambda}\right|_{\lambda=0},$
(3)

where we have assumed that the Hellmann-Feynman theorem holds in order to set
$\frac{\partial}{\partial{\lambda}}\langle\Psi_{\lambda}\left|\mathcal{H}_{\lambda}\right|\Psi_{\lambda}\rangle=\langle\Psi_{\lambda}\left|\frac{\partial{\mathcal{H}_{\lambda}}}{\partial{\lambda}}\right|\Psi_{\lambda}\rangle.$
The previous equation reproduces the standard expression found in the
literature Toulouse  _et al._ (2010).

The calculation of the correlation energy thus requires to evaluate the term
$\frac{dE_{\lambda}}{d\lambda}\equiv\dot{E}_{\lambda}$ along the AC path,
where the dependence on the coupling constant is retained by the ground state
wave function. The derivative of the ground state energy can be written as:

$\displaystyle\frac{dE_{\lambda}}{d\lambda}=$
$\displaystyle\frac{1}{2}\int\int
d\mathbf{x}d\mathbf{x}^{\prime}v(\mathbf{x},\mathbf{x}^{\prime})\langle\psi^{\dagger}(\mathbf{x})\psi^{\dagger}(\mathbf{x}^{\prime})\psi(\mathbf{x}^{\prime})\psi(\mathbf{x})\rangle_{\lambda}$
$\displaystyle-\int\int
d\mathbf{x}d\mathbf{x}^{\prime}\langle\mathbf{x}^{\prime}|\hat{V}^{\rm{xc}}_{0}|\mathbf{x}\rangle\langle\psi^{\dagger}(\mathbf{x}^{\prime})\psi(\mathbf{x})\rangle_{\lambda}.$
(4)

The last term in the previous equation represents contributions related to the
change of the one-particle Green’s function to the correlation energy Klimeš
_et al._ (2015); Ren  _et al._ (2011). In the present study, we disregard this
term and concentrate on the first term on the r.h.s. of Eq. (missing)4. In
fact, in second order, the last term is exactly zero for the HEG (see e.g.
chapter 22 in Ref.Gross  _et al._ , 1991), suggesting that its contribution
should be small.

It is shown in appendix A that the expression above can be recast in terms of
the four point polarisation propagator $P_{\lambda}$ Toyoda  _et al._ (1975),
fulfilling the BSE as detailed in the next section. The correlation energy
then reads:

$\begin{split}E_{c}&=-\frac{1}{2}\int_{0}^{1}d\lambda\int\frac{d\omega}{2\pi}\int\prod_{i=1}^{4}d\mathbf{x}_{i}v(\mathbf{x}_{1},\mathbf{x}_{4})\\\
&\delta(\mathbf{x}_{4}-\mathbf{x}_{2})\delta(\mathbf{x}_{3}-\mathbf{x}_{1})\left[P_{\lambda}(\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{3},\mathbf{x}_{4};\omega)-\right.\\\
&\left.-P_{0}(\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{3},\mathbf{x}_{4};\omega)\right].\\\
\end{split}$ (5)

With the aid of the polarisation propagator’s spectral representation obtained
in appendix B, it is possible to carry out the frequency integration
analytically, if the inter-particle interaction is not frequency dependent
Fetter and Walecka (1971).

###  II.2 Bethe-Salpeter equation along the AC path

The Bethe-Salpeter equation $P={P}_{0}+{P}_{0}IP$ relates the fully
interacting polarisation propagator $P$ to the independent particle-hole
propagator $P_{0}$, invoking an interaction kernel $I$. The kernel describes
all possible interactions between particles and holes Strinati (1988). Solving
this equation is generally a formidable task, which warrants the use of
approximations. A commonly used one is to make the interaction kernel
instantaneous and to include only Hartree and screened exchange interactions
in the kernel $I(1,2,3,4)=V(1,2,3,4)-W_{0}(2,1,3,4)$ Albrecht  _et al._
(1998); Rohlfing and Louie (2000); Hanke and Sham (1979, 1980). Using the
conventional space-time notation (where the superscript ’+’ denotes positive
infinitesimal time shifts), one can write:

$\displaystyle V(1,2,3,4)=v(1,4^{+})\delta(4,2^{+})\delta(3,1^{+})$ (6)
$\displaystyle W_{0}(1,2,3,4)=w(1,4^{+})\delta(4,2^{+})\delta(3,1^{+}).$ (7)

Here $V$ and $W_{0}$ are the four-point generalisations of the bare Coulomb
($v$) and screened ($w$) interactions, and the shorthand $(i)=({\bf
r}_{i},\sigma_{i},t_{i})$ refers to space, spin and time degrees of freedom.
The adiabatic switching of many-body interactions (presented in the previous
section) requires the irreducible kernel to be linearly scaled by the coupling
constant $\lambda$ and, correspondingly, the polarisation propagator
$P_{\lambda}$ to be evaluated for each point along the AC path. The resulting
BSE will then be parametrically dependent on two quantities: the coupling
constant $\lambda$ and the frequency variable $\omega$, provided that the
kernel is static. This approximation for the kernel implies that
$w(\mathbf{x}_{1},\mathbf{x}_{2})$ in Eq. (missing)7 is evaluated as
$w(\mathbf{x}_{1},\mathbf{x}_{2})=\\{v^{\frac{1}{2}}\cdot\epsilon^{-1}(\omega=0)\cdot
v^{\frac{1}{2}}\\}(\mathbf{x}_{1},\mathbf{x}_{2})$, where internal space
variables are integrated over and the dielectric function in this study is
computed at the RPA level.

It is convenient to introduce the ‘particle-hole’ basis:
$\Phi_{M}^{R}(\mathbf{x,x^{\prime}})=\psi_{a}(\mathbf{x})\psi^{*}_{i}(\mathbf{x}^{\prime})$
and
$\Phi_{M}^{A}(\mathbf{x,x^{\prime}})=\psi_{i}(\mathbf{x})\psi^{*}_{a}(\mathbf{x}^{\prime})$,
for the resonant $(R)$ and antiresonant $(A)$ component labelled by the
“superindex” $M\equiv(i,a)$. The BSE can then be recast in this basis to
yield:

$\displaystyle\mathbf{P}_{\lambda}(\omega)=\mathbf{P}_{0}(\omega)+\mathbf{P}_{0}(\omega)*\mathbf{I}_{\lambda}(\omega=0)*\mathbf{P}_{\lambda}(\omega).$
(8)

The symbol ‘$*$’ indicates the usual matrix product in the chosen orbital
basis. The independent particle-hole propagator is diagonal in this basis and
can be expressed as:

$\displaystyle P_{0}(\mathbf{x}_{1},$
$\displaystyle\mathbf{x}_{2},\mathbf{x}_{1}^{\prime},\mathbf{x}_{2}^{\prime};\omega)=\operatorname{tr}\left\\{\mathbf{P}_{0}(\omega)\right\\}=$
$\displaystyle\sum_{i,a}f_{i}\bar{f}_{a}\left(\dfrac{\psi_{a}(\mathbf{x}_{1})\psi^{*}_{i}(\mathbf{x}_{1}^{\prime})\psi_{i}(\mathbf{x}_{2})\psi^{*}_{a}(\mathbf{x}_{2}^{\prime})}{\omega-\Delta
E_{ia}}+\right.$
$\displaystyle\left.+\dfrac{\psi_{i}(\mathbf{x}_{1})\psi^{*}_{a}(\mathbf{x}_{1}^{\prime})\psi_{a}(\mathbf{x}_{2})\psi^{*}_{i}(\mathbf{x}_{2}^{\prime})}{-\omega-\Delta
E_{ia}}\right),$ (9)

where $f_{n}$ is the (fractional) occupancy for the orbital $n$ and
$\bar{f}_{n}=1-f_{n}$, the indices $i$ and $a$ in the summations above go over
the occupied and unoccupied single particle states respectively. $\Delta
E_{ia}=(\epsilon_{a}-\epsilon_{i})$ is the independent particle energy
difference.

Since the BSE kernel is static, it is possible to invert the matrix equation
Eq. (missing)8 for each frequency point. In appendix B it is shown that this
is, however, not necessary since the spectral decomposition for the
polarisation propagator $\mathbf{P}_{\lambda}(\omega)$ can be constructed
solving the non-Hermitian generalised eigenvalue problem (EVP):

$\displaystyle\begin{pmatrix}\mathbf{A}&\mathbf{B}\\\
\mathbf{B}^{*}&\mathbf{A}^{*}\end{pmatrix}\begin{pmatrix}\mathbf{X}\\\
\mathbf{Y}\end{pmatrix}=\begin{pmatrix}\mathbf{\Omega}&\mathbf{0}\\\
\mathbf{0}&-\mathbf{\Omega}\end{pmatrix}\begin{pmatrix}\mathbf{X}\\\
\mathbf{Y}\end{pmatrix},$ (10)

which does not depend on $\omega$. The matrix elements for $\mathbf{A}$ and
$\mathbf{B}$ are given by Fuchs  _et al._ (2008):

$\displaystyle A_{iajb}=$ $\displaystyle(\Delta E_{ia})+\lambda\langle
aj|V|ib\rangle-\lambda\langle aj|W|bi\rangle,$ (11) $\displaystyle B_{iajb}=$
$\displaystyle\lambda\langle ab|V|ij\rangle-\lambda\langle ab|W|ji\rangle.$
(12)

Their time-ordered diagrammatic representation is shown in Figure 1.

The dimension of this non-Hermitian eigenvalue problem is given by the sum of
particle-hole and hole-particle transitions. It is possible, however, to
reduce the matrix size to a single particle-hole subspace by exploiting the
time-reversal symmetry of the system. This statement will be proved in the
next section where we introduce a suitable basis for the matrix representation
as in Ref. Sander  _et al._ (2015).

Fig. 1: Goldstone diagrams for the particle-hole interaction included in the
BSE kernel: particle states are shown as right arrows $\rightarrow$ (labelled
by unoccupied ground state orbitals, $a,b$) and holes by left arrows
$\leftarrow$ (labelled by occupied orbitals $i,j$). The bare and screened
Coulomb interaction is represented by a simple wavy line and a double wavy
line respectively. Diagrams (A1) and (B1) represent the resonant and
antiresonant particle-hole annihilation and simultaneous creation, whereas
diagrams (A2) and (B2) represent the resonant and antiresonant particle-hole
scattering process.

###  II.3 Bloch representation

For extended systems it is convenient to switch to a reciprocal space
representation of the electronic problem, in order to capture more readily its
translational invariance. A fully homogeneous system has the highest
translational symmetry, which translates into a spherical Fermi surface. In
this study, we choose to discretise the infinite, translational invariant
system into a unit cell of given symmetry and subject to periodic boundary
conditions. This approach is akin to the QMC benchmarks against which we will
compare the outcome of our calculations. As a consequence, this discretisation
changes the Fermi surface, which is now replaced by a convex polyhedron: this
tends to the exact Fermi surface as the number of unit cell replicas (i.e.
$\mathbf{k}$-point sampling the reciprocal space) in our calculations is
increased. Each single particle quantum number $a,b,i,j$ can be mapped into a
set specifying the particles’ band (valence band for the hole states and
conduction band for particle states) as well as their Bloch wave-vectors,
$\mathbf{k},\mathbf{k^{\prime}}$. We have:

$\displaystyle\begin{cases}a\to(a,\mathbf{k^{\prime}});&i\to(i,\mathbf{k^{\prime}}\pm\mathbf{q})\\\
b\to(b,\mathbf{k});&j\to(j,\mathbf{k}\pm\mathbf{q})\end{cases}$ (13)

where $\mathbf{q}$ is the momentum of the photon impinging the system and the
$\pm$ sign applies to resonant and antiresonant electron-hole pairs
respectively. We denote the Fourier components of the Coulomb potential with
$v_{\mathbf{q}}(\mathbf{r},\mathbf{r^{\prime}})$, then the two electron
integrals in $\mathbf{A}$ are given by integrating over spin
$\sigma,\varsigma$ and space variables:

$\displaystyle\langle
aj|V|ib\rangle=\sum_{\begin{subarray}{c}\sigma,\sigma^{\prime}\\\
\varsigma,\varsigma^{\prime}\end{subarray}}\int$ $\displaystyle
d\mathbf{r}d\mathbf{r^{\prime}}\psi^{*}_{a,\mathbf{k^{\prime}}}(\mathbf{r},\sigma)\psi^{*}_{j,\mathbf{k+q}}(\mathbf{r^{\prime}},\varsigma^{\prime})v_{\mathbf{q}}(\mathbf{r},\mathbf{r^{\prime}})\times$
$\displaystyle\psi_{i,\mathbf{k^{\prime}+q}}(\mathbf{r},\varsigma)\psi_{b,\mathbf{k}}(\mathbf{r^{\prime}},\sigma^{\prime})\delta_{\sigma,\varsigma}\delta_{\sigma^{\prime},\varsigma^{\prime}}$
(14) $\displaystyle\langle
aj|W_{0}|bi\rangle=\sum_{\begin{subarray}{c}\sigma,\sigma^{\prime}\\\
\varsigma,\varsigma^{\prime}\end{subarray}}\int$ $\displaystyle
d\mathbf{r}d\mathbf{r^{\prime}}\psi^{*}_{a,\mathbf{k^{\prime}}}(\mathbf{r},\sigma)\psi^{*}_{j,\mathbf{k+q}}(\mathbf{r^{\prime}},\varsigma^{\prime})w_{\mathbf{q}}(\mathbf{r},\mathbf{r^{\prime}})\times$
$\displaystyle\psi_{b,\mathbf{k}}(\mathbf{r},\sigma^{\prime})\psi_{i,\mathbf{k^{\prime}+q}}(\mathbf{r^{\prime}},\varsigma)\delta_{\sigma,\sigma^{\prime}}\delta_{\varsigma,\varsigma^{\prime}}$
(15)

and similarly for the $\mathbf{B}$ matrix elements:

$\displaystyle\langle
ab|V|ij\rangle=\sum_{\begin{subarray}{c}\sigma,\sigma^{\prime}\\\
\varsigma,\varsigma^{\prime}\end{subarray}}\int$ $\displaystyle
d\mathbf{r}d\mathbf{r^{\prime}}\psi^{*}_{a,\mathbf{k^{\prime}}}(\mathbf{r},\sigma)\psi^{*}_{b,\mathbf{k}}(\mathbf{r^{\prime}},\sigma^{\prime})v_{\mathbf{q}}(\mathbf{r},\mathbf{r^{\prime}})\times$
$\displaystyle\psi_{i,\mathbf{k^{\prime}+q}}(\mathbf{r},\varsigma)\psi_{j,\mathbf{k-q}}(\mathbf{r^{\prime}},\varsigma^{\prime})\delta_{\sigma,\varsigma}\delta_{\sigma^{\prime},\varsigma^{\prime}}$
(16) $\displaystyle\langle
ab|W_{0}|ji\rangle=\sum_{\begin{subarray}{c}\sigma,\sigma^{\prime}\\\
\varsigma,\varsigma^{\prime}\end{subarray}}\int$ $\displaystyle
d\mathbf{r}d\mathbf{r^{\prime}}\psi^{*}_{a,\mathbf{k^{\prime}}}(\mathbf{r},\sigma)\psi^{*}_{b,\mathbf{k}}(\mathbf{r^{\prime}},\sigma^{\prime})w_{\mathbf{q}}(\mathbf{r},\mathbf{r^{\prime}})\times$
$\displaystyle\psi_{j,\mathbf{k-q}}(\mathbf{r},\varsigma^{\prime})\psi_{i,\mathbf{k^{\prime}+q}}(\mathbf{r^{\prime}},\varsigma)\delta_{\sigma,\varsigma^{\prime}}\delta_{\sigma^{\prime},\varsigma}.$
(17)

For the Hartree components in subsection II.3 and subsection II.3 the
interaction potential comprises of a single mode $\mathbf{q}$ equal to the
difference $\mathbf{k}-\mathbf{k^{\prime}}$. Whereas, the screened Coulomb
interaction in subsection II.3 and subsection II.3 can be decomposed in its
Fourier modes:

$\displaystyle w_{\mathbf{q}}(\mathbf{r,r^{\prime}})=4\pi
e^{2}\sum_{\mathbf{G,G^{\prime}}}\frac{e^{-\mathrm{i}(\mathbf{q+G})\cdot\mathbf{r^{\prime}}}\,\epsilon^{-1}_{\mathbf{q+G},\mathbf{q+G^{\prime}}}\,e^{\mathrm{i}(\mathbf{q+G^{\prime}})\cdot\mathbf{r}}}{|\mathbf{q+G}||\mathbf{q+G^{\prime}}|}.$
(18)

The spin structure of the matrix elements in Eqs. (II.3), (II.3) and (II.3)
spans the singlet subspace, where electron and hole states of each pair have
the same spin. The exchange part in the $\mathbf{A}$ matrix (subsection II.3)
is defined in principle also on the triplet subspace where the incoming
electron and hole have opposite spin orientations Rohlfing and Louie (2000).
If the Hamiltonian is time-reversal invariant (i.e. spin-orbit coupling is not
present) then the singlet and triplet solutions factorise Rödl  _et al._
(2008) and only the former is required to evaluate optical properties and the
correlation energy.

We can now move on to prove that, thanks to the time reversal symmetry of the
wave function, the original EVP can be downfolded into an Hermitian EVP. The
time reversal property of the Schrödinger equation implies that, for a single
particle state, the Bloch functions observe
$\psi_{n,\mathbf{-k}}(\mathbf{r})=\psi^{*}_{n,\mathbf{k}}(\mathbf{r})$ and the
one-electron eigenvalue are also the same. Since the basis functions spanning
the space of resonant and antiresonant components are linearly independent,
one can selectively invert the single particles’ wave-vectors in the
antiresonant component as follows:

$\begin{aligned} \begin{matrix}\\\
\scriptstyle{\mathbf{k}\rightarrow-\bar{\mathbf{k}}}\big{\\{}\\\
\end{matrix}\end{aligned}\begin{pmatrix}\mathbf{A}&\makebox[-6.0pt][c]{$\smash{\overbrace{\phantom{\begin{matrix}\mathbf{B}\end{matrix}}}^{\text{$\mathbf{k}\rightarrow-\bar{\mathbf{k}}$}}}$}\mathbf{B}\\\
\mathbf{B}^{*}&\mathbf{A}^{*}\end{pmatrix}\rightarrow\begin{pmatrix}\mathbf{A}&\mathbf{B}^{\prime}\\\
\mathbf{B}^{\prime}&\mathbf{A}\end{pmatrix}.$ (19)

It is well-known that the $\mathbf{A}$ matrix is Hermitian in this
representation Ring and Schuck (1980), hence the transformation above will
simply map $\mathbf{A}^{*}$ into $\mathbf{A}$. For the exchange part in
subsection II.3, the transformed matrix elements in $\mathbf{B}^{\prime}$ are
given by (with the replacements
$\psi_{j,\mathbf{-\bar{k}-q}}\to\psi^{*}_{j,\mathbf{\bar{k}+q}},\psi^{*}_{b,\mathbf{-\bar{k}}}\to\psi_{b,\mathbf{\bar{k}}}$):

$\displaystyle\int$ $\displaystyle
d\mathbf{r}d\mathbf{r^{\prime}}\psi^{*}_{a,\mathbf{k^{\prime}}}(\mathbf{r})\psi^{*}_{j,\mathbf{k+q}}(\mathbf{r})w_{\mathbf{q}}(\mathbf{r},\mathbf{r^{\prime}})\psi_{b,\mathbf{k}}(\mathbf{r^{\prime}})\psi_{i,\mathbf{k^{\prime}+q}}(\mathbf{r^{\prime}}).$
(20)

Here, we have dropped the bar over the $\mathbf{k}$ index and suppressed the
spin indices for simplicity. The Hartree component turns out to be identical
to the expression in subsection II.3, hence requiring no extra computational
step.

Furthermore, the transformation applied to $\mathbf{B}^{*}$ gives exactly the
same matrix elements in Eq. (missing)20 as we now show. Starting with the
matrix elements in subsection II.3 and applying the transformation
$\mathbf{k}^{\prime}\to-\mathbf{\bar{k}}^{\prime}$ to the antiresonant
electron-hole pair $(a,i)$, one obtains:

$\displaystyle\langle\bar{a}b|W_{0}|j\bar{i}\rangle^{*}=\int
d\mathbf{r}d\mathbf{r^{\prime}}$
$\displaystyle\psi^{*}_{j,\mathbf{k+q}}(\mathbf{r})\psi^{*}_{i,\mathbf{-\bar{k}^{\prime}-q}}(\mathbf{r^{\prime}})$
$\displaystyle\times
w^{*}_{\mathbf{q}}(\mathbf{r,r^{\prime}})\psi_{a,-\mathbf{\bar{k}^{\prime}}}(\mathbf{r})\psi_{b,\mathbf{k}}(\mathbf{r^{\prime}}).$
(21)

In the static approximation to the screening potential
$w_{\mathbf{q}}(\mathbf{r,r^{\prime}})$, the dielectric function at the RPA
level is an even function of $|\mathbf{q}|$ and has zero imaginary part, hence
it is possible to replace $w^{*}_{\mathbf{q}}(\mathbf{r,r^{\prime}})$ with
$w_{-\mathbf{q}}(\mathbf{r,r^{\prime}})$ in the equation above. Finally, by
swapping the integration variables, it is easy to see that the exchange
component in $\mathbf{B}^{*^{\prime}}$ reduces to Eq. (missing)20. The
reasoning above can be straightforwardly repeated for the Hartree component in
$\mathbf{B}^{*}$ to give $\langle aj|V|ib\rangle$. This result completes the
proof of Eq. (missing)19. In the following we will drop the prime when
referring to the $\mathbf{B}^{\prime}$ matrix.

We can now proceed to reduce the dimensionality of the EVP (10). Adding and
subtracting the individual equations contained in Eq. (missing)10 and then
solving one of the two equations, let’s say with respect to $(\mathbf{X-Y})$,
one obtains the squared problem Jørgensen and Simons (1981):

$\displaystyle(\mathbf{A}+\mathbf{B})(\mathbf{A}-\mathbf{B})(\mathbf{X}+\mathbf{Y})=\mathbf{\Omega}^{2}(\mathbf{X}+\mathbf{Y}).$

The expression above can be converted into a conventional EVP if the matrices
are positive definite. Then the matrix $(\mathbf{A}-\mathbf{B})^{\frac{1}{2}}$
is unique and the eigenvalue problem finally reads:

$\displaystyle(\mathbf{A}-\mathbf{B})^{\frac{1}{2}}(\mathbf{A}+\mathbf{B})$
$\displaystyle(\mathbf{A}-\mathbf{B})^{\frac{1}{2}}(\mathbf{A}-\mathbf{B})^{\frac{1}{2}}(\mathbf{X}+\mathbf{Y})=$
$\displaystyle=\mathbf{\Omega}^{2}(\mathbf{A}-\mathbf{B})^{\frac{1}{2}}(\mathbf{X}+\mathbf{Y}).$
(22)

The difference $\mathbf{P}_{\lambda}(\omega)-\mathbf{P}_{0}(\omega)$ in
LABEL:eq:EcP can be integrated over $\omega$ thanks to the spectral
decomposition in Eq. (missing)9 and in Eq. (missing)42. The correlation energy
can then be expressed using the correlation part of the two-particle density
matrix $\mathcal{P}_{\lambda}$:

$\displaystyle
E_{c}=\frac{1}{2}\int_{0}^{1}d\lambda\operatorname{tr}\\{\mathbf{V}\mathcal{P}_{\lambda}\\}$
(23)

with
$\mathcal{P}_{\lambda}=(\mathbf{X}+\mathbf{Y})_{\lambda}(\mathbf{X}+\mathbf{Y})_{\lambda}^{*}-\mathbf{1}$
and $\mathbf{V}$ the Coulomb matrix in subsection II.3 Furche (2001a); Ángyán
_et al._ (2011).

##  III Computational details

The HEG is in principle specified by the average electron density
$n=\frac{N_{e}}{V}$ (with $V$ the unit cell volume and $N_{e}$ number of
electrons in it) and the number of electrons in a given spin configuration
(either $N_{e,\uparrow}$ or $N_{e,\downarrow}$). These quantities are easily
translated into the Wigner-Seitz radius, $r_{s}$ given by $\frac{4}{3}\pi
r_{s}^{3}=\frac{V}{N_{e}}$ and the spin polarisation
$\zeta=\frac{|N_{e,\uparrow}-N_{e,\downarrow}|}{N_{e}}$. The discretisation of
the Bloch wave-vectors introduced in subsection II.3, however, requires also a
specific choice for the $\mathbf{k}$-point sampling and for the symmetry of
the simulation cell.

Calculations were performed on a $\Gamma-$point centred, uniform
$\mathbf{k}$-point mesh of dimensions $N_{k}\times N_{k}\times N_{k}$.
Furthermore, we applied simple cubic unit cells in all our simulations. As
already emphasized before results are independent of the choice of the unit
cell, once $\mathbf{k}$-point convergence is reached. We confirmed for
instance that a face-centred cubic cell and a simple cubic cell result in the
same RPA energy. Simple cubic cells are, however, often used for quantum
Monte-Carlo simulations. For instance, a $3\times 3\times 3$
$\mathbf{k}$-point grid with 2 electrons in the unit cell corresponds exactly
to a $3\times 3\times 3$ super cell with 54 electrons which was often used in
QMC simulations, because the resultant electronic configuration has only fully
occupied or entirely unoccupied one-electron orbitals Lin  _et al._ , 2001.

The VASP code Kresse and Furthmüller (1996); Kresse (1996) has been used for
all calculations, which required a three-stage computational procedure. The
one-body Green’s function variation was disregarded along the coupling
constant integral, as discussed above. The solidity of this approximation will
be checked by comparing the resulting correlation energy against the QMC
benchmark. For a given $(r_{s},\zeta)$ the procedure to evaluate the
correlation energy starts with a self-consistent calculation at the Kohn-Sham
level. The plane-wave basis set is specified by the energy cutoff $E_{\rm
cut}$ and includes plane waves with kinetic energy smaller than the given
cutoff value. This step is followed by a self-consistent evaluation of the
quasi-particle energies at the $QP-GW_{0}$ level Shishkin and Kresse (2007,
2006). The screening potential $W_{0}$ is evaluated at the RPA level with
orbitals and eigenvalues from the previous DFT calculation. To ensure
consistency of the results with previous calculations Grüneis  _et al._ (2009)
for the evaluation of the correlation energy we set the cutoff for the
response function’s basis set (encutgw) and densities in the two electron
integrals to the same value $E_{\rm cut}$. The frequency sampling has been
carried out on a linear grid, specified by its maximum value, set equal to 1.5
$E_{\rm cut}$ for all $(r_{s},\zeta)$. The grid point density has been chosen
as 0.1 for $r_{s}$=0.5, 0.8 and 0.2 for $r_{s}$=1.0 (collectively referred as
the high density region), then progressively increased to 0.6 and 1.0 in the
intermediate ($r_{s}$=2.0, 3.0) and low densities ($r_{s}$=4.0, 5.0),
respectively. The evaluation of quasi-particle energies, necessary to specify
the dressed propagator, has been performed iteratively. We set the number of
iterations to three; the resulting quasi-particle energies have a mean
accuracy of $\approx$ 1 mRy. These are the only quantities being updated
(plane waves being exact eigenfunctions for the HEG). Thus, the one-body
Green’s function retains a simple one electron form.

The BSE was recast as in subsection II.3 and has been implemented in VASP; all
virtual unoccupied orbitals spanned by the plane-wave basis were included. The
diagonalisation of subsection II.3 is required for each value of the coupling
constant, whereas the Coulomb kernel (constant along the AC-path) is computed
and stored once and for all at the beginning of the BSE step. The integral in
Eq. (missing)23 is evaluated numerically with a Gauss-Legendre quadrature
scheme. As few as 2 points are sufficient to produce results converged within
the order of meV.

Fig. 2:  Left panel: convergence for the BSE correlation energy as a function
of the total number of plane waves; the energy zero has been set equal to the
linearly extrapolated complete basis energy limit for each value of
$(r_{s},\zeta)$. Right panel: total energy convergence as a function of the
number of points included in the reciprocal space sampling. Here the shown
value is calculated by taking the difference between the total energy for any
given $\mathbf{k}$-mesh and the value at the preceding sampling.

##  IV Results and Discussion

Fig. 3: Direct (left) and exchange (right) diagrams contributing in second
order to the correlation energy.

We start by considering the correlation energy’s convergence properties for a
given $(r_{s},\zeta)$ against the total number of plane waves. The impact of
the basis set incompleteness is shown in Figure 2 (left panel): for the set of
densities considered, a linear convergence behaviour with the inverse of the
number of bands is observed Harl and Kresse (2008). The complete basis set
limit has been estimated for both the paramagnetic ($\zeta=0$) and
ferromagnetic ($\zeta=1$) case with a linear regression. For the ferromagnetic
system a particularly fast convergence is observed as a function of the number
of plane waves. To the best of our knowledge, an analytical wave-vector
analysis is not available for arbitrary spin polarisations. However, in Ref.
Gori-Giorgi and Perdew, 2002 the pair correlation function limit for large
wave-vectors is said to be determined by the kinks in the many electron wave
function. For fully spin polarized systems these are completely absent, this
then rationalises the fast convergence observed. We will comment on this point
in more detail below.

The assessment of the total energy convergence with respect to the
$\mathbf{k}$-point grid density is reported in the right panel of Figure 2.
Here we report the total energy for any given $\mathbf{k}$-point mesh minus
the value at the preceding $\mathbf{k}$-point sampling divided by the energy
value obtained with the most accurate $\mathbf{k}$-point grid for a given
level of theory. For each $\mathbf{k}$-point grid, the total energy has been
computed on the $GW_{0}$ reference state and includes the one-electron kinetic
energy, the exact exchange energy, a correction term related to the occurrence
of fractional occupancies (Ref. Harl  _et al._ , 2010), as well as the
correlation contribution. We include in the figure also the RPA energy
convergence that provides an upper bound on the relative error. The
corresponding error on the total energy is less than 1% for the ferromagnetic
system and less than 0.5% for the paramagnetic ground state. On the other
hand, convergence of the correlation energy alone proves more difficult. To
estimate the correlation energy convergence with respect to both the basis set
and the $\mathbf{k}$-point sampling we employ an extrapolation scheme adapted
from Ref. Klimeš  _et al._ , 2014:

$\displaystyle E^{\infty}_{c}(N_{k}\rightarrow\infty)\approxeq
E_{c}^{r}(N_{k}\rightarrow\infty)-E_{c}^{r}(N_{k})+E^{\infty}_{c}(N_{k}),$
(24)

where the superscript of the correlation energy $E_{c}$ indicates the
dependence on the basis set. The symbol r indicates that the calculations were
performed using a small number of bands (32 bands here), whereas $\infty$
indicates that the results have been extrapolated to the infinite basis set
limit. In practice, we found that it is sufficient to determine the basis set
correction ($-E_{c}^{r}(N_{k})+E^{\infty}_{c}(N_{k})$) using $N_{k}=3$, i.e.
$3\times 3\times 3$ $\mathbf{k}$-points.

To validate this somewhat involved approach, we first calculated the
correlation energy starting from the Kohn-Sham ground state, i.e. RPA @ LDA.
This should reproduce the previously published analytic results of Ref. Singwi
_et al._ (1968) for the paramagnetic system and of Ref. Perdew and Zunger
(1981) for the spin-polarised case. We tested the agreement for the case
$(r_{s}=1,\zeta=0)$ where we found it necessary to include up to $N_{k}=18$
points to reproduce the analytical results within 4 meV. Fewer
$\mathbf{k}$-points ($N_{k}=16$) are necessary for the spin-polarised case to
reproduce the data with a similar accuracy. We stress that the RPA results
reported below have been calculated on the $GW_{0}$ reference and hence can
not match the historical RPA correlation energies, owing to the
renormalisation of the propagator lines in the RPA response function.

Given the current computational limitations, we adapt the extrapolation scheme
above to the case of post-RPA calculations with the replacement:

$\displaystyle E_{c}^{r}(N_{k}\rightarrow\infty)\approxeq
E_{c}^{r}(N_{k}=8)+\Delta_{N_{k}},$ (25)

where the correlation energy is computed at the considered level of theory
whereas the $\Delta_{N_{k}}$ correction corrects for the $\mathbf{k}$-point
incompleteness error. This value is difficult to calculate accurately for the
BSE. In previous work Marsman  _et al._ (2009), it was found that the second
order exchange reduces the direct correlation energies by about 1/3. The
situation is similar here, with exchange contributions reducing the
correlation energy by roughly 30%. In line with this, the $\mathbf{k}$-point
errors are generally 1/3 smaller for BSE than for RPA (compare Fig. 2, but
note that the left panel presents relative errors). Overall, it therefore
seems sensible to determine the $\mathbf{k}$-point error using the
computationally efficient RPA, but to reduce the RPA $\mathbf{k}$-point
correction by a factor $\frac{2}{3}$ if exchange is included.

Fig. 4: Values of correlation energy for the full BSE (solid line) and for the
RPAsX approximation (dash-dotted line) as a function of the Wigner-Seitz
radius for the paramagnetic (left) and ferromagnetic (right) HEG. Results are
compared against the QMC simulations in Ref. Ceperley and Adler, 1980 (round
symbols) and in Ref. Ortiz and Ballone, 1994 (crosses). The random phase
approximation on the same reference state (broken line) is also included.

The correlation energies computed for a set of $(r_{s},\zeta)$ are shown in
Figure 4 by a solid line and compared against the Quantum Monte Carlo
estimates by Ceperly and Adler (CA in the following) Ceperley and Adler (1980)
(indicated by circles in the figure) and by Ortiz and Ballone (OB in the
following) Ortiz and Ballone (1994) (‘$\times$’ symbols in the figure). Finite
size effects impact also these benchmarks Fraser  _et al._ (1996) and
different extrapolation schemes, which do not necessarily fulfil the
variational principle, have been employed. Since we did not want to judge
which calculation is more accurate, we compare to both results. The RPA@GW0
estimate of the correlation energy (dotted line in Figure 4), is obtained as
usually by including the bubble diagrams only, i.e. excluding the ladder
diagrams (A2) and (B2) in Figure 1. The newly proposed RPA with screened
exchange (dash-dotted line), where only bubbles and diagrams (B2) are included
is also shown (see below). Starting our analysis from the full BSE
calculations (solid line in Figure 4), in the paramagnetic case the CA results
are reproduced with a mismatch of 0.06 eV for the highest density. The error
progressively increases in the high density region and becomes significant
with a deviation of $\approx$20% for $r_{s}\geq$4.0 a.u.. This pattern is also
observed for the ferromagnetic system, with the absolute deviation, however,
decreasing to $\approx$10% at low densities.

The general trend of the present BSE results can be rationalised as follows:
in the low density region the BSE makes a sizeable error, which is of unknown
origin, but could be related to the lack of particle-particle (or hole-hole)
ladder diagrams to describe short range interactions Yasuhara (1972). Indeed,
for the paramagnetic system (where short range interactions are most
important) we observe that the BSE and RPA results are close in energy and
both deviate from the QMC estimate. In the high density region the BSE
reproduces well the QMC estimates for both values of spin polarisation.

Short range interactions are small in the ferromagnetic case owing to the
Pauli exclusion principle and are incorrectly included by the RPA. This self-
correlation error degrades the agreement between the QMC estimates and the RPA
Paier  _et al._ (2012) (see left panel in Figure 4). The self-correlation
error is cancelled in leading order by the BSE, as we will discuss now. In the
low density region the large wave-vector contributions to the correlation
energy become predominant Giuliani and Vignale (2005). At large wave-vectors,
the screened interaction $W_{0}$ equals the bare interaction $V$, since
screening is weak. The leading contribution to the BSE becomes then equal to
the second order contribution in Møller-Plesset perturbation theory (MP2)
given in Ref. Jansen  _et al._ , 2010. This can be shown by determining the
leading term in the correlation part of the two body density matrix. For the
Hartree term and the RPA the equivalence has also been derived in Ref. Marsman
_et al._ , 2009. In second order, the correlation energy is then simply given
by the standard textbook equation:

$\displaystyle E^{(2)}_{c}=\frac{1}{2}\sum_{ijab}\frac{\langle
ij|ab\rangle}{\Delta E_{ia}+\Delta E_{jb}}\big{\\{}\langle
ab|ij\rangle-\langle ab|ji\rangle\big{\\}}.$ (26)

From this, two important observations follow.

(i) The Hartree contribution $\langle ab|ij\rangle$ and the exchange term
$-\langle ab|ji\rangle$ cancel each other for $j=i$ or for $a=b$ for any spin
orientation. This self-correlation error cancellation is obviously observed
for all values of $\zeta$ since only same-spin electrons are affected by the
self-correlation error. In other words, MP2 is self-correlation free and this
property is shared in second order by the BSE. This property is also observed
by the SOSEX approximation Grüneis  _et al._ (2009). However, for SOSEX the
corresponding second order diagram is introduced somewhat ad hoc by anti-
symmetrizing the direct-RPA coupled cluster amplitudes. For the BSE, the
diagrams are naturally included at the level of the two particle propagator
and two body density matrix. The BSE should therefore improve upon the SOSEX,
in particular, if static correlation effects are importantHenderson and
Scuseria (2010).

(ii) The expression for the correlation energy in second order also allows us
to understand why the correlation energy converges so quickly with respect to
the number of orbitals for the ferromagnetic case. For a ferromagnetic system,
only one spin component is present. For large wave-vectors, the dominant terms
to the correlation energy stem from $i=j$, and the direct term is exactly
cancelled by the exchange term. This explains the fast convergence of the
correlation energy with respect to the plane wave cutoff for the ferromagnetic
case. For a paramagnetic system, where both spin components are present, it is
possible to carry out the summation over the spin degrees of freedom
implicitly present in the spin-orbital representation above. This leads to a
factor $2^{l}$, with $l$ being the number of fermionic loops in the diagrams
in Figure 3. Hence, if Eq. (missing)26 is evaluated for spatial orbitals
Paldus and Čížek (1975), a factor four enters for the direct term, whereas a
factor two is present for the exchange term. Thus the two terms do not cancel
each other in the long wave-length limit, more precisely the direct term is
reduced by a factor two by the exchange term. This causes the slow convergence
for the non-spin polarized system.

Given the relevance of the HEG as a model for condensed matter, it is not
surprising that other research groups have recently evaluated its properties
as well. For instance, the use of the renormalised ALDA kernel with the
inclusion of correlation effects seems to reproduce exactly the HEG
correlation energy in the limit of very low densities ($r_{s}>$ 10) Patrick
and Thygesen (2015). However larger deviations are observed in the ”metallic”
density region ($r_{s}\in[0.5,5.0]$). Also, a recent study by the de
Gironcoli’s group Colonna  _et al._ (2014) has addressed the role of
correlations in the HEG by employing exact exchange (EXX) TDDFT. A feature of
this method is to include first order changes of the one-particle propagator
due to exact exchange, as well as particle-hole exchange interactions in the
two-point response function (due to adiabatically switching on the exact
exchange between electrons). This response function is then used to form the
exchange kernel as proposed by Görling Görling (1998). An issue that this
approach faces is the occurrence of imaginary frequencies when the response
function is diagonalised (see Eq. (10) in Colonna  _et al._ (2014)). Although
it is possible to circumvent this problem by including the ladder diagrams
(A2) and (B2) only to first order in the response function, this also spoils
the results to some extent Colonna  _et al._ (2014). Other difficulties in the
EXX-TDDFT method stem from the inversion of the non-interacting response
function to evaluate the exchange kernel (see Eq. (7) in Bleiziffer  _et al._
(2012)) and on the critical dependence on the basis set size Hirata  _et al._
(2001).

Also in the present BSE approach instabilities at lower electron densities can
be present (as in our case for $r_{s}\geq$5). These are witnessed as imaginary
frequencies $\Omega$ appearing in subsection II.3. The instabilities can cause
the matrices $\mathbf{A+B}$ or $\mathbf{A-B}$ to have negative eigenvalues. If
both matrices have negative eigenvalues, it has been argued that the ground
state is unstable with respect to particle-hole excitations Jørgensen and
Simons (1981). In the present case, however, we find that only the
$\mathbf{A+B}$ matrix has negative eigenvalues, for both the $GW_{0}$ and LDA
reference states. This finding is consistent with a recent theoretical
investigation Uimonen  _et al._ (2015) which has confirmed the presence of
negative-frequency modes in the BSE polarisation propagator.

In our case the instabilities originate from the particle-hole ladder
diagrams, as also witnessed in the case of the already mentioned TDDFT
calculations. To resolve this problem, we suggest to disregard all particle-
hole ladder diagrams in the $\mathbf{A}$ matrix— specifically those shown in
(A2) of Figure 1 —when solving the BSE equation. In the following, we will
refer to this approximation as RPA with screened exchange (RPAsX).

We believe this choice to be sensible for the following reasons. First, the
$\mathbf{A}$ matrix does not contribute to the correlation energy in the
perturbation expansion to second order Ángyán  _et al._ (2011); Heßelmann and
Görling (2011). Thus removing the diagrams (A2) will only change the
correlation energy in third order; in second order, for instance, we still
recover the exact MP2 energy (if the screened exchange kernel is replaced by
the bare one). The approximation is also closely related to the AC-SOSEX
method Grüneis  _et al._ (2009) but improves upon it. In the AC-SOSEX the
conventional RPA polarisation propagator is evaluated. As opposed to the
direct RPA, the polarisation propagator is then contracted against the
$\mathbf{B}$ matrix (containing the sum of the direct Coulomb interaction and
the bare exchange interaction, diagrams (B1) and (B2) in Eq. (missing)23)
Jansen  _et al._ (2010). In RPAsX, the propagator now includes, as it should,
the exchange term in the $\mathbf{B}$ matrix (B1 as well as B2). Furthermore,
the exchange terms (B2) now include a screened exchange interaction instead of
a bare interaction, which effectively mimics higher order ladder diagrams. It
is clear from Figure 4 that this approximation is particularly successful for
the HEG, as it completely prevents the occurrence of unstable solutions and
yields excellent agreement with the QMC data for both spin polarisation values
as shown in Figure 4 by a dash-dotted line. We observe a slight upward shift
of the correlation energy compared to the full BSE approach, which leads to
undercorrelation in comparison with CA by 0.03, 0.02 eV for $r_{s}$=0.5 and
0.8, 1.0 a.u. respectively. For $r_{s}\geq$2 the RPAsX estimate of the
correlation energy remains within 6% of the CA estimate for $\zeta=0$ (it
overcorrelates to a greater extent in comparison with OB). For the spin-
polarised case RPAsX lies within the range of values spanned by the different
QMC simulations in the high density region. As the density is lowered there is
a progressive increase of the correlation energy in RPAsX, which still
compares very favourably with CA values, reaching a maximum mismatch of 20 meV
for the lowest density considered. RPAsX correlation energy is mostly parallel
to the BSE results, but moves closer to the exact results. We admit that the
very good agreement with the QMC results must be fortuitous to some extend,
since the RPAsX still neglects particle-particle and hole-hole ladder
diagrams. But obviously the neglect of these diagrams cancels the also
neglected particle-hole ladder diagram at least for the HEG.

Finally, we consider a partially spin polarised system. The non-interacting
response function is diagonal in the spin basis and this property is conserved
by the BSE kernel, because opposite spin components interact only through the
direct diagrams (A1) and (B1) in Figure 1. As it was mentioned, increasing the
spin polarisation decreases the number of electrons able to interact at short
range. To assess the change introduced in the correlation energy $E_{c}$ by
modifying the fraction of short-ranged vs long-range interactions, the spin
enhancement function $\Upsilon(\zeta;r_{s})$ is commonly introduced:

$\displaystyle\Upsilon(\zeta;r_{s})=\frac{E_{c}(\zeta,r_{s})-E_{c}(0,r_{s})}{E_{c}(1,r_{s})-E_{c}(0,r_{s})}.$
(27)

Fig. 5: The spin enhancement function $\Upsilon$ as a function of the
fractional polarisation at density $r_{s}$=2. Calculations were performed on a
6$\times$6$\times$6 $\mathbf{k}$-point grid with complete basis set
extrapolation.

In Figure 5 we report the spin enhancement function computed at $r_{s}$=2
versus the spin polarisation $\zeta$ for various levels of theory. Results are
compared against the Perdew-Wang interpolation formula Perdew and Wang (1992)
with parameters estimated by Ortiz and Ballone Ortiz and Ballone (1994)
(dashed line). RPA results are evaluated on the DFT reference state
(represented in the figure by the ‘$\times$’ symbols): there is a noticeable
departure from the QMC estimate, with the RPA overestimating the short range
correlations, as expected. Changing reference state (RPA on $GW_{0}$ indicated
by ‘+’ symbols in the figure) actually increases the mismatch with the QMC
estimate. The BSE results, shown by squares in the figure, reproduce
exceedingly well the benchmark calculations for all values of $\zeta$
considered. Given the high computational cost required by these calculations
we do not report RPAsX for fractional spin polarisations. However, since the
BSE makes a more sizeable error than RPAsX (see Figure 4) and yet it
reproduces $\Upsilon$ very accurately, it is reasonable to assume that RPAsX
will also be very accurate for fractional values of $\zeta$.

##  V Conclusions

In this study we have evaluated the correlation energy for the homogeneous
electron gas by calculating the fluctuation contributions from the Bethe-
Salpeter equation. The implementation of the BSE has taken advantage of the
time inversion symmetry, converting it into an Hermitian quadratic eigenvalue
problem. Both the resonant and the antiresonant contributions to the particle-
hole interaction have been included in the BSE (going beyond the so-called
Tamm-Dancoff approximation).

The BSE kernel has been set up consistently with Hedin’s $GW_{0}$
approximation. This level of theory has been employed also to describe the
system’s electronic structure, taken as a starting point for the ensuing BSE
calculation. The occurrence of unstable solutions using the full Bethe-
Salpeter kernel has prompted us to seek an approximation scheme, still able to
compare favourably with QMC references. The approximation proposed, called
RPAsX, is by construction consistent with MP2 and reproduces the QMC
correlation energy very well for values of $r_{s}\geq$ 1.0 a.u.. We certainly
plan to test this approximation for a wider class of systems and materials.
There are however several obstacles that need to be solved before this
approach can be applied routinely. The most important one is that for the HEG
changes related to changes of the one-particle Green’s function and one
particle orbitals can be neglected (compared Eq. (missing)4). This is not the
case for real materials, where self-consistency of the orbitals van
Schilfgaarde  _et al._ (2006) or changes of the one particle Green’s function
when going to the interacting case are very relevant Klimeš  _et al._ (2015).
In order to apply the RPAsX to real systems it will be necessary to relax the
approximations made here. With this provision, we are confident that this
approach can be extended to realistic solid state and molecular systems.

## Acknowledgments

Financial support from the Austrian Science Found (FWF) for the SFB ViCoM
project no. F41 is kindly acknowledged.

##  Appendix A Derivation of Eq. 5

Fig. 6: Graphical representation for the two-body Green’s function. We
represent the time ordering $t_{1}=t_{2}$ which corresponds to having two
electron-hole pairs interacting instantaneously

In this appendix we derive the expression for the correlation energy in terms
of the polarisation propagator. The correlation energy involves the two-body
density matrix and, as shown in the main text in Eq. (missing)4, is given by

$\displaystyle\left(\frac{dE_{\lambda}}{d\lambda}\right)_{\rm{2B}}\equiv\dot{E}_{\lambda}^{c}=\frac{1}{2}\int
dt_{1}\int\prod_{i=1}^{4}d\mathbf{x}_{i}\delta(\mathbf{x}_{3}-\mathbf{x}_{1})\delta(\mathbf{x}_{4}-\mathbf{x}_{2})v(\mathbf{x}_{1},\mathbf{x}_{2})\langle\psi^{\dagger}(\mathbf{x}_{4},t^{+}_{1})\psi^{\dagger}(\mathbf{x}_{3},t_{1})\psi(\mathbf{x}_{1},t^{+}_{1})\psi(\mathbf{x}_{2},t_{1})\rangle_{\lambda}.$

The four field operators in this expression suggest that one can make a
connection to Green’s function theory, specifically to the two particle
Green’s function. We also note that the expression above only involves equal
time limits (all involved time points are equal). To make a connection to
Green’s function theory, we first need to introduce a time dependence in the
field operators by defining the time evolution operator in the interaction
picture (as it is commonly done in many body perturbation theory):
$\hat{\psi}(1)=e^{i\mathcal{H}_{0}t_{1}}\hat{\psi}(\mathbf{x}_{1})e^{-i\mathcal{H}_{0}t_{1}}$
and
$\hat{\psi}^{\dagger}(2)=e^{i\mathcal{H}_{0}t_{2}}\hat{\psi}^{\dagger}(\mathbf{x}_{2})e^{-i\mathcal{H}_{0}t_{2}}$.
The $n$-body Green’s operator is as usually defined as:

$\displaystyle\mathcal{G}_{n}$
$\displaystyle(1,2,...,n;1^{\prime},2^{\prime},...,n^{\prime})=$
$\displaystyle\frac{1}{i^{n}}\mathcal{T}\\{\hat{\psi}(1)\hat{\psi}(2)\dots\hat{\psi}(n)\hat{\psi}^{\dagger}(n^{\prime})\dots\hat{\psi}^{\dagger}(2^{\prime})\hat{\psi}^{\dagger}(1^{\prime})\\}$
(29)

where $\mathcal{T}$ is the time ordering operator. It is worth pointing out
that the operators $\mathcal{G}_{n}$ do not depend on the value of the
coupling constant $\lambda$. For $n=1$ the corresponding (greater) Green’s
function is obtained by evaluating the expectation value of the Green’s
operator for the ground-state wave function at coupling $\lambda$,
$\Psi_{\lambda}$. Since we are only interested in equal time limits here, we
can restrict the second time to fulfil $t^{\prime}=t+0^{+}=t^{+}$:

$\displaystyle\langle\mathcal{G}_{1}(\mathbf{x}_{1},t^{+};\mathbf{x}_{2},t)\rangle_{\lambda}$
$\displaystyle\equiv$ $\displaystyle
G_{1,\lambda}(\mathbf{x}_{1},t^{+};\mathbf{x}_{2},t)$
$\displaystyle=i\langle\psi^{\dagger}(\mathbf{x}_{2},t)\psi(\mathbf{x}_{1},t^{+})\rangle_{\lambda}$
(30)

At equal time limits, the Green’s function can be also related to the density
matrix $\langle\hat{n}\rangle_{\lambda}\equiv
n_{\lambda}(\mathbf{x},\mathbf{y})=-iG_{1,\lambda}(\mathbf{x},t;\mathbf{y},t)$.

For the two particle Green’s function, it is standard texbook material Gross
_et al._ (1991) (commutator relations or Wick’s theorem) to show that it
satisfies:

$\displaystyle G_{2,\lambda}$
$\displaystyle(\mathbf{x}_{1},t^{+}_{1},\mathbf{x}_{2},t_{2};\mathbf{x}_{3},t_{1},\mathbf{x}_{4},t^{+}_{2})$
$\displaystyle=$
$\displaystyle-\langle\psi^{\dagger}(\mathbf{x}_{3},t_{1})\psi(\mathbf{x}_{1},t^{+}_{1})\rangle_{\lambda}\langle\psi^{\dagger}(\mathbf{x}_{4},t^{+}_{2})\psi(\mathbf{x}_{2},t_{2})\rangle_{\lambda}$
$\displaystyle+\langle\psi^{\dagger}(\mathbf{x}_{4},t^{+}_{2})\psi(\mathbf{x}_{1},t^{+}_{1})\rangle_{\lambda}\langle\psi^{\dagger}(\mathbf{x}_{3},t_{1})\psi(\mathbf{x}_{2},t_{2})\rangle_{\lambda}$
$\displaystyle+\langle\psi^{\dagger}(\mathbf{x}_{4},t^{+}_{2})\psi^{\dagger}(\mathbf{x}_{3},t_{1})\psi(\mathbf{x}_{1},t^{+}_{1})\psi(\mathbf{x}_{2},t_{2})\rangle_{\lambda},$
(31)

where the last term on the r.h.s. is the correlation part of the two-body
Green’s function that is also present in the AC fluctuation-dissipation
theorem (ACFDT) correlation energy LABEL:eq:dEdl2B. The labelling convention
for four-point quantities is given in Figure 6 and it is consistent with Ref.
Starke and Kresse, 2012 up to a time inversion. Since the correlation energy
involves only equal time limits, we again set $t_{1}=t_{2}=t$. This allows us
to rewrite Eq. (missing)31 in terms of the ordinary density operators
$\hat{n}$ and their expectation values:

$\displaystyle G_{2,\lambda}$
$\displaystyle(\mathbf{x}_{1},t^{+},\mathbf{x}_{2},t;\mathbf{x}_{3},t,\mathbf{x}_{4},t^{+})$
$\displaystyle=$
$\displaystyle-\langle\psi^{\dagger}(\mathbf{x}_{3},t)\psi(\mathbf{x}_{1},t^{+})\rangle_{\lambda}\langle\psi^{\dagger}(\mathbf{x}_{4},t^{+})\psi(\mathbf{x}_{2},t)\rangle_{\lambda}$
$\displaystyle+n_{\lambda}(\mathbf{x}_{4},\mathbf{x}_{1})n_{\lambda}(\mathbf{x}_{3},\mathbf{x}_{2})-\langle\hat{n}(\mathbf{x}_{4},\mathbf{x}_{1})\hat{n}(\mathbf{x}_{3},\mathbf{x}_{2})\rangle_{\lambda}$
$\displaystyle=$
$\displaystyle-\langle\psi^{\dagger}(\mathbf{x}_{3},t)\psi(\mathbf{x}_{1},t^{+})\rangle_{\lambda}\langle\psi^{\dagger}(\mathbf{x}_{4},t^{+})\psi(\mathbf{x}_{2},t)\rangle_{\lambda}$
$\displaystyle-\langle\delta\hat{n}(\mathbf{x}_{4},\mathbf{x}_{1})\delta\hat{n}(\mathbf{x}_{3},\mathbf{x}_{2})\rangle_{\lambda},$
(32)

where in going to the last line, we have introduced the density fluctuation
operator
$\langle\delta\hat{n}\rangle_{\lambda}\equiv\langle\hat{n}\rangle_{\lambda}-n_{\lambda}$.
The last term in Eq. (missing)32 can be identified with the polarisation
propagator (also called density-fluctuation density-fluctuation response
function), i.e.
$P_{\lambda}\equiv\langle\delta\hat{n}\delta\hat{n}\rangle_{\lambda}$. One can
then rearrange Eq. (missing)32 into the usual relation connecting the two
particle propagator and the polarisation propagator Stefanucci and van Leeuwen
(2013):

$\displaystyle P_{\lambda}$
$\displaystyle(\mathbf{x}_{1},t^{+},\mathbf{x}_{2},t;\mathbf{x}_{3},t,\mathbf{x}_{4},t^{+})=$
$\displaystyle-
G_{2,\lambda}(\mathbf{x}_{1},t^{+},\mathbf{x}_{2},t;\mathbf{x}_{3},t,\mathbf{x}_{4},t^{+})$
$\displaystyle+G_{1,\lambda}(\mathbf{x}_{1},t^{+};\mathbf{x}_{3},t)G_{1,\lambda}(\mathbf{x}_{2},t;\mathbf{x}_{4},t^{+}).$
(33)

Using Eq. (missing)31 and Eq. (missing)33, one can rewrite the correlation
part of $G_{2,\lambda}$ in the ACFDT equation using the polarisation
propagator and obtain the following compact equation for the correlation
energy:

$\displaystyle\dot{E}^{\rm{c}}_{\lambda}=\frac{1}{2}\int\prod_{i=1}^{4}d\mathbf{x}_{i}\,V$
$\displaystyle\left[G_{1,\lambda}(\mathbf{x}_{1},t^{+};\mathbf{x}_{4},t^{+})G_{1,\lambda}(\mathbf{x}_{2},t;\mathbf{x}_{3},t)\right.$
$\displaystyle\left.-P_{\lambda}(\mathbf{x}_{1},t^{+},\mathbf{x}_{2},t;\mathbf{x}_{3},t,\mathbf{x}_{4},t^{+})\right]$
(34)

From Eq. (missing)34 it is easy to separate out the exchange contribution:

$\begin{split}\dot{E}^{\rm{x}}_{\lambda}&=-\frac{1}{2}\int\prod_{i=1}^{4}d\mathbf{x}_{i}\delta(\mathbf{x}_{3}-\mathbf{x}_{1})\delta(\mathbf{x}_{4}-\mathbf{x}_{2})\\\
&\qquad\qquad
v(\mathbf{x}_{1},\mathbf{x}_{2})n_{\lambda}(\mathbf{x}_{2},\mathbf{x}_{3})n_{\lambda}(\mathbf{x}_{1},\mathbf{x}_{4})\\\
&=-\frac{1}{2}\int
d\mathbf{x}_{1}d\mathbf{x}_{2}v(\mathbf{x}_{1},\mathbf{x}_{2})n_{\lambda}(\mathbf{x}_{2},\mathbf{x}_{1})n_{\lambda}(\mathbf{x}_{1},\mathbf{x}_{2}).\end{split}$

The remaining contribution in Eq. (missing)34 can be replaced into Eq.
(missing)3 to form the correlation energy assuming that the one body
contributions (contained in the exchange part) remain constant along the AC
path. LABEL:eq:EcP is then obtained by Fourier transforming with respect to
the infinitesimal time difference.

##  Appendix B Spectral representation for the polarisation propagator

The starting point to derive the spectral representation for the interacting
polarisation propagator $P_{\lambda}(\omega)$ is to invert Eq. (missing)8 in
the particle-hole basis:

$\displaystyle\mathbf{P}_{\lambda}(\omega)=\left[\mathbf{P}_{0}^{-1}(\omega)-\mathbf{I}_{\lambda}\right]^{-1}.$
(35)

We note that for real matrices an analogous but formally somewhat different
presentation can be found in Furche Furche (2001b). Poles in the polarisation
propagator will be located at those frequencies that fulfill the condition
$\operatorname{det}\\{\mathbf{P}_{0}^{-1}(\omega)-\mathbf{I}_{\lambda}\\}=0$.
Given the spectral representation for $\mathbf{P_{0}(\omega)}$ in Eq.
(missing)9 it is straighforward to construct its inverse, both being diagonal
matrices:

$\displaystyle\mathbf{P}_{0}^{-1}(\omega)$
$\displaystyle=\begin{pmatrix}-\Delta E_{ia}&\mathbf{0}\\\ \mathbf{0}&\Delta
E_{ia}\end{pmatrix}+\omega\begin{pmatrix}\mathbf{1}&\mathbf{0}\\\
\mathbf{0}&\mathbf{-1}\end{pmatrix}$
$\displaystyle=\mathbf{D}+\omega\mathbf{\Delta}.$ (36)

In going to the last line we have introduced the supermatrices $\mathbf{D}$
and $\mathbf{\Delta}$ for convenience. The polarisation propagator can then be
compactly expressed as:
$\mathbf{P}_{\lambda}(\omega)=\left[\mathbf{D}-\mathbf{I}_{\lambda}+\omega\mathbf{\Delta}\right]^{-1}=-\left[\mathbf{M}_{\lambda}-\omega\mathbf{\Delta}\right]^{-1}$,
where the square supermatrix $\mathbf{M}_{\lambda}$ corresponds to the one
introduced in Eq. (missing)10. Since the determinant is invariant under
similarity transformations, it is convenient to introduce the matrix:

$\displaystyle\mathbf{Z}_{\lambda}=\begin{pmatrix}\mathbf{X}_{\lambda}&\mathbf{Y}_{\lambda}^{*}\\\
\mathbf{Y}_{\lambda}&\mathbf{X}_{\lambda}^{*}\end{pmatrix}$ (37)

that fulfils the generalised eigenvalue equation

$\displaystyle\mathbf{M}_{\lambda}\mathbf{Z}_{\lambda}=\mathbf{\Delta}\mathbf{Z}_{\lambda}\mathbf{\Omega}_{\lambda}$
(38)

which is nothing but a re-statement of Eq. (missing)10 taking into account the
complex conjugate solution. The matrix $\mathbf{Z}_{\lambda}$ satisfies the
”symplectic” normalisation condition
$\mathbf{Z}_{\lambda}\mathbf{\Delta}\mathbf{Z}_{\lambda}^{\dagger}=\mathbf{\Delta}$;
this constraint, together with $\mathbf{\Delta}$ being idempotent, are
sufficient to show that $\mathbf{Z}_{\lambda}$ is unitary, and that the
condition

$\displaystyle\mathbf{\Delta}\mathbf{Z}_{\lambda}^{\dagger}\mathbf{\Delta}=\mathbf{Z}_{\lambda}^{-1}$
(39)

holds.

The relation in Eq. (missing)39 can be exploited to construct the spectral
representation of the polarisation propagator starting from Eq. (missing)38:

$\displaystyle\\{\mathbf{M}_{\lambda}-\omega\mathbf{\Delta}\\}\mathbf{Z}_{\lambda}$
$\displaystyle=\mathbf{\Delta}\mathbf{Z}_{\lambda}(\mathbf{\Omega}_{\lambda}-\omega)$
(40) $\displaystyle\Rightarrow\mathbf{M}_{\lambda}-\omega\mathbf{\Delta}$
$\displaystyle=\mathbf{\Delta}\mathbf{Z}_{\lambda}(\mathbf{\Omega}_{\lambda}-\omega)\mathbf{\Delta}\mathbf{Z}_{\lambda}^{\dagger}\mathbf{\Delta}.$
(41)

Inverting the right hand side of the previous equation one finally obtains:

$\displaystyle\mathbf{P}_{\lambda}(\omega)=\mathbf{Z}_{\lambda}\left(\omega-\mathbf{\Omega}_{\lambda}\right)^{-1}\mathbf{\Delta}\mathbf{Z}_{\lambda}^{\dagger}.$
(42)

The presence of a simple pole analytic structure for both polarisation
propagators entering LABEL:eq:EcP is of fundamental importance because it
allows the straightforward frequency integration necessary to evaluate the
correlation part of the two body density matrix $\mathcal{P}_{\lambda}$.

## References

  * Giuliani and Vignale (2005) G. Giuliani and G. Vignale,  _Quantum theory of the electron liquid_ (Cambridge University Press, Cambridge, UK; New York, 2005).
  * Endo  _et al._ (1999) T. Endo, M. Horiuchi, Y. Takada, and H. Yasuhara, Phys. Rev. B 59, 7367 (1999).
  * Hubbard (1958) J. Hubbard, Proc. R. Soc. A Math. Phys. Eng. Sci. 244, 199 (1958).
  * Gori-Giorgi  _et al._ (2000) P. Gori-Giorgi, F. Sacchetti, and G. Bachelet, Phys. Rev. B 61, 7353 (2000).
  * Sturm and Gusarov (2000) K. Sturm and A. Gusarov, Phys. Rev. B 62, 16474 (2000).
  * Caruso  _et al._ (2013) F. Caruso, D. Rohr, M. Hellgren, X. Ren, P. Rinke, A. Rubio, and M. Scheffler, Phys. Rev. Lett. 110, 146403 (2013).
  * Stefanucci and van Leeuwen (2013) G. Stefanucci and R. van Leeuwen,  _Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction_ (Cambridge University Press, Cambridge, UK; New York, 2013).
  * Furche (2001a) F. Furche, Phys. Rev. B 64, 195120 (2001a).
  * Miyake  _et al._ (2002) T. Miyake, F. Aryasetiawan, T. Kotani, M. van Schilfgaarde, M. Usuda, and K. Terakura, Phys. Rev. B 66, 245103 (2002).
  * Niquet  _et al._ (2003) Y. Niquet, M. Fuchs, and X. Gonze, Phys. Rev. A 68, 032507 (2003).
  * Jiang and Engel (2007) H. Jiang and E. Engel, J. Chem. Phys. 127, 184108 (2007), arXiv:0707.1615 .
  * Furche (2008) F. Furche, J. Chem. Phys. 129, 114105 (2008).
  * Scuseria  _et al._ (2008) G. E. Scuseria, T. M. Henderson, and D. C. Sorensen, J. Chem. Phys. 129, 231101 (2008).
  * Harl and Kresse (2008) J. Harl and G. Kresse, Phys. Rev. B 77, 045136 (2008).
  * Harl  _et al._ (2010) J. Harl, L. Schimka, and G. Kresse, Phys. Rev. B 81, 115126 (2010).
  * Toulouse  _et al._ (2009) J. Toulouse, I. Gerber, G. Jansen, A. Savin, and J. G. Ángyán, Phys. Rev. Lett. 102, 096404 (2009).
  * Toulouse  _et al._ (2010) J. Toulouse, W. Zhu, J. G. Ángyán, and A. Savin, Phys. Rev. A 82, 032502 (2010).
  * Jansen  _et al._ (2010) G. Jansen, R.-F. Liu, and J. G. Ángyán, J. Chem. Phys. 133, 154106 (2010).
  * Olsen and Thygesen (2013) T. Olsen and K. S. Thygesen, Phys. Rev. B 87, 075111 (2013).
  * Heßelmann and Görling (2010) A. Heßelmann and A. Görling, Mol. Phys. 108, 359 (2010).
  * Macher  _et al._ (2014) M. Macher, J. Klimeš, C. Franchini, and G. Kresse, J. Chem. Phys. 140, 084502 (2014).
  * Teale  _et al._ (2010) A. M. Teale, S. Coriani, and T. Helgaker, J. Chem. Phys. 132, 164115 (2010).
  * Bleiziffer  _et al._ (2012) P. Bleiziffer, A. Heßelmann, and A. Görling, J. Chem. Phys. 136, 134102 (2012).
  * Kaltak  _et al._ (2014) M. Kaltak, J. Klimeš, and G. Kresse, J. Chem. Theory Comput. 10, 2498 (2014).
  * Grüneis  _et al._ (2009) A. Grüneis, M. Marsman, J. Harl, L. Schimka, and G. Kresse, J. Chem. Phys. 131, 154115 (2009).
  * Schimka  _et al._ (2010) L. Schimka, J. Harl, A. Stroppa, A. Grüneis, M. Marsman, F. Mittendorfer, and G. Kresse, Nat. Mater. 9, 741 (2010).
  * García-González  _et al._ (2007) P. García-González, J. J. Fernández, A. Marini, and A. Rubio, J. Phys. Chem. A 111, 12458 (2007).
  * Harl and Kresse (2009) J. Harl and G. Kresse, Phys. Rev. Lett. 103, 4 (2009).
  * Ren  _et al._ (2011) X. Ren, A. Tkatchenko, P. Rinke, and M. Scheffler, Phys. Rev. Lett. 106, 153003 (2011).
  * Mori-Sánchez  _et al._ (2012) P. Mori-Sánchez, A. Cohen, and W. Yang, Phys. Rev. A 85, 042507 (2012), arXiv:arXiv:0903.4403v1 .
  * Henderson and Scuseria (2010) T. M. Henderson and G. E. Scuseria, Mol. Phys. 108, 2511 (2010).
  * Olsen and Thygesen (2014) T. Olsen and K. S. Thygesen, J. Chem. Phys. 140, 164116 (2014).
  * Albrecht  _et al._ (1998) S. Albrecht, L. Reining, R. Del Sole, and G. Onida, Phys. Rev. Lett. 80, 4510 (1998).
  * Rohlfing and Louie (2000) M. Rohlfing and S. G. Louie, Phys. Rev. B 62, 4927 (2000).
  * Yan  _et al._ (2012) J. Yan, K. W. Jacobsen, and K. S. Thygesen, Phys. Rev. B 86, 045208 (2012).
  * Gatti and Sottile (2013) M. Gatti and F. Sottile, Phys. Rev. B 88, 155113 (2013).
  * Del Puerto  _et al._ (2006) M. L. Del Puerto, M. L. Tiago, and J. Chelikowsky, Phys. Rev. Lett. 97, 1 (2006), arXiv:0608521 [cond-mat] .
  * Duchemin and Blase (2013) I. Duchemin and X. Blase, Phys. Rev. B 87, 245412 (2013), arXiv:1301.7230 .
  * Linderberg (1967) J. Linderberg, Int. J. Quantum Chem. 1, 719 (1967).
  * Gunnarsson and Lundqvist (1976) O. Gunnarsson and B. Lundqvist, Phys. Rev. B 13, 4274 (1976).
  * Langreth and Perdew (1977) D. C. Langreth and J. P. Perdew, Phys. Rev. B 15, 2884 (1977).
  * Harris and Jones (1974) J. Harris and R. O. Jones, J. Phys. F Met. Phys. 4, 1170 (1974).
  * Singwi  _et al._ (1968) K. S. Singwi, M. Tosi, R. Land, and A. Sjölander, Phys. Rev. 176, 589 (1968).
  * Klimeš  _et al._ (2015) J. Klimeš, M. Kaltak, E. Maggio, and G. Kresse, J. Chem. Phys. 143, 102816 (2015).
  * Harris (1984) J. Harris, Phys. Rev. A 29, 1648 (1984).
  * van Aggelen  _et al._ (2014) H. van Aggelen, Y. Yang, and W. Yang, J. Chem. Phys. 140, 18A511 (2014).
  * Gross  _et al._ (1991) E. K. U. Gross, E. Runge, and O. Heinonen,  _Many-Particle Theory_ (IOP Publishing, Bristol, Philadelphia and New York, 1991).
  * Toyoda  _et al._ (1975) K. Toyoda, T. Iwai, and O. Tanimoto, Prog. Theor. Phys. 54, 333 (1975).
  * Fetter and Walecka (1971) A. L. Fetter and J. D. Walecka,  _Quantum Theory of Many-Particle Systems_ (McGraw-Hill, Inc., New York, 1971) p. 601.
  * Strinati (1988) G. Strinati, La Riv. del Nuovo Cim. 11, 1 (1988).
  * Hanke and Sham (1979) W. Hanke and L. J. Sham, Phys. Rev. Lett. 43, 387 (1979).
  * Hanke and Sham (1980) W. Hanke and L. J. Sham, Phys. Rev. B 21, 4656 (1980).
  * Fuchs  _et al._ (2008) F. Fuchs, C. Rödl, a. Schleife, and F. Bechstedt, Phys. Rev. B 78, 085103 (2008).
  * Sander  _et al._ (2015) T. Sander, E. Maggio, and G. Kresse, Phys. Rev. B 92, 045209 (2015).
  * Rödl  _et al._ (2008) C. Rödl, F. Fuchs, J. Furthmüller, and F. Bechstedt, Phys. Rev. B 77, 184408 (2008).
  * Ring and Schuck (1980) P. Ring and P. Schuck,  _The Nuclear Many-Body Problem_ , 3rd ed. (Springer-Verlag, Berlin / Heidelberg, 1980).
  * Jørgensen and Simons (1981) P. Jørgensen and J. Simons,  _Second Quantization-Based Methods in Quantum Chemistry_ (Academic Press, New York, 1981).
  * Ángyán  _et al._ (2011) J. G. Ángyán, R.-F. Liu, J. Toulouse, and G. Jansen, J. Chem. Theory Comput. 7, 3116 (2011).
  * Lin  _et al._ (2001) C. Lin, F. H. Zong, and D. M. Ceperley, Phys. Rev. E 64, 016702 (2001), arXiv:0101339 [cond-mat] .
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996).
  * Kresse (1996) G. Kresse, Phys. Rev. B 54, 11169 (1996).
  * Shishkin and Kresse (2007) M. Shishkin and G. Kresse, Phys. Rev. B 75, 235102 (2007).
  * Shishkin and Kresse (2006) M. Shishkin and G. Kresse, Phys. Rev. B 74, 35101 (2006).
  * Gori-Giorgi and Perdew (2002) P. Gori-Giorgi and J. P. Perdew, Phys. Rev. B 66, 165118 (2002).
  * Klimeš  _et al._ (2014) J. Klimeš, M. Kaltak, and G. Kresse, Phys. Rev. B 90, 075125 (2014), arXiv:arXiv:1404.3101v1 .
  * Perdew and Zunger (1981) J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
  * Marsman  _et al._ (2009) M. Marsman, A. Grüneis, J. Paier, and G. Kresse, J. Chem. Phys. 130, 184103 (2009).
  * Ceperley and Adler (1980) D. M. Ceperley and B. J. Adler, Phys. Rev. Lett. 45, 566 (1980).
  * Ortiz and Ballone (1994) G. Ortiz and P. Ballone, Phys. Rev. B 50, 1391 (1994).
  * Fraser  _et al._ (1996) L. M. Fraser, W. M. C. Foulkes, G. Rajagopal, R. J. Needs, S. D. Kenny, and a. J. Williamson, Phys. Rev. B 53, 1814 (1996).
  * Yasuhara (1972) H. Yasuhara, Solid State Commun. 11, 1481 (1972).
  * Paier  _et al._ (2012) J. Paier, X. Ren, P. Rinke, G. E. Scuseria, A. Grüneis, G. Kresse, and M. Scheffler, New J. Phys. 14, 043002 (2012).
  * Paldus and Čížek (1975) J. Paldus and J. Čížek, Adv. Quantum Chem. 9, 105 (1975).
  * Patrick and Thygesen (2015) C. E. Patrick and K. S. Thygesen, J. Chem. Phys. 143, 102802 (2015).
  * Colonna  _et al._ (2014) N. Colonna, M. Hellgren, and S. de Gironcoli, Phys. Rev. B 90, 125150 (2014).
  * Görling (1998) A. Görling, Int. J. Quantum Chem. 69, 265 (1998).
  * Hirata  _et al._ (2001) S. Hirata, S. Ivanov, I. Grabowski, R. J. Bartlett, K. Burke, and J. D. Talman, J. Chem. Phys. 115, 1635 (2001).
  * Uimonen  _et al._ (2015) A.-M. Uimonen, G. Stefanucci, Y. Pavlyukh, and R. van Leeuwen, Phys. Rev. B 91, 115104 (2015).
  * Heßelmann and Görling (2011) A. Heßelmann and A. Görling, Mol. Phys. 109, 2473 (2011).
  * Perdew and Wang (1992) J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
  * van Schilfgaarde  _et al._ (2006) M. van Schilfgaarde, T. Kotani, and S. Faleev, Phys. Rev. Lett. 96, 226402 (2006).
  * Starke and Kresse (2012) R. Starke and G. Kresse, Phys. Rev. B 85, 75119 (2012).
  * Furche (2001b) F. Furche, J. Chem. Phys. 114, 5982 (2001b).