Linear response theory (Oddershede and Jorgensen, 1977; Casida, 1995; Petersilka et al., 1996) in MBPT is described by the Bethe-Salpeter equation. (Strinati, 1988). The standard BSE within the static GW approximation (referred to as BSE@GW in this work, which means the use of GW quasiparticle energies to build the independent-particle excitation energies and of the GW self-energy to build the static exchange-correlation kernel) can be recast, assuming a closed-shell reference state, as a non-Hermitian eigenvalue problem known as Casida equations: A λ B λ − B λ − A λ X m λ Y m λ = Ω m λ X m λ Y m λ , (8) where Ω m λ is the mth excitation energy with eigenvector ( X m λ Y m λ ) ⊺ at interaction strength λ, ^⊺ is the matrix transpose, and we have assumed real-valued spatial orbitals. The non-interacting and physical systems correspond to λ = 0 and 1, respectively. The matrices A ^λ and B ^λ are of size OV × OV, where O and V are the number of occupied and virtual orbitals, respectively, and O+ V is the total number of spatial orbitals. Introducing the so-called Mulliken notation for the bare two-electron integrals ( i j | k l ) = ∬ d r 1 d r 2 ϕ i ( r 1 ) ϕ j ( r 1 ) v c ( r 1 − r 2 ) ϕ k ( r 2 ) ϕ l ( r 2 ) , (9) and the corresponding (static) screened Coulomb potential matrix elements W i j , k l ( ω = 0 ) = ∬ d r 1 d r 2 ϕ i ( r 1 ) ϕ j ( r 1 ) W ( r 1 , r 2 ; ω = 0 ) ϕ k ( r 2 ) ϕ l ( r 2 ) , (10) the BSE matrix elements read (Maggio and Kresse, 2016). A a r , b s λ , σ σ ′ = δ a b δ r s ( ϵ r QP − ϵ a QP ) + λ α σ σ ′ ( a r | s b ) − W a b , s r ( ω = 0 ) , (11a) B a r , b s λ , σ σ ′ = λ α σ σ ′ ( a r | b s ) − W a s , b r ( ω = 0 ) , (11b) where ϵ i QP are the GW quasiparticle energies, and α _↑↓ = 2 and α _↑↑ = 0 for singlet (i.e., spin-conserved) and triplet (i.e., spin-flip) excitations, respectively.

In the absence of instabilities (i.e., when A ^λ − B ^λ is positive-definite), (Dreuw and Head-Gordon, 2005), Eq. 8 is usually transformed into an Hermitian eigenvalue problem of half the dimension ( A λ − B λ ) 1 / 2 ( A λ + B λ ) ( A λ − B λ ) 1 / 2 V m λ = ( Ω m λ ) 2 V m λ , (12) where the excitation amplitudes are ( X λ + Y λ ) m = ( Ω m λ ) − 1 / 2 ( A λ − B λ ) + 1 / 2 V m λ , (13a) ( X λ − Y λ ) m = ( Ω m λ ) + 1 / 2 ( A λ − B λ ) − 1 / 2 V m λ . (13b)

Singlet ( Ω m ↑ ↓ ≡ Ω m λ = 1 , ↑ ↓ ) and triplet ( Ω m ↑ ↑ ≡ Ω m λ = 1 , ↑ ↑ ) excitation energies are obtained by diagonalizing Eq. 8 at λ = 1.

Our goal here is to compare the BSE correlation energy E c BSE obtained using two formulas, namely the trace (or plasmon) formula (Rowe, 1968; Ring and Schuck, 1980) and the expression obtained using the adiabatic-connection fluctuation-dissipation theorem (ACFDT) formalism. (Furche and Van Voorhis, 2005; Toulouse et al., 2009, 2010; Hellgren and von Barth, 2010; Angyan et al., 2011; Heßelmann and Görling, 2011; Colonna et al., 2014; Maggio and Kresse, 2016; Holzer et al., 2018; Loos et al., 2020). The two approaches have been recently compared at the random-phase approximation (RPA) level for the case of Be₂, (Li et al., 2020), showing similar improved performances at the RPA@GW@PBE level with respect to the RPA@PBE level and an impressive accuracy by introducing BSE (BSE@GW@HF) correction in the trace formula. Here we would like to get more insights into the quality of these two approaches.

