A formally exact one-frequency-only Bethe-Salpeter-like equation. Similarities and differences between GW +BSE and self-consistent RPA

TL;DR

This paper introduces a new, exact Bethe-Salpeter-like equation that depends on a single frequency, simplifying the response function calculation and enabling more advanced approximations, demonstrated through the Hubbard molecule.

Contribution

A formally exact one-frequency-only Bethe-Salpeter-like equation is developed, contrasting with the standard multi-frequency BSE, and its advantages are illustrated with analytical solutions.

Findings

Exact analytical solution for the Hubbard molecule.

Simplification of the response function calculation.

Discussion of similarities and differences between GW+BSE and self-consistent RPA.

Abstract

A formally exact Bethe-Salpeter-like equation for the linear-response function is introduced with a kernel which depends only on the one frequency of the applied field. This is in contrast with the standard Bethe-Salpeter equation (BSE) which involves multiple-frequency integrals over the kernel and response functions. From the one-frequency kernel, known approximations are straightforwardly recovered. However, the present formalism lends itself to more powerful approximations. This is demonstrated with the exact analytical solution of the Hubbard molecule. Similarities and differences of the $GW$+BSE approach with the self-consistent random-phase approximation (RPA) is also discussed.

Equations

$$H=\sum _{k_1{k}_2}{ϵ}_{k_1{k}_2}{c}_{k_1}^{†}{c}_{k_2}+\frac{1}{4}\sum _{k_1{k}_2{k}_3{k}_4}{\bar{v}}_{k_1{k}_2{k}_3{k}_4}{c}_{k_1}^{†}{c}_{k_2}^{†}{c}_{k_4}{c}_{k_3},$$

$${\bar{v}}_{k_1{k}_2{k}_3{k}_4}=⟨k_1{k}_2|v|k_3{k}_4⟩-⟨k_1{k}_2|v|k_4{k}_3⟩,$$

$$⟨k_1{k}_2|v|k_3{k}_4⟩=\int drd{r}^{\prime }{\varphi }_{k_1}^{\ast }(r){\varphi }_{k_2}^{\ast }(r^{\prime })v(r,r^{\prime }){\varphi }_{k_3}(r){\varphi }_{k_4}(r^{\prime }).$$

$$G_{k{k}^{\prime }}(t-t^{\prime })=-i⟨0|\mathrm{T}\{c_k(t)c_{k^{\prime }}^{†}(t^{\prime })\}|0⟩,$$

$$G_{k{k}^{\prime }}^{0^{-1}}(t-t^{\prime })=\delta (t-t^{\prime })({\delta }_{k{k}^{\prime }}i{\partial }_t-{ϵ}_{k{k}^{\prime }}),$$

$$\sum _{k_1}\int d{t}_1$$

$$j_k=[c_k,V]=\frac{1}{2}\sum _{k_2{k}_3{k}_4}{\bar{v}}_{k{k}_2{k}_3{k}_4}{c}_{k_2}^{†}{c}_{k_4}{c}_{k_3}.$$

$$(i{\partial }_t-{ϵ}_k)G_k(t-t^{\prime })=\delta (t-t^{\prime })+\int d{t}_1{\Sigma }_k(t-t_1)G_k(t_1-t^{\prime }).$$

$${\Sigma }_k(t-t^{\prime })=-i\int d{t}_1⟨0|\mathrm{T}\{j_k(t)c_k^{†}(t_1)\}|0⟩G_k^{-1}(t_1-t^{\prime }),$$

$$G_k^{-1}(t-t^{\prime })=G_k^{0^{-1}}(t-t^{\prime })-{\Sigma }_k(t-t^{\prime }).$$

$${\Sigma }_k(t-t^{\prime })$$

$${\Sigma }_k(t-t^{\prime })=T_k(t-t^{\prime })-C_k^{\mathrm{red}}(t-t^{\prime }),$$

$$T_k(t-t^{\prime })=V_k^{\mathrm{MF}}\delta (t-t^{\prime })-i⟨0|\mathrm{T}\{j_k(t)j_k^{†}(t^{\prime })\}|0⟩,$$

$$C_k^{\mathrm{red}}(t-t^{\prime })$$

$$V_k^{\mathrm{MF}}=⟨0|\{[c_k,V],c_k^{†}\}|0⟩=\sum _{k^{\prime }}{\bar{v}}_{k{k}^{\prime }k{k}^{\prime }}{n}_{k^{\prime }},$$

$$n_k=⟨0|c_k^{†}{c}_k|0⟩$$

$${\Sigma }_k(t-t^{\prime })=V_k^{\mathrm{MF}}\delta (t-t^{\prime })-i⟨0|\mathrm{T}\{j_k(t)j_k^{†}(t^{\prime })\}|0{⟩}^{\mathrm{irr}},$$

$$G_k=G_k^0+G_k^0{T}_k{G}_k^0\mathrm{with}{T}_k={\Sigma }_k+{\Sigma }_k{G}_k^0{T}_k.$$

$$R_{k_1{k}_2{k}_1^{\prime }{k}_2^{\prime }}(t-t^{\prime })=-i⟨0|\mathrm{T}\{c_{k_2}^{†}(t)c_{k_1}(t)c_{k_1^{\prime }}^{†}(t^{\prime })c_{k_2^{\prime }}(t^{\prime })\}|0⟩$$

$$\int$$

$${\tilde{R}}_{k_1{k}_2}^{0^{-1}}(t-t^{\prime })=\delta (t-t^{\prime })(i{\partial }_t-{ϵ}_{k_1}+{ϵ}_{k_2}),$$

$$J_{k_1{k}_2}$$

$$N_{k_1{k}_2{k}_1^{\prime }{k}_2^{\prime }}$$

$$N_{k_1{k}_2}=n_{k_2}-n_{k_1}=|n_{k_2}-n_{k_1}|N_{k_1{k}_2}^0,$$

$$N_{k_1{k}_2}^0=1\mathrm{for}{k}_1>k_2\mathrm{and}-1\mathrm{for}{k}_1<k_2,$$

$${\tilde{R}}^{0^{-1}}R=N+F,$$

$${\tilde{R}}^{0^{-1}}\tilde{R}=1+\tilde{F}=1+\tilde{F}{\tilde{R}}^{-1}\tilde{R}\equiv 1+K\tilde{R},$$

$${\tilde{R}}^{-1}={\tilde{R}}^{0^{-1}}-K,$$

$$K=\tilde{F}{\tilde{R}}^{-1}=\tilde{F}[{\tilde{R}}^{0^{-1}}-K].$$

