Static and dynamic Bethe–Salpeter equations in the T-matrix approximation

The GW approximation¹ of many-body perturbation theory² is becoming a method of choice to target charged excitations (i.e., ionization potentials and electron affinities) in molecular systems.^3–7 These so-called quasiparticle energies can be experimentally measured from direct and inverse photoemission spectroscopies. From a more theoretical point of view, GW corresponds to an elegant resummation of all direct ring diagrams from the particle–hole (ph) channel, which is particularly justified in the high-density or weakly correlated regime.^8,9 Within the GW approximation, the self-energy—one of the key quantities of Hedin’s equations¹—reads as

where G is the one-body Green’s function, W is the dynamically screened Coulomb potential, and, e.g., 1 ≡ (σ₁, r₁, t₁) is a composite coordinate gathering spin, space, and time variables.

Alternatives to GW do exist. For example, the T-matrix (or Bethe–Goldstone) approximation, first introduced in nuclear physics,^10–14 then in condensed matter physics,^15–26,119 and more recently in quantum chemistry,^27,28 sums to infinity the ladder diagrams from the particle–particle (pp) channel and is justified in the low-density or strongly correlated regime.^13–15,29 While the two-point screened interaction W is the cornerstone of GW, the T-matrix approximation relies on a more complex (four-point) effective interaction—the so-called T matrix—yielding the following self-energy:

The natural idea of combining the ph and pp channels is also possible and has been explored, for example, in the Hubbard dimer within many-body perturbation theory (see Ref. 22 and references therein) and the uniform electron gas³⁰ within coupled-cluster theory.²⁹ 

One of the key features of the T-matrix approximation is its exactness up to the second order, thanks to the inclusion of second-order exchange diagrams. This class of diagrams, which are particularly important in few-electron molecular systems^31–34 and explain the improvement brought by the second-order screened exchange (SOSEX) correction applied to GW,^35–37 is well known to be missing in the GW approximation. Moreover, unlike W in the GW approximation, the T-matrix approximation also contains spin–flip terms; the spin structure of the T-matrix allows one to describe important processes such as the emission of spin waves in ferromagnetics.³⁸ 

In this work, we focus on neutral excitations and we explore how the T-matrix approximation performs within the Bethe–Salpeter equation (BSE) of many-body perturbation theory.^39–42 

Let us consider closed-shell electronic systems consisting of N electrons and K one-electron basis functions. The number of singly occupied and virtual (i.e., unoccupied) spinorbitals is O = N and V = K − O, respectively. Let us denote the pth spinorbital as ψ_p(x) and its one-electron energy as $εp$⁠. The composite variable x = (σ, r) gathers spin (σ) and spatial (r) variables. We assume real quantities throughout this paper: i, j, k, and l are occupied orbitals; a, b, c, and d are unoccupied orbitals; p, q, r, and s indicate arbitrary orbitals; and m labels single excitations, while n labels double electron attachments or double electron detachments.

By definition, in the quasiparticle approximation, the one-body Green’s function is² 

where η is a positive infinitesimal and its nature is completely defined by the set of orbitals and the corresponding energies that are used to build it. For example, G^HF(x₁, x₂;  ω) is the Hartree–Fock (HF) Green’s function built from HF spinorbitals $ψpHF(x)$ and energies $εpHF$⁠.

Contrary to the GW approximation that relies on the (two-point) dynamically screened Coulomb potential W computed from a ph-random-phase approximation (ph-RPA) problem to target charged excitations,^1–6 here we consider the GT approximation where one employs the (four-point) T matrix obtained from solving the pp-RPA equations.

The non-Hermitian pp-RPA problem reads^43–53 

where the elements of various matrices are defined as

and

are two-electron integrals in the spinorbital basis, i.e.,

In the absence of instabilities (which should not appear in Coulombic systems with repulsive interactions only⁴⁶), the pp-RPA problem yields V(V − 1)/2 positive eigenvalues $ΩnN+2$ and O(O − 1)/2 negative eigenvalues $ΩnN−2$⁠, which correspond to double attachments and double detachments, respectively. The pp-RPA correlation energy is given by^45,46

