Quasiparticles in Neon using the Faddeev Random Phase Approximation
C. Barbieri
Gesellschaft für Schwerionenforschung, Planckstr. 1, D-64291, Darmstadt, Germany
D. Van Neck
Laboratory of Theoretical Physics, Ghent University, Proeftuinstraat 86, B-9000 Gent, Belgium
W.H. Dickhoff
Department of Physics, Washington University, St. Louis, MO 63130, USA
(Dated: November 2, 2018)
The spectral function of the closed-shell Neon atom is computed by expanding the electron self-energy
through a set of Faddeev equations. This method describes the coupling of single-particle degrees of free-
dom with correlated two-electron, two-hole, and electron-hole pairs. The excitation spectra are obtained using
the Random Phase Approximation, rather than the Tamm-Dancoff framework employed in the third-order al-
gebraic diagrammatic contruction [ADC(3)] method. The difference between these two approaches is studied,
as well as the interplay between ladder and ring diagrams in the self-energy. Satisfactory results are obtained
for the ionization energies as well as the energy of the ground state with the Faddeev-RPA scheme that is also
appropriate for the high-density electron gas.
PACS numbers: 31.10.+z,31.15.Ar
I. INTRODUCTION
Ab initio treatments of electronic systems become unwork-
able for sufficiently complex systems. On the other hand,
the Kohn-Sham formulation [1] of density functional the-
ory (DFT) [2] incorporates many-body correlations (beyond
Hartree-Fock), while only single-particle (sp) equations must
be solved. Due to this simplicity DFT is the only feasible
approach in some modern applications of electronic structure
theory. There is therefore a continuing interest both in devel-
oping new and more accurate functionals and in studying con-
ceptual improvements and extensions to the DFT framework.
In particular it is found that DFT can handle short-range inter-
electronic correlations quite well, while there is room for im-
provements in the description of long-range (van der Waals)
forces and dissociation processes.
Microscopic theories offer some guidance in the devel-
opment of extensions to DFT. Orbital dependent function-
als can be constructed using many-body perturbation theory
(MBPT) [3, 4]. More recently, the development of general
ab initio DFT [5, 6] addressed the lack of a systematic im-
provement in DFT methods. In this approach one considers
an expansion of the exact ground-state wave function (e.g.,
MBPT or coupled cluster) from a chosen reference determi-
nant. Requiring that the correction to the density vanishes at a
certain level of perturbation theory allows one to construct the
corresponding approximation to the Kohn-Sham potential.
A different route has been proposed in Ref. [7] by devel-
oping a quasi-particle (QP)-DFT formalism. In the QP-DFT
approach the full spectral function is decomposed in the con-
tribution of the QP excitations, and a remainder or background
part. The latter part is complicated, but does not need to be
known accurately: it is sufficient to have a functional model
for the energy-averaged background part to set up a single-
electron selfconsistency problem that generates the QP exci-
tations. Such an approach is appealing since it contains the
well-developed standard Kohn-Sham formulation of DFT as
a special case, while at the same time emphasis is put on the
correct description of QPs, in the Landau-Migdal sense [8].
Hence, it can provide an improved description of the dynam-
ics at the Fermi surface. Given the close relation between QP-
DFT and the Green’s function (GF) formulation of many-body
theory [9, 10], it is natural to employ ab initio calculations
in the latter formalism to investigate the structure of possible
QP-DFT functionals. In this respect it is imperative to iden-
tify which classes of diagrams are responsible for the correct
description of the QP physics.
Some previous calculations, based on GF theory, have fo-
cused on a self-consistent treatment of the self-energy at the
second order [11, 12, 13] for simple atoms and molecules. For
the atomic binding energies it was found that the bulk of corre-
lations, beyond Hartree-Fock, are accounted for while signifi-
cant disagreement with experiment persists for QP properties
like ionization energies and electron affinities. The formal-
ism beyond the second-order approximation was taken up in
Ref. [14, 15, 16, 17, 18] by employing a self-energy of the
GW type [19]. In this approach, the random phase approx-
imation (RPA) in the particle-hole (ph) channel is adopted
to allow for possible collective effects on the atomic excited
states. The latter are coupled to the sp states by means of dia-
grams like the last two in Fig. 1(c). Two variants of the G0W0
formalism were employed in Ref. [14] (where the subscript
“0” indicates that non-dressed propagators are used). In the
first only the direct terms of the interelectron Coulomb po-
tential are taken into account. In the second version, also the
exchange terms are included when diagonalizing the ph space
[generalized RPA (GRPA)] and in constructing the self-energy
[generalized GW (GGW)]. Although the exchange terms are
known to be crucial in order to reproduce the experimentally
observed Rydberg sequence in the excitation spectrum of neu-
tral atoms, they were found to worsen the agreement between
the theoretical and experimental ionization energies [14].
http://arxiv.org/abs/0704.1542v2
In the GW approach the sp states are directly coupled
with the two-particle–one-hole (2p1h) and the two-hole–one-
particle (2h1p) spaces. However, only partial diagonalizations
(namely, in the ph subspaces) are performed. This procedure
unavoidably neglects Pauli correlations with the third particle
(or hole) outside the subspace. In the case of the GGW ap-
proach, this leads to a double counting of the second order
self-energy which must be corrected for explicitly [20, 21].
We note that simply subtracting the double counted diagram is
not completely satisfactory here, since it introduces poles with
negative residues in the self-energy. More important, the inter-
action between electrons in the two-particle (pp) and two-hole
(hh) subspaces are neglected altogether in (G)GW. Clearly, it
is necessary to identify which contributions, beyond GGW,
are needed to correctly reproduce the QP spectrum.
In this respect, it is known that highly accurate descrip-
tions of the QP properties in finite systems can be obtained
with the algebraic diagrammatic construction (ADC) method
of Schirmer and co-workers [22]. The most widely used third-
order version [ADC(3)] is equivalent to the so-called extended
2p1h Tamm-Dancoff (TDA) method [23] and allows to pre-
dict ionization energies with an accuracy of 10-20 mH in
atoms and small molecules. Upon inspection of its diagram-
matic content, the ADC(3) self-energy is seen to contain all
diagrams where TDA excitations are exchanged between the
three propagator lines of the intermediate 2p1h or 2h1p prop-
agation. The TDA excitations are constructed through a diag-
onalization in either 2p1h or 2h1p space, and neglect ground-
state correlations. However, it is clear that use of TDA leads
to difficulties for extended systems. In the high-density elec-
tron gas e.g., the correct plasmon spectrum requires the RPA
in the ph channel, rather than TDA.
In order to bridge the gap between the QP description in
finite and extended systems, it seems therefore necessary to
develop a formalism where the intermediate excitations in the
2p1h/2h1p propagator are described at the RPA level. This
can be achieved by a formalism based on employing a set of
Faddeev equations, as proposed in Ref. [24] and subsequently
applied to nuclear structure problems [25, 26, 27]. In this ap-
proach the GRPA equations are solved separately in the ph and
pp/hh subspaces. The resulting polarization and two-particle
propagators are then coupled through an all-order summation
that accounts completely for Pauli exchanges in the 2p1h/2h1p
spaces. This Faddeev-RPA (F-RPA) formalism is required if
one wants to couple propagators at the RPA level or beyond.
Apart from correctly incorporating Pauli exchange, F-RPA
takes the explicit inclusion of ground-state correlations into
account, and can therefore be expected to apply to both finite
and extended systems. The ADC(3) formalism is recovered
as an approximation by neglecting ground-state correlations
in the intermediate excitations (i.e. replacing RPA with TDA
phonons).
In this work we consider the Neon atom and apply the F-
RPA method to a nonrelativistic electronic problem for the
first time. The relevant features of the F-RPA formalism (also
extensively treated in Ref. [24]), are introduced in Sect. II.
The application to the Neon atom is discussed in Sec. III,
where we also investigate the separate effects of the ladder
and ring series on the self-energy, as well as the differences
between including TDA and RPA phonons. Our findings are
summarized in Sec. IV. Some more technical aspects are rele-
gated to the appendix, where the interested reader can find the
derivation of the Faddeev expansion for the 2p1h/2h1p propa-
gator, adapted from Ref. [24]. In particular, the approach used
to avoid the multiple-frequency dependence of the Green’s
functions is discussed in App. A 1, along with its basic as-
sumptions. The explicit expressions of the Faddeev kernels
are given in App. A 3. Together with Ref. [24], the appendix
provides sufficient information for an interested reader to ap-
ply the formalism.
II. FORMALISM
The theoretical framework of the present study is that of
propagator theory, where the object of interest is the sp prop-
agator, instead of the many-body wave function. In this paper
we will employ the convention of summing over repeated in-
dices, unless specified otherwise. Given a complete orthonor-
mal basis set of sp states, labeled by α,β,..., the sp propagator
can be written in its Lehmann representation as [9, 10]
gαβ(ω) =
ω − ε+n + iη
ω − ε−k − iη
, (1)
where Xnα = 〈Ψ
0 〉 (Y
α = 〈Ψ
|cα|Ψ
0 〉) are the
spectroscopic amplitudes, cα (c
) are the second quantization
destruction (creation) operators and ε+n = E
n − E
EN0 − E
). In these definitions, |ΨN+1n 〉, |Ψ
〉 are the
eigenstates, and EN+1n , E
k the eigenenergies of the (N ± 1)-
electron system. Therefore, the poles of the propagator reflect
the electron affinities and ionization energies.
The sp propagator solves the Dyson equation
gαβ(ω) = g
αβ(ω) + g
αγ(ω)Σ
γδ(ω)gδβ(ω) , (2)
which depends on the irreducible self-energy Σ⋆(ω). The lat-
ter can be written as the sum of two terms
αβ(ω) = Σ
Vαλ,µν Rµνλ,γδε(ω) Vγδ,βε , (3)
where ΣHF represents the Hartree-Fock diagram for the self-
energy. In Eqs. (2) and (3), g0(ω) is the sp propagator for
the system of noninteracting electrons, whose Hamiltonian
contains only the kinetic energy and the electron-nucleus at-
traction. The Vαβ,γδ represent the antisymmetrized matrix el-
ements of the interelectron (Coulomb) repulsion. Note that
in this work we only consider antisymmetrized elements of
the interaction, hence, our result for the ring summation al-
ways compare to the generalized GW approach. Equation (3)
introduces the 2p1h/2h1p-irreducible propagator R(ω), which
carries the information concerning the coupling of sp states to
more complex configurations. Both Σ⋆(ω) and R(ω) have a
perturbative expansion as a power series in the interelectron
FIG. 1: (Color online) a) Diagrammatic expansion of R(ω) in terms
of the (antisymmetrized) Coulomb interaction and undressed prop-
agators. b) R(ω) is related to the self-energy according to Eq. (3).
c) By substituting the diagrams a) in the latter equation, one finds the
perturbative expansion of the self-energy.
interaction V̂ . Some of the diagrams appearing in the expan-
sion of R(ω) are depicted in Fig. 1, together with the corre-
sponding contributions to the self-energy. Note that already at
zero order in R(ω) (three free lines with no mutual interaction)
the second order self-energy is generated.
Different approximations to the self-energy can be con-
structed by summing particular classes of diagrams. In this
work we are interested in the summation of rings and ladders,
through the (G)RPA equations. In order to include such effects
in R(ω), we first consider the polarization propagator describ-
ing excited states in the N-electron system
Παβ,γδ(ω) =
〈ΨN0 |c
n 〉 〈Ψ
γcδ|Ψ
ENn − E
〈ΨN0 |c
γcδ|Ψ
n 〉 〈Ψ
ENn − E
, (4)
and the two-particle propagator, that describes the addi-
tion/removal of two electrons
gIIαβ,γδ(ω) =
〈ΨN0 |cβcα|Ψ
n 〉 〈Ψ
|ΨN0 〉
EN+2n − E
〈ΨN0 |c
|ΨN−2
〉 〈ΨN−2
|cβcα|Ψ
EN0 − E
. (5)
We note that the expansion of R(ω) arises from applying the
equations of motion to the sp propagator (1), which is associ-
ated to the ground state |ΨN0 〉. Hence, all the Green’s functions
appearing in this expansion will also be ground state based, in-
cluding Eqs. (4) and (5). However the latter contain, in their
Lehmann representations, all the relevant information regard-
ing the excitation of ph and pp/hh collective modes. The ap-
proach of Ref. [24] consists in computing these quantities by
=g II
Π(ph)
(pp/hh)II
(pp/hh)
FIG. 2: (Color online) Diagrammatic equations for the polarization
(above) and the two-particle (below) propagators in the (G)RPA ap-
proach. Dashed lines are always antisymmetrized Coulomb matrix
elements and the full lines represent free (undressed) propagators.
solving the ring-GRPA and the ladder-RPA equations [10],
which are depicted for propagators in Fig. 2. In the more
general case of a self-consistent calculation, a fragmented in-
put propagator can be used and the corresponding dressed
(G)RPA [D(G)RPA] equations [10, 28] solved [see Eqs. (A2a)
and (A2b)]. Since the propagators (4) and (5) reflect two-body
correlations, they still have to be coupled to an additional sp
propagator in order to obtain the corresponding approximation
for the 2p1h and 2h1p components of R(ω). This is achieved
by solving two separate sets of Faddeev equations.
Taking the 2p1h case as an example, one can split R(2p1h)(ω)
in three different components R̄(i)(ω) (i = 1, 2, 3) that differ
from each other by the last pair of lines that interact in their
diagrammatic expansion,
R̄(2p1h)
αβγ,µνλ
(ω) =
αβγ,µνλ(ω) −G
βαγ,µνλ(ω)
i=1,2,3
R̄(i)
αβγ,µνλ
(ω) ,
where G0
(ω) is the 2p1h propagator for three freely propa-
gating lines. These components are solutions of the following
set Faddeev equations [29]
R̄(i)
αβγ,µνλ
(ω) = G0
αβγ,µ′ν′λ′(ω) Γ
µ′ν′λ′,µ′′ν′′λ′′
µ′′ν′′λ′′ ,µνλ
(ω) + R̄(k)
µ′′ν′′λ′′ ,µνλ
(ω) (7)
µ′′ν′′λ′′ ,µνλ(ω) −G
ν′′µ′′λ′′ ,µνλ(ω)
, i = 1, 2, 3
where (i, j, k) are cyclic permutations of (1, 2, 3). The inter-
action vertices Γ(i)(ω) contain the couplings of a ph or pp/hh
collective excitation and a freely propagating line. These are
given in the Appendix in terms of the polarization (4) and
two-particle (5) propagators. Equations. (7) include RPA-like
phonons and fully describe the resulting energy dependence
of R(ω). However, they still neglect energy-independent
contributions–even at low order in the interaction–that also
correspond to relevant ground-state correlations. The latter
can be systematically inserted according to
R(2p1h)
αβγ,µνλ
(ω) = Uµνλ,µ′ν′λ′ R̄
(2p1h)
µ′ν′λ′,µ′′ν′′λ′′
(ω) U†
µ′′ν′′λ′′ ,µνλ
, (8)
where R(ω) is the propagator we employ in Eq. (3), R̄(ω) is
the one obtained by solving Eqs. (7), U ≡ I + ∆U, and I is the
(pp/hh)
Π(ph)
g II (pp/hh)
Π(ph)
FIG. 3: (Color online) Example of one of the diagrams that are
summed to all orders by means of the Faddeev Eqs. (7) (left). The
corresponding contribution to the self-energy, obtained upon inser-
tion into Eq. (3), is also shown (right).
identity matrix. Following the algebraic diagrammatic con-
struction method [22, 23], the energy independent term ∆U
was determined by expanding Eq. (8) in terms of the inter-
action and imposing that it fulfills perturbation theory up to
first order (corresponding to third order in the self-energy).
The resulting ∆U, employed in this work, is the same as in
Ref. [23] and is reported in App. A 3 for completeness. It
has been shown that the additional diagrams introduced by
this correction are required to obtain accurate QP properties.
Equations. (7) and (8) are valid only in the case in which a
mean-field propagator is used to expand R(ω). This is the case
of the present work, which employs Hartree-Fock sp propaga-
tors as input. The derivation of these equations for the gen-
eral case of a fragmented propagator is given in the appendix.
More details about the actual implementation of the Faddeev
formalism to 2p1h/2h1p propagation have been presented in
Ref. [24]. The calculation of the 2h1p component of R(ω)
follows completely analogous steps.
It is important to note that the present formalism includes
the effects of ph and pp/hh motion to be included simulta-
neously, while allowing interferences between these modes.
These excitations are evaluated here at the RPA level and are
then coupled to each other by solving Eqs. (7). This generates
diagrams as the one displayed in Fig. 3, with the caveat that
two phonons are not allowed to propagate at the same time.
Equations. (7) also assure that Pauli correlations are properly
taken into account at the 2p1h/2h1p level. In addition, one can
in principle employ dressed sp propagators in these equations
to generate a self-consistent solution. If we neglect the lad-
der propagator gII (ω) (5) in this expansion, we are left with
the ring series alone and the analogous physics ingredients as
for the generalized GW approach. However, this differs from
GGW due to the fact that no double counting of the second-
order self-energy occurs, since the Pauli exchanges between
the polarization propagator and the third line are properly ac-
counted for (see Fig. 3). Alternatively, one can suppress the
polarization propagator to investigate the effects of pp/hh lad-
ders alone.
It is instructive to replace in the above equations all RPA
phonons with TDA ones; this amounts to allowing only
forward-propagating diagrams in Fig. 2, and is equivalent to
separate diagonalisations in the spaces of ph, pp and hh con-
l 0 1 2 3 4 5 6
rw 2.0 4.0 0.0 0.0 0.0 0.0 0.0
no 12 21 10 10 5 5 5
TABLE I: Parameters that define the sp basis: radius of the confining
wall rw (in atomic units) and number of orbits no used for different
partial waves l. The value of cw is always set to 5 a.u..
figurations, relative to the HF ground state. It can be shown
that using these TDA phonons to sum all diagrams of the type
in Fig. 3 reduces to one single diagonalization in the 2p1h or
2h1p spaces. Therefore, Eqs. (7) and (8) with TDA phonons
lead directly to the “extended” 2p1h TDA of Ref. [23], which
was later shown to be equivalent to ADC(3) in the general
ADC framework [22]. The Faddeev expansion formalism of
Ref. [24] creates the possibility to go beyond ADC(3) by in-
cluding RPA phonons. This is more satisfactory in the limit
of large systems. At the same time, the computational cost
remains modest since only diagonalizations in the 2p1h/2h1p
spaces are required.
Note that complete self-consistency requires the use of
fragmented (or dressed) propagators in the evaluation of all in-
gredients leading to the self-energy. This is outside the scope
of the present paper, but we included partial selfconsistency
by taking into account the modifications to the HF diagram
by employing the correlated one-body density matrix and it-
erating to convergence. This is relatively simple to achieve,
since the 2p1h/2h1p propagator is only evaluated once with
the input HF propagators. Below we will give results with and
without this partial selfconsistency at the HF level.
III. RESULTS
Calculations have been performed using two different
model spaces: (1) a standard quantumchemical Gaussian ba-
sis set, aug-cc-pVTZ for Neon [30], with Cartesian repre-
sentation of the d and f functions; (2) a numerical basis set
based on HF and subsequent discretization of the continuum,
to be detailed below. The aug-cc-pVTZ basis set was used
primarily to check our formalism with the ADC(3) result in
literature (i.e. [31], where this basis was employed). The
HF+continuum basis allows to approach, at least for the ion-
ization energies, the results for the full sp space (basis set
limit).
The HF+continuum is the same discrete model space em-
ployed previously in Refs. [11, 14]. It consists of: (1) Solv-
ing on a radial grid the HF problem for the neutral atom; (2)
Adding to this fixed nonlocal HF potential a parabolic poten-
tial wall of the type U(r) = θ(r − rw)cw(r − rw)
2, placed at a
distance rw of the nucleus. The latter eigenvalue problem has
a basis of discrete eigenstates. This basis is truncated by spec-
ifying some largest angular momentum lmax and the number
of virtual states for each value of l ≤ lmax. (3) Solve the HF
problem again, without the potential wall, in this truncated
discrete space. The resulting basis set is used for the subse-
F-TDA F-RPA F-TDAc F-RPAc Expt.
2p -0.799 -0.791 -0.803 -0.797 (0.94) -0.793 (0.92)
2s -1.796 -1.787 -1.802 -1.793 (0.90) -1.782 (0.85)
1s -32.126 -32.087 -32.140 -32.102 (0.86) -31.70
Etot -128.778 -128.772 -128.836 -128.840 -128.928
TABLE II: Results with the aug-cc-pVTZ basis. The first three rows
list the energies of the main sp fragments below the Fermi level, as
predicted by different self-energies. F-TDA/F-RPA refers to the Fad-
deev summation with TDA/RPA phonons, respectively. In all cases
the self-energy was corrected at third order through Eq. (8). The
suffix “c” refers to partial selfconsistency, when the static (HF-type)
self-energy is consistent with the correlated density matrix. Without
“c” the pure HF self-energy was taken. In the F-RPAc column the
strength of the fragment is indicated between brackets. The last row
is the total electronic binding energy. The experimental values are
taken from Refs. [32, 33]. All energies are in atomic units.
quent Green’s function calculations.
When a sufficiently large number of states is retained after
truncation, the final results should approach the basis set limit.
In particular the results should not depend on the choice of the
auxiliary confining potential. This was verified in Ref. [11] for
the second-order, and in Ref. [14] for the G0W0 self-energy; in
these cases the self-energy is sufficiently simple that extensive
convergence checks can be made for various choices of the
auxiliary potential. The parameters of the confining wall and
the number of sp states kept in the basis set was optimized in
Ref. [11], by requiring that the ionization energy is converged
to about 1 mH for the second-order self-energy. In Ref. [14]
the same choice of basis set was also seen to bring the ioniza-
tion energy for the G0W0 self-energy near convergence. For
completeness, the details of this basis are reported in Table I.
While the self-energy in the present paper is too complicated
to allow similar convergence checks, it seems safe to assume
that basis set effects will affect the calculated ionization ener-
gies by at most 5 mH.
In Table II we compare, for the aug-cc-pVTZ basis, the ion-
ization energies of the main single-hole configurations when
TDA or RPA phonons are employed in the Faddeev construc-
tion (this is labeled F-TDA and F-RPA, respectively, in the
table). Note that use of TDA phonons corresponds to the
usual ADC(3) self-energy. We find that the replacement of
TDA with RPA phonons provides more screening, leading to
slightly less bound poles which are shifted towards the exper-
imental values. This shift increases with binding energy. As
discussed at the end of Sec. I, one can include consistency of
the static part of the self-energy. About eight iterations are
needed for convergence. This is a nonnegligible correction,
providing about 5 mH more binding (i.e. larger ionization en-
ergies) for the valence/subvalence 2p and 2s, 15 mH for the
deeply bound 1s, and 60 mH to the total binding energy. Our
converged result for the Faddeev-TDA self-energy (labeled F-
TDAc in Table II) is in good agreement with the ADC(3) value
for the 2p ionization energy (-0.804 H) quoted in [31], as it
should be.
The analogous results obtained with the larger
F-TDA F-RPA F-TDAc F-RPAc Expt.
2p -0.807 -0.799 -0.808 -0.801 (0.94) -0.793 (0.92)
2s -1.802 -1.792 -1.804 -1.795 (0.91) -1.782 (0.85)
1s -32.136 -32.097 -32.142 -32.104 (0.81) -31.70
Etot -128.863 -128.857 -128.883 -128.888 -128.928
TABLE III: Results with the HF+continuum basis set from Table I.
See also the caption of Table II.
-40 -35 -30
 [a.u.]