$$\int d{t}_1{\tilde{R}}_{k_1{k}_2}^{0^{-1}}(t-t_1){\tilde{R}}_{k_1{k}_2{k}_1^{\prime }{k}_2^{\prime }}(t_1-t^{\prime })={\delta }_{k_1{k}_1^{\prime }}{\delta }_{k_2{k}_2^{\prime }}\delta (t-t^{\prime })$$

Full text

I Introduction

The development of efficient many-body approaches is an active research field in quantum chemistry and various branches of physics, such as condensed-matter, cold-atom, atomic, molecular, and nuclear physics. Originally developed in the framework of subnuclear and nuclear physics to describe bound states of systems of two interacting particles like the deuteron, the Bethe-Salpeter equation (BSE) Salpeter and Bethe (1951) has become an approach commonly used also in solid-state and condensed-matter physics Hanke and Sham (1974, 1975, 1979); Strinati (1982); Onida et al. (1995); Albrecht et al. (1998); Rohlfing and Louie (1998); Benedict et al. (1998); Olevano and Reining (2001), atomic physics Li et al. (2017), and quantum chemistry Jacquemin et al. (2015); Bruneval et al. (2015); Baumeier et al. (2012).

The fact that the standard BSE can be demonstrated Casalbuoni et al. (2010) to be equivalent to the Ward identities and the Hedin integral equation for the vertex Hedin (1965) enables a natural transfer of approximations, i.e. the Hedin $GW$ approximation Hedin (1965) on the self-energy toward the BSE kernel. The idea behind the $GW$ approach to tackle correlations simply by the introduction of only screening, i.e. the simple replacement of the bare two-body interaction $v$ by a screened interaction $W$, can be directly transferred to an approximation to the irreducible BSE kernel, which is hence written as $W$ instead of the time-dependent Hartree-Fock (TDHF) exchange kernel. In contrast with the TDHF exchange kernel, which for electronic systems is the opposite of the static Coulomb interaction $v(r,r^{\prime})=1/|r-r^{\prime}|$, a BSE kernel at the same level of the $GW$ approximation should in principle be frequency dependent since it relies on the dynamically screened Coulomb interaction $W(r,r^{\prime},\omega)$. This frequency dependence, which can be worked out, though with some difficulties, when calculating the $GW$ self-energy, and which is an important ingredient to have quasiparticle energies more in agreement with experiment, implies multiple-frequency integrals in the BSE and represented so far an insurmountable obstacle to the resolution of the full BSE. For this reason almost all BSE calculations were obliged to neglect the dynamical dependence of the BSE kernel and solve a static BSE. This approach often called $GW$+BSE Martin et al. (2016) which uses a dynamical $W(r,r^{\prime},\omega)$ in the self-energy and a static $W(r,r^{\prime},\omega=0)$ in the BSE kernel has nevertheless provided good results in agreement with experiment and exact solutions Onida et al. (1995); Albrecht et al. (1998); Rohlfing and Louie (1998); Benedict et al. (1998); Olevano and Reining (2001); Li et al. (2017).

It is difficult to estimate how important can be dynamical effects beyond the static BSE. Nevertheless, it is often conjectured that deviations of the static BSE solution from experiment can be solely due to dynamical BSE effects. A tentative list might include effects associated to double excitations in quantum chemistry Rebolini et al. (2013) or to electron-hole screened interaction in metals Marini and Del Sole (2003). Efforts to study dynamical BSE effects and introduce a real frequency dependence into the BSE have recently been attempted Romaniello et al. (2009); Sangalli et al. (2011); Zhang et al. (2013); Rebolini and Toulouse (2016). The standard BSE is an equation over two-body Green/correlation functions (kernel and response functions), i.e. functions of four space-time points. In systems with time-translation invariance, there is one-time degree-of-freedom less, which means always functions of three time differences, or their three Fourier transformed frequencies. The full dynamical BSE involves a so far numerically intractable integration over frequencies in the kernel and in the response function. Recent efforts Romaniello et al. (2009); Sangalli et al. (2011) have tried to redefine a kernel which incorporates frequency integration, to finally arrive at a more easily solvable one-frequency equation. Another approach Zhang et al. (2013); Rebolini and Toulouse (2016) has considered the coupling of the linear-response function to uncorrelated two-particle-two-hole (2p-2h) states. The coupling of the linear-response function to collective states plus free particle-hole (p-h) states to account for double excitations has been discussed in Ref. Romaniello et al. (2009).

All previous works followed the route which starts from the multi-frequency standard BSE and tries to reduce the number of involved frequencies, so as to end up with an equation with just only the one frequency of the external field. The purpose of this work is to follow a different route: we introduce from the beginning a formally exact one-frequency BSE-like equation, i.e. depending only on the frequency of the external field, for a linear-response function. In particular, this means that also the integral kernel $K$ depends only on the one frequency of the external field. Explicit expressions for $K$ will be elaborated in terms of well defined correlation functions and higher Green functions. For readers interested right away to see the final result, they may consult Eqs. (19), (LABEL:kernel2-total), (LABEL:kernel2-static), (LABEL:kernel2), and (37). Starting from these expressions we then rederive the approximate expressions given in the literature mentioned above. However, since our expressions are more general, they lend themselves to more far-reaching approximations without loosing the advantage of a one-frequency only approach. This is demonstrated with the exact solution of the Hubbard molecule. But we will also point out that the response function calculated in this way keeps all desirable qualities of the standard random-phase approximation (RPA), such as fulfillment of sum rules and conservation laws. Our derivation is based on the equation-of-motion (EOM) technique applied to an appropriately defined four-point one-frequency linear-response function.

The paper is organized as follows. We will use the example of the EOM technique for establishing the Dyson equation for the one-body Green function presented in Sec. II to introduce the key points of the derivation of the one-frequency-only BSE-like equation which will then be presented in Sec. III. In Sec. IV we will establish the connection of the present formalism to the previous approaches of Refs. Romaniello et al. (2009); Sangalli et al. (2011); Rebolini and Toulouse (2016) and with the standard $GW$+BSE approximation, making parallels also with the self-consistent random-phase approximation (SCRPA) Dukelsky and Schuck (1990); Jemai et al. (2005); Delion et al. (2005); Schuck and Tohyama (2016). We will present in Sec. V a short application of our formalism to the Hubbard molecule which in this way can be solved exactly. Finally, Sec. VI contains our conclusions and outlook.

Atomic units are used throughout this work.

II Rederivation of the one-body Dyson equation