The ground-state correlation energy within the trace formula is calculated as E c Tr@BSE = E c , ↑ ↓ Tr@BSE + E c , ↑ ↑ Tr@BSE = 1 2 ∑ m Ω m ↑ ↓ − Tr ( A ↑ ↓ ) + 1 2 ∑ m Ω m ↑ ↑ − Tr ( A ↑ ↑ ) , (14) where A ^σσ′ ≡ A ^λ=1,σσ′ is defined in Eq.11a and   Tr denotes the matrix trace. We note that the trace formula is an approximate expression of the correlation energy since it relies on the so-called quasi-boson approximation and on the killing condition on the zeroth-order Slater determinant ground state (Li et al., 2020 for more details). Note that here both sums in Eq. 14 run over all resonant (hence real- and complex-valued) excitation energies while they are usually restricted to the real-valued resonant BSE excitation energies. Thus, the Tr@BSE correlation energy is potentially a complex-valued function in the presence of singlet and/or triplet instabilities.

The ACFDT formalism, (Furche and Van Voorhis, 2005), instead, provides an in-principle exact expression for the correlation energy within time-dependent density-functional theory (TDDFT). (Runge and Gross, 1984; Petersilka et al., 1996; Ullrich, 2012). In practice, however, one always ends up with an approximate expression, which quality relies on the approximations to the exchange-correlation potential of the KS system and to the kernel of the TDDFT linear response equations. In this work, therefore, we use the ACFDT expression within the BSE formalism and we explore how well it performs and how it compares to the trace Eq. 14.

Within the ACFDT framework, only the singlet states do contribute for a closed-shell ground state, and the ground-state BSE correlation energy E c AC@BSE = 1 2 ∫ 0 1 d λ Tr ( K ↑ ↓ P λ , ↑ ↓ ) (15) is obtained via integration along the adiabatic connection path from the non-interacting system at λ = 0 to the physical system λ = 1, where K = A ̃ λ = 1 B λ = 1 B λ = 1 A ̃ λ = 1 (16) is the interaction kernel, (Angyan et al., 2011; Holzer et al., 2018; Loos et al., 2020) A ̃ a r , b s λ , σ σ ′ = α σ σ ′ λ ( a r | s b ) , and P λ = Y λ ( Y λ ) ⊺ Y λ ( X λ ) ⊺ X λ ( Y λ ) ⊺ X λ ( X λ ) ⊺ − 0 0 0 1 (17) is the correlation part of the two-body density matrix at interaction strength λ. Here again, the AC@BSE correlation energy might become complex-valued in the presence of singlet instabilities.

Note that the trace and ACFDT formulas yield, for any set of eigenstates, the same correlation energy at the RPA level. (Angyan et al., 2011). Moreover, in contrast to density-functional theory where the electron density is fixed along the adiabatic path, (Langreth and Perdew, 1979; Gunnarsson and Lundqvist, 1976; Zhang and Burke, 2004), at the BSE@GW level, the density is not maintained as λ varies. Therefore, an additional contribution to Eq. 15 originating from the variation of the Green’s function along the adiabatic connection should, in principle, be added. However, as commonly done within RPA (Toulouse et al., 2009, 2010; Angyan et al., 2011; Colonna et al., 2014) and BSE, (Holzer et al., 2018; Loos et al., 2020), we neglect this additional contribution.