F-RPA(ring)
F-RPA(ladder)
F-RPA
FIG. 4: (Color online) Spectral function for the s states in Ne
obtained with various self-energy approximations. From the top
down: the second-order (Σ(2)) self-energy, the F-RPA(ring), the F-
RPA(ladder), and the full F-RPA self-energy. The strength is given
relative to the Hartree-Fock occupation of each shell. Only fragments
with strength larger than Z > 0.005 are shown.
HF+continuum basis are given in Table III, which al-
lows to assess overall stability and basis set effects. We find
exactly the same trends as for aug-cc-pVTZ. In particular
the reduction of ionization energies from the replacement of
TDA with RPA phonons is almost independent of the basis
set used, while the effect of including partial consistency
is roughly halved. Overall, the ionization states are always
more bound with the larger basis set; while the basis set
limit could be still more bound than the present results with
the HF+continuum basis set, it is likely (based on the G0W0
extrapolation in Ref. [14]) that the difference does not exceed
5 mH.
As discussed in Sec. I, the the F-RPA self-energy contains
RPA excitations of both ph type (ring diagrams) and pp/hh
type (ladder diagrams). It is instructive to analyze their sepa-
rate contributions to the final ionization energies, in order to
1s 2s 2p
HF -32.77 (1.00) -1.931 (1.00) -0.850 (1.00)
Σ(2) -31.84 (0.74) -1.736 (0.88) -0.747 (0.91)
G0W0 -31.14 (0.85) -1.774 (0.91) -0.801 (0.94)
F-RPA (ring) -31.82 (0.73) -1.636 (0.56) -0.730 (0.80)
F-RPA (ladder) -32.04 (0.87) -1.802 (0.95) -0.781 (0.96)
F-RPA -32.10 (0.81) -1.792 (0.91) -0.799 (0.94)
Exp. -31.70 -1.782 (0.85) -0.793 (0.92)
TABLE IV: Energy (in a.u.) and strength (bracketed numbers) of the
main fragments in the spectral function of Neon, generated by differ-
ent self-energies. Results for the HF+continuum basis. Consecutive
rows refer to: (1) HF; (2) second-order self-energy; (3) G0W0 results
from Ref. [14]; (4) F-RPA self-energy with only ph rings retained;
(5) F-RPA self-energy with only pp/hh ladders retained; (6) Com-
plete F-RPA self-energy. In all F-RPA results the self-energy was
corrected at third order through Eq. (8). The static self-energy was
pure HF (no partial self-consistency). The experimental values are
taken from Refs. [32, 33].
understand how the F-RPA self-energy is related to the stan-
dard (G)GW self-energy. Table IV compares the results for
the ionization energies, obtained with the second-order self-
energy, to different approximations for including the ring sum-
mations. As one can see, the second-order self-energy gener-
ates an l=1 sp energy of -0.747 mH, which is 46 mH above
the empirical 2p ionization energy. The G0W0 self-energy,
which includes the ring summation with only direct Coulomb
matrix elements, improves this result and brings it close to ex-
periment. The 2s behaves in a similar way. Unfortunately,
including the exchange terms of the interelectron repulsion in
the GG0W0 method turns out to have the opposite effect (the
2p ionization energy becomes -0.712 H [14] [41]) and the
agreement with experiment is lost. Obviously, GG0W0 is too
simplistic to account for exchange in the ph channel.
With the F-RPA(ring) self-energy one can go one step fur-
ther and employ the Faddeev expansion to also force proper
Pauli exchange correlations in the 2p1h/2h1p spaces. As
shown in Table IV, this enhances the screening due to the
exchange interaction terms, leading to even less binding for
the 2s and 2p. The corrections relative to the second-order
self-energy can be large (100 mH for the 2s state) and in the
direction away from the experimental value. We also note that
the larger shift, in the 2s orbit, is accompanied by an increase
of the fragmentation (see Fig. 4 and Tab. IV). Similar observa-
tions were also made in Ref. [14] for other atoms: in general
ring summations in the direct channel alone bring the quasi-
hole peaks close to the experiment. This agreement is then
spoiled as soon as one includes proper exchange terms in the
self-energy. On the other hand, exchange in the ph channel
is required to reproduce the correct Rydberg sequence in the
excitation spectrum of neutral atoms. So further corrections
must arise from other diagrams, and obviously the summation
of ladder diagrams can play a relevant role, since these con-
tribute to the expansion of the self-energy at the same level as
that of the ring diagrams.
The result when only including ladder-type RPA phonons
in the F-RPA self-energy is also shown in Table IV. One can
see that pp/hh ladders do actually work in the opposite way
as the ph channel ring diagrams, and have the same order of
magnitude with, e.g., a shift of 66 mH for the 2s relative to
the second-order result. When combined with the ring dia-
grams in the full F-RPA self-energy, the agreement with ex-
periment is restored again. Note that the final result cannot be
obtained by adding the contributions of rings and ladders, but
depends nontrivially on the interplay between these classes of
diagrams thereby pointing to significant interference effects.
With the F-RPA(ring) self-energy, where only the contribu-
tions of the ph channel are included, the main peaks listed in
Table IV are not only considerably shifted but also strongly
depleted, e.g. a strength of only 0.56 for the main 2s peak.
The complete spectral function for the l = 0 strength in Fig. 4
shows that the depletion of the main fragment is accompanied
by strong fragmentation over several states. While correlation
effects are overestimated in F-RPA(ring), they are suppressed
in F-RPA(ladder), where only the pp/hh ladders are included
in the self-energy. In this case one finds a spectral distribu-
tion closer to the HF one, with a main 2s fragment of strength
0.95 and less fragmentation than the the second-order self-
energy. The spectral distribution generated by the complete
F-RPA self-energy is again a combination of the above ef-
fects. The strength of the deeply bound 1s orbital behaves in
an analogous way. The strength of the main peak is reduced
but several satellite levels appear due to the mixing with 2h1p
configurations. In all the calculations reported in Fig. 4 we
found a summed l = 0 strength exceeding 0.98 in the interval
[-40 H, -30 H] which can be associated with the 1s orbital,
and this remains true even in the presence of strong correla-
tions using the F-RPA(ring) self-energy. Of course, the mix-
ing with 3h2p configurations, not included in this work, may
further contribute to the fragmentation pattern in this energy
region.
IV. CONLUSIONS AND DISCUSSION
In conclusion, the electronic self-energy for the Ne atom
was computed by the F-RPA method which includes –
simultaneously– the effects of both ring and ladder diagrams.
This was accomplished by employing an expansion of the
self-energy based on a set of Faddeev equations. This tech-
nique was originally proposed for nuclear structure applica-
tions [24] and is described in the appendix. At the level of
the self-energy one sums all diagrams where the three propa-
gator lines of the intermediate 2p1h or 2h1p propagation are
connected by repeated exchange of RPA excitations in both
the ph and the pp/hh channel. This differs from the ADC(3)
formalism in the fact that the exchanged excitations are of
the RPA type, rather than the TDA type, and therefore take
ground-state correlations effects into account. The coupling
to the external points of the self-energy uses the same modi-
fied vertex as in ADC(3), which must be introduced to include
consistently all third-order perturbative contributions.
The resulting main ionization energies in the Neon atom
are at least of the same quality, and even somewhat improved,
compared to the ADC(3) result. Note that, numerically, F-
RPA can be implemented as a diagonalization in 2p1h/2h1p
space implying about the same cost as ADC(3). The present
study also shows that in localized electronic systems subtle
cancellations occur between the ring and ladder series. In par-
ticular, only a combination of the ring and ladder series leads
to sensible results, as the separate ring series tends to correct
the second-order result in the wrong direction.
Since the limit to extended systems requires an RPA treat-
ment of excitations, the F-RPA method holds promise to
bridge the gap between an accurate description of quasiparti-
cles in both finite and extended systems. In particular, the GW
treatment of the electron gas has been shown to yield excellent
binding energies, but poor quasiparticle properties [34, 35].
Further progress beyond GW theory requires a consistent in-
corporation of exchange in the ph channel. The F-RPA tech-
nique may be highly relevant in this respect. A common
framework for calculating accurate QP properties in both fi-
nite and extended systems, is also important for constrain-
ing functionals in quasiparticle density functional theory (QP-
DFT) [7].
Finally, complete self-consistency requires sizable compu-
tational efforts for bases as large as the HF+continuum basis
used here. It would nevertheless represent an important exten-
sion of the present work, since it is related to the fulfillment
of conservation laws [36, 37]. These issues will be addressed
in future work.
Acknowledgments
This work was supported by the U.S. National Science
Foundation under grant PHY-0652900.
APPENDIX A: FADDEEV EXPANSION OF THE 2P1H/2H1P
PROPAGATOR
Although only the one-energy (or two-time) part of the
2p1h/2h1p propagator enters the definition of the self energy,
Eq. (3), a full resummation of all its diagrammatic contribu-
tions would require to treat explicitly the dependence on three
separate frequencies, corresponding to the three final lines in
the expansion of R(ω). For example, inserting the RPA ring
(ladder) series in R(ω) implies the propagation of a ph (pp/hh)
pair of lines both forward and backward in time, while the
third line remains unaffected. A way out of this situation is
to solve the Bethe-Salpeter-like equations for the polarization
and ladder propagators separately and then to couple them to
the additional line. If it is assumed that different phonons do
not overlap in time, the three lines in between phonon struc-
tures will propagate only in one time direction [see figures (3)
and (5)]. In this situation the integration over several frequen-
cies can be circumvented following the prescription detailed
in the next subsection. This approach will be discussed in
the following for the general case of a fully fragmented prop-
gator, in order to derive a set of Faddeev equations capable
FIG. 5: (Color online) Diagrammatic representation of Eq. (A3).
Double lines represent fully dressed sp Green’s funcions which, how-
ever, are restricted to propagate only in one time direction [i.e., only
one of the two terms on the r.h.s. of Eq. (1) is retained]. The Faddeev
Eqs. (A9) and (7) allow for both forward and backward propagation
of the phonons Γ(π)(ω) and Γ(II)(ω) as long as these do not overlap
in time. For the propagators, time ordereing is asumed with forwad
propagation in the upward direction.
of dressing the sp propagator self-consistently. Since the for-
ward (2p1h) and the backward (2h1p) parts of R(ω) decouple
in two analogous sets of equations, it is sufficient to focus on
the first case alone.
1. Multiple frequencies integrals
We start by considering the effective interactions in the ph
and pp/hh channels that correspond to Eqs. (4) and (5) stripped
of the external legs. In the present work, these are the follow-
ing two-time objects:
αβ,γδ
(ω) = Vαδ,βγ + Vαν,βµ Π
µν,ρσ(ω) Vρδ,σγ (A1a)
= Vαδ,βγ +
ω − επn + iη
ω + επn′ − iη
αβ,γδ
(ω) = Vαβ,γδ + Vαβ,µν g
µν,ρσ(ω) Vρσ,γδ (A1b)
= Vαδ,βγ +
ω − εΓ+n + iη
ω − εΓ−
where the residues and poles for the ring series are Ωn
〈ΨNn |c
µcν|Ψ
0 〉Vµβ,να and ε
n = E
0 . For the ladders, ∆
〈ΨN+2n |c
0 〉Vµν,αβ and ∆
= Vαβ,µν〈Ψ
k |cµcν|Ψ
0 〉, with
poles εΓ+n = E
0 and ε
k = E
. Equations. (A1)
solve the ring and ladder RPA equations, respectively
αβ,γδ
(ω) = Vαδ,βγ (A2a)
αβ,µν
gµρ(ω + ω1)gσν(ω1) Vρδ,σγ ,
αβ,γδ
(ω) = Vαβ,γδ (A2b)
αβ,µν
gµρ(ω − ω1)gνσ(ω1) Vρσ,γδ .
To display how the phonons (A1a) and (A1b) enter the ex-
pansion of R(ω), we perform explicitly the frequency integrals
for the diagram of Fig. 5. Since it is assumed that the separate
propagators lines evolve only in one time direction, only the
forwardgoing (g>(ω)) or the backwardgoing (g<(ω)) part of
Eq. (1) must be included for particles and holes, respectively.
After some algebra, one obtains
∆Rαβγ,µνλ(ω) =
g>αα1 (ω −Ω) g
(ω1) g
(ω1 −Ω) Γ
β1γ1,σ1λ1
(Ω) g>σ1σ2 (s + Ω − ω)
α1σ2,µ1ν1
(s) g>µ1µ(s − ω2) g
(ω2) g
(s − ω)
ω − (ε+n1 + ε
) + iη
Vβ1λ1,γ1σ1 +
ω − (ε+n1 + ε
nπ) + iη
[ω − επ
− ε+n1 − ε
+ ε−k3
− ε+n4 + ε
− ε+n2 + ε
][−επ
− ε+n4 + ε
ω − (ε+n1 + ε
) + iη
Vα1σ2,µ1ν1 +
+,nII
+,nII
ω − (εΓ+nII − ε
) + iη
[ω + εΓ−kII − ε
− ε+n4 − ε
− ε+n6 + ε
] ∆−,kIIα1σ2
−,kII
− ε+n1 − ε
][εΓ−
− ε+n5 − ε
ω − (ε+n5 + ε
− ε−k7
) + iη
ω − (εΓ−
− ε+n4 − ε
) − iη
−,kII
−,kII
− ε+n2 + ε
] [−επ
− ε+n4 + ε
] [εΓ−
− ε+n1 − ε
] [εΓ−
− ε+n5 − ε
The last term in this expression contains an energy denom-
inator that involves the simultaneous propagation of two
phonons. Thus, it will be discarded in accordance with our
assumptions. It must be stressed that similar terms, with over-
lapping phonons, imply the explicit contribution of at least
3p2h/3h2p. A proper treatment of these would require a non
trivial externsion of the present formalism, which is beyond
the scope of this paper.
The remaining part in Eq. (A3) is the relevant contribution
for our purposes. This has the correct energy dependence of
a product of denominators that correspond to the intermediate
steps of propagation. All of these involve configurations that
have at most 2p1h character. Although, ground state corre-
lations are implicitely included by having already resummed
the RPA series. Still, this term does not factorize in a prod-
uct of separate Green’s functions due to the summations over
the fragmentation indices ni and ki [labeling the eigenstates of
the (N±1)-electron systems]. This is overcome if one defines
the matrices G0>(ω), Γ(1,2)(ω) and Γ(3)(ω), with elements (no
implicit summation used)
αnαβnβγkγ; µnµνnνλkλ
(ω) = δnα,nµ δnβ,nν δkγ,kλ
ω − (ε+nα + ε
− ε−kγ
) + iη
, (A4a)
αnαβnβγkγ; µnµνnνλkλ
(ω) = Γ(2)>
βnβαnαγkγ; νnνµnµλkλ
(ω) =
δα,µ δnα ,nµ
Vβλ,γν +
ω − (ε+nα + ε
nπ) + iη
[ω − επ
− ε+nα − ε
+ ε−kγ
− ε+nν + ε
− ε+nβ + ε
][−επ
− ε+nν + ε
, (A4b)
αnαβnβγkγ; µnµνnνλkλ
(ω) =
δγ,λ δkγ ,kλ
Vαβ,µν +
+,nII
+,nII
ω − (εΓ+nII − ε
) + iη
[ω + εΓ−kII − ε
− ε+nβ − ε
− ε+nν + ε
] ∆−,kII
−,kII
− ε+nα − ε
][εΓ−
− ε+nµ − ε
. (A4c)
In these definitions, the row and column indices are ordered to
represent at first two quasiparticle lines and then a quasihole.
The index ‘i’ in Γ(i)> refer to the line that propagates indepen-
dently along with the phonon. Using Eqs.(A4), the first term
on the r.h.s. of Eq. (A3) can be written as
∆R(2p1h)
αβγ,µνλ
(ω) = (A5)
nα nβ kγ
nµ nν kλ
G0>(ω)Γ(1)>(ω)G0>(ω)Γ(3)>(ω)G0>(ω)
αnαβnβγkγ; µnµνnνλkλ
Eq. (A5) generalizes to diagrams involving any number of
phonon insertions, as long as the terms involving two or more
simultaneous phonons are dropped. Based on this relation, we
use the following prescription to avoid performing integrals
over frequencies. One extends all the Green’s functions to ob-
jects depending not only on the sp basis’ indices (α, β, γ) but
also on the indices labeling quasi-particles and holes (ni and
ki). Whether a given argument represents a particle or an hole
depends on the type of line being propagated. At this point one
can perform calculations working with only two-time quanti-
ties. The standard propagator is recovered at the end by sum-
ming the “extended” one over the quasi-particle/hole indices.
2. Faddeev expansion
The 2p1h/2h1p propagator that includes the full resumma-
tion of both the ladder and ring diagrams at the (G)RPA level
is the solution of the following Bethe-Salpeter-like equation,
Rαβγ,µνλ(ω1, ω2, ω3) = (A6)
gαµ(ω1)gβν(ω2) − gβµ(ω2)gαν(ω1)
gλγ(−ω3) +
gββ1(ω2)gγ1γ(−ω3)Vβ1σ,γ1ρ
Rαρσ,µνλ(ω1, s, ω2 + ω3 − s)
+ gαα1 (ω1)gγ1γ(−ω3)Vα1σ,γ1ρ
Rρβσ,µνλ(s, ω2, ω1 + ω3 − s)
gαα1 (ω1)gββ1(ω2)Vα1β1,ρσ
Rρσγ,µνλ(s, ω1 + ω2 − s, ω3)
If this equation is solved, a double integration of R(ω1, ω2, ω3)
would yield the two-time propagator R(ω) contributing to
Eq. (3). However, the numerical solution of Eq. (A6) appears
beyond reach of the present day computers and one needs
to avoid dealing directly with multiple frequencies integrals.
The strategy used is to first solve the RPA equations (A2a)
and (A2b) separately. Once this is done it is necessary to re-
arrange the series (A6) in such a way that only the resummed
phonons appear. Following the formalism introduced by Fad-
deev [29, 38], we identify the components R(i)(ω) with the
three terms between curly brakets in Eq. (A6). By employ-
ing Eqs. (A2a) and (A2b) one is lead to the following set of
equations [42],
αβγ,µνλ
(ω1, ω2, ω3) = gαα1 (ω1)gββ1(ω2)gγ1γ(−ω3) (A7)
ds1 ds2 ds3
α1β1γ1,µ1ν1λ1
(ω1, ω2, ω3; s1, s2, s3)
gµ1µ(s1)gν1ν(s2) − gν1µ(s2)gµ1ν(s1)
gλλ1(−s3)
µ1ν1λ1,µνλ
(s1, s2, s3) + R
µ1ν1λ1,µνλ
(s1, s2, s3)
, i = 1, 2, 3 ,
where (i,j,k) are cyclic permutations of (1,2,3) and the inter-
action vertices Γ(i)(ω1, ω2, ω3) are given by
αβγ,µνλ
(ω1, ω2, ω3;ω4, ω5, ω6) = (A8a)
βαγ,νµλ
(ω2, ω1, ω3;ω5, ω4, ω6) =
= δ(ω1 − ω4)δ(ω2 + ω3 − ω5 − ω6)g
αµ(ω1)Γ
βγ,νλ
(ω2 + ω3) ,
αβγ,µνλ
(ω1, ω2, ω3;ω4, ω5, ω6) = (A8b)
δ(ω3 − ω6)δ(ω1 + ω2 − ω4 − ω5)g
λγ (−ω3)Γ
αβ,µν
(ω1 + ω2) .
Finally, we apply the prescription of Sec. A 1 and substitute
R(ω1, ω2, ω3) with its extended but two-time version R(ω).
This leads to the following set of Faddeev equations which
propagate 2p1h forward in time,
R̄(i)
αnαβnβγkγ; µnµνnνλkλ
(ω) =
αnαβnβγkγ; α′n′αβ′n
γ′k′γ
(ω) Γ(i)
α′n′αβ′n
γ′k′γ; µ′n
′n′νλ′k
αnαβnβγkγ; µnµνnνλkλ
(ω) −G0
αnαβnβγkγ; µnµνnνλkλ
× R̄( j)
µ′n′µν′n
; µnµνnνλkλ
(ω) + R̄(k)
µ′n′µν′n
; µnµνnνλkλ
i = 1, 2, 3 . (A9)
Since the full energy dependence is retained in Eq. (A7),
the self-energy corresponding to its solution, R(ω1, ω2, ω3), is
complete up to third order [see Eq. (3)]. This is no longer the
case after the reduction to a two-time propagator. In particu-
lar, the approximation that only forward 2p1h propagation is
allowed between different phonons implies that all diagrams
with different time propagation of their external lines are ne-
glected in Eqs. (A9). However, these terms are not energy
dependent and can be can be reinserted in a systematic way a
posteriori as in Eq. (8). In the general case,
Rαβγ,µνλ(ω) = (A10)
U (2p1h)
αβγ; α′n′αβ′n
γ′k′γ
R̄(2p1h)
α′n′αβ′n
γ′k′γ; µ′n
′n′νλ′k
(ω) U (2p1h) †
µ′n′µν′n
; µνλ
U (2p1h)
αβγ; µnµνnνλkλ
= δαµδβνδγλ + ∆U
(2p1h)
αβγ; µnµνnνλkλ
, (A11)
where the correction ∆U can be determined by comparison
with perturbation theory.
The vertices (A4), that appear in Eqs. (A9), and U (2p1h) are
expressed in terms of the fully fragmented propagator. There-
fore, this approach allows to obtain self-consistent solutions
of the sp Green’s function [25]. Whenever, like in this work,
only a mean-field propagator is employed as input there ex-
ist a one-to-one correspondence between the fragmentation
indices and the sp basis. This is expressed by the relations
Xnα = δn,α(1 − δα∈F) and Y
α = δk,αδα∈F , where F represents
the set of occupied orbits. In this case, it is possible to drop
one set of indices so that Eqs. (A9) and (A10) simplify into
the form (7) and (8).
3. Faddeev vertices
In practical applications, it is worth to note that the poles
of the free propagator G0(ω), Eq. (A4a), do not contribute to
the kernel of Eqs. (A9). This can be proven by employing the
closure relations for the RPA problem, in the form obtained
by extracting the free poles in Eqs. (A2). As an example, for
the forward poles of the ladder propagator these are
ω→ε+n1
+ε+n2
[(ω − ε+n1 − ε
) × (Eq. A2b)] =⇒ (A12)
µν,γδ
(ω = ε+n1 + ε
) = 0 , ∀n1, n2 ,
and similarly for other cases. Making use of these relations
one can derive the following working expression of the ker-
nels of the 2p1h Faddeev equations (no implicit summations
used)
G0>(ω)Γ(1)>(ω)
αnαβnβγkγ; µnµνnνλkλ
G0>(ω)Γ(2)>(ω)
βnβαnαγkγ; νnνµnµλkλ
= (A13a)
= δnα,nµ
[επnπ − ε
][ω − (ε+nα + ε
nπ) + iη]
− ε+nβ + ε
][−επ
− ε+nν + ε
G0>(ω)Γ(3)>(ω)
αnαβnβγkγ; µnµνnνλkλ
= (A13b)
= δkγ,kλ
+,nII
+,nII
[εΓ+nII − ε
− ε+nβ][ω − (ε
) + iη]
−,kII
−,kII
− ε+nα − ε
][εΓ−
− ε+nµ − ε
After substituting Eq. (A10) into (3), one needs the working expression for the matrix product V U(2p1h) (where V is the inter-
electron interaction). The minimum correction that guaranties to reproduce all third order self-energy diagrams is
V U(2p1h)
α; µnµνnνλkλ
= Vαλ,µν +
Vαλ,γ1δ1 Y
Vγ2δ2,µν
2 [ε−
− ε+nµ − ε
(A14)
Vαδ1,µγ1 Y
Vγ2λ,δ2ν
− ε+nδ − ε
Vαδ1,νγ1 Y
Vγ2λ,δ2µ
− ε+nδ − ε
The case of 2h1p is handled in a completely analogous way along the steps of Secs. (A 1) and (A 2). After extending
R(ω1, ω2, ω3) to depend on the fragmentation indices (k1,k2,n), the 2h1p equivalent of Eq. (A9) is obtained with the follow-
ing definitions of the kernels,
G0>(ω)Γ(1)>(ω)
αkαβkβγnγ; µkµνkνλnλ
G0>(ω)Γ(2)>(ω)
βkβαkαγnγ; νkνµkµλnλ
= (A15a)
= δkα,kµ
+ ε+nγ ][ω − (ε
) − iη]
[επnπ − ε
+ ε+nγ ][ε
nπ − ε
+ ε+nλ]
G0>(ω)Γ(3)>(ω)
αkαβkβγnγ; µkµνkνλnλ
= (A15b)
= δnγ,nλ
−,nII
−,nII
][ω − (εΓ−
− ε+nγ ) − iη]
+,kII
+,kII
[εΓ+nII − ε
− ε−kβ
][εΓ+nII − ε
− ε−kν
and correction to the external legs,
V U(2h1p)
α; µkµνkνλnλ
= Vαλ,µν +
Vαλ,γ1δ1 X
Vγ2δ2,µν
2 [ε−kµ + ε
− ε+nγ − ε
(A16)
Vαδ1,µγ1 X
Vγ2λ,δ2ν
[ε−kδ + ε
− ε+nγ − ε
Vαδ1,νγ1 X
Vγ2λ,δ2µ
[ε−kδ + ε
− ε+nγ − ε
It should be pointed out that while the prescription of
Sec. A 1 allows sp lines to propagate only in one time di-
rection, it allows for backward propagation of the phonons.
These contributions translate directly into the energy indepen-
dent terms of Eqs. (A13) and (A15) and are a direct conse-
quence of the inversion pattern typical of RPA theory. These
terms have normally a weaker impact than the direct ones on
the solutions of Eqs. (A9). However, it is show in Ref. [24]
that they are crucial to guarantee the exact separation of the
spurious solutions—always introduced by the Faddeev for-
malism [39, 40]—if RPA phonons are used. For the same rea-
sons, the last terms in curly brackets of Eqs. (A13) and (A15)
should be dropped whenever Tamm-Dancoff (TDA) phonons
are propagated.
The approach followed in this work for solving Eqs. (A9)
is to transform them into a matrix representations [24].
Once this is done, one is left with an eigenvalue prob-
lem that depends only on the 2p1h (2h1p) configurations
(n, n′, k) [(k, k′, n)]. The spurious states are known ex-
actly [24] and can be projected out analytically to reduce the
computational load. In any case, they would give vanishing
contributions to Eq. (3).
[1] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
[2] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
[3] A. Görling and M. Levy, Phys. Rev. A 50, 196 (1994).
[4] A. Görling, J. Chem. Phys. 123, 062203 (2005).
[5] R. J. Bartlett et al., J. Chem. Phys. 122, 034104 (2005); J.
Chem. Phys. 123, 062205 (2005).
[6] P. Mori-Sanchez, Q. Wu, and W. T. Yang, J. Chem. Phys. 123,
062204 (2005).
[7] D. Van Neck, S. Verdonck, G. Bonny, P. W. Ayers, and M.
Waroquier, Phys. Rev. A 74, 042501(2006).
[8] A. B. Migdal, Theory of Finite Fermi Systems and Applications
to Atomic Nuclei (John Wiley and Sons, New York, 1967).
[9] A. L. Fetter and J. D. Walecka, Quantum Theory of Many-
Particle Physics (McGraw-Hill, New York, 1971).
[10] W. H. Dickhoff and D. Van Neck, Many-Body Theory Exposed!
(World Scientific, Singapore, 2005).
[11] D. Van Neck, K. Peirs, and M. Waroquier, J. Chem. Phys. 115,
15 (2001).
[12] K. Peirs, D. van Neck, and M. Waroquier, J. Chem. Phys. 117,
4095 (2002).
[13] N. E. Dahlen and R. van Leeuwen, J. Chem. Phys. 122, 164102
(2005).
[14] S. Verdonck, D. Van Neck, P. W. Ayers, M. Waroquier, Phys.
Rev. A 74, 062503 (2006).
[15] E.L. Shirley and R.M. Martin, Phys. Rev. B 47, 15404 (1993).
[16] N.E. Dahlen, R. van Leeuwen, and U. von Barth, Phys. Rev. A
73, 012511 (2006).
[17] N.E. Dahlen and U. von Barth, Phys. Rev. B 69, 195102 (2004).
[18] A. Stan, N.E. Dahlen and R. van Leeuwen, Europhys. Lett. 76,
298 (2006).
[19] L. Hedin, Phys. Rev. 139, A796 (1965).
[20] N. Fukuda, F. Iwamoto, and K. Sawada, Phys. Rev. 135, A932
(1964).
[21] C. Barbieri, N. Paar, R. Roth, and P. Papakostantinou,
arXiv:nucl-th/0608011v1.
[22] J. Schirmer, L.S. Cederbaum, and O. Walter, Phys. Rev. A 28,
1237 (1983)
[23] O. Walter and J. Schirmer, J. Phys. B:At. Mol. Phys. 14, 3805
(1981).
[24] C. Barbieri and W. H. Dickhoff, Phys. Rev. C 63, 034313
(2001).
[25] C. Barbieri and W. H. Dickhoff, Phys. Rev. C 65, 064313
http://arxiv.org/abs/nucl-th/0608011
(2002).
[26] W. H. Dickhoff and C. Barbieri, Prog. Part. Nucl. Phys. 52, 337
(2004).
[27] C. Barbieri, Phys. Lett. B 643, 268 (2006).
[28] W. J. W. Geurts, K. Allaart, and W. H. Dickhoff, Phys. Rev. C
50, 514 (1994).
[29] L. D. Faddeev, Zh. Éksp. Teor. Fiz. 39 1459 (1961) [Sov. Phys.
JETP 12, 1014 (1961)].
[30] T.H. Dunning Jr., J. Chem. Phys. 90, 1007 (1989); D.E. Woon
and T.H. Dunning Jr, J. Chem. Phys. 98, 1358 (1993)
[31] A.B. Trofimov and J. Schirmer, J. Chem. Phys. 123, 144115
(2005)
[32] NIST Atomic Spectra Database, NIST
Standard Reference Database #78,
http://physics.nist.gov/PhysRefData/ASD/in-dex.html
[33] A. Thompson et al., X-ray Data Booklet (Lawrence Berkeley
National Laboratory, Berkeley, CA, 2001), and references cited
therein.
[34] U. von Barth and B. Holm, Phys. Rev. B 54, 8411 (1996); B.
Holm and U. von Barth, Phys. Rev. B 57, 2108 (1998); B.
Holm, Phys. Rev. Lett. 83, 788 (1999).
[35] P. Garcia-Gonzales and R.W. Godby, Phys. Rev. B 63, 075112
(2001)
[36] G. Baym and L.P. Kadanoff, Phys. Rev. 124, 287 (1961).
[37] G. Baym, Phys. Rev. 127, 1391 (1962).
[38] C. J. Joachain, Quantum Collision Theory, (North-Holland,
Amsterdam, 1975).
[39] S. K. Adhikari and W. Glöckle, Phys. Rev. C 19, 616 (1979).
[40] J. W. Evans and D. K. Hoffman, J. Math. Phys. 22, 2858 (1981).
[41] Note that the G0W0 and the GG0W0 results of Ref. [14] were
obtained by retaining only the diagonal part of the electron self-
energy, Σ⋆αα(ω) in the HF+continuum basis. This approximation
was not made in the present work. The error in the ionization
energies by retaining the diagonal approximation is quite small
(about ∼2 mH [14] for the Ne atom), but larger effects are pos-
sible for the total binding energy.
[42] Note that the present definitions of the R(i) differ from the ones
of Ref. [24] which contain the additional term 1
[G0 − G0,ex].
The two different forms of the Faddeev equations that result can
be easily related into each other and are completely equivalent.
The definition used here agrees with the standard literature on
the subject [29, 38].
http://physics.nist.gov/PhysRefData/ASD/in-dex.html

  The spectral function of the closed-shell Neon atom is computed by expanding
the electron self-energy through a set of Faddeev equations. This method
describes the coupling of single-particle degrees of freedom with correlated
two-electron, two-hole, and electron-hole pairs. The excitation spectra are
obtained using the Random Phase Approximation, rather than the Tamm-Dancoff
framework employed in the third-order algebraic diagrammatic contruction
[ADC(3)] method. The difference between these two approaches is studied, as
well as the interplay between ladder and ring diagrams in the self-energy.
Satisfactory results are obtained for the ionization energies as well as the
energy of the ground state with the Faddeev-RPA scheme that is also appropriate
for the high-density electron gas.

Quasiparticles in Neon using the Faddeev Random Phase Approximation
C. Barbieri
Gesellschaft für Schwerionenforschung, Planckstr. 1, D-64291, Darmstadt, Germany
D. Van Neck
Laboratory of Theoretical Physics, Ghent University, Proeftuinstraat 86, B-9000 Gent, Belgium
W.H. Dickhoff
Department of Physics, Washington University, St. Louis, MO 63130, USA
(Dated: November 2, 2018)
The spectral function of the closed-shell Neon atom is computed by expanding the electron self-energy
through a set of Faddeev equations. This method describes the coupling of single-particle degrees of free-
dom with correlated two-electron, two-hole, and electron-hole pairs. The excitation spectra are obtained using
the Random Phase Approximation, rather than the Tamm-Dancoff framework employed in the third-order al-
gebraic diagrammatic contruction [ADC(3)] method. The difference between these two approaches is studied,
as well as the interplay between ladder and ring diagrams in the self-energy. Satisfactory results are obtained
for the ionization energies as well as the energy of the ground state with the Faddeev-RPA scheme that is also
appropriate for the high-density electron gas.
PACS numbers: 31.10.+z,31.15.Ar
I. INTRODUCTION
Ab initio treatments of electronic systems become unwork-
able for sufficiently complex systems. On the other hand,
the Kohn-Sham formulation [1] of density functional the-
ory (DFT) [2] incorporates many-body correlations (beyond
Hartree-Fock), while only single-particle (sp) equations must
be solved. Due to this simplicity DFT is the only feasible
approach in some modern applications of electronic structure
theory. There is therefore a continuing interest both in devel-
oping new and more accurate functionals and in studying con-
ceptual improvements and extensions to the DFT framework.
In particular it is found that DFT can handle short-range inter-
electronic correlations quite well, while there is room for im-
provements in the description of long-range (van der Waals)
forces and dissociation processes.
Microscopic theories offer some guidance in the devel-
opment of extensions to DFT. Orbital dependent function-
als can be constructed using many-body perturbation theory
(MBPT) [3, 4]. More recently, the development of general
ab initio DFT [5, 6] addressed the lack of a systematic im-
provement in DFT methods. In this approach one considers
an expansion of the exact ground-state wave function (e.g.,
MBPT or coupled cluster) from a chosen reference determi-
nant. Requiring that the correction to the density vanishes at a
certain level of perturbation theory allows one to construct the
corresponding approximation to the Kohn-Sham potential.
A different route has been proposed in Ref. [7] by devel-
oping a quasi-particle (QP)-DFT formalism. In the QP-DFT
approach the full spectral function is decomposed in the con-
tribution of the QP excitations, and a remainder or background
part. The latter part is complicated, but does not need to be
known accurately: it is sufficient to have a functional model
for the energy-averaged background part to set up a single-
electron selfconsistency problem that generates the QP exci-
tations. Such an approach is appealing since it contains the
well-developed standard Kohn-Sham formulation of DFT as
a special case, while at the same time emphasis is put on the
correct description of QPs, in the Landau-Migdal sense [8].
Hence, it can provide an improved description of the dynam-
ics at the Fermi surface. Given the close relation between QP-
DFT and the Green’s function (GF) formulation of many-body
theory [9, 10], it is natural to employ ab initio calculations
in the latter formalism to investigate the structure of possible
QP-DFT functionals. In this respect it is imperative to iden-
tify which classes of diagrams are responsible for the correct
description of the QP physics.
Some previous calculations, based on GF theory, have fo-
cused on a self-consistent treatment of the self-energy at the
second order [11, 12, 13] for simple atoms and molecules. For
the atomic binding energies it was found that the bulk of corre-
lations, beyond Hartree-Fock, are accounted for while signifi-
cant disagreement with experiment persists for QP properties
like ionization energies and electron affinities. The formal-
ism beyond the second-order approximation was taken up in
Ref. [14, 15, 16, 17, 18] by employing a self-energy of the
GW type [19]. In this approach, the random phase approx-
imation (RPA) in the particle-hole (ph) channel is adopted
to allow for possible collective effects on the atomic excited
states. The latter are coupled to the sp states by means of dia-
grams like the last two in Fig. 1(c). Two variants of the G0W0
formalism were employed in Ref. [14] (where the subscript
“0” indicates that non-dressed propagators are used). In the
first only the direct terms of the interelectron Coulomb po-
tential are taken into account. In the second version, also the
exchange terms are included when diagonalizing the ph space
[generalized RPA (GRPA)] and in constructing the self-energy
[generalized GW (GGW)]. Although the exchange terms are
known to be crucial in order to reproduce the experimentally
observed Rydberg sequence in the excitation spectrum of neu-
tral atoms, they were found to worsen the agreement between
the theoretical and experimental ionization energies [14].
http://arxiv.org/abs/0704.1542v2
In the GW approach the sp states are directly coupled
with the two-particle–one-hole (2p1h) and the two-hole–one-
particle (2h1p) spaces. However, only partial diagonalizations
(namely, in the ph subspaces) are performed. This procedure
unavoidably neglects Pauli correlations with the third particle
(or hole) outside the subspace. In the case of the GGW ap-
proach, this leads to a double counting of the second order
self-energy which must be corrected for explicitly [20, 21].
We note that simply subtracting the double counted diagram is
not completely satisfactory here, since it introduces poles with
negative residues in the self-energy. More important, the inter-
action between electrons in the two-particle (pp) and two-hole
(hh) subspaces are neglected altogether in (G)GW. Clearly, it
is necessary to identify which contributions, beyond GGW,
are needed to correctly reproduce the QP spectrum.
In this respect, it is known that highly accurate descrip-
tions of the QP properties in finite systems can be obtained
with the algebraic diagrammatic construction (ADC) method
of Schirmer and co-workers [22]. The most widely used third-
order version [ADC(3)] is equivalent to the so-called extended
2p1h Tamm-Dancoff (TDA) method [23] and allows to pre-
dict ionization energies with an accuracy of 10-20 mH in
atoms and small molecules. Upon inspection of its diagram-
matic content, the ADC(3) self-energy is seen to contain all
diagrams where TDA excitations are exchanged between the
three propagator lines of the intermediate 2p1h or 2h1p prop-
agation. The TDA excitations are constructed through a diag-
onalization in either 2p1h or 2h1p space, and neglect ground-
state correlations. However, it is clear that use of TDA leads
to difficulties for extended systems. In the high-density elec-
tron gas e.g., the correct plasmon spectrum requires the RPA
in the ph channel, rather than TDA.
In order to bridge the gap between the QP description in
finite and extended systems, it seems therefore necessary to
develop a formalism where the intermediate excitations in the
2p1h/2h1p propagator are described at the RPA level. This
can be achieved by a formalism based on employing a set of
Faddeev equations, as proposed in Ref. [24] and subsequently
applied to nuclear structure problems [25, 26, 27]. In this ap-
proach the GRPA equations are solved separately in the ph and
pp/hh subspaces. The resulting polarization and two-particle
propagators are then coupled through an all-order summation
that accounts completely for Pauli exchanges in the 2p1h/2h1p
spaces. This Faddeev-RPA (F-RPA) formalism is required if
one wants to couple propagators at the RPA level or beyond.
Apart from correctly incorporating Pauli exchange, F-RPA
takes the explicit inclusion of ground-state correlations into
account, and can therefore be expected to apply to both finite
and extended systems. The ADC(3) formalism is recovered
as an approximation by neglecting ground-state correlations
in the intermediate excitations (i.e. replacing RPA with TDA
phonons).
In this work we consider the Neon atom and apply the F-
RPA method to a nonrelativistic electronic problem for the
first time. The relevant features of the F-RPA formalism (also
extensively treated in Ref. [24]), are introduced in Sect. II.
The application to the Neon atom is discussed in Sec. III,
where we also investigate the separate effects of the ladder
and ring series on the self-energy, as well as the differences
between including TDA and RPA phonons. Our findings are
summarized in Sec. IV. Some more technical aspects are rele-
gated to the appendix, where the interested reader can find the
derivation of the Faddeev expansion for the 2p1h/2h1p propa-
gator, adapted from Ref. [24]. In particular, the approach used
to avoid the multiple-frequency dependence of the Green’s
functions is discussed in App. A 1, along with its basic as-
sumptions. The explicit expressions of the Faddeev kernels
are given in App. A 3. Together with Ref. [24], the appendix
provides sufficient information for an interested reader to ap-
ply the formalism.
II. FORMALISM
The theoretical framework of the present study is that of
propagator theory, where the object of interest is the sp prop-
agator, instead of the many-body wave function. In this paper
we will employ the convention of summing over repeated in-
dices, unless specified otherwise. Given a complete orthonor-
mal basis set of sp states, labeled by α,β,..., the sp propagator
can be written in its Lehmann representation as [9, 10]
gαβ(ω) =
ω − ε+n + iη
ω − ε−k − iη
, (1)
where Xnα = 〈Ψ
0 〉 (Y
α = 〈Ψ
|cα|Ψ
0 〉) are the
spectroscopic amplitudes, cα (c
) are the second quantization
destruction (creation) operators and ε+n = E
n − E
EN0 − E
). In these definitions, |ΨN+1n 〉, |Ψ
〉 are the
eigenstates, and EN+1n , E
k the eigenenergies of the (N ± 1)-
electron system. Therefore, the poles of the propagator reflect
the electron affinities and ionization energies.
The sp propagator solves the Dyson equation
gαβ(ω) = g
αβ(ω) + g
αγ(ω)Σ
γδ(ω)gδβ(ω) , (2)
which depends on the irreducible self-energy Σ⋆(ω). The lat-
ter can be written as the sum of two terms
αβ(ω) = Σ
Vαλ,µν Rµνλ,γδε(ω) Vγδ,βε , (3)
where ΣHF represents the Hartree-Fock diagram for the self-
energy. In Eqs. (2) and (3), g0(ω) is the sp propagator for
the system of noninteracting electrons, whose Hamiltonian
contains only the kinetic energy and the electron-nucleus at-
traction. The Vαβ,γδ represent the antisymmetrized matrix el-
ements of the interelectron (Coulomb) repulsion. Note that
in this work we only consider antisymmetrized elements of
the interaction, hence, our result for the ring summation al-
ways compare to the generalized GW approach. Equation (3)
introduces the 2p1h/2h1p-irreducible propagator R(ω), which
carries the information concerning the coupling of sp states to
more complex configurations. Both Σ⋆(ω) and R(ω) have a
perturbative expansion as a power series in the interelectron
FIG. 1: (Color online) a) Diagrammatic expansion of R(ω) in terms
of the (antisymmetrized) Coulomb interaction and undressed prop-
agators. b) R(ω) is related to the self-energy according to Eq. (3).
c) By substituting the diagrams a) in the latter equation, one finds the
perturbative expansion of the self-energy.
interaction V̂ . Some of the diagrams appearing in the expan-
sion of R(ω) are depicted in Fig. 1, together with the corre-
sponding contributions to the self-energy. Note that already at
zero order in R(ω) (three free lines with no mutual interaction)
the second order self-energy is generated.
Different approximations to the self-energy can be con-
structed by summing particular classes of diagrams. In this
work we are interested in the summation of rings and ladders,
through the (G)RPA equations. In order to include such effects
in R(ω), we first consider the polarization propagator describ-
ing excited states in the N-electron system
Παβ,γδ(ω) =
〈ΨN0 |c
n 〉 〈Ψ
γcδ|Ψ
ENn − E
〈ΨN0 |c
γcδ|Ψ
n 〉 〈Ψ
ENn − E
, (4)
and the two-particle propagator, that describes the addi-
tion/removal of two electrons
gIIαβ,γδ(ω) =
〈ΨN0 |cβcα|Ψ
n 〉 〈Ψ
|ΨN0 〉
EN+2n − E
〈ΨN0 |c
|ΨN−2
〉 〈ΨN−2
|cβcα|Ψ
EN0 − E
. (5)
We note that the expansion of R(ω) arises from applying the
equations of motion to the sp propagator (1), which is associ-
ated to the ground state |ΨN0 〉. Hence, all the Green’s functions
appearing in this expansion will also be ground state based, in-
cluding Eqs. (4) and (5). However the latter contain, in their
Lehmann representations, all the relevant information regard-
ing the excitation of ph and pp/hh collective modes. The ap-
proach of Ref. [24] consists in computing these quantities by
=g II
Π(ph)
(pp/hh)II
(pp/hh)
FIG. 2: (Color online) Diagrammatic equations for the polarization
(above) and the two-particle (below) propagators in the (G)RPA ap-
proach. Dashed lines are always antisymmetrized Coulomb matrix
elements and the full lines represent free (undressed) propagators.
solving the ring-GRPA and the ladder-RPA equations [10],
which are depicted for propagators in Fig. 2. In the more
general case of a self-consistent calculation, a fragmented in-
put propagator can be used and the corresponding dressed
(G)RPA [D(G)RPA] equations [10, 28] solved [see Eqs. (A2a)
and (A2b)]. Since the propagators (4) and (5) reflect two-body
correlations, they still have to be coupled to an additional sp
propagator in order to obtain the corresponding approximation
for the 2p1h and 2h1p components of R(ω). This is achieved
by solving two separate sets of Faddeev equations.
Taking the 2p1h case as an example, one can split R(2p1h)(ω)
in three different components R̄(i)(ω) (i = 1, 2, 3) that differ
from each other by the last pair of lines that interact in their
diagrammatic expansion,
R̄(2p1h)
αβγ,µνλ
(ω) =
αβγ,µνλ(ω) −G
βαγ,µνλ(ω)
i=1,2,3
R̄(i)
αβγ,µνλ
(ω) ,
where G0
(ω) is the 2p1h propagator for three freely propa-
gating lines. These components are solutions of the following
set Faddeev equations [29]
R̄(i)
αβγ,µνλ
(ω) = G0
αβγ,µ′ν′λ′(ω) Γ
µ′ν′λ′,µ′′ν′′λ′′
µ′′ν′′λ′′ ,µνλ
(ω) + R̄(k)
µ′′ν′′λ′′ ,µνλ
(ω) (7)
µ′′ν′′λ′′ ,µνλ(ω) −G
ν′′µ′′λ′′ ,µνλ(ω)
, i = 1, 2, 3
where (i, j, k) are cyclic permutations of (1, 2, 3). The inter-
action vertices Γ(i)(ω) contain the couplings of a ph or pp/hh
collective excitation and a freely propagating line. These are
given in the Appendix in terms of the polarization (4) and
two-particle (5) propagators. Equations. (7) include RPA-like
phonons and fully describe the resulting energy dependence
of R(ω). However, they still neglect energy-independent
contributions–even at low order in the interaction–that also
correspond to relevant ground-state correlations. The latter
can be systematically inserted according to
R(2p1h)
αβγ,µνλ
(ω) = Uµνλ,µ′ν′λ′ R̄
(2p1h)
µ′ν′λ′,µ′′ν′′λ′′
(ω) U†
µ′′ν′′λ′′ ,µνλ
, (8)
where R(ω) is the propagator we employ in Eq. (3), R̄(ω) is
the one obtained by solving Eqs. (7), U ≡ I + ∆U, and I is the
(pp/hh)
Π(ph)
g II (pp/hh)
Π(ph)
FIG. 3: (Color online) Example of one of the diagrams that are
summed to all orders by means of the Faddeev Eqs. (7) (left). The
corresponding contribution to the self-energy, obtained upon inser-
tion into Eq. (3), is also shown (right).
identity matrix. Following the algebraic diagrammatic con-
struction method [22, 23], the energy independent term ∆U
was determined by expanding Eq. (8) in terms of the inter-
action and imposing that it fulfills perturbation theory up to
first order (corresponding to third order in the self-energy).
The resulting ∆U, employed in this work, is the same as in
Ref. [23] and is reported in App. A 3 for completeness. It
has been shown that the additional diagrams introduced by
this correction are required to obtain accurate QP properties.
Equations. (7) and (8) are valid only in the case in which a
mean-field propagator is used to expand R(ω). This is the case
of the present work, which employs Hartree-Fock sp propaga-
tors as input. The derivation of these equations for the gen-
eral case of a fragmented propagator is given in the appendix.
More details about the actual implementation of the Faddeev
formalism to 2p1h/2h1p propagation have been presented in
Ref. [24]. The calculation of the 2h1p component of R(ω)
follows completely analogous steps.
It is important to note that the present formalism includes
the effects of ph and pp/hh motion to be included simulta-
neously, while allowing interferences between these modes.
These excitations are evaluated here at the RPA level and are
then coupled to each other by solving Eqs. (7). This generates
diagrams as the one displayed in Fig. 3, with the caveat that
two phonons are not allowed to propagate at the same time.
Equations. (7) also assure that Pauli correlations are properly
taken into account at the 2p1h/2h1p level. In addition, one can
in principle employ dressed sp propagators in these equations
to generate a self-consistent solution. If we neglect the lad-
der propagator gII (ω) (5) in this expansion, we are left with
the ring series alone and the analogous physics ingredients as
for the generalized GW approach. However, this differs from
GGW due to the fact that no double counting of the second-
order self-energy occurs, since the Pauli exchanges between
the polarization propagator and the third line are properly ac-
counted for (see Fig. 3). Alternatively, one can suppress the
polarization propagator to investigate the effects of pp/hh lad-
ders alone.
It is instructive to replace in the above equations all RPA
phonons with TDA ones; this amounts to allowing only
forward-propagating diagrams in Fig. 2, and is equivalent to
separate diagonalisations in the spaces of ph, pp and hh con-
l 0 1 2 3 4 5 6
rw 2.0 4.0 0.0 0.0 0.0 0.0 0.0
no 12 21 10 10 5 5 5
TABLE I: Parameters that define the sp basis: radius of the confining
wall rw (in atomic units) and number of orbits no used for different
partial waves l. The value of cw is always set to 5 a.u..
figurations, relative to the HF ground state. It can be shown
that using these TDA phonons to sum all diagrams of the type
in Fig. 3 reduces to one single diagonalization in the 2p1h or
2h1p spaces. Therefore, Eqs. (7) and (8) with TDA phonons
lead directly to the “extended” 2p1h TDA of Ref. [23], which
was later shown to be equivalent to ADC(3) in the general
ADC framework [22]. The Faddeev expansion formalism of
Ref. [24] creates the possibility to go beyond ADC(3) by in-
cluding RPA phonons. This is more satisfactory in the limit
of large systems. At the same time, the computational cost
remains modest since only diagonalizations in the 2p1h/2h1p
spaces are required.
Note that complete self-consistency requires the use of
fragmented (or dressed) propagators in the evaluation of all in-
gredients leading to the self-energy. This is outside the scope
of the present paper, but we included partial selfconsistency
by taking into account the modifications to the HF diagram
by employing the correlated one-body density matrix and it-
erating to convergence. This is relatively simple to achieve,
since the 2p1h/2h1p propagator is only evaluated once with
the input HF propagators. Below we will give results with and
without this partial selfconsistency at the HF level.
III. RESULTS
Calculations have been performed using two different
model spaces: (1) a standard quantumchemical Gaussian ba-
sis set, aug-cc-pVTZ for Neon [30], with Cartesian repre-
sentation of the d and f functions; (2) a numerical basis set
based on HF and subsequent discretization of the continuum,
to be detailed below. The aug-cc-pVTZ basis set was used
primarily to check our formalism with the ADC(3) result in
literature (i.e. [31], where this basis was employed). The
HF+continuum basis allows to approach, at least for the ion-
ization energies, the results for the full sp space (basis set
limit).
The HF+continuum is the same discrete model space em-
ployed previously in Refs. [11, 14]. It consists of: (1) Solv-
ing on a radial grid the HF problem for the neutral atom; (2)
Adding to this fixed nonlocal HF potential a parabolic poten-
tial wall of the type U(r) = θ(r − rw)cw(r − rw)
2, placed at a
distance rw of the nucleus. The latter eigenvalue problem has
a basis of discrete eigenstates. This basis is truncated by spec-
ifying some largest angular momentum lmax and the number
of virtual states for each value of l ≤ lmax. (3) Solve the HF
problem again, without the potential wall, in this truncated
discrete space. The resulting basis set is used for the subse-
F-TDA F-RPA F-TDAc F-RPAc Expt.
2p -0.799 -0.791 -0.803 -0.797 (0.94) -0.793 (0.92)
2s -1.796 -1.787 -1.802 -1.793 (0.90) -1.782 (0.85)
1s -32.126 -32.087 -32.140 -32.102 (0.86) -31.70
Etot -128.778 -128.772 -128.836 -128.840 -128.928
TABLE II: Results with the aug-cc-pVTZ basis. The first three rows
list the energies of the main sp fragments below the Fermi level, as
predicted by different self-energies. F-TDA/F-RPA refers to the Fad-
deev summation with TDA/RPA phonons, respectively. In all cases
the self-energy was corrected at third order through Eq. (8). The
suffix “c” refers to partial selfconsistency, when the static (HF-type)
self-energy is consistent with the correlated density matrix. Without
“c” the pure HF self-energy was taken. In the F-RPAc column the
strength of the fragment is indicated between brackets. The last row
is the total electronic binding energy. The experimental values are
taken from Refs. [32, 33]. All energies are in atomic units.
quent Green’s function calculations.
When a sufficiently large number of states is retained after
truncation, the final results should approach the basis set limit.
In particular the results should not depend on the choice of the
auxiliary confining potential. This was verified in Ref. [11] for
the second-order, and in Ref. [14] for the G0W0 self-energy; in
these cases the self-energy is sufficiently simple that extensive
convergence checks can be made for various choices of the
auxiliary potential. The parameters of the confining wall and
the number of sp states kept in the basis set was optimized in
Ref. [11], by requiring that the ionization energy is converged
to about 1 mH for the second-order self-energy. In Ref. [14]
the same choice of basis set was also seen to bring the ioniza-
tion energy for the G0W0 self-energy near convergence. For
completeness, the details of this basis are reported in Table I.
While the self-energy in the present paper is too complicated
to allow similar convergence checks, it seems safe to assume
that basis set effects will affect the calculated ionization ener-
gies by at most 5 mH.
In Table II we compare, for the aug-cc-pVTZ basis, the ion-
ization energies of the main single-hole configurations when
TDA or RPA phonons are employed in the Faddeev construc-
tion (this is labeled F-TDA and F-RPA, respectively, in the
table). Note that use of TDA phonons corresponds to the
usual ADC(3) self-energy. We find that the replacement of
TDA with RPA phonons provides more screening, leading to
slightly less bound poles which are shifted towards the exper-
imental values. This shift increases with binding energy. As
discussed at the end of Sec. I, one can include consistency of
the static part of the self-energy. About eight iterations are
needed for convergence. This is a nonnegligible correction,
providing about 5 mH more binding (i.e. larger ionization en-
ergies) for the valence/subvalence 2p and 2s, 15 mH for the
deeply bound 1s, and 60 mH to the total binding energy. Our
converged result for the Faddeev-TDA self-energy (labeled F-
TDAc in Table II) is in good agreement with the ADC(3) value
for the 2p ionization energy (-0.804 H) quoted in [31], as it
should be.
The analogous results obtained with the larger
F-TDA F-RPA F-TDAc F-RPAc Expt.
2p -0.807 -0.799 -0.808 -0.801 (0.94) -0.793 (0.92)
2s -1.802 -1.792 -1.804 -1.795 (0.91) -1.782 (0.85)
1s -32.136 -32.097 -32.142 -32.104 (0.81) -31.70
Etot -128.863 -128.857 -128.883 -128.888 -128.928
TABLE III: Results with the HF+continuum basis set from Table I.
See also the caption of Table II.
-40 -35 -30
 [a.u.]