To set the stage, we first present a short derivation of the Dyson equation for the standard one-body Green function by the EOM technique, highlighting the points over which we will base the derivation of the one-frequency-only BSE-like equation in the next section.

We consider the most generic Hamiltonian, $H=H_{0}+V$, composed by the non-interacting (kinetic plus external potential) Hamiltonian $H_{0}$ and the two-body interaction operator $V$, which we write in terms of creation/annihilation operators $c^{\dagger}_{k}$ and $c_{k}$ on an arbitrary orthonormal basis set of orbitals $\{\phi_{k}(r)\}$ as

[TABLE]

where $\epsilon_{k_{1}k_{2}}$ are the matrix elements of the non-interacting Hamiltonian $H_{0}$ over the orthonormal basis set, and

[TABLE]

are the antisymmetrized matrix elements of the Coulomb interaction $v(r,r^{\prime})=1/|r-r^{\prime}|$, or more precisely, detailing the notation,

[TABLE]

We work in a three-dimensional space and $r$ and $k$ are meant as 3D vectors, but we can generalize to 1D and 2D; the spin degree of freedom $\sigma$ is always implied and can be included for spin-polarized cases in $k$ and in $r$ and summed over whenever $r$ is integrated out.

We remind the definition of the one-body Green function

[TABLE]

where $\mathrm{T}\{o(t)o^{\prime}(t^{\prime})\}=\theta(t-t^{\prime})o(t)o^{\prime}(t^{\prime})-\theta(t^{\prime}-t)o^{\prime}(t^{\prime})o(t)$ is the time-ordering product between fermion operators $o$, $c_{k}(t)=e^{iHt}c_{k}e^{-iHt}$ is the time-dependent annihilation operator in the Heisenberg formalism (and similarly for the time-dependent creation operator $c^{\dagger}_{k^{\prime}}(t^{\prime})$), and $|0\rangle$ is the ground state. We can then introduce the non-interacting Green function $G^{0}$, associated to the non-interacting Hamiltonian $H^{0}$, and its inverse

[TABLE]

by which we can write out a first EOM for $G$

[TABLE]

where we have introduced the operator

[TABLE]

The term containing $j_{k}$ is a two-body Green function with a particular time ordering.

For simplicity and without loss of generality, we will henceforth write the equations for the case of homogeneous systems, where $k$ stands for momentum (and spin) and $\epsilon_{k}$ is the kinetic energy. This is very similar to work in the natural spin-orbital basis, also sometimes called canonical basis, obtained from the diagonalization of the one-body density matrix in the case of inhomogeneous or finite systems. Let us now write the well-known Dyson equation Fetter and Walecka (1970)

[TABLE]

Using Eq. (6), the self-energy is then formally given by

[TABLE]

where we introduced the inverse of the Green function defined by $G^{-1}G=1$ in short-hand notation. From the Dyson equation (8), this inverse can be expressed as

[TABLE]

The self-energy can therefore be written as

[TABLE]

In this equation, $G_{k}^{0^{-1}}(t_{1}-t^{\prime})$ can be applied on the left using a second EOM (we should realize that $\partial_{t_{1}}$ contained in $G_{k}^{0^{-1}}(t_{1}-t^{\prime})$ normally acts to the right and, thus, one has to perform an integration par parts over $t_{1}$ to make it act to the left, which changes $i\partial_{t_{1}}$ into $-i\partial_{t_{1}}$), we then arrive at

[TABLE]

where

[TABLE]

which is a kind of one-body T-matrix, and

[TABLE]

The usual mean-field potential is given by

[TABLE]

where $\{..\}$ stands for the anticommutator and

[TABLE]

are the occupation numbers. We should mention that the mean-field potential in Eq. (15) is only diagonal in a homogeneous system. In a finite system, this is not necessarily the case in spite of the fact that in the natural spin-orbital (canonical) basis the density matrix is diagonal. However, to avoid heavy formulas, we will always from now on assume that the mean-field is also diagonal. It is easily recognized that the second term of the expression for the above one-body $T$-matrix is expressed as a three-body propagator of the two-particle-one-hole (2p-1h) plus two-hole-one-particle (2h-1p) type. This 3-body Green function contains the so-called one-line reducible Feynman graphs, which is easily verified by perturbation theory. By definition, a self-energy should not contain such contributions which can be “cut” into two pieces by cutting a single fermion line at a given time. It is again easily verified by perturbation theory that the second term on the right-hand side in Eq. (12) just does nothing else than taking out of the $T$-matrix all reducible terms. Therefore, in short, we can write the self-energy as

[TABLE]

where the index “irr” indicates that the corresponding correlation function should be one-line irreducible. Expression (17) is therefore a formally exact and compact expression for the self-energy. Please also note that the expression is very symmetric, which is well suited for introducing approximate forms of the self-energy. For completeness, let us also write an expression for the one-body Green function in the following way (with $G^{0^{-1}}G^{0}=1$)

[TABLE]

Notice that in Eq. (17) we have an index “irr”, so that the single particle $T$-matrix of Eq. (18) is different from $\Sigma$. Here, we did not write out the time dependencies and integrals. In frequency space there are no integrations and it becomes an algebraic equation as, by the way, the Dyson equation itself. Please note that Eq. (18) has the usual form connecting a Green function to the scattering $T$-matrix. However, here the $T$-matrix is defined for a many-body system. Taking out of $T$ the one-line reducible contributions changes $TG^{0}$ into $\Sigma G$, that is we also have the relation $T=\Sigma+\Sigma G^{0}T$ as indicated in Eq. (18).

After this hopefully pedagogic and relatively elaborate presentation of well-known many-body relations on the one-body Green function, let us now turn to the two-body case and response function.

III Response function and Bethe-Salpeter-like equation

III.1 Derivation of the one-frequency Bethe-Salpeter-like equation

We will derive for the two-time response function defined by, with $k_{1}\neq k_{2}$ and $k_{1}^{\prime}\neq k_{2}^{\prime}$,

[TABLE]