In the following we will illustrate how the starting point can influence the resulting quasiparticle energies. The Green’s function obtained from the one-shot G ₀ W ₀ does not satisfy particle-hole symmetry, the fundamental gap is underestimated (top right panel of Figure 2) yet more accurate than G _HF W _HF (top left panel of Figure 2), the quasiparticle intensity relative to the bonding component is close to the exact result up to U/t ≈ 16 (center right panel of Figure 2), while overestimated for the antibonding components. Moreover, we note that the intensities of the two poles of the bonding component crosses at U/t = 24. This means that if we sort the quasiparticle and the satellite according to their intensity at a given U/t, the nature of the two poles is interchanged when one increases U/t, which results in a discontinuity in the QP energy. Meanwhile, the total number of particle is not conserved (N < 2). For G ₀ W ₀ we found a small deviation from N = 2 for small U/t (e.g. N = 1.98828 at U = 1), which becomes larger by increasing the interaction (e.g. N = 1.55485 for U/t = 10). Instead, starting from G _HF the particle number is always conserved. We checked that for the self-consistent calculations the total particle number is conserved, as it should.

Considering G ₀ as starting point in evGW, we encounter the problem described in Ref. (Véril et al., 2018), namely the discontinuity of various key properties (such as the fundamental gap in the top right panel of Figure 2) with respect to the interaction strength U/t. This issue is solved, for the Hubbard dimer, by considering a better starting point or using the fully self-consistent scheme scGW. Note, however, that improving the starting point does not always cure the discontinuity problem as this issue stems from the quasiparticle approximation itself. Full self-consistency, instead, avoids systematically discontinuities since no distinction is made between quasiparticle and satellites. Unfortunately, full self-consistency is much more involved from a computational point of view and, moreover, it does not give an overall improvement of the various properties of interest, at least for the Hubbard dimer, for which G _HF W _HF is to be preferred. For more realistic (molecular) systems, it was shown in Ref. (Berger et al., 2020). that the computationally cheaper self-consistent COHSEX scheme solves the problem of multiple quasiparticle solutions.

In Figure 3, we report the real part of the singlet and triplet excitation energies obtained from the solution of Eq. 8 for λ = 1. For comparison, we report also the exact excitation energies obtained as differences of the excited- and ground-state total energies of the Hubbard dimer obtained by diagonalizing the Hamiltonian (18) in the Slater determinant basis { 1 ↑ , 1 ↓ , 1 ↑ , 2 ↓ , 1 ↓ , 2 ↑ , 2 ↑ , 2 ↓ } built from the sites [Ref. (Romaniello et al., 2009a) for the exact total energies]. For the singlet manifold, this yields, for the single excitation Ω 1 ↑ ↓ and double excitation Ω 2 ↑ ↓ , the following expressions: Ω 1 ↑ ↓ = 1 2 U + 16 t 2 + U 2 , Ω 2 ↑ ↓ = 16 t 2 + U 2 , (27) while the unique triplet transition energy is Ω 1 ↑ ↑ = 1 2 − U + 16 t 2 + U 2 . (28)

Of course, one cannot access the double excitation within the static approximation of BSE, (Strinati, 1988; Romaniello et al., 2009b; Loos and Blase, 2020), so only the lowest singlet and triplet excitations, Ω 1 ↑ ↓ and Ω 1 ↑ ↑ , are studied below.

Using one-shot G _HF W _HF quasiparticle energies (BSE@G _HF W _HF) produces complex excitation energies (see right panel of Figure 3). We find the same scenario also with other flavors of GW (not reported in the figure), such as scGW. The occurrence of complex poles and singlet/triplet instabilities at the BSE level are well documented (Holzer et al., 2018; Blase et al., 2020; Loos et al., 2020) and is not specific to the Hubbard dimer. For example, one finds complex poles also for H₂ along its dissociation path, (Li and Olevano, 2021), but also for larger diatomic molecules. (Loos et al., 2020). For U/t > 12.4794, the singlet energy becomes pure imaginary, the same is observed for the triplet energy for 7.3524 < U/t < 12.4794. These two points corresponds to discontinuities in the first derivative of the excitation energies with respect to U/t (Figure 3). The BSE excitation energies are good approximations to their exact analogs only for U/t ≲ 2 for the singlet and U/t ≲ 6 for the triplet. Using exact quasiparticle energies instead produces real excitation energies, with the singlet energy in very good agreement with the exact result; the triplet energy, instead, largely overestimates the exact value. This seems to suggest that complex poles are caused by the approximate nature of the GW quasiparticle energies, although, of course, the quality of the kernel also plays a role. Indeed, setting W = 0 but using GW QP energies, BSE yields real-valued excitation energies. It would be interesting to further investigate this issue by using the exact kernel together with GW QP energies. This is left for future work.