Considering the time structure of the T-matrix approximation as $T(1,3;2,4)=−δ(t1−t3)δ(t2−t4)T(x1,x2;x4,x3;t1−t2)$⁠,² the frequency-dependent T-matrix self-energy can be obtained from the Fourier transform of Eq. (2) as

The correlation part of the T matrix can be constructed from the knowledge of the pp-RPA eigenvalues and eigenvectors. In the spinorbital basis, it is defined as $Tpq,rsc=Tpq,rs−pq‖rs$ and it has the following form:⁵² 

with

Combining Eqs. (3) and (10), the correlation part of the T-matrix self-energy reads^2,22,27,28

While the dynamical GW self-energy corresponds to the downfolding of the 2h1p and 2p1h configurations on the 1h and 1p configurations via their coupling with the 1h1p configurations, respectively,⁵⁴ Eq. (12) shows that, in the case of the T-matrix approximation, the same 2h1p and 2p1h configurations are downfolded on the 1p and 1h configurations via their coupling with the 2h and 2p configurations, respectively.

Within the (perturbative) one-shot GT scheme (labeled G₀T₀ in the following), the quasiparticle energies are obtained via the linearization of the quasiparticle equation,^55–63 i.e.,

where we have assumed a HF starting point, and

is the renormalization factor or weight of the quasiparticle solution. Other levels of (partial) self-consistency can be considered like the “eigenvalue” self-consistent GT (evGT)^57,63–67 or the quasiparticle self-consistent GT (qsGT)^68–72 schemes.

As the one-body Green’s function is the pillar of the GW and GT approximations, similarly the two-body Green’s function G₂ is the central quantity of the BSE formalism of many-body perturbation theory^39–42 via its link with the two-body correlation function L, which satisfies the following Dyson equation:

where

and

is the so-called BSE kernel that takes into account the variation of Σ with respect to the variation of G. By taking into account the interaction of the excited electron and its hole left behind (the infamous excitonic effect), the BSE is able to reliably model (neutral) optical excitations as measured by absorption spectroscopy. The moderate cost of the BSE [which scales as $O(K4)$ in its standard implementation] and its all-round accuracy are the main reasons behind its growing popularity in the molecular electronic structure community.^{41,42,67,73–91}

In order to target neutral (singly) excited states, we first consider the static version of the BSE employing the GT quasiparticle energies [see Eq. (13)] as well as the T-matrix kernel [i.e., Σ^GT in Eq. (17)]. In this case, the BSE@GT linear eigenvalue problem simply reads

with

The eigenvalues $ΩmBSE$ of Eq. (18) provide OV singlet (i.e., spin-conserved) and OV triplet (i.e., spin–flip) single excitations. Note that the spin structure of the BSE@GT equations is analogous to the BSE@GW version,⁹² and one can separately compute singlet and triplet excitation energies. Neglecting the coupling between excitations and de-excitations, i.e., $BBSE=0$⁠, is known as the Tamm–Dancoff approximation (TDA).

Due to the frequency-independent nature of the static BSE, it is well known that one cannot access double (and higher) excitations.^92–97 In order to go beyond the static approximation, it is possible to consider, within the dynamical TDA (dTDA) that neglects the frequency dependence of the coupling block B, the dynamical version of the BSE (dBSE).^40,94,96 In this case, one must solve the (non-linear) dynamical eigenvalue problem,

with

where, by following Strinati’s seminal work,⁴⁰ one can derive the following expression for the elements of the dynamical T matrix:

from which one can check that we recover the static expression (10) in the limit $ΩnN±2→∞$⁠. Equation (22) highlights the interesting dynamical structure of the T matrix, where, similar to the dBSE@GW scheme,^40,94,96 the 2h2p configurations are downfolded on the 1h1p configurations.⁹⁸ Additional details about the derivation of Eq. (22) are reported in the  Appendix.

Because solving a non-linear eigenvalue problem is computationally challenging, here we rely on the perturbative scheme developed in Ref. 96 in order to access dynamically corrected single excitations for which additional relaxation effects coming from higher excitations are taken into account.^{92,94,96,97,99–104} Below, we quickly recap this dynamical perturbative scheme.