an exact equation which has the same structure as the Dyson equation for the one-body Green function derived above. The inequalities $k_{1}\neq k_{2}$ and $k_{1}^{\prime}\neq k_{2}^{\prime}$ are not independent of the one-body basis: for homogeneous matter the indices stand for momenta and spin and then the inequalities concern the momenta. For finite systems the indices correspond to the canonical basis. With this definition, we have $\langle 0|c^{{\dagger}}_{k_{2}}c_{k_{1}}|0\rangle=0$ and $\langle 0|c^{{\dagger}}_{k_{1}^{\prime}}c_{k_{2}^{\prime}}|0\rangle=0$, so that the quantity $R$ is the same as the linear-response function often denoted by $\chi$ in condensed-matter physics or quantum chemistry Rebolini et al. (2013). Let us further note that in Eq. (19) we have chosen a definite ordering of the fermion operators. This stems from the fact that we are considering a one-body-density-matrix/one-body-density-matrix correlation function. Notably there will appear an integral kernel which also depends on only one time difference or on one frequency. For people used to multi-time Green functions, this may seem surprising. However, this is not unknown in the literature Dukelsky et al. (1998). There also exists, e.g., the Mori-Zwanzig formalism for correlation functions of statistical physics Mori (1965); Hedin (1961). Furthermore, in nuclear physics, the EOM formalism developed by Rowe Hedin (1968), and further developed in Refs. Dukelsky et al. (1998); Dukelsky and Schuck (1990); Jemai et al. (2005); Delion et al. (2005); Schuck and Tohyama (2016) (with more references therein), is closely related to what we will present here. However, these facts seem to be very little known in the condensed-matter and chemical physics communities where one often struggles to get rid of eventually superfluous frequency dependencies of the integral kernel of the response function which are inherent to the so-called Hedin equations Hedin (1965). Introducing a single frequency integral kernel from the start and not a posteriori will turn out to have several advantages. For example, though we will recover, e.g., certain aspects of the $W$ kernel of the BSE as used in the $GW$ approach, we will also see more clearly what kind of approximations are involved with the use of static and dynamic forms of $W$ in the BSE and how eventually to go beyond in a systematic way.

So, let us start as before by writing down the first EOM for the response function

[TABLE]

where

[TABLE]

which is a straightforward extension of the one-body case. We have also introduced

[TABLE]

and the so-called norm kernel

[TABLE]

with

[TABLE]

where the sign factor $N^{0}$ is given by

[TABLE]

and therefore $N^{0}_{k_{1}k_{2}}N^{0}_{k_{1}k_{2}}=1$. Please note that the one-body density matrix $\langle 0|c^{\dagger}_{k_{2}}c_{k_{1}}|0\rangle$ is diagonal for our assumed homogeneous system (or in the canonical basis) and we suppose that it is also diagonal in spin. One recognizes in Eqs. (23)-(24) the phase-space factors from the standard RPA when the occupation numbers $n_{k}$ are replaced by their step function form, $n^{0}_{k}$, when using the Hartree-Fock (HF) ground state. In general, however, the occupation numbers are the correlated ones, different from zero and one. It is remarked that this norm factor is a different feature with respect to the one-body Green-function case. Note also that, contrary to the one-body case, the quantity ${\tilde{R}}^{0^{-1}}$ introduced in Eq. (21) is not exactly the inverse of the non-interacting response function $R^{0}$, but instead we have in short-hand notation ${\tilde{R}}^{0^{-1}}R^{0}=N$ where $N$ is the norm matrix.

We now proceed exactly in analogy with the one-body case. Because of the presence of the norm matrix $N$ in Eq. (20), we first have to divide it out by multiplying Eq. (20) by the inverse of $N$. Writing Eq. (20) schematically as

[TABLE]

we obtain by division with $N$

[TABLE]

with $\tilde{R}=RN^{-1}$ and $\tilde{F}=FN^{-1}$. So we arrive at a BSE-like equation of the form

[TABLE]

with the kernel $K$ given by

[TABLE]

With explicit notations, the BSE-like equation with a one-frequency kernel can thus be written as

[TABLE]

with

[TABLE]

We apply then the EOM a second time as in the one-body case and obtain the final expression of the kernel

[TABLE]

with a purely static contribution

[TABLE]

and a dynamic contribution

[TABLE]

Please note the complete analogy of this expression with Eq. (17). At this point some discussion is in order: we realize that the right-hand side of Eq. (LABEL:kernel2) corresponds to a four-body Green function of the 2p-2h and 2h-2p type. It contains therefore double p-h excitations. The index “irr” may seem less evident than in the one-body case. One may, however, verify again by perturbation theory that everything works exactly as in the one-body case and that the subtraction of the matrix $K$ in ${\tilde{R}}^{0^{-1}}-K$ exactly eliminates all p-h reducible contributions of the 2p-2h/2h-2p Green function. The p-h (two-line) irreducibility is just the analog of the one-line irreducibility in the one-body case. Up to some technical details to be discussed below, we thus have derived, as announced, a BSE-like equation with a one-frequency kernel obtained by Fourier transforming into frequency space the time dependence of the kernel in Eq. (LABEL:kernel2), which at equilibrium depends only on the time difference $t-t^{\prime}$. The frequency-space BSE-like equation that we have obtained is thus

[TABLE]

where

[TABLE]

are the one-body energies with mean-field shifts included. Please notice that in Eq. (37) the kernel $K^{0}$ is now without the mean-field contribution, i.e. in Eq. (LABEL:kernel2-static) $V$ has been replaced by $V-V^{\text{MF}}$ where $V^{\text{MF}}$ is the mean-field potential operator. Not to introduce new symbols, from now on, $K$ should always be understood in this way.

This needs, however, further elaboration and discussions. Actually the existence of the kernel $K$ hinges entirely on the existence of the inverse of the one-frequency response function $\tilde{R}$, via $K={\tilde{R}}^{{0}^{-1}}-\tilde{R}^{-1}$. Again, this is in complete analogy to the case of the Dyson equation: $\Sigma=G^{0^{-1}}-G^{-1}$. For readers who may doubt about the existence of $\tilde{R}^{-1}$, we announce that below we will find approximate expressions for $K$ which reproduce known expressions from the literature which have been derived from the Hedin equations. We note that, as this was the case with the one-body self-energy, also here the single-frequency kernel in Eq. (LABEL:kernel2-total) splits into a purely static (instantaneous) and a dynamic (time-dependent) part. It is quite suggestive to interpret the purely static term $K^{0}$ as some kind of higher mean field. Below, we will give an explicit expression for it and will see that it contains static p-h correlation functions. Viewing the ground state as containing a gas of p-h quantum fluctuations, one can then interpret the purely static term as the (frequency-independent) mean field of those fluctuations. We will refer to $K^{0}$ as a “particle-hole mean field” and shall show below in which way it is related to a specific form of $W$ in the $GW$+BSE approach.

However, before that, let us transform the BSE-like equation by returning from $\tilde{R}$ to the original linear-response function $R$. It is straightforward to show that the latter then obeys the following equation

[TABLE]