For the Hubbard dimer, we have E ^HF = − 2t + U/2, and the correlation energy given in Eq. 15 can be calculated analytically. After a lengthy but simple derivation, one gets E c AC@BSE = − U 2 + t 2 − 2 U 2 2 U ( 2 t + 3 U ) { Δ ϵ G W − 1 2 ( t + U ) [ − U 2 + 2 ( t + U ) Δ ϵ G W ] U ( 2 t + 3 U ) + 2 ( t + U ) Δ ϵ G W } − t + 2 U 2 U ( 2 t + 3 U ) 3 t + 4 U 2 t + 3 U + t U Δ ϵ G W a t a n − U U ( 2 t + 3 U ) 2 Δ ϵ G W ( t + U ) + [ − U 2 + 2 ( t + U ) Δ ϵ G W ] [ U ( 2 t + 3 U ) + 2 ( t + U ) Δ ϵ G W ] .

Results are reported in Figure 4 and are compared with the exact correlation energy (Romaniello et al., 2009a) E c = − 16 t 2 + U 2 2 + 2 t . (29)

The AC@BSE correlation energy does not possess the correct asymptotic behavior for small U, as Taylor expanding Eq. 29 for small U, we obtain E c A C @ B S E = − U 2 32 t − 5 U 3 96 t 2 + 323 U 4 6144 t 3 + O ( U 4 ) , (30) while the exact correlation energy behaves as E c = − U 2 16 t + U 4 1024 t 3 + O ( U 6 ) . (31)

Moreover, we found that the radius of convergence of the small-U/t expansion of E c A C @ B S E is very small due to a square-root branch point for U/t ≈ − 2/3.

In the case of the trace formula Eq. 14, the singlet and triplet contributions behave as E c , ↑ ↓ Tr@BSE = − U 2 32 t − 7 U 3 128 t 2 + 99 U 4 2048 t 3 + O ( U 5 ) , (32a) E c , ↑ ↑ Tr@BSE = − U 2 32 t + 7 U 3 128 t 2 − 157 U 4 2048 t 3 + O ( U 5 ) , (32b) which guarantees the correct asymptotic behavior for the total Tr@BSE correlation energy E c Tr@BSE = − U 2 16 t − 29 U 4 1024 t 3 + O ( U 5 ) , (33) and cancels the cubic term (as it should).

The trace formula is strongly affected by the appearance of the imaginary excitation energies: as shown in Figure 4 where we plot the real and complex components of the BSE@G _HF W _HF correlation energy as functions of U/t at various levels of theory, irregularities (i.e., discontinuities in the first derivative of the energy) appear at the values of U/t for which the triplet and singlet energies become purely imaginary. The ACFDT expression, instead, is more stable over the range of U/t considered here with only a small cusp on the energy surface at the singlet instability point after which the real part of E _c ^AC@BSE behaves linearly with respect to U/t. Overall, however, the correlation energy obtained by the trace formula is almost on top of its exact counterpart over a wide range of U/t, with a rather small contribution from the triplet component, i.e., E c , ↑ ↑ Tr@BSE ≪ E c , ↑ ↓ Tr@BSE . For comparison purposes, the RPA correlation energy, which is obtained from the trace or ACDFT formula using BSE@G _HF W _HF with W = 0 in the BSE kernel, is also reported in Figure 4. Both formulas yield the same correlation energies as expected, and they show no irregularities thanks to the fact that BSE excitation energies are real-valued at the RPA level. Also correlation energies obtained using BSE@exact (also shown in Figure 4) do not show irregularities for the same reason. Moreover, they show a visible upshift with respect to the corresponding AC@BSE@G _HF W _HF and Tr@BSE@G _HF W _HF results, which worsens the agreement with the exact correlation energy. Finally, we observe that both expressions for the correlation energy (at BSE@GW level) produce better results than the Galitskii-Migdal Eq. 4, as one can see from Figure 4, in particular at large U/t.