Based on Rayleigh–Schrödinger perturbation theory, the non-linear eigenproblem (20) can be split as a zeroth-order static reference and a first-order dynamic perturbation such that

with

As usual, one can naturally expand the Sth BSE excitation energy and its corresponding eigenvector as

Solving the static BSE [see Eq. (18)] yields the (zeroth-order) static $ΩSBSE$ excitation energies and their corresponding eigenvectors $XSBSE$ and $YSBSE$⁠. The first-order correction to the Sth excitation energy is, within the dTDA,

This correction can be renormalized by computing, at no extra cost, the renormalization factor, which reads

This yields our final expression for the dynamically corrected BSE excitation energies,

Note again that the present perturbative scheme does not allow one to access double excitations as only excitations calculated within the static approach can be dynamically corrected.

The present formalism has been implemented in the electronic structure package QuAcK,¹⁰⁵ which is freely available at https://github.com/pfloos/QuAcK. We consider here only systems with closed-shell singlet ground states. Thus, the GW and GT calculations are performed by considering a (restricted) HF starting point and standard Gaussian basis sets (defined with cartesian functions) are employed. Note that all quasiparticle energies, which are obtained via Eq. (13), are corrected in the same way. Finally, the infinitesimal η is set to zero for all calculations. The evGT and qsGT schemes have been also implemented but are not considered here, mainly because, for the small molecular systems studied here (see below), we have observed very small differences between one-shot and self-consistent quasiparticle energies. Although the dynamical correction is computed throughout in the dTDA, the zeroth-order excitonic Hamiltonian [see Eq. (18)] is always the “full” BSE static Hamiltonian, i.e., without TDA. Reference full configuration interaction (FCI) calculations have been performed using quantum package.¹⁰⁶ 

In terms of computational cost, the overall scaling of BSE@GT is equivalent to BSE@GW as they both correspond to seeking the lowest eigenvalues of a matrix of size (2OV × 2OV). Searching iteratively for the lowest eigenstates can be routinely performed via Davidson’s algorithm with a $O(K4)$ computational cost.¹⁰⁷ The cost of the dynamical correction, which is thoroughly discussed in Ref. 96, is more expensive but is again equivalent in both formalisms. Both the computational cost associated with the computation of the T-matrix and screening W scale as $O(K6)$ in their standard implementation as one must obtain all the eigenvalues and eigenvectors of the pp-RPA and ph-RPA problems, respectively.¹⁰⁸ However, the prefactor of the pp-RPA calculation is significantly larger than its ph-RPA counterpart due to the larger size of the pp-RPA matrices and its non-Hermitian nature.^{45–48,50–52} Moreover, within the T-matrix formalism, one must compute both the singlet and triplet contributions of the T-matrix, while for singlet states, only the singlet part of W is required. Although similar approaches remain to be developed for the T-matrix formalism, contour deformation and density fitting techniques can be efficiently implemented in the case of GW to reduce the scaling to $O(K3)$⁠.^109–111 

As a first didactical example, we consider the lowest singlet and triplet excited states of the hydrogen molecule H₂ and the variation of the respective vertical transition energies upon dissociation. The excitation energies associated with these low-lying excited states are represented in Fig. 1 as the function of the internuclear distance R_H–H at the FCI (black), BSE@G₀W₀ (blue), and BSE@G₀T₀ (red) levels with the cc-pVTZ basis. The variation of the HOMO and LUMO quasiparticle energies, as well as the HOMO–LUMO gap computed at the G₀W₀ and G₀T₀ levels, is depicted in Fig. 2. This shows that, as already observed in the Hubbard dimer²² and in molecular systems,²⁷ the G₀W₀ and G₀T₀ quasiparticle energies are similar near the Fermi level.