The reader may be worried that we did not get rid of the possibly troublesome norm factor $N_{k_{1}k_{2}}=n_{k_{2}}-n_{k_{1}}$ in the denominator in Eqs. (LABEL:kernel2-static) and (LABEL:kernel2) implying that there may be numerical troubles for situations where $n_{k_{1}}\simeq n_{k_{2}}$. Actually, there are good reasons for this division. It is analogous to, e.g., what happens with the generator coordinate method (GCM) where also the norm kernel has to be diagonalized and configurations corresponding to zero eigenvalues be eliminated Ring and Schuck (1980). This always happens when expanding the quantity of interest into a non-orthogonal basis set (here the products $c^{\dagger}_{k_{2}}c_{k_{1}}$), a feature which is also underpinning our approach. Actually the present approach is practically equivalent to the EOM of Rowe Hedin (1968) (see also Ref. Schuck and Tohyama (2016)), where one expands an excited state into a series of components where higher and higher many-body operators act on the formally exact ground state. Exactly the same type of norm factors $N$ as here appear on the right-hand-side of an eigenvalue problem as in Eq. (38) below. We will not further discuss this very general case here as we will henceforth work in the p-h/h-p subspace (for definition, see below), where this problem does not appear, and replace the norm kernel in the denominator by its HF expression. Taking higher-order corrections of the norm in the denominator into account has probably little influence on the accuracy of the results as suggested by some explicit examples Delion et al. (2005).

Before going on, let us transform our BSE-like equation into an eigenvalue problem. As just mentioned, this will be done in the p-h/h-p subspace

[TABLE]

with $h,h^{\prime}$ referring to hole states ($h,h^{\prime}\leq k_{\rm F}$, where $k_{\rm F}$ is the Fermi momentum) and $p,p^{\prime}$ referring to particle states ($p,p^{\prime}>k_{\rm F}$). This equation is of the typical RPA form as described, e.g., in Ref. Ring and Schuck (1980). The present generic equation is, however, potentially much more general because in principle the $A$ and $B$ matrices depend on the eigenvalues $\Omega_{\nu}$ and on the amplitudes $(X^{\nu},Y^{\nu})$, the latters being related to the ground state $|0\rangle$ and excited state $|\nu\rangle$ by $X^{\nu}_{ph}=\langle 0|c^{\dagger}_{h}c_{p}|\nu\rangle$ and $Y^{\nu}_{ph}=\langle 0|c^{\dagger}_{p}c_{h}|\nu\rangle$. The $A$ and $B$ matrices are related to the one-frequency kernel $K$ in Eq. (LABEL:kernel2-total) by

[TABLE]

For example, to first order in the interaction this gives

[TABLE]

where the occupation factors $n_{i}$ have been replaced by their uncorrelated form $n^{0}_{i}$, and Eq. (38) reduces to the standard RPA equation (with exchange) or TDHF. We will not further elaborate on the eigenvalue form of our approach and rather continue investigating the one-frequency kernel $K$.

III.2 The purely static part of the kernel $K$

Let us now discuss the $K^{0}$ term of the kernel and see how far it is related to the static $W$ kernel of the $GW$ approach. To establish an explicit form for $K^{0}$, we have to evaluate the double commutator contained in the particle-hole mean-field part of Eq. (LABEL:kernel2-static). One finds

[TABLE]

where

[TABLE]

is the fully correlated part (i.e. the cumulant) of the two-body density matrix. We see that $K^{0}$ is divided into four parts: the first term on the right-hand side is the usual RPA antisymmetrized interaction term. We should realize that in this first term the norm factor on the right of the interaction has been divided out [see Eq. (LABEL:kernel2-static)] and that, contrary to standard RPA, the occupation factors are in principle not the HF ones but the correlated ones. Neglecting all the terms involving $C$ in Eq. (41) but keeping correlations in the occupancies leads to the so-called renormalized RPA (r-RPA) briefly explained further in the Appendix. The next two terms are the one-body self-energy contributions (either the hole or the particle is not connected to the interaction). The remaining two-body correlation terms connect particles and holes. They can be qualified as screening terms and we want to investigate them further. The screening terms can be divided into two groups: the first two terms correspond to an exchange of p-h fluctuations between the particle and hole and are, therefore, responsible for the screening of the long-range part of the interaction. This can be seen from the ordering of the indices $k_{1}$ and $k_{3}$ in the matrix element of the interaction. Clearly a creator and a destructor are correlated. The second two terms correspond to an exchange of p-p/h-h fluctuations, that is they sum p-p/h-h ladder diagrams. They take care of the short-range correlations. Let us mention that neglecting the dynamic kernel, a self-consistent scheme for the two-body correlation function can be established, since it is given by integrating $R(\omega)$ over the frequency in the upper/lower half complex plane. This self-consistent scheme is referred to as SCRPA. It has the nice quality that all desirable properties of standard RPA, such as the fulfillment of the sum rule and conservation laws are maintained. This is explicitly shown in Ref. Delion et al. (2016). In the past, it has produced encouraging results for several non-trivial model cases Dukelsky and Schuck (1990); Jemai et al. (2005); Delion et al. (2005); Schuck and Tohyama (2016). Let us mention that Eq. (41) has been given earlier Dukelsky et al. (1998) and that it has recently also been derived by Chatterjee and Pernal Chatterjee and Pernal (2012) for applications in chemical physics including, however, also diagonal configurations.

In order to establish a connection with the static screened interaction $W$ of the $GW$+BSE approach, we consider in more detail the p-h fluctuation terms. As an example, let us consider the fourth term on the right-hand side of Eq. (41) and evaluate it first to second order in the interaction. Since $C$ is at least of first order, we will elaborate this and get the corresponding $K^{0}$ to second order. First let us give the relation between $C$ and the linear-response function $R$

[TABLE]

where $\bar{n}_{k}=1-n_{k}$ and

[TABLE]

The reader may wonder why there is the last term on the right-hand side of above Eq. (44). The point is that since $R$ is the solution of the BSE-like equation, it does not contain such disconnected terms where the time $t$ does not communicate with time $t^{\prime}$, see the definition of $R$ in Eq. (19). This is also easily seen in solving, e.g., Eq. (37) to lowest order, that is without the kernel $K$ and using the HF form of norm $N$, which leads to Eq. (46) below. However, on the left-hand side, in the expectation value of the two-body-density-matrix operators, such terms are contained and, therefore, we have to add them on the right-hand side. A good way to see this is to evaluate Eq. (43) in the HF approximation where $C=0$ by definition. Then the right-hand side must also be zero which is only the case if the extra term is added. Let us now expand the response function in Eq. (44) up to first order

[TABLE]

where $R^{0}$ is the non-interacting HF linear-response function

