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

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{\o}ller-Plesset perturbation theory, it is self-interaction free in second order. Results for the correlatione 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 occurence 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.


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 [1]. As the Coulomb repulsion among electrons is turned on, correlations between particles (prevalently of opposite spin) [2] set in. This mechanism produces even stronger deviations from the freeelectron 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 [3]. 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 [4], its complete characterisation is still lacking [5]. 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 * Electronic address: georg.kresse@univie.ac.at 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 = − 1 i δΣ δG . Hence, the degree of sophistication in the BSE approach can be systematically improved by including more selfenergy 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 [6]. In both cases, a sub-class of self-energy diagrams, the so-called "bubble" diagrams are summed up to infinite order [7]. While the RPA approach has witnessed alternate fortune over the years, there has been a strong resurgence of interest in the field recently [8][9][10][11][12][13][14][15][16][17][18]. Application of the RPA method to molecules [8,[19][20][21][22][23] and solid state systems [9,15,[24][25][26][27] 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) [14]. On the other hand, systematic underestimation of binding energies [8,28,29] and non-physical features for dissociation curves in diatomic systems [30,31] are some shortcomings of the RPA. To improve on this, the BSE has been exploited to describe hydrogen dissociation [32], to gain access to optical properties of semiconductors and insulators [33][34][35] and to study excitonic effects in extended [36] and molecular systems [37,38]. 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 [39] or Kohn-Sham [40][41][42] 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 [43] and between spin-parallel electrons [2] 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 section II B 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.

A. Exchange and Correlation Energy evaluated along the AC path
The adiabatic connection formalism was originally introduced in the DFT framework [40][41][42] 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 oneand two-body contributions. The two-body part of the correlation energy can be estimated solving the BSE, as described in section II B. The impact of one-body contributions has been the subject of recent studies [29,44]. 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: The field operatorsψ(x) and their adjoint are given in the ordinary Schrödinger picture; the Coulomb potential is defined in the spin-space basis as v(x, x ) = δ σ,σ |r−r | with the usual space-time notation (x, t) = (r, σ, t) for the space (r), spin (σ) and time (t) variables. The oneparticle operatorV xc λ 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 Σ xc (ω).
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 potentialV xc 0 . The latter also varies with varying coupling constant and different switching can be realised forV xc λ . These are either designed to keep the electron density fixed along the AC path [41], or assume a linear dependence on λ [45,46]. Here we use the linear switchingV xc 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: If the symmetry of the ground state does not change along the AC path, then |Ψ 0 evolves into |Ψ 1 for non-degenerate ground states thanks to the Gell-Mann and Low theorem. It is easy to see that with Then the correlation energy can be expressed as (see also [17,44]): where we have assumed that the Hellmann-Feynman theorem holds in order to set ∂ ∂λ Ψ λ |H λ | Ψ λ = Ψ λ ∂H λ ∂λ Ψ λ . The previous equation reproduces the standard expression found in the literature [17].
The calculation of the correlation energy thus requires to evaluate the term dE λ dλ ≡Ė λ 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: The last term in the previous equation represents contributions related to the change of the one-particle Green's function to the correlation energy [29,44]. In the present study, we disregard this term and concentrate on the first term on the r.h.s. of Eq. (4). In fact, in second order, the last term is exactly zero for the HEG (see e.g. chapter 22 in Ref. 47), 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 λ [48], fulfilling the BSE as detailed in the next section. The correlation energy then reads: 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 [49].
Here V and W 0 are the four-point generalisations of the bare Coulomb (v) and screened (w) interactions, and the shorthand (i) = (r i , σ 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 λ and, correspondingly, the polarisation propagator P λ to be evaluated for each point along the AC path. The resulting BSE will then be parametrically dependent on two quantities: the coupling constant λ and the frequency variable ω, provided that the kernel is static. This approximation for the kernel implies that w(x 1 , x 2 ) in Eq. (7) is evaluated as , 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: , for the resonant (R) and antiresonant (A) component labelled by the "superindex" M ≡ (i, a). The BSE can then be recast in this basis to yield: 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: where f n is the (fractional) occupancy for the orbital n andf n = 1 − f n , the indices i and a in the summations above go over the occupied and unoccupied single particle states respectively. ∆E ia = ( a − i ) is the independent particle energy difference.
Since the BSE kernel is static, it is possible to invert the matrix equation Eq. (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 P λ (ω) can be constructed solving the non-Hermitian generalised eigenvalue problem (EVP): which does not depend on ω. The matrix elements for A and B are given by [53]: 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. [54]. : Goldstone diagrams for the particle-hole interaction included in the BSE kernel: particle states are shown as right arrows → (labelled by unoccupied ground state orbitals, a, b) and holes by left arrows ← (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.

C. 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. 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, k, k . We have: where q is the momentum of the photon impinging the system and the ± sign applies to resonant and antiresonant electron-hole pairs respectively. We denote the Fourier components of the Coulomb potential with v q (r, r ), then the two electron integrals in A are given by integrating over spin σ, ς and space variables: and similarly for the B matrix elements: For the Hartree components in Eq. (14) and Eq. (16) the interaction potential comprises of a single mode q equal to the difference k − k . Whereas, the screened Coulomb interaction in Eq. (15) and Eq. (17) can be decomposed in its Fourier modes:

5
The spin structure of the matrix elements in Eqs. (14), (16) and (17) spans the singlet subspace, where electron and hole states of each pair have the same spin. The exchange part in the A matrix (Eq. (15)) is defined in principle also on the triplet subspace where the incoming electron and hole have opposite spin orientations [34]. If the Hamiltonian is time-reversal invariant (i.e. spinorbit coupling is not present) then the singlet and triplet solutions factorise [55] 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 ψ n,−k (r) = ψ * n,k (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: It is well-known that the A matrix is Hermitian in this representation [56], hence the transformation above will simply map A * into A. For the exchange part in Eq. (17), the transformed matrix elements in B are given by (with Here, we have dropped the bar over the k index and suppressed the spin indices for simplicity. The Hartree component turns out to be identical to the expression in Eq. (14), hence requiring no extra computational step. Furthermore, the transformation applied to B * gives exactly the same matrix elements in Eq. (20) as we now show. Starting with the matrix elements in Eq. (17) and applying the transformation k → −k to the antiresonant electron-hole pair (a, i), one obtains: In the static approximation to the screening potential w q (r, r ), the dielectric function at the RPA level is an even function of |q| and has zero imaginary part, hence it is possible to replace w * q (r, r ) with w −q (r, r ) in the equation above. Finally, by swapping the integration variables, it is easy to see that the exchange component in B * reduces to Eq. (20). The reasoning above can be straightforwardly repeated for the Hartree component in B * to give aj|V |ib . This result completes the proof of Eq. (19). In the following we will drop the prime when referring to the B matrix.
We can now proceed to reduce the dimensionality of the EVP (10). Adding and subtracting the individual equations contained in Eq. (10) and then solving one of the two equations, let's say with respect to (X − Y), one obtains the squared problem [57]: The expression above can be converted into a conventional EVP if the matrices are positive definite. Then the matrix (A − B) 1 2 is unique and the eigenvalue problem finally reads: The difference P λ (ω) − P 0 (ω) in Eq. (5) can be integrated over ω thanks to the spectral decomposition in Eq. (9) and in Eq. (B8). The correlation energy can then be expressed using the correlation part of the two-particle density matrix P λ : [8,58].

III. COMPUTATIONAL DETAILS
The HEG is in principle specified by the average electron density n = Ne 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,↑ or N e,↓ ). These quantities are easily translated into the Wigner-Seitz radius, r s given by 4 3 πr 3 s = V Ne and the spin polarisation ζ = |N e,↑ −N e,↓ | Ne . The discretisation of the Bloch wavevectors introduced in section II C, however, requires also a specific choice for the k-point sampling and for the symmetry of the simulation cell.
Calculations were performed on a Γ−point centred, uniform k-point mesh of dimensions N k × N k × 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 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×3×3 k-point grid with 2 electrons in the unit cell corresponds exactly to a 3 × 3 × 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 59.
The VASP code [60,61] 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 , ζ) 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 cut and includes plane waves with kinetic energy smaller than the given cutoff value. This step is followed by a selfconsistent evaluation of the quasi-particle energies at the QP − GW 0 level [62,63]. 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 [25] 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 cut . The frequency sampling has been carried out on a linear grid, specified by its maximum value, set equal to 1.5 E cut for all (r s , ζ). 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 ≈ 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 Eq. (22) and has been implemented in VASP; all virtual unoccupied orbitals spanned by the plane-wave basis were included. The diagonalisation of Eq. (22) 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. (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.

IV. RESULTS AND DISCUSSION
We start by considering the correlation energy's convergence properties for a given (r s , ζ) 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 [14]. The complete basis set limit has been estimated for both the paramagnetic (ζ = 0) and ferromagnetic (ζ = 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. 64 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 k-point grid density is reported in the right panel of Figure 2. Here we report the total energy for any given k-point mesh minus the value at the preceding k-point sampling divided by the energy value obtained with the most accurate k-point grid for a given level of theory. For each 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. 15), 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 k-point sampling we employ an extrapolation scheme adapted from Ref. 65: 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 ∞ 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 r c (N k ) + E ∞ c (N k )) using N k = 3, i.e. 3 × 3 × 3 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. [43] for the paramagnetic system and of Ref. [66] for the spin-polarised case. We tested the agreement for the case (r s = 1, ζ = 0) where we found it necessary to include up to N k = 18 points to reproduce the analytical results within 4 meV. Fewer 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 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 (rs, ζ). 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 k-mesh and the value at the preceding sampling. calculations with the replacement: where the correlation energy is computed at the considered level of theory whereas the ∆ N k correction corrects for the k-point incompleteness error. This value is difficult to calculate accurately for the BSE. In previous work [67], 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 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 k-point error using the computationally efficient RPA, but to reduce the RPA k-point correction by a factor 2 3 if exchange is included. The correlation energies computed for a set of (r s , ζ) 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) [68] (indicated by circles in the figure) and by Ortiz and Ballone (OB in the following) [69] ('×' symbols in the figure). Finite size effects impact also these benchmarks [70] 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@GW 0 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 (dashdotted 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 ≈20% for r s ≥4.0 a.u.. This pattern is also observed for the ferromagnetic system, with the absolute deviation, however, decreasing to ≈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 holehole) ladder diagrams to describe short range interactions [71]. 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  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 [1]. 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. 18. 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. 67. In second order, the correlation energy is then simply given by the standard textbook equation: From this, two important observations follow.
(i) The Hartree contribution ab|ij and the exchange term − ab|ji 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 ζ 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 [25]. 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 important [31].
(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 spinorbital 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. (26) is evaluated for spatial orbitals [73], 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) [74]. However larger deviations are observed in the "metallic" density region (r s ∈ [0.5, 5.0]). Also, a recent study by the de Gironcoli's group [75] 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 oneparticle 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 [76]. An issue that this approach faces is the occurrence of imaginary frequencies when the response function is diagonalised (see Eq. (10) in [75]). 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 [75]. 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 [23]) and on the critical dependence on the basis set size [77].
Also in the present BSE approach instabilities at lower electron densities can be present (as in our case for r s ≥5). These are witnessed as imaginary frequencies Ω appearing in Eq. (22). The instabilities can cause the matrices A + B or 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 [57]. In the present case, however, we find that only the A + B matrix has negative eigenvalues, for both the GW 0 and LDA reference states. This finding is consistent with a recent theoretical investigation [78] which has confirmed the presence of negativefrequency modes in the BSE polarisation propagator.
In our case the instabilities originate from the particlehole 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 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 A matrix does not contribute to the correlation energy in the perturbation expansion to second order [58,79]. 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 [25] 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 B matrix (containing the sum of the direct Coulomb inter-action and the bare exchange interaction, diagrams (B1) and (B2) in Eq. (23)) [18]. In RPAsX, the propagator now includes, as it should, the exchange term in the 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 ≥2 the RPAsX estimate of the correlation energy remains within 6% of the CA estimate for ζ = 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 particleparticle 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 Υ(ζ; r s ) is commonly introduced: In Figure 5 we report the spin enhancement function computed at r s =2 versus the spin polarisation ζ for various levels of theory. Results are compared against the Perdew-Wang interpolation formula [80]   shown by squares in the figure, reproduce exceedingly well the benchmark calculations for all values of ζ 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 Υ very accurately, it is reasonable to assume that RPAsX will also be very accurate for fractional values of ζ.

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 ≥ 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. (4)). This is not the case for real materials, where self-consistency of the orbitals [81] or changes of the one particle Green's function when going to the interacting case are very relevant [44]. 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.  6: Graphical representation for the two-body Green's function. We represent the time ordering t1 = t2 which corresponds to having two electron-hole pairs interacting instantaneously action picture (as it is commonly done in many body perturbation theory):ψ(1) = e iH0t1ψ (x 1 )e −iH0t1 and ψ † (2) = e iH0t2ψ † (x 2 )e −iH0t2 . The n-body Green's operator is as usually defined as: where T is the time ordering operator. It is worth pointing out that the operators G n do not depend on the value of the coupling constant λ. For n = 1 the corresponding (greater) Green's function is obtained by evaluating the expectation value of the Green's operator for the groundstate wave function at coupling λ, Ψ λ . Since we are only interested in equal time limits here, we can restrict the second time to fulfil t = t + 0 + = t + : At equal time limits, the Green's function can be also related to the density matrix n λ ≡ n λ (x, y) = −iG 1,λ (x, t; y, t).
For the two particle Green's function, it is standard texbook material [47] (commutator relations or Wick's theorem) to show that it satisfies: 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 Eq. (A1). The labelling convention for fourpoint quantities is given in Figure 6 and it is consistent with Ref. 82 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. (A4) in terms of the ordinary density operatorsn and their expectation values: where in going to the last line, we have introduced the density fluctuation operator δn λ ≡ n λ − n λ . The last term in Eq. (A5) can be identified with the polarisation propagator (also called density-fluctuation densityfluctuation response function), i.e. P λ ≡ δnδn λ . One can then rearrange Eq. (A5) into the usual relation connecting the two particle propagator and the polarisation propagator [7]: P λ (x 1 , t + , x 2 , t; x 3 , t, x 4 , t + ) = − G 2,λ (x 1 , t + , x 2 , t; x 3 , t, x 4 , t + ) + G 1,λ (x 1 , t + ; x 3 , t)G 1,λ (x 2 , t; x 4 , t + ).
Using Eq. (A4) and Eq. (A6), one can rewrite the correlation part of G 2,λ in the ACFDT equation using the polarisation propagator and obtain the following compact equation for the correlation energy: dx i V G 1,λ (x 1 , t + ; x 4 , t + )G 1,λ (x 2 , t; x 3 , t) −P λ (x 1 , t + , x 2 , t; x 3 , t, x 4 , t + ) (A7) From Eq. (A7) it is easy to separate out the exchange contribution: The remaining contribution in Eq. (A7) can be replaced into Eq. (3) to form the correlation energy assuming that the one body contributions (contained in the exchange part) remain constant along the AC path. Eq. (5) is then obtained by Fourier transforming with respect to the infinitesimal time difference. The starting point to derive the spectral representation for the interacting polarisation propagator P λ (ω) is to invert Eq. (8) in the particle-hole basis: We note that for real matrices an analogous but formally somewhat different presentation can be found in Furche [83]. Poles in the polarisation propagator will be located at those frequencies that fulfill the condition det{P −1 0 (ω) − I λ } = 0. Given the spectral representation for P 0 (ω) in Eq. (9) it is straighforward to construct its inverse, both being diagonal matrices: In going to the last line we have introduced the supermatrices D and ∆ for convenience. The polarisation propagator can then be compactly expressed as: where the square supermatrix M λ corresponds to the one introduced in Eq. (10). Since the determinant is invariant under similarity transformations, it is convenient to introduce the matrix: that fulfils the generalised eigenvalue equation which is nothing but a re-statement of Eq. (10) taking into account the complex conjugate solution. The matrix Z λ satisfies the "symplectic" normalisation condition Z λ ∆Z † λ = ∆; this constraint, together with ∆ being idempotent, are sufficient to show that Z λ is unitary, and that the condition holds.
The relation in Eq. (B5) can be exploited to construct the spectral representation of the polarisation propagator starting from Eq. (B4): Inverting the right hand side of the previous equation one finally obtains: The presence of a simple pole analytic structure for both polarisation propagators entering Eq. (5) is of fundamental importance because it allows the straightforward frequency integration necessary to evaluate the correlation part of the two body density matrix P λ .