Overall, as evidenced by Fig. 1, the performance of BSE@G₀W₀ and BSE@G₀T₀ is analogous for this system. For the lowest singlet excited state of $B1Σu+$ symmetry, the T-matrix-based formalism is slightly better when R_H–H increases but ultimately fails to reproduce the FCI results. For the $E1Σg+$ state, BSE@G₀T₀ is more accurate than BSE@G₀W₀ for small bond length and the scenario is reversed after the avoided crossing with the doubly excited state of $F1Σg+$ symmetry. Of course, both formalisms cannot “see” the $F1Σg+$ states as the static BSE formalism is blind to double excitations. Therefore, it cannot model properly the avoided crossing between $E1Σg+$ and $F1Σg+$ states. For the $B’1Σu+$ and C¹Π_u states, BSE@G₀W₀ and BSE@G₀T₀ reproduces fairly well the FCI potential energy curves with a modest preference for the latter.

Similar observations can be made for the triplet states, the GW- and GT-based formalisms yielding very similar excitation energies, except for the c³Π_u state for which BSE@G₀W₀ has clearly the edge. Moreover, triplet instabilities seem to affect BSE@G₀T₀ slightly earlier than BSE@G₀W₀.

In Fig. 3, we show the energy shift provided by the dynamical correction for the lowest singlet and lowest triplet excited states of H₂ as the function of R_H–H. These dynamically corrected schemes are labeled dBSE@G₀W₀ and dBSE@G₀T₀. For the singlet state of $B1Σu+$ symmetry, the dynamical correction slightly improves the excitation energies at small internuclear distances for both schemes, while, for larger bond lengths, an improvement is only visible at the T-matrix level. Note that, for this system with few electrons, the dynamical corrections are quite small in magnitude. In the case of the triplet state of $b3Σu+$ symmetry, the dynamical correction worsens the results compared to FCI, especially in the case of BSE@G₀W₀.

As a second example, we consider the symmetric dissociation of linear molecule beryllium hydride (BeH₂), a system for which one can assume that screening plays a more important role than in the previous example. The variation of the lowest singlet and triplet excitation energies as the function of the distance R_Be–H is shown in Fig. 4, while the quasiparticle energies of the frontier orbitals and the associated (fundamental) gap computed at the G₀W₀ and G₀T₀ levels are depicted in Fig. 5. All calculations are performed with the cc-pVDZ basis. Again, one notes that the G₀W₀ and G₀T₀ quasiparticle energies are very similar near the Fermi level. Therefore, one can safely assume that any significant variation of the excitation energies computed within the GW- and GT-based formalisms originates mainly from their distinct kernel. The excitation energies computed with the dynamical schemes, dBSE@G₀W₀ and dBSE@G₀T₀, are reported as thin solid lines. Here, one can show that dynamical corrections improve agreement between BSE and FCI in most cases.

For the four lowest singlet excited states (left panel of Fig. 4), dBSE@G₀T₀ is clearly better than dBSE@G₀W₀, while the opposite trend is observed for the four lowest triplet states (right panel of Fig. 4). Note that, for large R_Be–H, the two BSE-based schemes provide only a qualitative description of the excited states with errors of several eVs. Nonetheless, the overall ordering of the excited states is globally respected.

As a third and final example, we compute the excitation energies associated with the two lowest singlet and two lowest triplet excited states of water at equilibrium geometry (see Fig. 6). Note that all these excited states are of Rydberg nature and correspond to n → 3s and n → 3p transitions for the B₁ and A₂ states, respectively.¹¹³ In addition to the BSE-based models studied in the present paper, we have selected well-known wave function methods,^114–117 namely, configuration interaction with singles (CIS), CIS with perturbative doubles [CIS(D)], time-dependent Hartree–Fock (TDHF), and FCI (taken as reference), and computed the excitation energies of these transitions. It is worth mentioning here that the TDHF (or RPAx) equations within the TDA are strictly equivalent to the CIS equation¹¹⁷ and that CIS(D) is a simple perturbative double correction to CIS and can be considered as an excited-state analog of second-order Møller–Plesset perturbation theory.^115,116

Two key observations can be made: (i) BSE@G₀W₀ is by far the best performer with a slight overestimation on the order of 0.1 eV (as compared to FCI) and (ii) BSE@G₀T₀ systematically underestimates the excitation energies [similarly to CIS(D)] and outperforms CIS, CIS(D), and TDHF for the singlet states only. These general trends are also observed for other systems, and they well evidence the crucial role of screening in GW, hinting that a screened version of the T-matrix formalism as proposed in Ref. 22 might be a promising way for improvement.