[TABLE]

Inserting Eq. (LABEL:v-1) into Eq. (44) and then Eq. (44) into Eq. (43), and using $N^{0}_{k_{1}k_{2}}=N^{0^{-1}}_{k_{1}k_{2}}$, one obtains for the fourth term in the $K^{0}$ kernel

[TABLE]

Please note that the lowest-order term in Eq. (LABEL:v-1) is cancelled by the second term on the right-hand side of Eq. (44). The expressions (47) and (48) (see below for the latter) are the only ones which contribute to $K^{0}$ at second order with a p-h bubble exchange. We see this from the fact that the index $k_{2}$ in Eq. (47) is a hole, then $k_{1}$ must be a particle and, since $k_{4}$ is a particle, $k_{3}$ must be a hole because our convention is that the index couple $k_{1}k_{2}$ or $k_{3}k_{4}$ can only be p-h or h-p. As we see, the term in Eq. (47) contributes to the $B$ matrix in Eq. (39). In analogy, we obtain for the fifth term in Eq. (41)

[TABLE]

Again this term only contributes to the $B$ matrix of Eq. (39). Both terms correspond to the first two terms in Eq. (34) of Ref. Rebolini and Toulouse (2016). If we treat the last two (p-p/h-h) terms of our Eq. (41) in the same way as the p-h terms, we also reproduce the other two terms in Eq. (34) of Ref. Rebolini and Toulouse (2016)

[TABLE]

As before, these terms only contribute to the $B$ matrix of Eq. (39).

From Eq. (41), it is clear that in Eq. (LABEL:v-1) we can replace $R^{0}_{ll^{\prime}}\delta_{ll_{1}}\delta_{l^{\prime}l_{3}}$ by the full linear-response function $R_{ll^{\prime}l_{1}l_{3}}$ what leads to a better approximation where the exchange p-h bubble $ll^{\prime}$ is replaced, e.g., by the RPA or even higher approximations. In general we have for $R$ in frequency space

[TABLE]

where $\eta\to 0^{+}$. Actually this can be done also in Eq. (49) where one can resum the pp ladders taking care of the short-range correlations. We will not further dwell on those extensions of our formalism for the moment.

Let us now consider the self-energy corrections in Eq. (41). For instance, let us extract a further interaction. For example, we obtain (indicating the time variables as subscripts for compactness)

[TABLE]

and an analogous expression for the renormalization of the particle line. Evaluating, as before, $C$ to first order, one obtains for the second and third term

[TABLE]

We realize that this expression contributes only to the $A$ matrix and that the second-order contribution to $K^{0}$ is now complete.

Of course, as in the case of the screening terms, we also can sum the p-h bubbles to a full linear-response function. For this, we should factorize the three-body propagator in Eq. (51) into a product of a response function and a one-body propagator

[TABLE]

As indicated, there are four different ways to do this factorization and, thus, we multiply the final expression with this factor to obtain

[TABLE]

Proceeding with the second self-energy correction in the same way, one obtains an analogous expression.

We will discuss the relation with the static form of $W$ below in Sec. IV. The well-known problem coming from the approximation of Eq. (53) is that if the response function is replaced by its uncorrelated p-h response in Eq. (54), one does not recover the correct second-order expression of the kernel Danielewicz and Schuck (1994). The result is by a factor of two too large and, therefore, one usually subtracts the second-order contribution in order to obtain the correct lowest-order contribution of $\Sigma$ to the kernel and also of the ensuing RPA correlation energy Danielewicz and Schuck (1994). In principle this subtraction procedure is not completely correct, since the corresponding imaginary part of the self-energy has no definite sign. How this can be fixed is explained in Refs. Danielewicz and Schuck (1994); Adachi and Schuck (1989). We will not dwell on this here and ignore this subtlety in the remainder of the paper, supposing that the uncorrelated terms are small, and we will concentrate on the comparison with approximations given in the literature Romaniello et al. (2009); Sangalli et al. (2011); Rebolini and Toulouse (2016) where this problem is also not addressed.

In our approach, the one-body self-energy corrections appear directly in the purely static part of the integral kernel. So the self-energy corrections are treated consistently with the screening terms. This may be important because it is known that often there are significant cancellations between both contributions. Since screening and self-energy corrections in the purely static part $K^{0}$ of the kernel involve again the response function $R$, as explained above, a self-consistent cycle can be established.

III.3 The dynamic part of the one-frequency kernel $K$

Let us now discuss the time-dependent (dynamic) part $K^{\mathrm{dyn}}$ of the interaction kernel in Eq. (LABEL:kernel2) which is a p-h irreducible four-body propagator of the 2p-2h type. It is straightforward to evaluate it to lowest order. Since the dynamic kernel is already explicitly of second order in the interaction, it is sufficient to evaluate the 2p-2h propagator to lowest (HF) order. Using $J_{k_{1}k_{2}}=c_{k_{2}}^{\dagger}j_{k_{1}}-j^{\dagger}_{k_{2}}c_{k_{1}}$ and dropping for now, for simplicity, the factor $N_{k_{1}^{\prime}k_{2}^{\prime}}^{{0^{-1}}}$ which just gives a minus sign for the contribution to the $B$ matrix, the full expression of the second-order dynamic kernel is then

[TABLE]

where the subscript “0” indicates that this term is evaluated to lowest order. The first two terms are self-energy corrections recognizable by the index pair $k_{1},k^{\prime}_{1}$ or $k_{2},k^{\prime}_{2}$ whereas the other two terms have “mixed” indices. These expressions describe the decay of a p-h mode into uncorrelated (incoherent) 2p-2h states. The four terms have different meanings. The two terms with either $jj^{\dagger}$ or $j^{\dagger}j$ describe, as just mentioned, dynamic self-energy corrections to the particle and the hole states, respectively. The other two terms with $jj$ or $j^{\dagger}j^{\dagger}$ describe a p-h exchange between the particle and the hole. Such incoherent processes have already been considered a long time ago by Landau Landau (1957) in his study of the damping of zero sound in a Fermi liquid. A detailed study of this is given in Ref. Adachi and Schuck (1989).

To obtain the spectral representation of $K^{{\rm dyn},(2)}$, we just need to consider the generic propagator

[TABLE]

calculate its Fourier transform

[TABLE]

with $\delta_{k_{1}k_{2}k^{\prime}_{1}k^{\prime}_{2}}=\delta_{k_{1}k^{\prime}_{1}}\delta_{k_{2}k^{\prime}_{2}}-\delta_{k_{1}k^{\prime}_{2}}\delta_{k_{2}k^{\prime}_{1}}$, and use this in Eq. (LABEL:K-dyn-0). The obtained expression is well known in the nuclear physics literature Wambach (1988). More recently such expressions have been derived by Rebolini and Toulouse Rebolini and Toulouse (2016), but starting from the Hedin equations and without including the self-energy corrections.