F-RPA(ring)
F-RPA(ladder)
F-RPA
FIG. 4: (Color online) Spectral function for the s states in Ne
obtained with various self-energy approximations. From the top
down: the second-order (Σ(2)) self-energy, the F-RPA(ring), the F-
RPA(ladder), and the full F-RPA self-energy. The strength is given
relative to the Hartree-Fock occupation of each shell. Only fragments
with strength larger than Z > 0.005 are shown.
HF+continuum basis are given in Table III, which al-
lows to assess overall stability and basis set effects. We find
exactly the same trends as for aug-cc-pVTZ. In particular
the reduction of ionization energies from the replacement of
TDA with RPA phonons is almost independent of the basis
set used, while the effect of including partial consistency
is roughly halved. Overall, the ionization states are always
more bound with the larger basis set; while the basis set
limit could be still more bound than the present results with
the HF+continuum basis set, it is likely (based on the G0W0
extrapolation in Ref. [14]) that the difference does not exceed
5 mH.
As discussed in Sec. I, the the F-RPA self-energy contains
RPA excitations of both ph type (ring diagrams) and pp/hh
type (ladder diagrams). It is instructive to analyze their sepa-
rate contributions to the final ionization energies, in order to
1s 2s 2p
HF -32.77 (1.00) -1.931 (1.00) -0.850 (1.00)
Σ(2) -31.84 (0.74) -1.736 (0.88) -0.747 (0.91)
G0W0 -31.14 (0.85) -1.774 (0.91) -0.801 (0.94)
F-RPA (ring) -31.82 (0.73) -1.636 (0.56) -0.730 (0.80)
F-RPA (ladder) -32.04 (0.87) -1.802 (0.95) -0.781 (0.96)
F-RPA -32.10 (0.81) -1.792 (0.91) -0.799 (0.94)
Exp. -31.70 -1.782 (0.85) -0.793 (0.92)
TABLE IV: Energy (in a.u.) and strength (bracketed numbers) of the
main fragments in the spectral function of Neon, generated by differ-
ent self-energies. Results for the HF+continuum basis. Consecutive
rows refer to: (1) HF; (2) second-order self-energy; (3) G0W0 results
from Ref. [14]; (4) F-RPA self-energy with only ph rings retained;
(5) F-RPA self-energy with only pp/hh ladders retained; (6) Com-
plete F-RPA self-energy. In all F-RPA results the self-energy was
corrected at third order through Eq. (8). The static self-energy was
pure HF (no partial self-consistency). The experimental values are
taken from Refs. [32, 33].
understand how the F-RPA self-energy is related to the stan-
dard (G)GW self-energy. Table IV compares the results for
the ionization energies, obtained with the second-order self-
energy, to different approximations for including the ring sum-
mations. As one can see, the second-order self-energy gener-
ates an l=1 sp energy of -0.747 mH, which is 46 mH above
the empirical 2p ionization energy. The G0W0 self-energy,
which includes the ring summation with only direct Coulomb
matrix elements, improves this result and brings it close to ex-
periment. The 2s behaves in a similar way. Unfortunately,
including the exchange terms of the interelectron repulsion in
the GG0W0 method turns out to have the opposite effect (the
2p ionization energy becomes -0.712 H [14] [41]) and the
agreement with experiment is lost. Obviously, GG0W0 is too
simplistic to account for exchange in the ph channel.
With the F-RPA(ring) self-energy one can go one step fur-
ther and employ the Faddeev expansion to also force proper
Pauli exchange correlations in the 2p1h/2h1p spaces. As
shown in Table IV, this enhances the screening due to the
exchange interaction terms, leading to even less binding for
the 2s and 2p. The corrections relative to the second-order
self-energy can be large (100 mH for the 2s state) and in the
direction away from the experimental value. We also note that
the larger shift, in the 2s orbit, is accompanied by an increase
of the fragmentation (see Fig. 4 and Tab. IV). Similar observa-
tions were also made in Ref. [14] for other atoms: in general
ring summations in the direct channel alone bring the quasi-
hole peaks close to the experiment. This agreement is then
spoiled as soon as one includes proper exchange terms in the
self-energy. On the other hand, exchange in the ph channel
is required to reproduce the correct Rydberg sequence in the
excitation spectrum of neutral atoms. So further corrections
must arise from other diagrams, and obviously the summation
of ladder diagrams can play a relevant role, since these con-
tribute to the expansion of the self-energy at the same level as
that of the ring diagrams.
The result when only including ladder-type RPA phonons
in the F-RPA self-energy is also shown in Table IV. One can
see that pp/hh ladders do actually work in the opposite way
as the ph channel ring diagrams, and have the same order of
magnitude with, e.g., a shift of 66 mH for the 2s relative to
the second-order result. When combined with the ring dia-
grams in the full F-RPA self-energy, the agreement with ex-
periment is restored again. Note that the final result cannot be
obtained by adding the contributions of rings and ladders, but
depends nontrivially on the interplay between these classes of
diagrams thereby pointing to significant interference effects.
With the F-RPA(ring) self-energy, where only the contribu-
tions of the ph channel are included, the main peaks listed in
Table IV are not only considerably shifted but also strongly
depleted, e.g. a strength of only 0.56 for the main 2s peak.
The complete spectral function for the l = 0 strength in Fig. 4
shows that the depletion of the main fragment is accompanied
by strong fragmentation over several states. While correlation
effects are overestimated in F-RPA(ring), they are suppressed
in F-RPA(ladder), where only the pp/hh ladders are included
in the self-energy. In this case one finds a spectral distribu-
tion closer to the HF one, with a main 2s fragment of strength
0.95 and less fragmentation than the the second-order self-
energy. The spectral distribution generated by the complete
F-RPA self-energy is again a combination of the above ef-
fects. The strength of the deeply bound 1s orbital behaves in
an analogous way. The strength of the main peak is reduced
but several satellite levels appear due to the mixing with 2h1p
configurations. In all the calculations reported in Fig. 4 we
found a summed l = 0 strength exceeding 0.98 in the interval
[-40 H, -30 H] which can be associated with the 1s orbital,
and this remains true even in the presence of strong correla-
tions using the F-RPA(ring) self-energy. Of course, the mix-
ing with 3h2p configurations, not included in this work, may
further contribute to the fragmentation pattern in this energy
region.
IV. CONLUSIONS AND DISCUSSION
In conclusion, the electronic self-energy for the Ne atom
was computed by the F-RPA method which includes –
simultaneously– the effects of both ring and ladder diagrams.
This was accomplished by employing an expansion of the
self-energy based on a set of Faddeev equations. This tech-
nique was originally proposed for nuclear structure applica-
tions [24] and is described in the appendix. At the level of
the self-energy one sums all diagrams where the three propa-
gator lines of the intermediate 2p1h or 2h1p propagation are
connected by repeated exchange of RPA excitations in both
the ph and the pp/hh channel. This differs from the ADC(3)
formalism in the fact that the exchanged excitations are of
the RPA type, rather than the TDA type, and therefore take
ground-state correlations effects into account. The coupling
to the external points of the self-energy uses the same modi-
fied vertex as in ADC(3), which must be introduced to include
consistently all third-order perturbative contributions.
The resulting main ionization energies in the Neon atom
are at least of the same quality, and even somewhat improved,
compared to the ADC(3) result. Note that, numerically, F-
RPA can be implemented as a diagonalization in 2p1h/2h1p
space implying about the same cost as ADC(3). The present
study also shows that in localized electronic systems subtle
cancellations occur between the ring and ladder series. In par-
ticular, only a combination of the ring and ladder series leads
to sensible results, as the separate ring series tends to correct
the second-order result in the wrong direction.
Since the limit to extended systems requires an RPA treat-
ment of excitations, the F-RPA method holds promise to
bridge the gap between an accurate description of quasiparti-
cles in both finite and extended systems. In particular, the GW
treatment of the electron gas has been shown to yield excellent
binding energies, but poor quasiparticle properties [34, 35].
Further progress beyond GW theory requires a consistent in-
corporation of exchange in the ph channel. The F-RPA tech-
nique may be highly relevant in this respect. A common
framework for calculating accurate QP properties in both fi-
nite and extended systems, is also important for constrain-
ing functionals in quasiparticle density functional theory (QP-
DFT) [7].
Finally, complete self-consistency requires sizable compu-
tational efforts for bases as large as the HF+continuum basis
used here. It would nevertheless represent an important exten-
sion of the present work, since it is related to the fulfillment
of conservation laws [36, 37]. These issues will be addressed
in future work.
Acknowledgments
This work was supported by the U.S. National Science
Foundation under grant PHY-0652900.
APPENDIX A: FADDEEV EXPANSION OF THE 2P1H/2H1P
PROPAGATOR
Although only the one-energy (or two-time) part of the
2p1h/2h1p propagator enters the definition of the self energy,
Eq. (3), a full resummation of all its diagrammatic contribu-
tions would require to treat explicitly the dependence on three
separate frequencies, corresponding to the three final lines in
the expansion of R(ω). For example, inserting the RPA ring
(ladder) series in R(ω) implies the propagation of a ph (pp/hh)
pair of lines both forward and backward in time, while the
third line remains unaffected. A way out of this situation is
to solve the Bethe-Salpeter-like equations for the polarization
and ladder propagators separately and then to couple them to
the additional line. If it is assumed that different phonons do
not overlap in time, the three lines in between phonon struc-
tures will propagate only in one time direction [see figures (3)
and (5)]. In this situation the integration over several frequen-
cies can be circumvented following the prescription detailed
in the next subsection. This approach will be discussed in
the following for the general case of a fully fragmented prop-
gator, in order to derive a set of Faddeev equations capable
FIG. 5: (Color online) Diagrammatic representation of Eq. (A3).
Double lines represent fully dressed sp Green’s funcions which, how-
ever, are restricted to propagate only in one time direction [i.e., only
one of the two terms on the r.h.s. of Eq. (1) is retained]. The Faddeev
Eqs. (A9) and (7) allow for both forward and backward propagation
of the phonons Γ(π)(ω) and Γ(II)(ω) as long as these do not overlap
in time. For the propagators, time ordereing is asumed with forwad
propagation in the upward direction.
of dressing the sp propagator self-consistently. Since the for-
ward (2p1h) and the backward (2h1p) parts of R(ω) decouple
in two analogous sets of equations, it is sufficient to focus on
the first case alone.
1. Multiple frequencies integrals
We start by considering the effective interactions in the ph
and pp/hh channels that correspond to Eqs. (4) and (5) stripped
of the external legs. In the present work, these are the follow-
ing two-time objects:
αβ,γδ
(ω) = Vαδ,βγ + Vαν,βµ Π
µν,ρσ(ω) Vρδ,σγ (A1a)
= Vαδ,βγ +
ω − επn + iη
ω + επn′ − iη
αβ,γδ
(ω) = Vαβ,γδ + Vαβ,µν g
µν,ρσ(ω) Vρσ,γδ (A1b)
= Vαδ,βγ +
ω − εΓ+n + iη
ω − εΓ−
where the residues and poles for the ring series are Ωn
〈ΨNn |c
µcν|Ψ
0 〉Vµβ,να and ε
n = E
0 . For the ladders, ∆
〈ΨN+2n |c
0 〉Vµν,αβ and ∆
= Vαβ,µν〈Ψ
k |cµcν|Ψ
0 〉, with
poles εΓ+n = E
0 and ε
k = E
. Equations. (A1)
solve the ring and ladder RPA equations, respectively
αβ,γδ
(ω) = Vαδ,βγ (A2a)
αβ,µν
gµρ(ω + ω1)gσν(ω1) Vρδ,σγ ,
αβ,γδ
(ω) = Vαβ,γδ (A2b)
αβ,µν
gµρ(ω − ω1)gνσ(ω1) Vρσ,γδ .
To display how the phonons (A1a) and (A1b) enter the ex-
pansion of R(ω), we perform explicitly the frequency integrals
for the diagram of Fig. 5. Since it is assumed that the separate
propagators lines evolve only in one time direction, only the
forwardgoing (g>(ω)) or the backwardgoing (g<(ω)) part of
Eq. (1) must be included for particles and holes, respectively.
After some algebra, one obtains
∆Rαβγ,µνλ(ω) =
g>αα1 (ω −Ω) g
(ω1) g
(ω1 −Ω) Γ
β1γ1,σ1λ1
(Ω) g>σ1σ2 (s + Ω − ω)
α1σ2,µ1ν1
(s) g>µ1µ(s − ω2) g
(ω2) g
(s − ω)
ω − (ε+n1 + ε
) + iη
Vβ1λ1,γ1σ1 +
ω − (ε+n1 + ε
nπ) + iη
[ω − επ
− ε+n1 − ε
+ ε−k3
− ε+n4 + ε
− ε+n2 + ε
][−επ
− ε+n4 + ε
ω − (ε+n1 + ε
) + iη
Vα1σ2,µ1ν1 +
+,nII
+,nII
ω − (εΓ+nII − ε
) + iη
[ω + εΓ−kII − ε
− ε+n4 − ε
− ε+n6 + ε
] ∆−,kIIα1σ2
−,kII
− ε+n1 − ε
][εΓ−
− ε+n5 − ε
ω − (ε+n5 + ε
− ε−k7
) + iη
ω − (εΓ−
− ε+n4 − ε
) − iη
−,kII
−,kII
− ε+n2 + ε
] [−επ
− ε+n4 + ε
] [εΓ−
− ε+n1 − ε
] [εΓ−
− ε+n5 − ε
The last term in this expression contains an energy denom-
inator that involves the simultaneous propagation of two
phonons. Thus, it will be discarded in accordance with our
assumptions. It must be stressed that similar terms, with over-
lapping phonons, imply the explicit contribution of at least
3p2h/3h2p. A proper treatment of these would require a non
trivial externsion of the present formalism, which is beyond
the scope of this paper.
The remaining part in Eq. (A3) is the relevant contribution
for our purposes. This has the correct energy dependence of
a product of denominators that correspond to the intermediate
steps of propagation. All of these involve configurations that
have at most 2p1h character. Although, ground state corre-
lations are implicitely included by having already resummed
the RPA series. Still, this term does not factorize in a prod-
uct of separate Green’s functions due to the summations over
the fragmentation indices ni and ki [labeling the eigenstates of
the (N±1)-electron systems]. This is overcome if one defines
the matrices G0>(ω), Γ(1,2)(ω) and Γ(3)(ω), with elements (no
implicit summation used)
αnαβnβγkγ; µnµνnνλkλ
(ω) = δnα,nµ δnβ,nν δkγ,kλ
ω − (ε+nα + ε
− ε−kγ
) + iη
, (A4a)
αnαβnβγkγ; µnµνnνλkλ
(ω) = Γ(2)>
βnβαnαγkγ; νnνµnµλkλ
(ω) =
δα,µ δnα ,nµ
Vβλ,γν +
ω − (ε+nα + ε
nπ) + iη
[ω − επ
− ε+nα − ε
+ ε−kγ
− ε+nν + ε
− ε+nβ + ε
][−επ
− ε+nν + ε
, (A4b)
αnαβnβγkγ; µnµνnνλkλ
(ω) =
δγ,λ δkγ ,kλ
Vαβ,µν +
+,nII
+,nII
ω − (εΓ+nII − ε
) + iη
[ω + εΓ−kII − ε
− ε+nβ − ε
− ε+nν + ε
] ∆−,kII
−,kII
− ε+nα − ε
][εΓ−
− ε+nµ − ε
. (A4c)
In these definitions, the row and column indices are ordered to
represent at first two quasiparticle lines and then a quasihole.
The index ‘i’ in Γ(i)> refer to the line that propagates indepen-
dently along with the phonon. Using Eqs.(A4), the first term
on the r.h.s. of Eq. (A3) can be written as
∆R(2p1h)
αβγ,µνλ
(ω) = (A5)
nα nβ kγ
nµ nν kλ
G0>(ω)Γ(1)>(ω)G0>(ω)Γ(3)>(ω)G0>(ω)
αnαβnβγkγ; µnµνnνλkλ
Eq. (A5) generalizes to diagrams involving any number of
phonon insertions, as long as the terms involving two or more
simultaneous phonons are dropped. Based on this relation, we
use the following prescription to avoid performing integrals
over frequencies. One extends all the Green’s functions to ob-
jects depending not only on the sp basis’ indices (α, β, γ) but
also on the indices labeling quasi-particles and holes (ni and
ki). Whether a given argument represents a particle or an hole
depends on the type of line being propagated. At this point one
can perform calculations working with only two-time quanti-
ties. The standard propagator is recovered at the end by sum-
ming the “extended” one over the quasi-particle/hole indices.
2. Faddeev expansion
The 2p1h/2h1p propagator that includes the full resumma-
tion of both the ladder and ring diagrams at the (G)RPA level
is the solution of the following Bethe-Salpeter-like equation,
Rαβγ,µνλ(ω1, ω2, ω3) = (A6)
gαµ(ω1)gβν(ω2) − gβµ(ω2)gαν(ω1)
gλγ(−ω3) +
gββ1(ω2)gγ1γ(−ω3)Vβ1σ,γ1ρ
Rαρσ,µνλ(ω1, s, ω2 + ω3 − s)
+ gαα1 (ω1)gγ1γ(−ω3)Vα1σ,γ1ρ
Rρβσ,µνλ(s, ω2, ω1 + ω3 − s)
gαα1 (ω1)gββ1(ω2)Vα1β1,ρσ
Rρσγ,µνλ(s, ω1 + ω2 − s, ω3)
If this equation is solved, a double integration of R(ω1, ω2, ω3)
would yield the two-time propagator R(ω) contributing to
Eq. (3). However, the numerical solution of Eq. (A6) appears
beyond reach of the present day computers and one needs
to avoid dealing directly with multiple frequencies integrals.
The strategy used is to first solve the RPA equations (A2a)
and (A2b) separately. Once this is done it is necessary to re-
arrange the series (A6) in such a way that only the resummed
phonons appear. Following the formalism introduced by Fad-
deev [29, 38], we identify the components R(i)(ω) with the
three terms between curly brakets in Eq. (A6). By employ-
ing Eqs. (A2a) and (A2b) one is lead to the following set of
equations [42],
αβγ,µνλ
(ω1, ω2, ω3) = gαα1 (ω1)gββ1(ω2)gγ1γ(−ω3) (A7)
ds1 ds2 ds3
α1β1γ1,µ1ν1λ1
(ω1, ω2, ω3; s1, s2, s3)
gµ1µ(s1)gν1ν(s2) − gν1µ(s2)gµ1ν(s1)
gλλ1(−s3)
µ1ν1λ1,µνλ
(s1, s2, s3) + R
µ1ν1λ1,µνλ
(s1, s2, s3)
, i = 1, 2, 3 ,
where (i,j,k) are cyclic permutations of (1,2,3) and the inter-
action vertices Γ(i)(ω1, ω2, ω3) are given by
αβγ,µνλ
(ω1, ω2, ω3;ω4, ω5, ω6) = (A8a)
βαγ,νµλ
(ω2, ω1, ω3;ω5, ω4, ω6) =
= δ(ω1 − ω4)δ(ω2 + ω3 − ω5 − ω6)g
αµ(ω1)Γ
βγ,νλ
(ω2 + ω3) ,
αβγ,µνλ
(ω1, ω2, ω3;ω4, ω5, ω6) = (A8b)
δ(ω3 − ω6)δ(ω1 + ω2 − ω4 − ω5)g
λγ (−ω3)Γ
αβ,µν
(ω1 + ω2) .
Finally, we apply the prescription of Sec. A 1 and substitute
R(ω1, ω2, ω3) with its extended but two-time version R(ω).
This leads to the following set of Faddeev equations which
propagate 2p1h forward in time,
R̄(i)
αnαβnβγkγ; µnµνnνλkλ
(ω) =
αnαβnβγkγ; α′n′αβ′n
γ′k′γ
(ω) Γ(i)
α′n′αβ′n
γ′k′γ; µ′n
′n′νλ′k
αnαβnβγkγ; µnµνnνλkλ
(ω) −G0
αnαβnβγkγ; µnµνnνλkλ
× R̄( j)
µ′n′µν′n
; µnµνnνλkλ
(ω) + R̄(k)
µ′n′µν′n
; µnµνnνλkλ
i = 1, 2, 3 . (A9)
Since the full energy dependence is retained in Eq. (A7),
the self-energy corresponding to its solution, R(ω1, ω2, ω3), is
complete up to third order [see Eq. (3)]. This is no longer the
case after the reduction to a two-time propagator. In particu-
lar, the approximation that only forward 2p1h propagation is
allowed between different phonons implies that all diagrams
with different time propagation of their external lines are ne-
glected in Eqs. (A9). However, these terms are not energy
dependent and can be can be reinserted in a systematic way a
posteriori as in Eq. (8). In the general case,
Rαβγ,µνλ(ω) = (A10)
U (2p1h)
αβγ; α′n′αβ′n
γ′k′γ
R̄(2p1h)
α′n′αβ′n
γ′k′γ; µ′n
′n′νλ′k
(ω) U (2p1h) †
µ′n′µν′n
; µνλ
U (2p1h)
αβγ; µnµνnνλkλ
= δαµδβνδγλ + ∆U
(2p1h)
αβγ; µnµνnνλkλ
, (A11)
where the correction ∆U can be determined by comparison
with perturbation theory.
The vertices (A4), that appear in Eqs. (A9), and U (2p1h) are
expressed in terms of the fully fragmented propagator. There-
fore, this approach allows to obtain self-consistent solutions
of the sp Green’s function [25]. Whenever, like in this work,
only a mean-field propagator is employed as input there ex-
ist a one-to-one correspondence between the fragmentation
indices and the sp basis. This is expressed by the relations
Xnα = δn,α(1 − δα∈F) and Y
α = δk,αδα∈F , where F represents
the set of occupied orbits. In this case, it is possible to drop
one set of indices so that Eqs. (A9) and (A10) simplify into
the form (7) and (8).
3. Faddeev vertices
In practical applications, it is worth to note that the poles
of the free propagator G0(ω), Eq. (A4a), do not contribute to
the kernel of Eqs. (A9). This can be proven by employing the
closure relations for the RPA problem, in the form obtained
by extracting the free poles in Eqs. (A2). As an example, for
the forward poles of the ladder propagator these are
ω→ε+n1
+ε+n2
[(ω − ε+n1 − ε
) × (Eq. A2b)] =⇒ (A12)
µν,γδ
(ω = ε+n1 + ε
) = 0 , ∀n1, n2 ,
and similarly for other cases. Making use of these relations
one can derive the following working expression of the ker-
nels of the 2p1h Faddeev equations (no implicit summations
used)
G0>(ω)Γ(1)>(ω)
αnαβnβγkγ; µnµνnνλkλ
G0>(ω)Γ(2)>(ω)
βnβαnαγkγ; νnνµnµλkλ
= (A13a)
= δnα,nµ
[επnπ − ε
][ω − (ε+nα + ε
nπ) + iη]
− ε+nβ + ε
][−επ
− ε+nν + ε
G0>(ω)Γ(3)>(ω)
αnαβnβγkγ; µnµνnνλkλ
= (A13b)
= δkγ,kλ
+,nII
+,nII
[εΓ+nII − ε
− ε+nβ][ω − (ε
) + iη]
−,kII
−,kII
− ε+nα − ε
][εΓ−
− ε+nµ − ε
After substituting Eq. (A10) into (3), one needs the working expression for the matrix product V U(2p1h) (where V is the inter-
electron interaction). The minimum correction that guaranties to reproduce all third order self-energy diagrams is
V U(2p1h)
α; µnµνnνλkλ
= Vαλ,µν +
Vαλ,γ1δ1 Y
Vγ2δ2,µν
2 [ε−
− ε+nµ − ε
(A14)
Vαδ1,µγ1 Y
Vγ2λ,δ2ν
− ε+nδ − ε
Vαδ1,νγ1 Y
Vγ2λ,δ2µ
− ε+nδ − ε
The case of 2h1p is handled in a completely analogous way along the steps of Secs. (A 1) and (A 2). After extending
R(ω1, ω2, ω3) to depend on the fragmentation indices (k1,k2,n), the 2h1p equivalent of Eq. (A9) is obtained with the follow-
ing definitions of the kernels,
G0>(ω)Γ(1)>(ω)
αkαβkβγnγ; µkµνkνλnλ
G0>(ω)Γ(2)>(ω)
βkβαkαγnγ; νkνµkµλnλ
= (A15a)
= δkα,kµ
+ ε+nγ ][ω − (ε
) − iη]
[επnπ − ε
+ ε+nγ ][ε
nπ − ε
+ ε+nλ]
G0>(ω)Γ(3)>(ω)
αkαβkβγnγ; µkµνkνλnλ
= (A15b)
= δnγ,nλ
−,nII
−,nII
][ω − (εΓ−
− ε+nγ ) − iη]
+,kII
+,kII
[εΓ+nII − ε
− ε−kβ
][εΓ+nII − ε
− ε−kν
and correction to the external legs,
V U(2h1p)
α; µkµνkνλnλ
= Vαλ,µν +
Vαλ,γ1δ1 X
Vγ2δ2,µν
2 [ε−kµ + ε
− ε+nγ − ε
(A16)
Vαδ1,µγ1 X
Vγ2λ,δ2ν
[ε−kδ + ε
− ε+nγ − ε
Vαδ1,νγ1 X
Vγ2λ,δ2µ
[ε−kδ + ε
− ε+nγ − ε
It should be pointed out that while the prescription of
Sec. A 1 allows sp lines to propagate only in one time di-
rection, it allows for backward propagation of the phonons.
These contributions translate directly into the energy indepen-
dent terms of Eqs. (A13) and (A15) and are a direct conse-
quence of the inversion pattern typical of RPA theory. These
terms have normally a weaker impact than the direct ones on
the solutions of Eqs. (A9). However, it is show in Ref. [24]
that they are crucial to guarantee the exact separation of the
spurious solutions—always introduced by the Faddeev for-
malism [39, 40]—if RPA phonons are used. For the same rea-
sons, the last terms in curly brackets of Eqs. (A13) and (A15)
should be dropped whenever Tamm-Dancoff (TDA) phonons
are propagated.
The approach followed in this work for solving Eqs. (A9)
is to transform them into a matrix representations [24].
Once this is done, one is left with an eigenvalue prob-
lem that depends only on the 2p1h (2h1p) configurations
(n, n′, k) [(k, k′, n)]. The spurious states are known ex-
actly [24] and can be projected out analytically to reduce the
computational load. In any case, they would give vanishing
contributions to Eq. (3).
[1] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
[2] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
[3] A. Görling and M. Levy, Phys. Rev. A 50, 196 (1994).
[4] A. Görling, J. Chem. Phys. 123, 062203 (2005).
[5] R. J. Bartlett et al., J. Chem. Phys. 122, 034104 (2005); J.
Chem. Phys. 123, 062205 (2005).
[6] P. Mori-Sanchez, Q. Wu, and W. T. Yang, J. Chem. Phys. 123,
062204 (2005).
[7] D. Van Neck, S. Verdonck, G. Bonny, P. W. Ayers, and M.
Waroquier, Phys. Rev. A 74, 042501(2006).
[8] A. B. Migdal, Theory of Finite Fermi Systems and Applications
to Atomic Nuclei (John Wiley and Sons, New York, 1967).
[9] A. L. Fetter and J. D. Walecka, Quantum Theory of Many-
Particle Physics (McGraw-Hill, New York, 1971).
[10] W. H. Dickhoff and D. Van Neck, Many-Body Theory Exposed!
(World Scientific, Singapore, 2005).
[11] D. Van Neck, K. Peirs, and M. Waroquier, J. Chem. Phys. 115,
15 (2001).
[12] K. Peirs, D. van Neck, and M. Waroquier, J. Chem. Phys. 117,
4095 (2002).
[13] N. E. Dahlen and R. van Leeuwen, J. Chem. Phys. 122, 164102
(2005).
[14] S. Verdonck, D. Van Neck, P. W. Ayers, M. Waroquier, Phys.
Rev. A 74, 062503 (2006).
[15] E.L. Shirley and R.M. Martin, Phys. Rev. B 47, 15404 (1993).
[16] N.E. Dahlen, R. van Leeuwen, and U. von Barth, Phys. Rev. A
73, 012511 (2006).
[17] N.E. Dahlen and U. von Barth, Phys. Rev. B 69, 195102 (2004).
[18] A. Stan, N.E. Dahlen and R. van Leeuwen, Europhys. Lett. 76,
298 (2006).
[19] L. Hedin, Phys. Rev. 139, A796 (1965).
[20] N. Fukuda, F. Iwamoto, and K. Sawada, Phys. Rev. 135, A932
(1964).
[21] C. Barbieri, N. Paar, R. Roth, and P. Papakostantinou,
arXiv:nucl-th/0608011v1.
[22] J. Schirmer, L.S. Cederbaum, and O. Walter, Phys. Rev. A 28,
1237 (1983)
[23] O. Walter and J. Schirmer, J. Phys. B:At. Mol. Phys. 14, 3805
(1981).
[24] C. Barbieri and W. H. Dickhoff, Phys. Rev. C 63, 034313
(2001).
[25] C. Barbieri and W. H. Dickhoff, Phys. Rev. C 65, 064313
http://arxiv.org/abs/nucl-th/0608011
(2002).
[26] W. H. Dickhoff and C. Barbieri, Prog. Part. Nucl. Phys. 52, 337
(2004).
[27] C. Barbieri, Phys. Lett. B 643, 268 (2006).
[28] W. J. W. Geurts, K. Allaart, and W. H. Dickhoff, Phys. Rev. C
50, 514 (1994).
[29] L. D. Faddeev, Zh. Éksp. Teor. Fiz. 39 1459 (1961) [Sov. Phys.
JETP 12, 1014 (1961)].
[30] T.H. Dunning Jr., J. Chem. Phys. 90, 1007 (1989); D.E. Woon
and T.H. Dunning Jr, J. Chem. Phys. 98, 1358 (1993)
[31] A.B. Trofimov and J. Schirmer, J. Chem. Phys. 123, 144115
(2005)
[32] NIST Atomic Spectra Database, NIST
Standard Reference Database #78,
http://physics.nist.gov/PhysRefData/ASD/in-dex.html
[33] A. Thompson et al., X-ray Data Booklet (Lawrence Berkeley
National Laboratory, Berkeley, CA, 2001), and references cited
therein.
[34] U. von Barth and B. Holm, Phys. Rev. B 54, 8411 (1996); B.
Holm and U. von Barth, Phys. Rev. B 57, 2108 (1998); B.
Holm, Phys. Rev. Lett. 83, 788 (1999).
[35] P. Garcia-Gonzales and R.W. Godby, Phys. Rev. B 63, 075112
(2001)
[36] G. Baym and L.P. Kadanoff, Phys. Rev. 124, 287 (1961).
[37] G. Baym, Phys. Rev. 127, 1391 (1962).
[38] C. J. Joachain, Quantum Collision Theory, (North-Holland,
Amsterdam, 1975).
[39] S. K. Adhikari and W. Glöckle, Phys. Rev. C 19, 616 (1979).
[40] J. W. Evans and D. K. Hoffman, J. Math. Phys. 22, 2858 (1981).
[41] Note that the G0W0 and the GG0W0 results of Ref. [14] were
obtained by retaining only the diagonal part of the electron self-
energy, Σ⋆αα(ω) in the HF+continuum basis. This approximation
was not made in the present work. The error in the ionization
energies by retaining the diagonal approximation is quite small
(about ∼2 mH [14] for the Ne atom), but larger effects are pos-
sible for the total binding energy.
[42] Note that the present definitions of the R(i) differ from the ones
of Ref. [24] which contain the additional term 1
[G0 − G0,ex].
The two different forms of the Faddeev equations that result can
be easily related into each other and are completely equivalent.
The definition used here agrees with the standard literature on
the subject [29, 38].
http://physics.nist.gov/PhysRefData/ASD/in-dex.html