We have derived and implemented, for the first time, the static and dynamic Bethe–Salpeter equations when one considers T-matrix quasiparticle energies and a T-matrix-based kernel. The performance of the static scheme and its perturbative dynamical correction have been assessed by computing the neutral excited states of several molecular systems. Our results suggest that, in the context of the computation of molecular excitation energies, the BSE@GT formalism performs best in few-electron systems where the electron density remains low. The overall accuracy of the present scheme still needs to be assessed for larger systems (where screening is known to be more important). For such purposes, a comprehensive benchmark study would be required and we are planning to do so in the future.

It would be interesting to investigate its performance for the computation of ground-state correlation energies within the adiabatic connection fluctuation dissipation formalism, where BSE@GW has been shown to be particularly outstanding.^89,90,118 The combination of GT and GW via the range separation of the Coulomb operator to avoid double counting of the low-order diagrams is also a promising avenue. The work along these lines is currently in progress. Finally, the unrestricted and spin–flip extensions of the present formalism are currently being developed.

See the supplementary material for the raw data of each figure.

The authors would like to thank Roberto Orlando for useful discussions. P.-F.L. thanks the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. 863481) for funding. This project has also received financial support from the Centre National de la Recherche Scientifique (CNRS) through the 80|PRIME program and has been supported through the EUR grant NanoX under Grant No. ANR-17-EURE-0009 in the framework of the “Programme des Investissements d’Avenir.”

The authors have no conflicts to disclose.

The data that support the findings of this study are available within the article and its supplementary material.

In order to derive the dynamical kernel $T̃ib,ajc$ given in Eq. (22), we follow Ref. 55 (see also Ref. 90) and start from the equation for the BSE amplitude,

where L₀ is given by Eq. (16a) and the T-matrix kernel is

Equation (A1) is derived by assuming that (i) the (resonant) pole ω_S = E_S − E₀ > 0 of L is isolated from the other poles (which is usually the case for neutral excitations in finite systems) and (ii) the poles of L₀ are different from ω_S (which is also generally the case). The so-called T-matrix self-energy Σ is given by Eq. (2) with^2,22

where v is the bare Coulomb operator, and we neglect the functional derivative δT/δG in the kernel Ξ. The first two terms in the right-hand side of Eq. (A3) are the Hartree and exchange contributions to the T-matrix, whereas the last term is the correlation contribution. Making the time dependence of Eq. (A1) explicit and defining $T(3,5,4,6)=−δ(τ35+)δ(τ64+)T(x3,x5,x4,x6;τ34)$⁠, one obtains

[where τ_ij = t_i − t_j and $τij+=ti+−tj$ with $ti+=ti+η$ (η → 0⁺)] and

Using the Fourier transform $G(τ)=∫dω2πG(ω)e−iωτ$⁠, changing the variable from t₃ to τ₃₄ and taking the limit $t1′=t1+$⁠, we have

Using the Lehman representation of the one-body Green’s function in the quasiparticle approximation given by Eq. (3), and multiply the left- and right-hand sides by $(εa−εi−ωS)∫dx1dx1′ψa(x1)ψi(x1′)$⁠, we obtain

(where Θ is the Heaviside step function) using the fact that

For the resonant case, i.e., ω_S > 0, we have

where $ĉp†$ and $ĉp$ are the usual creation and annihilation operators, respectively, and $N$ and $N,S$ are the ground state and the Sth excited state, respectively, of the N-electron system. After some algebraic steps, one obtains

where we have defined

Using the definition $Xia,S=N|ĉi†ĉa|N,S$⁠, we arrive at

where the spectral representation of the dynamical T-matrix is

with $Tib,ajc=Tib,aj−ib‖aj$ being the correlation part of $T$⁠. Equation (A12) represents a non-linear eigenvalue equation to calculate the positive excitation energies of a system, which can be rewritten as

with

If one drops the dynamical part $T̃c$⁠, one ends up with the usual time-dependent Hartree–Fock (TDHF) equations.¹¹⁷ To calculate the correlation contribution, one can employ Eq. (10) in Eq. (A13), and, after integration over the frequency, one obtains Eq. (22).