Instead of approximating the 2p-2h propagator by its uncorrelated expression, we can include higher-order effects. For example one can factorize it into a product of a response function and an uncorrelated p-h propagator. Or, one can factorize it into a product of two linear-response functions. The choice will depend on the physical situation. Such approximations have been considered in Ref. Schuck (1976). As in the case of the self-energy of the one-body Dyson equation, those factorizations do not give, however, the correct lowest-order limit of the kernel. If important, one has to correct for it. How this can be done consistently is explained, as already mentioned, in Refs. Danielewicz and Schuck (1994); Adachi and Schuck (1989). Let us give explicit expressions for the spectral representation of $K^{\rm dyn}$ for the case where we approximate the 2p-2h propagator into a product of a linear-response function times an uncorrelated p-h propagator. Typically, one will evaluate the response function with the RPA method. It is then easy to get the spectral form of $K^{\rm dyn}$ (skipping now the p-p/h-h contributions)

[TABLE]

We note that the first two terms on the right-hand side of the above equation correspond again to self-energy corrections, whereas the last two terms are contributions where p-h modes are exchanged between the particle and the hole. We also realize that these exchange contributions correspond to Eq. (31) of Ref. Romaniello et al. (2009). The “backward going” terms do not contribute as easily realized. If we consider the static limit ($\omega=0$), they should be considered together with Eqs. (47) and (48). This shall be discussed in more detail in the next section.

Before doing so, it may be worth showing how to include further p-h correlations in summing up the free p-h propagators, contained in Eq. (58), to extra RPA modes. This is most easily done by factorizing the 2p-2h propagator into a fully antisymmetrized product of two p-h response functions

[TABLE]

where

[TABLE]

is the Fourier transform into time space of Eq. (50). The Fourier transform of Eq. (59) into frequency space is then easily performed with Eq. (60)

[TABLE]

where “exchange terms” means that all exchange terms present in Eq. (59) should be included also here. Inserting Eq. (61) into Eq. (LABEL:kernel2) yields an expression equivalent to Eq. (23) of Ref. Sangalli et al. (2011) (see also Ref. Schuck (1976)). Notably only the first term with $+i\eta$ will survive, that is, it enters only the $A$ matrix, as also pointed out in Ref. Sangalli et al. (2011). Since it is fully antisymmetric between the two-particle states and two-hole states in entrance and exit channels, the approximation gives a conserving approximation for the response function Kadanoff and Baym (1962); Delion et al. (2016).

IV Comparison with $GW$+BSE

Let us now consider similarities and differences of the present approach to the response function and the $GW$+BSE scheme as commonly used in condensed-matter and chemical physics.

A first point consists in the fact that in the present formalism all Coulomb matrix elements are antisymmetrized [see Eq. (2)] whereas in the $GW$+BSE scheme all exchange matrix elements are usually absent besides the one contained in the first order of the screening term. This also concerns the $W$ used within the RPA in condensed-matter physics: only the bubble diagrams are resummed, as it was done in the original work of Bohm and Pines Bohm and Pines (1953). Including then the static limit of $K^{\rm dyn}$ (i.e. at $\omega=0$) in Eq. (58) to $K^{0}$ of Eq. (41) yields an expression very similar to the “excitonic” Hamiltonian $H^{2p,exc}$ in Eqs. (16) and (21) of Ref. Romaniello et al. (2009). However, there are also substantial differences and, for a detailed comparison, let us give our full static expression here (summing the p-h bubble exchange to a full response function and skipping the self-energy and p-p/h-h contributions for easier comparison)

[TABLE]

We see that the first two terms belong to the $B$ matrix and the last two terms to the $A$ matrix of Eq. (39). In the $GW$+BSE scheme all antisymmetrized matrix elements $\bar{v}_{k_{1}k_{2}k_{3}k_{4}}$ in Eq. (62) are replaced by only the direct term $v_{k_{1}k_{2}k_{3}k_{4}}$. In addition, in the denominators the differences of orbital energies are absent, so that only the RPA roots $\Omega_{\nu}$ remain, which corresponds to the static $W(0)$ of the $GW$+BSE kernel (see, e.g., Ref. Romaniello et al. (2009)). It is difficult to judge the combined effect of the two differences of $GW$+BSE with respect to the above expression in Eq. (62). The extra orbital energies in the denominators in our expressions have, however, certainly a reduction effect. A detailed numerical evaluation is out of the scope of the present work but shall eventually be presented in the future. It seems to us that the appearance of the orbital energies in the denominators has its justification. In the work of Romaniello et al. Romaniello et al. (2009) they also appear as an extra static contribution from their $\widetilde{W}$ expression in Eq. (27) in Ref. Romaniello et al. (2009) in putting therein $\omega_{\lambda}=0$. It is clear that, although the terms in Eq. (62) are instantaneous, there can never be an exact equal time process when an RPA mode crosses between the particle and the hole lines. There is always an infinitesimal time difference allowing for the orbital energies to appear in the denominators of Eq. (62).

A further difference of our EOM approach is that the self-energy contributions appear directly in the kernel. It is possible to resum them separately, which would lead to dressed quasi-particles (and quasi-holes), quite similarly to the $GW$+BSE scheme. In this respect we do not see any significant difference between the two approaches.

In our scheme we obtain the same (approximate) dynamic contributions to the kernel as obtained by Sangalli et al. Sangalli et al. (2011) [see their Eqs. (22) and (23)]. They also contain the self-energy contributions. It is also clear that those dynamic contributions only renormalize the $A$ matrix and give no contribution to the $B$ matrix. On the other hand in Ref. Sangalli et al. (2011), the $B$ matrix is not renormalized, not containing the additional correlations which are summed up in Eq. (41). In Ref. Rebolini and Toulouse (2016), the renormalization of the $B$ matrix is given only to lowest order. Let us also point out that the so-called time-blocking approximation (TBA) Litvinova and Tselyaev (2007) invented recently in nuclear physics to derive a kernel depending only on one frequency certainly has a very close relation with the procedures employed in Refs. Rebolini and Toulouse (2016) and Romaniello et al. (2009); Sangalli et al. (2011). It may be relevant to realize that the first two terms in Eq. (62) which derive from Eq. (41) and renormalize the $B$ matrix are an approximation to Eq. (41) [see the Appendix]. As we will see in the next section in some two-body problems it may be important to keep the full expression of Eq. (41).

V Illustration on the Hubbard molecule

The Hubbard model describes electrons on a lattice with the Coulomb interaction replaced by an on-site constant $U$. The well-known Hamiltonian is given by

[TABLE]

where $c^{\dagger}_{i\sigma}$ and $c_{i\sigma}$ are the electron creation and destruction operators at site $i$ with spin projection $\sigma$ and the $\hat{n}_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$ are the number operators for electrons at site $i$ with spin projection $\sigma$. As usual $t$ is the nearest-neighbor hopping integral. For demonstration purposes, in this work, we will limit ourselves to the simplest non-trivial case which is the one of two sites ($N_{\text{s}}=2$) with two electrons, the so-called Hubbard molecule. As the problem has already been solved exactly with the SCRPA method Jemai et al. (2005) derived from Rowe’s Hedin (1968) EOM, we only will outline the basic principle here using, however, the present approach. It is advantageous to write the Hamiltonian in momentum space (we consider periodic boundary conditions)

[TABLE]

where $\hat{n}_{{\bf k},\sigma}=c^{\dagger}_{{\bf k},\sigma}c_{{\bf k},\sigma}$ is the occupation number operator of the momentum-spin mode $({\bf k},\sigma)$ and $\epsilon_{k}=-2t\cos(k)$ are the one-body energies with the lattice spacing set to unity. Because of having only two electrons and the periodic boundary conditions, the only allowed momenta are $k_{1}=0$ and $k_{2}=-\pi$. Accordingly, we only have two types of p-h operators: $J_{\sigma}=c^{\dagger}_{k_{1},\sigma}c_{k_{2},\sigma}$ with $\sigma=\pm 1/2$. Let us introduce the “charge” and “spin” operators

[TABLE]

and consider the charge and spin linear-response functions

[TABLE]

We therefore have to consider two $2\times 2$ matrix response functions. For this very simple example it so happens that the dynamic part $K^{\rm dyn}$ of the one-frequency kernel decouples from the purely static part $K^{0}$, and only $K^{0}$ contributes in the p-h/h-p space. As seen from Eq. (41), the purely static kernel $K^{0}$ only contains static two-body correlation functions. They can be calculated from integrating ${R}^{(\pm)}(\omega)$ over the frequency in the upper/lower half complex plane. Since additionally the occupation numbers can also be expressed via the static two-body correlation functions as (see Ref. Jemai et al. (2005))

[TABLE]

we have a closed system of equations which can be solved. It turns out that the exact solution is obtained. This is explained in detail in Ref. Jemai et al. (2005) starting, however, with the equivalent EOM for RPA operators and not with the Green functions and we will not repeat the whole procedure here.

The fact that the two-body problem is solved exactly by the SCRPA in the Hubbard model is also found in several other models, like the Lipkin model Schuck and Tohyama (2016) and the pairing model Schuck and Tohyama (2016); Hirsch et al. (2002). However, it is not a general feature of SCRPA that it solves any two-body problem exactly. Generally there exist specific 2p-2h configurations which have to be taken into account when solving a two-body problem. It should also be pointed out that the two-body correlation functions in Eq. (41) cannot be further approximated if the exact solution for, e.g., the Hubbard molecule shall be obtained. Already the forms in Eqs. (47) and (48) are approximations to Eq. (41) even if the exchange bubble is resummed to a full linear-response function and it is likely that they will not maintain the exact solution. It is thus seen from this example that the present approach leads in a systematic way to manageable expressions which, if necessary, sum higher correlations than is the case with the $GW$+BSE approach.

Finally, a reduced version of SCRPA called r-RPA is presented in the Appendix. In a separate paper J. Li, N. Drummond, P. Schuck, and V. Olevano (2018) this approximation has been applied to the case of a real system: the helium atom. In that paper both the r-RPA and the $GW$+BSE solutions are compared to the exact Hylleraas solution evidencing their respective performances, similarities and differences, on a real system.

VI Conclusions

The objective of this work was three-fold. First, we derived a formally exact one-frequency-only BSE-like equation for the linear-response function whose integral kernel only depends on the single frequency of the applied field. Explicit expressions of this kernel in terms of higher Green functions are presented. They lend themselves very naturally to physically motivated approximations. Second, in this way, known approximations of a single frequency kernel derived from Hedin’s equations are straightforwardly recovered. It is shown that with our approach not only the second-order expressions for the static ($B$ matrix) and dynamic BSE kernel given in Ref. Rebolini and Toulouse (2016) can be recovered but that these second-order terms can naturally be resummed to full linear-response functions. This has also been shown in Refs. Romaniello et al. (2009); Sangalli et al. (2011) for the dynamic part but the renormalization of the static part ($B$ matrix), which is of the same order as the dynamic one, is missing there. Taking the static limit ($\omega=0$) of the dynamic part, we obtain a complete expression for the static limit of our kernel. Third, this then allows us to make a detailed comparison with the static limit of the kernel of the well known $GW$+BSE approach. Between both static approaches there exist, besides quite some similarities, also substantial differences which may be interesting to study further in future work with numerical examples. At the end of the paper, we also show that for the so-called Hubbard molecule the exact solution can be recovered from our approach. This is only possible with a consistent and fully resummed static kernel as presented here. Let us finally mention that the present formalism is very much related to the EOM introduced by Rowe and further elaborated in Ref. Schuck and Tohyama (2016).

VII Acknowledgements

P.S. is grateful to D. Delion, J. Dukelsky, M. Jemai, and M. Tohyama for a fruitful collaboration on the EOM method for correlation functions. Useful comments on the manuscript from M. Holzmann, K. Pernal, and P. Romaniello are appreciated.

Appendix: The renormalized RPA

A very much simplified version of SCRPA consists in neglecting in $K^{0}$ all the $C$ correlation function terms. Then, one obtains for the one-frequency BSE-like equation

[TABLE]

with

[TABLE]

We see that this renormalized RPA (r-RPA) equation is like the standard RPA besides the fact that the occupation numbers are the correlated ones and not the HF ones. We thus have to give an expression for the $n_{k}$’s which couple back to the RPA. Such an approximation for the occupation numbers $n_{k}$ has, e.g., been derived by Catara et al. F. Catara, G. Piccitto, M. Sambataro, and N. Van Giai (1996). The expressions are given by

[TABLE]

and

[TABLE]

where the two-body density matrix can directly be obtained from the linear-response function. We, therefore, have established a minimal self-consistent system of equations where the occupation numbers are calculated from the response function.