A thorough analytical and numerical characterization of the whole perturbation series of one-particle many-body Green’s function (MBGF) theory is presented in a pedagogical manner. Three distinct but equivalent algebraic (first-quantized) recursive definitions of the perturbation series of the Green’s function are derived, which can be combined with the well-known recursion for the self-energy. Six general-order algorithms of MBGF are developed, each implementing one of the three recursions, the ΔMPn method (where n is the perturbation order) [S. Hirata et al., J. Chem. Theory Comput. 11, 1595 (2015)], the automatic generation and interpretation of diagrams, or the numerical differentiation of the exact Green’s function with a perturbation-scaled Hamiltonian. They all display the identical, nondivergent perturbation series except ΔMPn, which agrees with MBGF in the diagonal and frequency-independent approximations at 1n3 but converges at the full-configuration-interaction (FCI) limit at n= (unless it diverges). Numerical data of the perturbation series are presented for Koopmans and non-Koopmans states to quantify the rate of convergence towards the FCI limit and the impact of the diagonal, frequency-independent, or ΔMPn approximation. The diagrammatic linkedness and thus size-consistency of the one-particle Green’s function and self-energy are demonstrated at any perturbation order on the basis of the algebraic recursions in an entirely time-independent (frequency-domain) framework. The trimming of external lines in a one-particle Green’s function to expose a self-energy diagram and the removal of reducible diagrams are also justified mathematically using the factorization theorem of Frantz and Mills. Equivalence of ΔMPn and MBGF in the diagonal and frequency-independent approximations at 1n3 is algebraically proven, also ascribing the differences at n = 4 to the so-called semi-reducible and linked-disconnected diagrams.

One-particle many-body Green’s function (MBGF) theory,1–30 also known as electron propagator theory, is one of the four pillars of ab initio electron-correlation theory,30 whose approximations are systematically improvable. Compared with the other three pillars, i.e., configuration-interaction (CI),30–33 coupled-cluster (CC),30,31,34,35 and many-body perturbation theories (MBPT),30,31,34,36 MBGF is mysterious as its Feynman–Dyson perturbation series converges at the exact basis-set solutions of a many-electron Schrödinger equation within the framework of a one-particle theory. Although the solution space of a one-particle theory is apparently much smaller than that of the many-electron Schrödinger equation, MBGF attains exactness by rendering the one-particle operator, specifically, its Dyson self-energy part, frequency dependent and thereby making each one-particle equation have multiple roots. Hence, MBGF bears similarity with equally mysterious Kohn–Sham (KS) density-functional theory (DFT),37,38 which is also a formally exact one-particle theory for a many-electron Schrödinger problem.39,40 In fact, Sham and Schlüter proposed mapping41 of a self-energy operator onto an exchange-correlation potential of KS DFT, of which the lowest-order (exchange-only) incarnation is the optimized effective potential of Talman and Shadwick.42–45 It can be further extended to include frequency dependence46–48 and electron-correlation effects.49–51 

What may add to the mysterious appearance of MBGF is the fact that the highest perturbation order of its methods developed so far is only three15,52–63 or four,20,54,64,65 the latter case using simplifying approximations. This is in contrast with CI,66,67 CC,68–72 or MBPT,73,74 all of which can be carried out at any arbitrary high order using the determinant-based algorithm of full configuration interaction (FCI).66 Even when limiting ourselves to efficient implementations, MBPT was extended to sixth order [MBPT(6)] in 1985 by Laidig et al.75 and fifth-order CC was reported in 2002 by Musiał et al.,76 in contrast to full MBGF(3) being the highest today. This means that the mathematical structure of MBGF may not be as fully understood as those of MBPT or CC, making a reliable and easy-to-understand algebraic method of derivation of higher-order MBGF presently unavailable. This may be one of the reasons why this theory has not been as fully embraced by molecular scientists as MBPT or CC is today. This situation is unfortunate because not only does MBGF directly compute electron binding energies, key parameters of most chemical processes, but also there is evidence that higher-order perturbation corrections are relatively more important in MBGF than in MBPT. Nonetheless, MBGF has recently enjoyed a surge of interest in molecular applications,39,40,77–84 sometimes in relation to CC theory.85,86

The persistent difficulty in fully understanding and advancing MBGF to higher orders can at least partly be traced to its expositions, which greatly differ from those of CI, CC, or MBPT. The latter are usually based on the time-independent Schrödinger equation with well-defined approximations made to their wave functions and energies. For example, MBPT can be defined algebraically by recursion equations derivable systematically by the Rayleigh–Schrödinger perturbation theory.31,34 A typical exposition of MBGF follows a different path.87–89 Instead of the Schrödinger equation, MBGF is based on the Dyson equation, the central entity of which is the self-energy. The derivation of MBGF typically starts by accepting the full equivalence of this equation with the Schrödinger equation and proceeds by postulating the diagrammatic Feynman–Dyson perturbation expansion9,17,18 of the self-energy. Whereas the Gell-Mann–Low theorem87,89,90 offers a time-dependent perturbation-theoretical justification of this diagrammatic derivation, it is an oblique91 argument for our time-independent (or frequency-domain) theory of stationary states.

Also, the diagrams of the self-energy are said to be limited to the linked and irreducible types only. These diagrammatic rules are often introduced on the basis of intuitive arguments and justified by somewhat imprecise mathematical logic. For instance, as shown in this study, the removal of reducible diagrams takes place as a result of systematic factorization of denominators,31,34,91,92 much like the cancellation of unlinked diagrams in MBPT.11,31,34,91–95 None of these delicate details seems to have been fully discussed. Furthermore, although diagrams are a powerful tool of derivation, they have pitfalls: The first report of MBGF(3) included twelve self-energy diagrams (of the Goldstone type) that are obtained by opening MBPT(3) diagrams by cutting one of its lines in all topologically distinct ways.53 Two years later, it was pointed out55 that there were six additional third-order self-energy diagrams, the so-called Σ() diagrams,96 obtained by a vertex insertion (see below) to the MBPT(2) diagrams. It is, therefore, not impossible to overlook an entire class of diagrams, which may start to appear only at some high order (although, more precisely, they are legitimately absent in the equation-of-motion formalism53 of MBGF).

The present authors have, therefore, been searching for an algebraic (first-quantized) definition of MBGF in the style of the Rayleigh–Schrödinger recursion equations of MBPT, with a view to justifying the diagrammatic derivation and to implementing a general-order MBGF method.

To this end, two of the present authors with two coauthors previously implemented97 what we call the ΔMPn method, originated by Pickup and Goscinski7 and extended by Chong et al.,98–100 in a determinant-based, general-order algorithm. In this method, the nth-order perturbation correction to the electron binding energy of a Koopmans state is defined as the difference in the MBPT(n) correlation correction between the N- and (N±1)-electron systems using the same N-electron Hartree–Fock (HF) reference. The MBPT(n) corrections are, in turn, defined unambiguously by the Rayleigh–Schrödinger recursion equations for any number of electrons to an arbitrary high perturbation order.

We numerically confirmed97 that ΔMPn reduces to MBGF(n) with the self-energy in the diagonal, frequency-independent approximation at 1n3 but converges at the nonapproximated MBGF(n) at n=, whose self-energy is nondiagonal and frequency dependent. This method, therefore, switches from the most approximate form of the self-energy at low orders to the nonapproximated one at an infinite order. This intriguing behavior occurs because it contains two unexpected and bizarre classes of diagrams, which we call the semi-reducible and linked-disconnected diagrams.97 They act to recuperate the effects of off-diagonal elements and frequency dependence of self-energy, respectively, in a perturbative manner. This finding exacerbated our concern that we may be missing an unknown number of important diagrams in higher-order MBGF.

In this study, we return to the core definition88 of a perturbation theory, which equates the nth-order perturbation correction of a given quantity to the nth derivative with respect to λ of the same quantity computed exactly with a perturbation-scaled Hamiltonian, Ĥ=H^0+λV^, where H^0 is the zeroth-order Hamiltonian and V^ is the fluctuation potential. This implies a universal strategy of deriving an algebraic definition of virtually any perturbation theory, which consists in analytically differentiating the FCI expression of a target quantity. Computationally, it also suggests the λ-variation method, which evaluates the low-order perturbation corrections by numerical differentiation of the FCI results of the perturbation-scaled Hamiltonian, allowing a general-order implementation of any perturbation theory, including perturbation series of the Green’s function and self-energy.

Adopting this strategy, we obtain three sets of Rayleigh–Schrödinger-like recursion equations defining a perturbation expansion of the exact one-particle Green’s function. They can be combined with the well-known recursion equation for the nondiagonal, frequency-dependent (i.e., nonapproximated) self-energy9 to determine its perturbation corrections to an arbitrarily high order. Using these algebraic definitions, we show, in an entirely time-independent framework, that unlinked-diagrammatic contributions in the Green’s functions cancel with one another, leaving only linked ones, after systematic factorization of their denominators.31,34,91,92 In other words, the factorization and linked-diagram theorems11,31,34,91–95 of MBPT in the time-independent form of Manne91 hold for the one-particle Green’s function, justifying its linked-diagrammatic perturbation expansion and proving its size-consistency.101,102 The factorization theorem31,34,91,92 also holds for the self-energy, demonstrating the mutual cancellation of reducible diagrams and thus justifying their irreducible-diagrammatic expansion. The most significant work in this area thus far is due to Kutzelnigg and Mukherjee,25 whose time-independent formulation differs from ours in that they adopt the second-quantized superoperator algebra and switch to diagrammatics early, not providing an algebraic recursion.

We also present six algorithms of general-order MBGF, five of which give identical perturbation series of electron binding energies of both Koopmans and non-Koopmans (shake-up or satellite) states. The first three implement the three recursions for the one-particle Green’s function in a determinant-based, general-order algorithm. The fourth is based on the λ-variation method, using only slightly modified FCI code and finite-difference approximations for the λ-differentiation. The fifth is the one that automatically generates all irreducible and linked self-energy diagrams order-by-order and synthesizes codes that evaluate them, precisely following three basic rules (see below) of the diagrammatic MBGF. The fact that all these five methods agree identically up to fifth order (within the precision of the finite-difference approximations of the λ-variation method) means that the existing diagrammatic rules are complete and correct, which is also implied by the linked- and irreducible-diagram theorems for the self-energy presented in this article. The last (sixth) general-order algorithm is the determinant-based ΔMPn method97 extended in this study to non-Koopmans states using the Hirschfelder–Certain degenerate perturbation theory.103 It gives a different perturbation series than MBGF but still converges at the FCI results for both Koopmans and non-Koopmans states. We also quantify the effect of the diagonal and/or frequency-independent approximation on the self-energy.

Furthermore, we illustrate a purely algebraic derivation of the first- and second-order self-energies, also highlighting its relationship with a diagrammatic derivation. The algebraic derivation is far more tedious than diagrammatic or superoperator derivation but is thoroughly systematic, involving only simple arithmetic operations and an equivalent of the Slater–Condon rules, which can thus place the diagrammatic technique on a firm mathematical footing. It also serves as a pedagogical example explaining some of the basic algebraic and diagrammatic mechanisms by which unlinked and reducible diagrams are erased. Finally, one of the three recursions is used to elucidate the precise relationship between ΔMPn and MBGF at several low perturbation orders.

The exact one-particle Green’s function9,18,30 of a closed-shell molecule is an m-by-m matrix (where m is the number of spinorbitals) whose elements are dependent on frequency (ω), namely,

Gpq(ω)=ΨN,0|p^(ωEN,0+H^)1q^|ΨN,0ΨN,0|ΨN,0+  ΨN,0|q^(ωH^+EN,0)1p^|ΨN,0ΨN,0|ΨN,0,
(1)

where ΨN,μ and EN,μ are the exact (FCI) N-electron (unnormalized) wave function and energy for the μth state (μ=0 being the ground state) and p^ (q^) is a creation (annihilation) operator of an electron in the pth (qth) spinorbital. Throughout this article, we adhere to the convention in which i, j, k, l, etc., refer to a canonical spinorbital occupied in the N-electron, ground-state HF wave function, a, b, c, d, etc., refer to a virtual canonical spinorbital, and p, q, r, s, etc., refer to either. N is an even number.

Inserting the resolution-of-the-identity, we can rewrite this as

Gpq(ω)=μΨN,0|p^|ΨN1,μΨN1,μ|q^|ΨN,0ΨN,0|ΨN,0(ωEN,0+EN1,μ)+μΨN,0|q^|ΨN+1,μΨN+1,μ|p^|ΨN,0ΨN,0|ΨN,0(ωEN+1,μ+EN,0),
(2)

where μ runs over all exact (N±1)-electron states. This indicates that the exact Green’s function diverges whenever ω agrees with an exact electron-detachment (EN,0EN1,μ) or electron-attachment (EN+1,μEN,0) energy. Therefore, the Green’s function matrix, which is as small as the Fock matrix, contains an entire spectrum of electron binding energies, which are exponentially many, by virtue of being frequency dependent. The ingenuity of MBGF as an exact one-particle theory can be traced to this definition.

There is hardly any difficulty constructing this exact Green’s function numerically from FCI wave functions and energies for N- and (N±1)-electron states and characterizing its behavior.72,97 What is at issue is the formulation and numerical calculation of a whole series of size-consistent perturbation approximations to this exact Green’s function. To define such a series, we adopt the Møller–Plesset partitioning104 of the exact Hamiltonian Ĥ,

Ĥ=H^(0)+λV^(1),
(3)

where the zeroth-order Hamiltonian H^(0) is the Fock operator. At this stage, λ, which is equal to unity, is merely there to facilitate the determination of the overall perturbation rank of each term (it will play a more active role in the λ-variation method). The zeroth-order or HF Green’s function is defined as

Gpq(0)(ω)=ΦN,0(0)|p^(ωEN,0(0)+H^(0))1q^|ΦN,0(0)ΦN,0(0)|ΦN,0(0)
+ΦN,0(0)|q^(ωH^(0)+EN,0(0))1p^|ΦN,0(0)ΦN,0(0)|ΦN,0(0)
(4)
=δpqΔωp,
(5)

where Δωp=ω𝜖p and 𝜖p is the energy of the pth canonical HF molecular spinorbital (MO). The zeroth-order wave function, ΦN,0(0), is the HF wave function for the ground state, and EN,0(0) is the sum of occupied orbital energies (different from the HF energy). Henceforth, we use the letter “Ψ” for a (potentially) multi-Slater-determinant wave function and the letter “Φ” for a single determinant wave function, and the perturbation rank is always indicated in the parenthesized superscript and those without such a qualifier are usually exact (FCI) counterparts. Equation (5) indicates that the zeroth-order Green’s function diverges whenever ω agrees with an electron binding energy in the Koopmans approximation.30 

In a typical exposition of MBGF,9,18,30 the exact and zeroth-order Green’s functions are said to be related to each other via the self-energy, Σ(ω), which together satisfy the so-called Dyson equation

G(ω)=G(0)(ω)+G(0)(ω)Σ(ω)G(ω),
(6)

where every factor in this equation is an m-by-m matrix whose elements are frequency dependent. It is, however, more accurate and instructive to characterize this equation as merely a definition of Σ(ω) in terms of G(ω) and G(0)(ω). In other words, Σ(ω) does not contain any new information which is not already encoded in G(ω), and we need to define G(ω) or its perturbative approximation before we can know Σ(ω) (not the other way around).

Multiplying from the left the inverse of G(0)(ω) and multiplying from the right the inverse of G(ω) with Eq. (6), we obtain the inverse Dyson equation

{G(0)(ω)}1={G(ω)}1+Σ(ω).
(7)

Recall that the usefulness of the exact Green’s function lies in the fact that it diverges when ω agrees with an exact electron binding energy. We seek such poles of G(ω). Since G(ω) does not have an inverse at these poles, the determinant of the left-hand side vanishes there. The working equation for ω that reports an exact electron binding energy, therefore, reads as

{G(0)(ω)}1Σ(ω)=0.
(8)

Furthermore, since

G(0)(ω)=(ω1𝜖)1,
(9)

where 1 is an m-by-m unit matrix and 𝜖 is the diagonal Fock matrix (𝜖pq=δpq𝜖p) in the canonical HF MO basis, Eq. (8) is equivalent to an m-by-m matrix eigenvalue equation of the form

{𝜖+Σ(ω)}U=Uω
(10)

or

U{𝜖+Σ(ω)}U=ω,
(11)

where U is a unitary matrix and ω is a diagonal matrix. This has a similar form and the same dimension as the Hartree–Fock–Roothaan equation and may be viewed as a correlation-corrected one-particle equation in the HF MO basis. The self-energy, Σ(ω), is then regarded as the matrix representation of the one-particle correlation operator, not unlike a correlation potential in KS DFT. The column vector of U whose corresponding eigenvalue is a self-consistent solution of Eq. (10) defines a new orbital that includes correlation effects (again, not unlike a KS orbital), which is known as the Dyson orbital. However, MBGF differs from KS DFT in that the exact limit of MBGF is known and a converging, systematic series of approximations is also available for Σ(ω), which is furthermore non-multiplicative, usually perturbative (non-variational), and frequency dependent. The beauty of the Dyson equation, in our view, lies in the fact that Σ(ω) defined by it has this (DFT-like) compelling physical meaning, while lending itself to a systematic perturbation expansion.

Equation (10) is a recursion equation in which unknown ω appears both in the left- and right-hand sides, which is solved typically by repeated diagonalizations in an iterative root search for ω. There are three simplifying approximations. In the diagonal approximation, this diagonalization step is avoided by the neglect of all off-diagonal elements of the self-energy,

𝜖p+Σpp(ω)=ω,   (diagonal),
(12)

which still requires a root search. In the frequency-independent approximation, on the other hand, the recursive structure of the Dyson equation is broken by substituting ω=𝜖p in the left-hand side,

U{𝜖+Σ(𝜖p)}U=ω,   (ω-independent),
(13)

which can be solved by one diagonalization for each state for which ω𝜖p (valid only for Koopmans states). In the diagonal, frequency-independent approximation, ω is obtained by one-shot evaluation (without diagonalization or root search) of the left-hand side of

𝜖p+Σpp(𝜖p)=ω,   (diagonal,ω-independent).
(14)

This is also appropriate only for Koopmans states.

However, the most important approximation is the perturbation expansion of Σ(ω). Since Σ(ω) is defined by G(ω) through the Dyson equation, the perturbation series of Σ(ω) is derived from the same of G(ω). Setting aside for the moment how the latter series is obtained (which is the focus of this work and will be thoroughly discussed below), here we first derive the algebraic equations connecting the two series.

First, we assume that the exact Green’s function is expanded in a perturbation series with the same Møller–Plesset partitioning of the Hamiltonian [Eq. (3)] as

G(ω)=G(0)(ω)+λG(1)(ω)+λ2G(2)(ω)+.
(15)

We also write the self-energy in a perturbation series as

Σ(ω)=λΣ(1)(ω)+λ2Σ(2)(ω)+λ3Σ(3)(ω)+,
(16)

in which the term proportional to λ0 is not present by definition (or it is zero). Substituting these series into the Dyson equation (6) and collecting terms with the same power of λ, we obtain

G(0)=G(0),
(17)
G(1)=G(0)Σ(1)G(0),
(18)
G(2)=G(0)Σ(2)G(0)+G(0)Σ(1)G(1),
(19)
G(3)=G(0)Σ(3)G(0)+G(0)Σ(2)G(1)+G(0)Σ(1)G(2),
(20)

etc., where the frequency argument is omitted for brevity. Generally, for n1, we have

G(n)=G(0)i=1nΣ(i)G(ni).
(21)

This can be inverted to yield an algebraic recursive definition of the self-energy in terms of lower-order self-energies and Green’s functions of the equal and lower orders,

Σ(n)=(G(0))1G(n)(G(0))1i=1n1Σ(i)G(ni)(G(0))1
(22)
=(G(0))1G(n)(G(0))1i=1n1Σ(i)G(0)Σ(ni)i=1n2j=1ni1Σ(i)G(0)Σ(j)G(0)Σ(nij),
(23)

where it should be understood (henceforth) that the terms with null range of summation are to be neglected. This recursion equation is considered well known.9 

Adhering to the conventional exposition of MBGF, we now switch to diagrammatic techniques. A perturbation correction to the self-energy is defined as a sum of Feynman–Goldstone diagrams, which are enumerated and algebraically interpreted following some sets of rules.9,11,18,87,89 These rules are mathematically justified by the Gell-Mann–Low time-dependent perturbation theory.90 However, this time-dependent logic is rather incongruous with the time-independent (frequency-domain) formulation employed up to this point and does not offer a concrete strategy for deriving programmable working equations of even low-order self-energies, let alone algebraic recursive definitions that can be used to implement a general-order algorithm. The primary objective of this work is, therefore, to fill this gap by presenting exactly the latter, i.e., easy-to-understand (if tedious) algebraic recursive definitions of both Green’s function and self-energy, which explains (if not prove) the diagrammatic rules rather than rely on them. Here, we delay this discussion but document the established diagrammatic rules first.

There are three basic rules for enumerating all nth-order self-energy diagrams. They are illustrated in Figs. 1–3, using only skeleton diagrams in which hole/particle distinctions are suppressed. Diagrammatic terminology is summarized in Table I.

TABLE I.

Diagrammatic terminology.

Closed diagram A diagram with no external line 
Open diagram A diagram with two external lines 
Disconnected diagram A diagram that consists of at least two 
 (either closed or open) parts, between 
 which no connecting path exists 
Connected diagram A diagram that is not disconnected 
Unlinked diagram A disconnected diagram with at least one 
 closed disconnected part 
Linked diagram A diagram that is not unlinked, 
 which is connected or disconnected 
 with open parts only 
Reducible diagram A self-energy diagram with at least one 
 pair of vertexes connected by a single line 
Irreducible diagram A self-energy diagram that is not reducible 
Insertion diagram A diagram product created by inserting 
 one or more inner diagrams into the outer 
 (principal) diagram. Its value is the product 
 of the values of disconnected parts whose 
 resolvent lines span only one part. 
 The outer diagram has two identical resolvent 
 lines above and below the point of insertion 
V vertex A filled circle, which represents V^ 
p vertex The terminus of the incoming dangling line, 
 which represents p^ 
q vertex The terminus of the outgoing dangling line, 
 which represents q^ 
Resolvent line A wiggly line denoting a denominator 
 consisting of orbital energies of the external 
 and internal lines intersected and possibly ω 
Internal line (edge) A line connecting two vertexes 
Long external line A line attached to only one vertex 
(dangling line)  
Short external line A dangling line trimmed at the vertex, 
(stub) which has no hole/particle distinction but 
 has in/out attribute 
Articulation line A single line connecting two vertexes, 
 the elimination of which leaves 
 the diagram disconnected 
Closed diagram A diagram with no external line 
Open diagram A diagram with two external lines 
Disconnected diagram A diagram that consists of at least two 
 (either closed or open) parts, between 
 which no connecting path exists 
Connected diagram A diagram that is not disconnected 
Unlinked diagram A disconnected diagram with at least one 
 closed disconnected part 
Linked diagram A diagram that is not unlinked, 
 which is connected or disconnected 
 with open parts only 
Reducible diagram A self-energy diagram with at least one 
 pair of vertexes connected by a single line 
Irreducible diagram A self-energy diagram that is not reducible 
Insertion diagram A diagram product created by inserting 
 one or more inner diagrams into the outer 
 (principal) diagram. Its value is the product 
 of the values of disconnected parts whose 
 resolvent lines span only one part. 
 The outer diagram has two identical resolvent 
 lines above and below the point of insertion 
V vertex A filled circle, which represents V^ 
p vertex The terminus of the incoming dangling line, 
 which represents p^ 
q vertex The terminus of the outgoing dangling line, 
 which represents q^ 
Resolvent line A wiggly line denoting a denominator 
 consisting of orbital energies of the external 
 and internal lines intersected and possibly ω 
Internal line (edge) A line connecting two vertexes 
Long external line A line attached to only one vertex 
(dangling line)  
Short external line A dangling line trimmed at the vertex, 
(stub) which has no hole/particle distinction but 
 has in/out attribute 
Articulation line A single line connecting two vertexes, 
 the elimination of which leaves 
 the diagram disconnected 

Rule 1 (“linked only”) states that, starting with a linked closed diagram of the MBPT(n) correction to the energy, we “cut” one of its lines and create an open diagram with two long external (dangling) lines. We then “trim” these dangling lines to obtain an nth-order self-energy diagram, containing n interaction vertexes (filled circles), n − 1 resolvent lines (wiggly lines), and two short external lines (or stubs). This process is illustrated in Fig. 1. Since the initial closed energy diagram is linked (as per the linked-diagram theorem11,31,34,91–95 of MBPT), the resulting open self-energy diagram is also always linked, ensuring its size-consistency. One can loosely associate diagram 2 of Fig. 1 with G(n) in the first term of Eq. (23). The multiplications of (G(0))1 from both sides of it then correspond to the trimming of the dangling lines (note that the physical meaning of a line is a zeroth-order Green’s function), exposing Σ(n) as diagram 3. However, as we show in  Appendix C, this is too simplistic an argument mathematically speaking, although it captures the correct physics.

FIG. 1.

Rule 1 (“linked only”). Linked, closed energy diagram 1 of MBPT(3) is cut open at the red wedge symbol, creating linked, open, third-order Green’s-function diagram 2 with two dangling lines. The dangling lines are trimmed at the wedges (red) to expose linked, open, third-order self-energy diagram 3 with two stubs (invisible). Hole/particle distinctions are suppressed.

FIG. 1.

Rule 1 (“linked only”). Linked, closed energy diagram 1 of MBPT(3) is cut open at the red wedge symbol, creating linked, open, third-order Green’s-function diagram 2 with two dangling lines. The dangling lines are trimmed at the wedges (red) to expose linked, open, third-order self-energy diagram 3 with two stubs (invisible). Hole/particle distinctions are suppressed.

Close modal

If we used only rule 1, we would be double-counting the so-called reducible self-energy diagrams. A reducible diagram is the one with at least one pair of vertexes (filled circles) connected by just a single line (“articulation line”), the elimination of which leaves the diagram disconnected. Rule 2 (“irreducible only”) requires the deletion of all such diagrams, an example of which is given in Fig. 2. The second and subsequent terms in Eq. (23) are a rough physical justification of this rule,9 but this is again not entirely accurate mathematically. A more rigorous justification of this rule is given in  Appendix C as the irreducible-diagram theorem.

FIG. 2.

Rule 2 (“irreducible only”). Linked, closed energy diagram 4 of MBPT(4) is cut open at the wedge (red), and then the two dangling lines (diagram 5) are trimmed, giving rise to reducible self-energy diagram 6, characterized by an articulation line connecting two second-order self-energy subdiagrams. Diagram 6 should be deleted. Hole/particle distinctions are suppressed.

FIG. 2.

Rule 2 (“irreducible only”). Linked, closed energy diagram 4 of MBPT(4) is cut open at the wedge (red), and then the two dangling lines (diagram 5) are trimmed, giving rise to reducible self-energy diagram 6, characterized by an articulation line connecting two second-order self-energy subdiagrams. Diagram 6 should be deleted. Hole/particle distinctions are suppressed.

Close modal

If we used only rules 1 and 2, we could be overlooking a whole class of diagrams. Rule 3 (“vertex insertion”) dictates that additional nth-order self-energy diagrams are obtained by inserting a vertex with a bubble (or a tadpole) into an (n − 1)th-order energy diagram, cutting open the bubble, and then trimming the dangling lines, an example of which is shown in Fig. 3. As noted in the Introduction, such diagrams were not included in the initial report of MBGF(3). They do not usually emerge from rule 1 because the parent closed diagrams are zero in MBPT(3) with the HF reference and are, therefore, usually unlisted to begin with. The corresponding self-energy diagrams are numerically significant. Rules 1 and 3 can, however, be consolidated into one requiring to connect n vertexes in all topologically distinct ways to produce open diagrams with two stubs. Alternatively, rule 3 can be absorbed into rule 1 modified to first enumerate all non-HF MBPT(n) diagrams34 and cut open a bubble implicit in each of the Fock operator vertexes.105 

FIG. 3.

Rule 3 (“vertex insertion”). A vertex with a bubble is inserted into a line of closed energy diagram 7 of MBPT(2) at the red wedge (with a concomitant addition of a new resolvent line). The bubble is cut open (diagram 8), creating open diagram 9. The two dangling lines are trimmed, producing open, third-order self-energy diagram 10. Alternatively, consider diagram 8 as a non-HF energy diagram of MBPT(3) and follow rule 1. Hole/particle distinctions are suppressed.

FIG. 3.

Rule 3 (“vertex insertion”). A vertex with a bubble is inserted into a line of closed energy diagram 7 of MBPT(2) at the red wedge (with a concomitant addition of a new resolvent line). The bubble is cut open (diagram 8), creating open diagram 9. The two dangling lines are trimmed, producing open, third-order self-energy diagram 10. Alternatively, consider diagram 8 as a non-HF energy diagram of MBPT(3) and follow rule 1. Hole/particle distinctions are suppressed.

Close modal

The resulting diagrams are algebraically interpreted by the rules given in Table II. For instance, second-order self-energy diagrams 11 and 12 in Fig. 4 are interpreted as the first and second terms, respectively, of the right-hand side,

Σpq(2)=(1)1+112i,a,biq||abab||ipΔωiab+(1)2+112i,j,aaq||ijij||apΔijωa
(24)

with

Δωiab=ω+𝜖i𝜖a𝜖b,
(25)
Δijωa=𝜖i+𝜖jω𝜖a.
(26)
TABLE II.

Algebraic interpretation rules of self-energy diagrams in the Hugenholtz style.9 

(1) Label the incoming short external line with general orbital 
 index p and the outgoing short external line with q 
(2) Label downgoing internal lines with hole indices i, j, k, etc., 
 and upgoing internal lines with particle indices a, b, c, etc. 
(3) Associate each vertex with left out, right out||left in, right in 
(4) Between each adjacent pair of vertexes, draw a resolvent line 
(5) Associate each resolvent line with 1Δijkabc, if the number of 
 intersections is even; 1Δijkωab if the number of hole 
 intersections exceeds the number of particle intersections; 1Δωijabc 
 otherwise. Here, i,j,k, (a,b,c,) are the labels 
 of the hole (particle) lines intersecting the resolvent line 
(6) Sum over all internal lines 
(7) Multiply (−1)h+l to the diagram with h hole lines and l loops 
 A fictitious loop need not be considered 
(8) Multiply 1/n! for each set of n equivalent lines. Two lines are 
 equivalent when they start from the same vertex and end at 
 the same vertex as well as have the identical hole/particle attribute 
(1) Label the incoming short external line with general orbital 
 index p and the outgoing short external line with q 
(2) Label downgoing internal lines with hole indices i, j, k, etc., 
 and upgoing internal lines with particle indices a, b, c, etc. 
(3) Associate each vertex with left out, right out||left in, right in 
(4) Between each adjacent pair of vertexes, draw a resolvent line 
(5) Associate each resolvent line with 1Δijkabc, if the number of 
 intersections is even; 1Δijkωab if the number of hole 
 intersections exceeds the number of particle intersections; 1Δωijabc 
 otherwise. Here, i,j,k, (a,b,c,) are the labels 
 of the hole (particle) lines intersecting the resolvent line 
(6) Sum over all internal lines 
(7) Multiply (−1)h+l to the diagram with h hole lines and l loops 
 A fictitious loop need not be considered 
(8) Multiply 1/n! for each set of n equivalent lines. Two lines are 
 equivalent when they start from the same vertex and end at 
 the same vertex as well as have the identical hole/particle attribute 
FIG. 4.

The second-order self-energy diagrams.

FIG. 4.

The second-order self-energy diagrams.

Close modal

Two of the present authors with two coauthors generalized97 the ΔMPn method of Pickup and Goscinski7 and of Chong et al.98–100 to higher orders. It defines the nth-order perturbation correction to the μth electron binding energy, denoted by Σ¯μ(n), as the MBPT(n) correction difference between the N- and (N±1)-electron states using the same N-electron HF reference wave function,

Σ¯μ(n)=EN,0(n)EN1,μ(n),
(27)
Σ¯μ(n)=EN+1,μ(n)EN,0(n),
(28)

for n1, where μ labels an (N±1)-electron state with the MBPT(n) correction of EN±1,μ(n). It is well known7 that this difference reduces (at the formalism level) to the nth-order self-energy in the diagonal, frequency-independent approximation at n = 2, the fact used for pedagogical purposes by Szabo and Ostlund30 to explain the relationship between MBPT and MBGF. Nobes et al.106 and Beste et al.107 proposed similar methods, which account for orbital-relaxation effects and are more practically useful.

The MBPT(n) corrections are, in turn, defined unambiguously by the algebraic recursion equations derivable straightforwardly from the Rayleigh–Schrödinger perturbation theory for the ground state of the N-electron system and for all Koopmans states of the (N±1)-electron systems,

EM,μ(0)=ΨM,μ(0)|H^(0)|ΨM,μ(0),
(29)
EM,μ(n)=ΨM,μ(0)|V^(1)|ΨM,μ(n1),   (n1).
(30)

Here, the nth-order correction to the wave function of the μth M-electron state, ΨM,μ(n), is given by the recursion equation,

(EM,μ(0)H^(0))ΨM,μ(n)=V^(1)ΨM,μ(n1)i=1nEM,μ(i)ΨM,μ(ni),
(31)

which can be initiated with the zeroth-order wave function of the form

ΨM,μ(0)=ΦM,μ(0),
(32)

where ΦM,μ(0) is a single Koopmans-type Slater determinant.

In our previous study,97 the above recursions were implemented in a determinant-based algorithm to realize ΔMPn at any arbitrary order (n) for Koopmans states only. Using this method, we numerically confirmed that the ΔMPn correction (Σ¯μ(n)) agrees identically with the nth-order self-energy in the diagonal, frequency-independent approximation for 1n3,

Σ¯μ(n)=Σpp(n)(𝜖p),
(33)

where p labels the spinorbital that is vacated or filled in the μth Koopmans state. As n, ΔMPn approaches MBGF(n) with the nondiagonal, frequency-dependent self-energy. We speculated97 that this switch occurs because two unexpected classes of diagrams (the semi-reducible diagrams and linked-disconnected diagrams) are included in ΔMPn, which, respectively, correct the diagonal and frequency-independent approximations systematically. Examples of such diagrams are given in Fig. 5.

FIG. 5.

A fourth-order semi-reducible diagram (13) and linked-disconnected diagram (14). Diagram 15 is an alternative way of drawing diagram 14. In these ΔMP4 diagrams, the resolvent lines are evaluated at ω=𝜖p (i.e., in the frequency-independent approximation). Hole/particle distinctions are suppressed.

FIG. 5.

A fourth-order semi-reducible diagram (13) and linked-disconnected diagram (14). Diagram 15 is an alternative way of drawing diagram 14. In these ΔMP4 diagrams, the resolvent lines are evaluated at ω=𝜖p (i.e., in the frequency-independent approximation). Hole/particle distinctions are suppressed.

Close modal

One may intuitively understand that diagram 13 recuperates the effect of off-diagonal Σpq(2)(𝜖p) at the fourth order, which is already included in MBGF(2) but not in ΔMP2. Diagrams 14 and 15 have a disconnected or an irreducible part with two resolvent lines. The second resolvent line arises from the differentiation of the denominator with respect to ω because

ω1Δωiab=1(Δωiab)2.
(34)

Hence, these diagrams account for the effect of frequency dependence in the self-energy in a perturbative (Taylor-series) manner without explicit frequency dependence. In  Appendix D, we algebraically show the equivalence of ΔMPn and approximate MBGF(n) at 1n3 and that the differences between them at n = 4 indeed correspond to these semi-reducible and linked-disconnected diagrams.

In this study, we further extend the ΔMPn method to non-Koopmans states, motivated by its need in two of the three recursions of MBGF(n) to be described below. The perturbation corrections to the wave function and energy of a non-Koopmans state in an (N±1)-electron system satisfy the identical recursions from the Rayleigh–Schrödinger perturbation theory [Eqs. (29)–(31)]. However, they alone cannot determine the zeroth-order wave functions for degenerate non-Koopmans states to initiate the recursion. Such zeroth-order wave functions are not single determinants but their linear combinations,

ΨM,μ(0)=νUM,μν(0)ΦM,ν(0),
(35)

where ν runs over all states having the same zeroth-order energy with the μth state and UM,μν(0) is a unitary matrix element. These zeroth-order wave functions are determined by the recursion plus some additional “consistency” conditions103 on wave functions, which satisfy

ΨM,μ(0)|ΨM,ν(0)=δμν,
(36)
ΨM,μ(0)|ΨM,μ(n)=δn0.
(37)

A procedure to determine these zeroth-order wave functions (and higher-order ones) is exceedingly complex but is well established and unambiguously documented by Hirschfelder and Certain.103 We will not reproduce here the 21 pages of this degenerate Hirschfelder–Certain perturbation theory (HCPT), except to note that they can be implemented in a determinant-based algorithm to enable ΔMPn for non-Koopmans as well as Koopmans states at any arbitrary order. We also emphasize that, once determined, the resulting perturbation corrections to the wave function and energy satisfy the same Rayleigh–Schrödinger recursion equations given above and correspond to the terms in the usual expansions in λ,

ΨM,μ=ΨM,μ(0)+λΨM,μ(1)+λ2ΨM,μ(2)+λ3ΨM,μ(3)+,
(38)
EM,μ=EM,μ(0)+λEM,μ(1)+λ2EM,μ(2)+λ3EM,μ(3)+,
(39)

where M can be N, N + 1, or N − 1. It should be noted that MBPT is a special case of the degenerate HCPT in that the latter applied to a non-degenerate zeroth-order state reduces exactly to the former. The HCPT is said to be equivalent to several other degenerate perturbation theories variously named.108 

The λ-variation method88 can evaluate virtually any perturbation series numerically up to high orders. It computes the nth-order perturbation correction of any given target quantity as the nth λ-derivative of the same quantity calculated by the FCI method with a scaled Hamiltonian, H^(0)+λV^(1).

For instance, the perturbation corrections to the Green’s function and self-energy are evaluated as

Gpq(n)=1n!nGpqλnλ=0,
(40)
Σpq(n)=1n!nΣpqλnλ=0,
(41)

where the exact (FCI) G and Σ are defined by Eqs. (1) and (7), respectively, and evaluated accordingly. The derivatives are obtained by a finite-difference numerical differentiation109,110 of these quantities calculated with a FCI program slightly modified so as to permit any given value for λ. There is a minimal risk of introducing formulation or programming errors.

An algebraic recursive definition of virtually any perturbation series can be obtained by analytically differentiating with λ the exact (FCI) definition of a quantity with a scaled Hamiltonian, H^(0)+λV^(1). The nth derivative is identified as the nth-order perturbation correction.

Using the resolution-of-the-identity with exact (N±1)-electron wave functions, the definition of the m-by-m matrix of the exact one-particle Green’s function (where m is the number of orbitals) becomes

Gpq(ω)=μ,νΨN,0|p^|ΨN1,μΨN1,μ|(ωEN,0+H^)1|ΨN1,νΨN1,ν|q^|ΨN,0ΨN,0|ΨN,0+μ,νΨN,0|q^|ΨN+1,μΨN+1,μ|(ωH^+EN,0)1|ΨN+1,νΨN+1,ν|p^|ΨN,0ΨN,0|ΨN,0.
(42)

Every factor except p^, q^, and ω is dependent on λ. The λ-derivatives of the wave functions and energies of the N- and (N±1)-electron states are already known; they are the perturbation corrections of the respective quantities given by MBPT or HCPT [Eqs. (38) and (39)].

Let us first consider the λ-derivative of the second factor in the numerator of each term,

ΨN±1,μ|(ω±EN,0H^)1|ΨN±1,ν  =GN±1,μν(0)+λGN±1,μν(1)+λ2GN±1,μν(2)+,
(43)

where GN±1,μν(n) (the frequency argument is omitted for brevity) is the nth-order perturbation correction to a many-particle Green’s function, which is an M-by-M matrix where M is the number of all (N±1)-electron states. Hereafter, a “Green’s function” without a qualifier refers to the one-particle Green’s function, G or Gpq.

Since the bra and ket wave functions are the eigenfunctions of Ĥ, these factors can be simplified, reducing Eq. (42) into Eq. (2). However, we must not make such simplification here because once these wave functions are expanded in a perturbation series, they are neither eigenfunctions nor orthogonal to one another across different states. In other words, whereas the left-hand side of Eq. (43) may be diagonal, the individual terms in the right-hand side are not. For a similar reason, we need to retain the normalizing denominator, ΨN,0|ΨN,0, in each term of Eq. (42); while the exact N-electron wave function may be presumed normalized, its perturbation expansions are not. On the other hand, the (N±1)-electron wave functions used in the resolution of the identity are normalized by the second factor in each numerator of Eq. (42), which involves an inverse operator.

Using the matrix identity,9 

(AB)1=A1+(AB)1BA1
(44)
=A1+A1B(AB)1
(45)
=A1+A1BA1+A1BA1BA1+A1BA1BA1BA1+
(46)

with

Aμν=ΨN±1,μ(0)|ω±EN,0(0)H^(0)|ΨN±1,ν(0),
(47)
Bμν=n=1i=0nΨN±1,μ(i)|ωH^(0)|ΨN±1,ν(ni)  ±n=1i=0n1ΨN±1,μ(i)|V^(1)|ΨN±1,ν(ni1)  n=1i=0nj=0niEN,0(j)ΨN±1,μ(i)|ΨN±1,ν(nij),
(48)

we obtain the recursion for GN±1,μν(n) as

GN±1,μν(n)=i=1nκGN±1,μκ(ni)VN±1,κν(i)GN±1,νν(0)
(49)

with

VN±1,μν(n)=i=0nΨN±1,μ(i)|ωH^(0)|ΨN±1,ν(ni)±i=0n1ΨN±1,μ(i)|V^(1)|ΨN±1,ν(ni1)i=0nj=0niEN,0(j)ΨN±1,μ(i)|ΨN±1,ν(nij),
(50)

which is initiated with the zeroth-order quantities given by

GN±1,μν(0)=δμν(ω±EN,0(0)EN±1,μ(0))1,
(51)

where EN,0(0) or EN±1,μ(0) is the sum of the HF energies of the spinorbitals occupied in the corresponding Slater determinant.

With the recursion of the many-particle Green’s function, the same for the one-particle Green’s function can be obtained as

Gpq(n)=i=0nj=0niμ,νzN1,μp(i)*GN1,μν(j)zN1,νq(nij)+i=0nj=0niμ,νzN+1,μq(i)*GN+1,μν(j)zN+1,νp(nij)i=1nD(i)Gpq(ni).
(52)

The reader is reminded that the above equation is valid for n = 0 with the last term with null summation range omitted and can, therefore, be used to initiate the recursion. Here, z’s are the perturbation corrections to the first and third factors in each numerator of Eq. (42) and are given by

zN1,νq(n)=i=0nΨN1,ν(i)|q^|ΨN,0(ni),
(53)
zN+1,νp(n)=i=0nΨN+1,ν(i)|p^|ΨN,0(ni).
(54)

In Eq. (52), D(i) originates from the normalizing denominators of Eq. (42). Using the scalar version of Eq. (44), the reciprocal of the denominator is expanded as

ΨN,0|ΨN,01=1λD(1)λ2D(2)+λ2D(1)D(1)λ3D(3)+λ3D(2)D(1)+λ3D(1)D(2)λ3D(1)D(1)D(1)+
(55)

with

D(1)=ΨN,0(0)|ΨN,0(1)+ΨN,0(1)|ΨN,0(0),
(56)
D(2)=ΨN,0(0)|ΨN,0(2)+ΨN,0(1)|ΨN,0(1)+ΨN,0(2)|ΨN,0(0),
(57)

and so forth, or more generally,

D(n)=i=0nΨN,0(i)|ΨN,0(ni).
(58)

This can be further simplified with ΨN,0(n)|ΨN,0(0)=δn0, but we delay this simplification until  Appendix B.

We have so far found it exceedingly difficult to reach this or any of the subsequent recursions using the superoperator algebra.5,7,9,18,111 In its present form, this language seems intractable because of the ambiguity of the choice of bra and ket wave functions when evaluating superoperator expectation values, their lack of Hermiticity, and other missing details needed for reliable derivations at all orders. However, it should be possible to sharpen this language into an even more useful form25 perhaps with the aid of the algebraic recursions presented here.

We have also been unsuccessful in using Löwdin’s perturbation theory based on the inner-projection technique112 to derive these recursions. As Öhrn and Born17 implied and Kutzelnigg and Mukherjee25 pointed out, this tends to treat correlation in the Green’s-function operator and ground-state wave function differently, unlike in the diagrammatic Feynman–Dyson perturbation theory or the foregoing λ-derivative derivation, where these two types of correlations are included on an equal footing. The inner-projection technique may instead be used to define a different but equally valid (perhaps even superior) perturbation series.

The same can be said about a CC expansion to the Green’s function,85,86 leading to frequency-independent diagrams, which also decouple (N + 1)- and (N − 1)-electron sectors. The time- or frequency-dependence and coupling of the two sectors is inevitable in the Feynman–Dyson perturbation expansion of the Green’s function, given its origin as a propagator in quantum electrodynamics,113–115 which “measures the local response of the field at a given point at a later time to a local disturbance of the field at another given point at an earlier time.”116 The CC expansion offers an alternative, systematic approximation series,72,117–123 which is potentially superior to an MBPT expansion underlying the Feynman–Dyson MBGF. See also Sec. 4.3 of Ref. 97, in which the relationship between the ΔMPn and CC theories for electron detachment and attachment is elucidated.

The resolution-of-the-identity operator inserted in Eq. (1) may be composed of any complete set of (N±1)-electron wave functions. Using the one made up of single Slater determinants, we can write the exact Green’s function as

Gpq(ω)=μ,νΨN,0|p^|ΦN1,μ(0)ΦN1,μ(0)|(ωEN,0+H^)1|ΦN1,ν(0)ΦN1,ν(0)|q^|ΨN,0ΨN,0|ΨN,0+μ,νΨN,0|q^|ΦN+1,μ(0)ΦN+1,μ(0)|(ωH^+EN,0)1|ΦN+1,ν(0)ΦN+1,ν(0)|p^|ΨN,0ΨN,0|ΨN,0.
(59)

Differentiation with respect to λ of this definition leads to an alternative, but equivalent, recursion for G’s. This definition has an advantage over the previous one, Eq. (42), in that the single Slater determinants do not depend on λ and their λ-derivatives need not be taken. The resulting recursion equations are, therefore, simpler and a general-order algorithm implementing them tends to be more stable and efficient, not requiring HCPT energies or wave functions.

Perturbation corrections to the many-particle Green’s functions are then defined as

ΦN±1,μ(0)|(ω±EN,0H^)1|ΦN±1,ν(0)  =G̃N±1,μν(0)+λG̃N±1,μν(1)+λ2G̃N±1,μν(2)+,
(60)

which are shown to obey the following recursion:

G̃N±1,μν(n)=±κG̃N±1,μκ(n1)ṼN±1,κν(1)G̃N±1,νν(0)i=1nG̃N±1,μν(ni)EN,0(i)G̃N±1,νν(0)
(61)

with

ṼN±1,μν(1)=ΦN±1,μ(0)|V^(1)|ΦN±1,ν(0).
(62)

These are obtained by differentiating the left-hand side of Eq. (60) with λ or by substituting in Eq. (44),

Aμν=ΦN±1,μ(0)|ω±EN,0(0)H^(0)|ΦN±1,ν(0),
(63)
Bμν=±ΦN±1,μ(0)|V^(1)|ΦN±1,ν(0)n=1EN,0(n)ΦN±1,μ(0)|ΦN±1,ν(0).
(64)

The recursion can be initiated with the zeroth-order conditions,

G̃N±1,μν(0)=δμν(ω±EN,0(0)EN±1,μ(0))1.
(65)

The many-particle Green’s function, G̃(n), of this subsection differs from G(n) in Sec. II E for n1.

With the many-particle Green’s function and the perturbation corrections to the first and third factors in each numerator of Eq. (59) given by

z̃N1,νq(n)=ΦN1,ν(0)|q^|ΨN,0(n),
(66)
z̃N+1,νp(n)=ΦN+1,ν(0)|p^|ΨN,0(n),
(67)

we can write the recursive definition of the one-particle Green’s function as

Gpq(n)=i=0nj=0niμ,νz̃N1,μp(i)*G̃N1,μν(j)z̃N1,νq(nij)+i=0nj=0niμ,νz̃N+1,μq(i)*G̃N+1,μν(j)z̃N+1,νp(nij)i=1nD(i)Gpq(ni),
(68)

for n1, where D(n) is already defined by Eq. (58). This is also valid for n = 0 with the last term omitted.

For its simplicity, this recursion is used for most numerical calculations in this article as well as in algebraic derivations of Σ(1) and Σ(2) in  Appendix A. It also serves as the basis of the linked-diagram and irreducible-diagram theorems of MBGF discussed in  Appendixes B and  C. We can have some glimpse of the outlines of these arguments in the recursion equations. The second term in Eq. (61) is the so-called renormalization term and consists of simple products (i.e., not tensor contractions with at least one common index) of lower-order MBPT energies and many-particle Green’s functions. They, therefore, correspond to unlinked diagrams, violating size-consistency if they persist. That these are subtracted from the parent (first) term suggests that the parent term also contains the same unlinked diagrams, which are canceled exactly by the second term, after some systematic factorization of the denominators in the former (see  Appendix B). Likewise, the last term in Eq. (68) is another renormalization term and is unlinked, which is expected to annihilate the same in the first two terms, after systematic denominator factorization. The subtraction of these manifestly unlinked terms should leave only the linked contributions in the one-particle Green’s function. The details are given in  Appendix B.

Yet another recursive definition of the perturbation series of the one-particle Green’s function can be obtained by inserting the resolution-of-the-identity in the basis of two complete sets of the zeroth-order and exact (FCI) wave functions of the (N±1)-electron states,

Gpq(ω)=μ,νΨN,0|p^|ΨN1,μΨN1,μ(0)|(ωEN,0+H^)1|ΨN1,νΨN1,ν(0)|q^|ΨN,0ΨN,0|ΨN,0+μ,νΨN,0|q^|ΨN+1,μΨN+1,μ(0)|(ωH^+EN,0)1|ΨN+1,νΨN+1,ν(0)|p^|ΨN,0ΨN,0|ΨN,0,
(69)

where ΨN±1,μ(0) is the zeroth-order HCPT wave function for the μth state, as given by Eq. (35). For a Koopmans state, it is a single Slater determinant and may be denoted by ΦN±1,μ(0), but in a more general case, it is a linear combination of degenerate Slater determinants, whose expansion coefficients are unknown a priori and need to be determined by the HCPT.103 The two sets are individually complete, justifying Eq. (69). We consider this expansion because the resulting recursion equations contain the ΔMPn expressions, thus clarifying the relationship between ΔMPn and MBGF, which is to be fully elucidated in  Appendix D.

As before, the following numerator factors can be expanded in a perturbation series as

ΨN±1,μ(0)|(ω±EN,0H^)1|ΨN±1,ν  =G¯N±1,μν(0)+λG¯N±1,μν(1)+λ2G¯N±1,μν(2)+,
(70)

where the many-particle Green’s function, G¯(n) (which differs from either G(n) or G̃(n)), is defined recursively as

G¯N±1,μν(n)=i=1nκG¯N±1,μκ(ni)V¯N±1,κν(i)G¯N±1,νν(0)
(71)

with

V¯N±1,μν(n)=(ω±EN,0(0)EN±1,μ(0))ΨN±1,μ(0)|ΨN±1,ν(n)+i=1nΨN±1,μ(0)|ΨN±1,ν(ni)Σ¯ν(i),
(72)

where Σ¯ν(n) is the ΔMPn correction for the νth state, defined by Eq. (27) or (28). This equation can be obtained by simply setting i = 0 in the summations of Eq. (50) and making use of the recursion [Eq. (31)] of the Rayleigh–Schrödinger perturbation theory, which is satisfied by a HCPT wave function, ΨN±1,ν(n). It is initiated with the zeroth-order quantities

G¯N±1,μν(0)=δμν(ω±EN,0(0)EN±1,μ(0))1.
(73)

The recursive definition of the one-particle Green’s function is then given by

Gpq(n)=i=0nj=0niμ,νxN1,μp(i)G¯N1,μν(j)yN1,νq(nij)+i=0nj=0niμ,νxN+1,μq(i)G¯N+1,μν(j)yN+1,νp(nij)i=1nD(i)Gpq(ni),
(74)

for n1 (with the last term omitted when n = 0) with

xN1,μp(n)=i=0nΨN,0(i)|p^|ΨN1,μ(ni),
(75)
yN1,νq(n)=ΨN1,ν(0)|q^|ΨN,0(n),
(76)
xN+1,μq(n)=i=0nΨN,0(i)|q^|ΨN+1,μ(ni),
(77)
yN+1,νp(n)=ΨN+1,ν(0)|p^|ΨN,0(n).
(78)

We have developed a symbolic computing program124 that generates self-energy diagrams according to the three diagram-enumeration rules illustrated in Figs. 1–3. It then synthesizes codes that evaluate the algebraic formulas of these diagrams obtained with the diagram-interpretation rules in Table II. This can be used to give order-by-order perturbation corrections to the self-energy up to high orders, although it may not be considered as a general-order algorithm, strictly speaking.

Table III gives the numerical values of the self-energy matrix of the BH molecule calculated in this way. The details of the calculation are given in the table caption. From here on, p and q label spatial orbitals and only the αα spin block of the self-energy matrix is shown. The αβ and βα blocks are zero, and the ββ block is the same as the αα block.

TABLE III.

The nth-order perturbation correction (in Eh) to the self-energy, Σ(n), of BH (rB-H=1.232 Å) at ω=0.2Eh obtained by the automatic generation and evaluation of diagrams. The minimal (STO-3G) basis set and the frozen core approximation were used. The highest-occupied MO (HOMO) corresponds to p = q = 3. The HF and FCI energies are −24.752 788 and 24.809629Eh, respectively.

Σpq(2)q = 23456
p = 2 −0.010 239 −0.000 984 0.000 000 −0.000 000 0.011 268 
−0.000 984 0.001 304 0.000 000 −0.000 000 0.007 050 
0.000 000 0.000 000 0.003 328 −0.000 000 0.000 000 
−0.000 000 −0.000 000 −0.000 000 0.003 328 −0.000 000 
0.011 268 0.007 050 0.000 000 −0.000 000 0.015 802 
Σpq(2)q = 23456
p = 2 −0.010 239 −0.000 984 0.000 000 −0.000 000 0.011 268 
−0.000 984 0.001 304 0.000 000 −0.000 000 0.007 050 
0.000 000 0.000 000 0.003 328 −0.000 000 0.000 000 
−0.000 000 −0.000 000 −0.000 000 0.003 328 −0.000 000 
0.011 268 0.007 050 0.000 000 −0.000 000 0.015 802 
Σpq(3)q = 23456
p = 2 −0.003 632 −0.004 702 0.000 000 −0.000 000 0.008 025 
−0.004 702 −0.004 586 0.000 000 −0.000 000 0.005 288 
0.000 000 0.000 000 0.012 711 −0.000 000 0.000 000 
−0.000 000 −0.000 000 −0.000 000 0.012 711 −0.000 000 
0.008 025 0.005 288 0.000 000 −0.000 000 0.009 849 
Σpq(3)q = 23456
p = 2 −0.003 632 −0.004 702 0.000 000 −0.000 000 0.008 025 
−0.004 702 −0.004 586 0.000 000 −0.000 000 0.005 288 
0.000 000 0.000 000 0.012 711 −0.000 000 0.000 000 
−0.000 000 −0.000 000 −0.000 000 0.012 711 −0.000 000 
0.008 025 0.005 288 0.000 000 −0.000 000 0.009 849 
Σpq(4)q = 23456
p = 2 −0.000 858 −0.004 281 0.000 000 −0.000 000 0.004 964 
−0.004 281 −0.004 399 0.000 000 −0.000 000 0.003 296 
0.000 000 0.000 000 0.011 417 −0.000 000 0.000 000 
−0.000 000 −0.000 000 −0.000 000 0.011 417 −0.000 000 
0.004 964 0.003 296 0.000 000 −0.000 000 0.005 991 
Σpq(4)q = 23456
p = 2 −0.000 858 −0.004 281 0.000 000 −0.000 000 0.004 964 
−0.004 281 −0.004 399 0.000 000 −0.000 000 0.003 296 
0.000 000 0.000 000 0.011 417 −0.000 000 0.000 000 
−0.000 000 −0.000 000 −0.000 000 0.011 417 −0.000 000 
0.004 964 0.003 296 0.000 000 −0.000 000 0.005 991 
Σpq(5)q = 23456
p = 2 0.000 316 −0.003 118 0.000 000 −0.000 000 0.003 026 
−0.003 118 −0.003 129 0.000 000 −0.000 000 0.001 957 
0.000 000 0.000 000 0.009 122 −0.000 000 0.000 000 
−0.000 000 −0.000 000 −0.000 000 0.009 122 −0.000 000 
0.003 026 0.001 957 0.000 000 −0.000 000 0.003 596 
Σpq(5)q = 23456
p = 2 0.000 316 −0.003 118 0.000 000 −0.000 000 0.003 026 
−0.003 118 −0.003 129 0.000 000 −0.000 000 0.001 957 
0.000 000 0.000 000 0.009 122 −0.000 000 0.000 000 
−0.000 000 −0.000 000 −0.000 000 0.009 122 −0.000 000 
0.003 026 0.001 957 0.000 000 −0.000 000 0.003 596 

Table IV lists perturbation corrections to the self-energy of the BH molecule, the same system as Table III but obtained with the λ-variation method. The symmetric seven-point finite-difference formulas109,110 were used for the first through fifth derivatives (corresponding to the first- through fifth-order perturbation corrections to the self-energy) with an evenly spaced grid centered at λ=0 with the grid spacing of Δλ=0.01. Very tight convergence of the FCI roots is essential for more precise results.

TABLE IV.

Same as Table III but obtained by the λ-variation method.

Σpq(2)q = 23456
p = 2 −0.010 239 0.000 984 0.000 000 0.000 000 −0.011 268 
0.000 984 0.001 304 0.000 000 0.000 000 0.007 050 
0.000 000 0.000 000 0.003 328 −0.000 000 0.000 000 
0.000 000 0.000 000 −0.000 000 0.003 328 −0.000 000 
−0.011 268 0.007 050 0.000 000 −0.000 000 0.015 802 
Σpq(2)q = 23456
p = 2 −0.010 239 0.000 984 0.000 000 0.000 000 −0.011 268 
0.000 984 0.001 304 0.000 000 0.000 000 0.007 050 
0.000 000 0.000 000 0.003 328 −0.000 000 0.000 000 
0.000 000 0.000 000 −0.000 000 0.003 328 −0.000 000 
−0.011 268 0.007 050 0.000 000 −0.000 000 0.015 802 
Σpq(3)q = 23456
p = 2 −0.003 632 0.004 702 −0.000 000 −0.000 000 −0.008 025 
0.004 702 −0.004 586 −0.000 000 −0.000 000 0.005 288 
−0.000 000 −0.000 000 0.012 711 0.000 000 −0.000 000 
−0.000 000 −0.000 000 0.000 000 0.012 711 −0.000 000 
−0.008 025 0.005 288 −0.000 000 −0.000 000 0.009 849 
Σpq(3)q = 23456
p = 2 −0.003 632 0.004 702 −0.000 000 −0.000 000 −0.008 025 
0.004 702 −0.004 586 −0.000 000 −0.000 000 0.005 288 
−0.000 000 −0.000 000 0.012 711 0.000 000 −0.000 000 
−0.000 000 −0.000 000 0.000 000 0.012 711 −0.000 000 
−0.008 025 0.005 288 −0.000 000 −0.000 000 0.009 849 
Σpq(4)q = 23456
p = 2 −0.000 858 0.004 281 −0.000 000 −0.000 000 −0.004 964 
0.004 281 −0.004 398 0.000 000 −0.000 000 0.003 296 
−0.000 000 0.000 000 0.011 416 0.000 000 0.000 000 
−0.000 000 −0.000 000 0.000 000 0.011 417 0.000 000 
−0.004 964 0.003 296 0.000 000 0.000 000 0.005 994 
Σpq(4)q = 23456
p = 2 −0.000 858 0.004 281 −0.000 000 −0.000 000 −0.004 964 
0.004 281 −0.004 398 0.000 000 −0.000 000 0.003 296 
−0.000 000 0.000 000 0.011 416 0.000 000 0.000 000 
−0.000 000 −0.000 000 0.000 000 0.011 417 0.000 000 
−0.004 964 0.003 296 0.000 000 0.000 000 0.005 994 
Σpq(5)q = 23456
p = 2 0.000 312 0.003 120 0.000 000 0.000 000 −0.003 028 
0.003 120 −0.003 131 0.000 000 0.000 000 0.001 958 
0.000 000 0.000 000 0.009 137 0.000 000 0.000 000 
0.000 000 0.000 000 0.000 000 0.009 126 0.000 004 
−0.003 028 0.001 958 0.000 000 0.000 004 0.003 593 
Σpq(5)q = 23456
p = 2 0.000 312 0.003 120 0.000 000 0.000 000 −0.003 028 
0.003 120 −0.003 131 0.000 000 0.000 000 0.001 958 
0.000 000 0.000 000 0.009 137 0.000 000 0.000 000 
0.000 000 0.000 000 0.000 000 0.009 126 0.000 004 
−0.003 028 0.001 958 0.000 000 0.000 004 0.003 593 

The fact that the two data sets in Tables III and IV agree with each other within the precision of finite-difference approximations (which seems no worse than 10−5Eh at the fifth order) strongly suggests (but never proves) that the diagram-enumeration and diagram-interpretation rules are complete and correct to all orders. Some off-diagonal elements of the self-energy have their signs reversed between the two data sets simply because of the arbitrariness of MO phases (the data in Table III were obtained with a different HF program than the rest of the results).

The ΔMPn method97 was extended to non-Koopmans states by implementing a general-order HCPT method103 using the determinant-based algorithm.66 Our implementation was verified by the λ-variation method, which computed up to the fourth-order HCPT corrections to the energies of non-Koopmans states. They were evaluated as the λ-derivatives of the FCI energies with the forward seven- or nine-point finite difference formulas with λ=0 as the lower end (not the center) of an evenly spaced grid with Δλ=0.01. This asymmetic grid choice was to safeguard against the states that are degenerate at λ=0 but become non-degenerate and can be ordered differently between λ<0 and λ>0. These derivatives were found to agree (within the precision of the finite-difference approximations) with our general-order HCPT implementation at several low perturbation orders (not shown). The HCPT wave functions and energies were also used in implementing two of the three algebraic recursions.

Figure 6 plots the ΔMPn electron binding energies of the electron-detached states [Eq. (27)] of the BH molecule. They converge rapidly at the FCI results (open circles) with increasing n for both Koopmans (red lines) and non-Koopmans states (blue lines) except for several non-Koopmans states (green lines) where they display clear signs of divergence. The most striking observation, however, is the massive first-order corrections to the energies of non-Koopmans states, bringing the latter in line with the FCI results. This is in sharp contrast with null first-order corrections for Koopmans states. This underscores the extremely large orbital-relaxation effect on non-Koopmans states, the majority of which seems to be captured by a mixing of degenerate single determinants [Eq. (35)] at the first order.

FIG. 6.

The ΔMPn electron binding energies of the (N − 1)-electron states of the BH molecule (the same as in Table III) as a function of the perturbation order (n).

FIG. 6.

The ΔMPn electron binding energies of the (N − 1)-electron states of the BH molecule (the same as in Table III) as a function of the perturbation order (n).

Close modal

We defer the comparison of the ΔMPn data with those from various approximations of MBGF until Sec. III D.

Recursion I, II, or III was invoked to determine high-order perturbation corrections to the one-particle Green’s function and self-energy. Taking recursion I as an example, the computational procedure is outlined as follows: (i) High-order perturbation corrections to the ground-state wave function and energy of the N-electron system are obtained with the determinant-based, general-order MBPT algorithm;73 (ii) the same for the ground and all excited states of the (N±1)-electron system are then determined by the determinant-based, general-order HCPT algorithm103 described in Sec. III C; (iii) with these, z’s and D’s are generated68,72 up to sufficiently high orders; (iv) for a given ω, recursion equations (49) and (50) are evaluated68,72 to construct higher-order perturbation corrections to the many-particle Green’s function (G) in the form of M-by-M matrices; (v) next, the recursion equation [Eq. (52)] is used to generate higher-order perturbation corrections to the one-particle Green’s function (G) in the form of m-by-m matrices; (vi) the perturbation corrections to the self-energy matrix (Σ) can be obtained recursively using Eq. (23).

Recursion II renders step (ii) unnecessary, making the whole calculation more stable and faster (relatively speaking because at any order these calculations cost more than a FCI calculation). All numerical data in this section are, therefore, obtained with recursion II.

Table V lists the perturbation corrections to the self-energy of the BH molecule obtained by the general-order algorithm outlined above using recursion II. We verified that the other two recursions give the identical results to all shown digits. They also agree with the data obtained by the automatic generation of diagrams (Table III) apart from nonessential sign change and also with the results of the λ-variation calculations within the numerical precision of the latter. This observation underscores the full mathematical equivalence of the three algebraic recursions, λ-variation, and diagrammatics.

TABLE V.

Same as Table III but obtained by a general-order MBGF algorithm based on the algebraic recursion II. The results obtained with recursion I or III are identical.

Σpq(2)q = 23456
p = 2 −0.010 239 0.000 984 0.000 000 0.000 000 −0.011 268 
0.000 984 0.001 304 0.000 000 −0.000 000 0.007 050 
0.000 000 0.000 000 0.003 328 0.000 000 −0.000 000 
0.000 000 −0.000 000 −0.000 000 0.003 328 −0.000 000 
−0.011 268 0.007 050 −0.000 000 −0.000 000 0.015 802 
Σpq(2)q = 23456
p = 2 −0.010 239 0.000 984 0.000 000 0.000 000 −0.011 268 
0.000 984 0.001 304 0.000 000 −0.000 000 0.007 050 
0.000 000 0.000 000 0.003 328 0.000 000 −0.000 000 
0.000 000 −0.000 000 −0.000 000 0.003 328 −0.000 000 
−0.011 268 0.007 050 −0.000 000 −0.000 000 0.015 802 
Σpq(3)q = 23456
p = 2 −0.003 632 0.004 702 0.000 000 0.000 000 −0.008 025 
0.004 702 −0.004 586 0.000 000 −0.000 000 0.005 288 
0.000 000 0.000 000 0.012 711 0.000 000 −0.000 000 
0.000 000 −0.000 000 0.000 000 0.012 711 −0.000 000 
−0.008 025 0.005 288 −0.000 000 −0.000 000 0.009 849 
Σpq(3)q = 23456
p = 2 −0.003 632 0.004 702 0.000 000 0.000 000 −0.008 025 
0.004 702 −0.004 586 0.000 000 −0.000 000 0.005 288 
0.000 000 0.000 000 0.012 711 0.000 000 −0.000 000 
0.000 000 −0.000 000 0.000 000 0.012 711 −0.000 000 
−0.008 025 0.005 288 −0.000 000 −0.000 000 0.009 849 
Σpq(4)q = 23456
p = 2 −0.000 858 0.004 281 0.000 000 0.000 000 −0.004 964 
0.004 281 −0.004 399 0.000 000 −0.000 000 0.003 296 
0.000 000 0.000 000 0.011 417 0.000 000 −0.000 000 
0.000 000 −0.000 000 0.000 000 0.011 417 −0.000 000 
−0.004 964 0.003 296 −0.000 000 −0.000 000 0.005 991 
Σpq(4)q = 23456
p = 2 −0.000 858 0.004 281 0.000 000 0.000 000 −0.004 964 
0.004 281 −0.004 399 0.000 000 −0.000 000 0.003 296 
0.000 000 0.000 000 0.011 417 0.000 000 −0.000 000 
0.000 000 −0.000 000 0.000 000 0.011 417 −0.000 000 
−0.004 964 0.003 296 −0.000 000 −0.000 000 0.005 991 
Σpq(5)q = 23456
p = 2 0.000 316 0.003 118 0.000 000 0.000 000 −0.003 026 
0.003 118 −0.003 129 0.000 000 −0.000 000 0.001 957 
0.000 000 0.000 000 0.009 122 0.000 000 −0.000 000 
0.000 000 −0.000 000 0.000 000 0.009 122 −0.000 000 
−0.003 026 0.001 957 −0.000 000 −0.000 000 0.003 596 
Σpq(5)q = 23456
p = 2 0.000 316 0.003 118 0.000 000 0.000 000 −0.003 026 
0.003 118 −0.003 129 0.000 000 −0.000 000 0.001 957 
0.000 000 0.000 000 0.009 122 0.000 000 −0.000 000 
0.000 000 −0.000 000 0.000 000 0.009 122 −0.000 000 
−0.003 026 0.001 957 −0.000 000 −0.000 000 0.003 596 
Σpq(6)q = 23456
p = 2 0.000 709 0.002 080 0.000 000 0.000 000 −0.001 836 
0.002 080 −0.002 031 0.000 000 −0.000 000 0.001 124 
0.000 000 0.000 000 0.007 036 0.000 000 −0.000 000 
0.000 000 −0.000 000 0.000 000 0.007 036 −0.000 000 
−0.001 836 0.001 124 −0.000 000 −0.000 000 0.002 117 
Σpq(6)q = 23456
p = 2 0.000 709 0.002 080 0.000 000 0.000 000 −0.001 836 
0.002 080 −0.002 031 0.000 000 −0.000 000 0.001 124 
0.000 000 0.000 000 0.007 036 0.000 000 −0.000 000 
0.000 000 −0.000 000 0.000 000 0.007 036 −0.000 000 
−0.001 836 0.001 124 −0.000 000 −0.000 000 0.002 117 
Σpq(7)q = 23456
p = 2 0.000 725 0.001 318 0.000 000 0.000 000 −0.001 089 
0.001 318 −0.001 282 0.000 000 −0.000 000 0.000 628 
0.000 000 0.000 000 0.005 286 0.000 000 −0.000 000 
0.000 000 −0.000 000 0.000 000 0.005 286 −0.000 000 
−0.001 089 0.000 628 −0.000 000 −0.000 000 0.001 207 
Σpq(7)q = 23456
p = 2 0.000 725 0.001 318 0.000 000 0.000 000 −0.001 089 
0.001 318 −0.001 282 0.000 000 −0.000 000 0.000 628 
0.000 000 0.000 000 0.005 286 0.000 000 −0.000 000 
0.000 000 −0.000 000 0.000 000 0.005 286 −0.000 000 
−0.001 089 0.000 628 −0.000 000 −0.000 000 0.001 207 

Table VI compares the MBGF(n) electron binding energies of the HOMO (highest-occupied MO) of the same system in various (diagonal, frequency-independent, and both) approximations to the self-energy against the one without any such approximations (“full”). The ΔMPn electron binding energy of the same orbital is also included.

TABLE VI.

The electron binding energy (in Eh) of the HOMO (p = 3 and 𝜖p=0.246538Eh) of the BH molecule (the same as in Table III) calculated by MBGF(n) using the self-energy in the diagonal, ω-independent, both, or no approximation (“full”). The electron binding energy of the same orbital in the ΔMPn approximation is also shown.

nFullDiagonalω-indep.Diagonal, ω-indep.ΔMPn
−0.246 54 −0.246 54a −0.246 54 −0.246 54a −0.246 54a 
−0.246 54 −0.246 54a −0.246 54 −0.246 54a −0.246 54a 
−0.244 11 −0.244 07a −0.244 05 −0.244 00a −0.244 00a 
−0.247 69 −0.247 61a −0.247 74 −0.247 66a −0.247 66a 
−0.251 13 −0.251 13 −0.251 40 −0.251 40 −0.251 52 
−0.253 45 −0.253 56 −0.253 92 −0.254 04 −0.254 18 
−0.254 86 −0.255 06 −0.255 47 −0.255 71 −0.255 71 
−0.255 69 −0.255 96 −0.256 41 −0.256 73 −0.256 51 
20 −0.257 00 −0.257 37 −0.257 92 −0.258 37 −0.257 00 
 −0.257 00a −0.257 37a −0.257 92 −0.258 37a −0.257 00a 
nFullDiagonalω-indep.Diagonal, ω-indep.ΔMPn
−0.246 54 −0.246 54a −0.246 54 −0.246 54a −0.246 54a 
−0.246 54 −0.246 54a −0.246 54 −0.246 54a −0.246 54a 
−0.244 11 −0.244 07a −0.244 05 −0.244 00a −0.244 00a 
−0.247 69 −0.247 61a −0.247 74 −0.247 66a −0.247 66a 
−0.251 13 −0.251 13 −0.251 40 −0.251 40 −0.251 52 
−0.253 45 −0.253 56 −0.253 92 −0.254 04 −0.254 18 
−0.254 86 −0.255 06 −0.255 47 −0.255 71 −0.255 71 
−0.255 69 −0.255 96 −0.256 41 −0.256 73 −0.256 51 
20 −0.257 00 −0.257 37 −0.257 92 −0.258 37 −0.257 00 
 −0.257 00a −0.257 37a −0.257 92 −0.258 37a −0.257 00a 
a

See also Ref. 97.

As described above, the self-energies were derived from the corresponding one-particle and many-particle Green’s functions evaluated first. The latter are always divergent in the frequency-independent approximation at ω=𝜖p because they contain dangling lines (and/or articulation lines) with a divergent resolvent factor of 1Δpω (see  Appendix C). Therefore, the self-energies in the frequency-independent approximation at ω=𝜖p were determined by a linear interpolation from those on an evenly spaced grid of ω’s avoiding this pole. In this table, we only show digits that are safely correct.

The MBGF(n) data with the three approximations to the self-energy are all different from one another at n2. These differences persist in the n= limit. The effects of the diagonal and frequency-independent approximations are roughly additive. This can be understood graphically from Fig. 7, in which the left- and right-hand sides of the Dyson equation [Eq. (11) or (12)] are plotted as a function of ω. Let us focus on the blue curve, which plots the left-hand side of the Dyson equation [Eq. (11)] with the self-energy summing up to the twentieth-order perturbation correction. The intersections (only one is visible, indicated by the blue circle) of this curve with the diagonal line, ω, are the roots of the Dyson equation with no approximation; the blue circle corresponds to the MBGF(20) electron binding energy with no approximation, which is practically exact (FCI). The green curve is a locus of the left-hand side of the Dyson equation [Eq. (12)] with the self-energy in the diagonal approximation. Its intersections (again only one is visible, which is the green cross) with the diagonal ω-line are the MBGF(20) electron binding energies in the diagonal approximation. On the other hand, the values of these curves at ω=𝜖p (indicated by arrows, where p labels the HOMO) are the electron binding energies in the frequency-independent approximation. Therefore, the red diamond occurs at the electron binding energy in the frequency-independent approximation, whereas the orange square in the diagonal and frequency-independent approximations. Because of the near-perfect parallelism of the two (blue and green) curves, the difference in the frequency-independent electron binding energy with and without the diagonal approximation is nearly the same in the frequency-dependent case. This explains the near-additivity of the effects of the two approximations.

FIG. 7.

The left-hand side of the Dyson equation [Eq. (11) or (12)] as a function of ω for the HOMO of the BH molecule (the same as in Table III) including up to the nth-order self-energy (n = 1, 2, 3, 4, 20). The intersections of these curves with the diagonal line, ω, are the roots (electron binding energies) of the Dyson equation. The vertical arrows occur at ω=𝜖p, where p corresponds to the HOMO, indicating the roots in the frequency-independent approximation.

FIG. 7.

The left-hand side of the Dyson equation [Eq. (11) or (12)] as a function of ω for the HOMO of the BH molecule (the same as in Table III) including up to the nth-order self-energy (n = 1, 2, 3, 4, 20). The intersections of these curves with the diagonal line, ω, are the roots (electron binding energies) of the Dyson equation. The vertical arrows occur at ω=𝜖p, where p corresponds to the HOMO, indicating the roots in the frequency-independent approximation.

Close modal

Furthermore, the small gradients of the curves mean that the frequency-independent approximation is accurate. Small gradients, in turn, are the result of the absence of nearby poles in the self-energy and may be expected to hold true for the HOMO and LUMO (lowest-unoccupied MO) in most any systems. The errors caused by the diagonal approximation are even smaller in this example, but it may be unsafe to generalize this observation obtained using the tiny basis set with extremely few off-diagonal self-energy matrix elements. It is, however, probably safe to state that the perturbation expansion of the self-energy tends to be the greatest source of errors.

As already established,97 the ΔMPn self-energy agrees with the one in the diagonal, frequency-independent approximation at 1n3 but converges at the full self-energy in the FCI limit. These data substantiate our speculation97 that ΔMPn does not agree with MBGF(n) or any of its approximation in the range of 4n< and, therefore, constitutes a distinct perturbation series.

Figures 8 and 9 plot the electron binding energies of the HOMO and LUMO, respectively, of the BH molecule as a function of the perturbation order calculated by MBGF(n) and ΔMPn. For this tiny system, the diagonal approximation has no effect on the LUMO electron binding energy.

FIG. 8.

The electron binding energies of the HOMO of the BH molecule (the same as in Table III) as a function of the perturbation order.

FIG. 8.

The electron binding energies of the HOMO of the BH molecule (the same as in Table III) as a function of the perturbation order.

Close modal
FIG. 9.

The electron binding energies of the LUMO of the BH molecule (the same as in Table III) as a function of the perturbation order. The effect of the diagonal approximation to the self-energy is null in this tiny example.

FIG. 9.

The electron binding energies of the LUMO of the BH molecule (the same as in Table III) as a function of the perturbation order. The effect of the diagonal approximation to the self-energy is null in this tiny example.

Close modal

The following observations can be made from these figures: (1) The second-order corrections can have the wrong sign, sometimes causing the second-order electron binding energies to be surprisingly poor; (2) at the third and higher orders, convergence is monotonic and uniform (as opposed to staircase-like as in MBPT) but not particularly rapid, making higher-order corrections relatively more important and worthwhile. It takes third- or even fourth-order corrections to undo the errors introduced at the second order; (3) the errors caused by the diagonal and/or frequency-independent approximations increase with the perturbation order but always stay much less than the variation with the perturbation order. These errors may be underestimated by low-order MBGF studies and will not reach their full magnitudes until at the tenth order or so; (4) for the HOMO, surprisingly, ΔMPn converges noticeably more rapidly than full MBGF, despite the fact that the former includes the off-diagonal elements and frequency dependence of the self-energy more approximately than the latter especially at low orders. This is, however, likely the result of cancellation of errors between the diagonal and frequency-independent approximations (which underestimate the electron binding energies) and low-order perturbation approximations (which overestimate the same). In fact, ΔMPn converges more slowly than full MBGF for the LUMO.

An important property of MBGF is the presence of multiple roots to a one-particle Dyson equation owing to the frequency dependence of the self-energy. If it were not for this, no one-particle equation could describe an exact theory as it would have very few roots as compared with far more electron-detached and attached many-electron states.

Figure 10 is the same as Fig. 7 with the displayed range of ω expanded to include the other roots of the Dyson equation than the Koopmans state. Here, the colorful curves (including the flat black line) are the left-hand side of the Dyson equation in the diagonal approximation [Eq. (12)], whereas the black diagonal line is ω itself. The intersections between them seen in ω<1.2Eh are the roots corresponding to non-Koopmans states, whereas those around ω0.25Eh are the Koopmans state (the area enlarged in Fig. 7).

FIG. 10.

Same as Fig. 7 in a wider range of ω.

FIG. 10.

Same as Fig. 7 in a wider range of ω.

Close modal

Figure 11 plots the same but for the LUMO. In this figure, non-Koopmans solutions are visible both in ω<0.7Eh and ω>1.2Eh, while the Koopmans roots are concentrated in a small area around ω0.27Eh. Although the perturbative self-energies beyond the first order are invariably divergent at many ω’s, the intersections (roots) can never be divergent for the simple reason that they must also be on the diagonal ω-line, which is analytic everywhere. This is how MBGF can describe a strongly correlated wave function, while being a perturbation theory (recall that non-Koopmans states are multi-determinantal even in their zeroth-order wave functions and are deemed strongly correlated). Put another way, since MBGF(n) is exact at n=, it cannot be divergent even at n<.

FIG. 11.

Same as Fig. 7 or 10 but for the LUMO.

FIG. 11.

Same as Fig. 7 or 10 but for the LUMO.

Close modal

Interestingly, because of the symmetric and antisymmetric forms of the poles in the odd- and even-order self-energies, respectively, the intersections corresponding to non-Koopmans roots seem to exist only in even-order MBGF, which is not too surprising considering that MBGF(0) and MBGF(1) clearly cannot have any non-Koopmans roots. Figure 12 plots the positions of these intersections around ω0.8Eh in Fig. 11 as a function of the perturbation order, together with the ΔMPn electron binding energies of what is likely the same state. For the reason stated above, we have the MBGF(n) data only at even n’s (data for n>8 were unavailable owing to too sharp a rise of the poles). ΔMPn, in contrast, reports electron binding energies of non-Koopmans states at any order, including at the zeroth, first, and odd orders. Furthermore, for this particular example, it is more rapidly convergent at the FCI limit. These observations may be taken to support ΔMPn over MBGF(n), but it should be recalled that the former is frequently divergent, while the latter (without the frequency-independent approximation) is not.

FIG. 12.

A non-Koopmans solution of the Dyson equation in the diagonal approximation to the self-energy for the LUMO of the BH molecule (the same as in Table III) as a function of the perturbation order. The ΔMPn electron binding energy for the same state is also plotted.

FIG. 12.

A non-Koopmans solution of the Dyson equation in the diagonal approximation to the self-energy for the LUMO of the BH molecule (the same as in Table III) as a function of the perturbation order. The ΔMPn electron binding energy for the same state is also plotted.

Close modal

Figures 13 and 14 plot the electron binding energies of the HOMO−1 and LUMO, respectively, of the H2O molecule as a function of the perturbation order. They generally support the conclusions drawn from the BH molecular data. MBGF(2) tends to overcorrect the correlation in electron binding energies much more than MBPT(2) does in total energies, sometimes performing as poorly as MBGF(0) (i.e., Koopmans’ theorem). This makes MBGF(3) more important than MBPT(3); we do not observe the staircase-like convergence in MBGF but a more monotonic and uniform convergence after the second order. The errors incurred by the diagonal and/or frequency-independent approximation are less significant than those from low-order perturbation approximations but can be comparable to the fourth-order corrections. In these two examples, ΔMPn converges less rapidly than full MBGF, reaffirming that the opposite case found in Fig. 8 is accidental.

FIG. 13.

The electron binding energies of the HOMO−1 of the H2O molecule (rO–H=0.967 Å, aH–O–H=107.6°) as a function of the perturbation order. The minimal (STO-3G) basis set and the frozen core approximation were used. The HOMO−1 corresponds to p = q = 4. The HF and FCI energies are −74.962 663 and 75.012842Eh, respectively.

FIG. 13.

The electron binding energies of the HOMO−1 of the H2O molecule (rO–H=0.967 Å, aH–O–H=107.6°) as a function of the perturbation order. The minimal (STO-3G) basis set and the frozen core approximation were used. The HOMO−1 corresponds to p = q = 4. The HF and FCI energies are −74.962 663 and 75.012842Eh, respectively.

Close modal
FIG. 14.

Same as Fig. 13 but for the LUMO.

FIG. 14.

Same as Fig. 13 but for the LUMO.

Close modal

Much of modern many-body physics depends on quantum field theory,113–115 in which a time- or frequency-dependent Green’s function plays a crucial role,116,125 which is represented by a (dressed) line in a linked Feynman diagram expressing a size-consistent perturbation correction. This article reports advances in the analytical, algorithmic, and numerical aspects of high-order MBGF for nonrelativistic molecular quantum mechanics.

The advances in the analytical aspect include the derivation of three distinct algebraic recursive definitions of perturbative corrections to the one-particle Green’s function. It is to be combined with the well-known recursion for the perturbative self-energy. On the basis of one of the three recursions, we have derived, purely algebraically, the spinorbital expressions of the first- and second-order self-energies. The derivation is tedious but first-principles and extensible to any order.

We have also shown that the one-particle Green’s function and self-energy are diagrammatically linked and size-consistent at any order, by applying the factorization theorem with insertion to unlinked diagrams, which are then found to cancel exactly the renormalization terms. The linked-diagram theorem in a time-independent picture just like the one by Manne91 for MBPT also holds for MBGF.

The sum of the one-particle Green’s-function diagrams that differ from one another in the time orderings of the p or q vertex is shown to be the one in which the spans of the resolvent lines are broken, facilitating the “trimming” of those dangling lines and exposing a self-energy diagram. The same systematic factorization of the resolvent lines intersecting articulation lines justifies the removal of all reducible self-energy diagrams. These are the main assertions of the irreducible-diagram theorem.

We have also presented an algebraic proof of the equivalence of ΔMPn and MBGF(n) in the diagonal, frequency-independent approximation at 1n3 and the assignment of the difference at n = 4 to the semi-reducible and linked-disconnected diagrams.

On the algorithmic aspect, we have implemented six pieces of computer code (many sharing the same subroutines) that can evaluate the perturbation corrections to the one-particle Green’s function and self-energy up to an arbitrary high order. They implement the three recursions, automatic generation and interpretation of diagrams, λ-variation, and ΔMPn, all but the last yielding the identical, nondivergent perturbation series, mutually verifying the formulations and implementations. In fact, the λ-variation method not only furnishes numerical benchmark data of virtually any perturbation series with a minimal risk of formulation or programming errors but also suggests a universal strategy of deriving an algebraic recursive definition of a given perturbation theory, which has been adopted in this work.

Numerical results of these general-order MBGF calculations for the tiniest systems affordable have led to the following conclusions: MBGF(2) can be rather poor; MBGF(n) at n3 shows monotonic and uniform (as opposed to staircase-like) convergence; MBGF(3) is relatively more important for electron binding energies than MBPT(3) for total energies; the effects of the diagonal and frequency-independent approximations are roughly additive for Koopmans states and they do not manifest fully until high orders (such as at the tenth order); ΔMPn can sometimes converge more rapidly than full MBGF towards the FCI limit for either a Koopmans or non-Koopmans state owing to error cancellation in the former; MBGF(n) may lack solutions for non-Koopmans states at certain values of n.

S.H. has been supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Grant No. DE-FG02-11ER16211 and by CREST, Japan Science and Technology Agency. A.E.D. has been supported by the SciDAC program of the Department of Energy, Office of Science, Office of Basic Energy Sciences under Grant No. DE-FG02-12ER46875. J.V.O. has been supported by the National Science Foundation through Grant No. CHE-1565760 to Auburn University. Dr. Rodney J. Bartlett, Dr. Karol Kowalski, Dr. Marcel Nooijen, Dr. Yngve Öhrn, and Dr. Takeshi Yanai are thanked for illuminating discussions and advice.

1. Algebraic derivation

Here, we present a purely algebraic derivation of the first- and second-order self-energy expressions in the canonical HF reference, which is distinct from the usual diagrammatic derivation or the one based on the superoperator algebra. It instead relies on recursion II given in Sec. II F and is carried out with a straightforward, albeit tedious, application of simple arithmetics and an equivalent of the Slater–Condon rules (i.e., the normal-ordered second quantization and Wick’s theorem34 facilitated by a symbolic computing program126).

Let us first consider the expressions of z̃’s defined by Eqs. (66) and (67). According to Eq. (31), the first-order perturbation correction to the ground-state N-electron wave function is given by

|ΨN,0(1)=R^V^(1)|ΦN,0(0)R^EN,0(1)|ΦN,0(0)
(A1)
=12!2!(T2(1))ijab|Φijab,
(A2)

where

R^=(EN,0(0)H^(0))1,
(A3)
EN,0(1)=ΦN,0(0)|V^(1)|ΦN,0(0),
(A4)

and

(T2(1))ijab=ab||ijΔijab,
(A5)

where Φijab denotes a doubly excited N-electron determinant. Einstein’s convention of implied summations over repeated indices (that do not appear in the left-hand sides) is used when summation ranges are obvious and unimportant for the discussion. At the second order, we have

|ΨN,0(2)=R^V^(1)|ΨN,0(1)R^EN,0(1)|ΨN,0(1)R^EN,0(2)|ΦN,0(0)
(A6)
=(T1(2))ia|Φia+12!2!(T2(2))ijab|Φijab+13!3!(T3(2))ijkabc|Φijkabc+14!4!(T4(2))ijklabcd|Φijklabcd,
(A7)

where

EN,0(2)=ΦN,0(0)|V^(1)|ΨN,0(1)=12!2!ij||ab(T2(1))ijab,
(A8)

and

(T1(2))ia=12jk||ib(T2(1))jkbaΔia+12ja||bc(T2(1))jibcΔia,
(A9)
(T2(2))ijab=12kl||ij(T2(1))klabΔijab+P(ab)P(ij)kb||ic(T2(1))kjcaΔijab+12ab||cd(T2(1))ijcdΔijab,
(A10)
(T3(2))ijkabc=P(ab|c)P(ij|k)lc||ij(T2(1))lkabΔijkabc+P(a|bc)P(i|jk)bc||id(T2(1))jkdaΔijkabc,
(A11)
(T4(2))ijklabcd=P(ab|cd)P(ij|kl)cd||ij(T2(1))klabΔijklabcd,
(A12)

where P(ab\cd), etc., are index-permutation operators, whose precise definition can be found in pp. 284 and 285 of Shavitt and Bartlett.34 

We can then find the formulas of z̃’s contracted with appropriate determinants as follows:

|ΦN1,ν(0)z̃N1,νp(0)=δip|Φi,
(A13)
|ΦN+1,ν(0)z̃N+1,νp(0)=δap|Φa,
(A14)

at the zeroth order, where Φi and Φa are, respectively, the 1h and 1p Koopmans single determinants. At the first order, we have

|ΦN1,ν(0)z̃N1,νp(1)=12!δbp|Φija(T2(1))ijab+12!2!δkp|Φijkab(T2(1))ijab,
(A15)
|ΦN+1,ν(0)z̃N+1,νp(1)=12!δjp|Φiab(T2(1))ijab+12!2!δcp|Φijabc(T2(1))ijab.
(A16)

Here, for instance, Φija is a 2h1p non-Koopmans single determinant. At the second order, we obtain

|ΦN1,ν(0)z̃N1,νp(2)  =δap|Φi(T1(2))iaδjp|Φija(T1(2))ia   12!δbp|Φija(T2(2))ijab+12!2!δkp|Φijkab(T2(2))ijab   +12!3!δcp|Φijkab(T3(2))ijkabc13!3!δlp|Φijklabc(T3(2))ijkabc   13!4!δdp|Φijklabc(T4(2))ijklabcd   +14!4!δmp|Φijklmabcd(T4(2))ijklabcd
(A17)

and

|ΦN+1,ν(0)z̃N+1,νp(2)  =δip|Φa(T1(2))iaδbp|Φiab(T1(2))ia   +12!δjp|Φiab(T2(2))ijab+12!2!δcp|Φijabc(T2(2))ijab   13!2!δkp|Φijabc(T3(2))ijkabc13!3!δdp|Φijkabcd(T3(2))ijkabc   +14!3!δlp|Φijkabcd(T4(2))ijklabcd   +14!4!δep|Φijklabcde(T4(2))ijklabcd.
(A18)

Next, we construct the order-by-order algebraic formula of G̃ (the many-particle Green’s function) using Eq. (61). We then form G (the one-particle Green’s function) with Eq. (68) followed by the evaluation of Σ (the self-energy) via Eq. (23).

The zeroth-order many-particle Green’s functions, G̃N±1,μν(0), are already given by Eq. (65) and are diagonal. This immediately leads to

Gpq(0)=δpqΔωp.
(A19)

According to Eq. (61), the first-order many-particle Green’s functions are defined by

G̃N±1,μν(1)=±G̃N±1,μμ(0)ṼN±1,μν(1)G̃N±1,νν(0)δμνG̃N±1,μμ(0)EN,0(1)G̃N±1,μμ(0),
(A20)

which are further expanded into an operator form as

G^N1(1)=|ΦN1,μ(0)G̃N1,μν(1)ΦN1,ν(0)|=12|Φijaka||ijΦk|ΔωaijΔωk+12|Φkij||kaΦija|ΔωkΔωaij14|Φklaij||klΦija|ΔωaklΔωaij+|Φkjaia||jbΦkib|ΔωajkΔωaik14|Φijkabab||ijΦk|ΔωabijkΔωk14|Φkij||abΦijkab|ΔωkΔωabijk+12|Φijkabal||jkΦilb|ΔωabijkΔωbil+14|Φijkabab||icΦjkc|ΔωabijkΔωcjk+12|Φilbjk||alΦijkab|ΔωbilΔωabijk+14|Φjkcic||abΦijkab|ΔωcjkΔωabijk+
(A21)

and

G^N+1(1)=|ΦN+1,μ(0)G̃N+1,μν(1)ΦN+1,ν(0)|=12|Φiabab||icΦc|ΔωiabΔωc12|Φcic||abΦiab|ΔωcΔωiab+14|Φiabab||cdΦicd|ΔωiabΔωicd|Φjacia||jbΦibc|ΔωjacΔωibc+14|Φijabcab||ijΦc|ΔωijabcΔωc+14|Φcij||abΦijabc|ΔωcΔωijabc12|Φijabcab||diΦjcd|ΔωijabcΔωjcd+14|Φijabcak||ijΦkbc|ΔωijabcΔωkbc12|Φjcddi||abΦijabc|ΔωjcdΔωijabc+14|Φkbcij||akΦijabc|ΔωkbcΔωijabc+.
(A22)

These sums do not terminate until all terms involving up to the nh(n − 1)p and nh(n + 1)p sectors are reached, where n is the number of electrons. Note that the EN,0(1) factor is absent in the above expressions owing to mutual cancellation of terms containing it.

From Eq. (68), we then find

Gpq(1)=z̃N1,μp(0)*G̃N1,μν(1)z̃N1,νq(0)+z̃N1,μp(0)*G̃N1,μν(0)z̃N1,νq(1)+z̃N1,μp(1)*G̃N1,μν(0)z̃N1,νq(0)+z̃N+1,μq(0)*G̃N+1,μν(1)z̃N+1,νp(0)+z̃N+1,μq(0)*G̃N+1,μν(0)z̃N+1,νp(1)+z̃N+1,μq(1)*G̃N+1,μν(0)z̃N+1,νp(0)D(1)Gpq(0)=0.
(A23)

In the last equality, we have used D(1) = 0,

z̃N1,μp(0)*G̃N1,μν(0)z̃N1,νq(1)=0,
(A24)
z̃N1,μp(1)*G̃N1,μν(0)z̃N1,νq(0)=0,
(A25)

and the corresponding identities in the (N + 1)-electron sector. These can be understood by noting that z̃N1(0) spans only 1h determinants [Eq. (A13)], whereas z̃N1(1) spans the disjoint space of 2h1p determinants [Eq. (A15)]. The use is also made of the identities,

z̃N1,μp(0)*G̃N1,μν(1)z̃N1,νq(0)=0,
(A26)
z̃N+1,μq(0)*G̃N+1,μν(1)z̃N+1,νp(0)=0,
(A27)

which can be verified using similar logic applied to the structure of G̃(1) in Eqs. (A21) and (A22). We then deduce

Σpq(1)={(G(0))1G(1)(G(0))1}pq=0.
(A28)

Moving on to the second order, according to recursion II [Eq. (61)], the many-particle Green’s functions are given by

G̃N±1,μν(2)=±G̃N±1,μκ(1)ṼN±1,κν(1)G̃N±1,νν(0)G̃N±1,μν(1)EN,0(1)G̃N±1,νν(0)δμνG̃N±1,μμ(0)EN,0(2)G̃N±1,μμ(0).
(A29)

The second-order one-particle Green’s function is then defined by its elements as

Gpq(2)=z̃N1,μp(0)*G̃N1,μν(2)z̃N1,νq(0)+z̃N1,μp(0)*G̃N1,μν(1)z̃N1,νq(1)+z̃N1,μp(1)*G̃N1,μν(1)z̃N1,νq(0)+z̃N1,μp(0)*G̃N1,μν(0)z̃N1,νq(2)+z̃N1,μp(2)*G̃N1,μν(0)z̃N1,νq(0)+z̃N1,μp(1)*G̃N1,μν(0)z̃N1,νq(1)+z̃N+1,μq(0)*G̃N+1,μν(2)z̃N+1,νp(0)+z̃N+1,μq(0)*G̃N+1,μν(1)z̃N+1,νp(1)+z̃N+1,μq(1)*G̃N+1,μν(1)z̃N+1,νp(0)+z̃N+1,μq(0)*G̃N+1,μν(0)z̃N+1,νp(1)+z̃N+1,μq(1)*G̃N+1,μν(0)z̃N+1,νp(0)+z̃N+1,μq(1)*G̃N+1,μν(0)z̃N+1,νp(1)D(1)Gpq(1)D(2)Gpq(0).
(A30)

Using Eqs. (A21) and (A22) for G̃(1) and Eqs. (A13)–(A18) for z̃’s, we can evaluate each term of G(2) in the (N − 1)-electron sector as follows:

z̃N1,μp(0)*G̃N1,μν(2)z̃N1,νq(0)=z̃N1,μp(0)*G̃N1,μκ(1)ṼN1,κν(1)G̃N1,νν(0)z̃N1,νq(0)  +z̃N1,μp(0)*G̃N1,μν(1)EN,0(1)G̃N1,νν(0)z̃N1,νq(0)  +δμνz̃N1,μp(0)*G̃N1,μμ(0)EN,0(2)G̃N1,μμ(0)z̃N1,μq(0)=12δkpδlqij||akal||ijΔωkΔωaijΔωl+14δkpδkqij||abab||ij(Δωk)2Δωabijk  12δkpδlqil||abab||ikΔωkΔωabiklΔωl+14δkpδkqij||abab||ij(Δωk)2Δijab.
(A31)

The last expression does not contain EN,0(1) for the same reason why G^N±1(1) [Eqs. (A21) and (A22)] do not, either, after mutual cancellation. Likewise,

z̃N1,μp(0)*G̃N1,μν(1)z̃N1,νq(1)=12δkpδcqij||ak(T2(1))ijacΔωkΔωaij+12δkpδlqil||ab(T2(1))ikabΔωkΔωabikl  14δkpδkqij||ab(T2(1))ijabΔωkΔωabijk=12δkpδcqij||akac||ijΔωkΔωaijΔijac+12δkpδlqil||abab||ikΔωkΔωabiklΔikab  14δkpδkqij||abab||ijΔωkΔωabijkΔijab,
(A32)
z̃N1,μp(0)*G̃N1,μν(0)z̃N1,νq(2)=δkpδcq(T1(2))kcΔωk=12δkpδcqij||ak(T2(1))ijacΔωkΔkc+12δkpδcqic||ab(T2(1))ikabΔωkΔkc=12δkpδcqij||akac||ijΔωkΔkcΔijac  +12δkpδcqic||abab||ikΔωkΔkcΔikab,
(A33)

and

z̃N1,μp(1)*G̃N1,μν(0)z̃N1,νq(1)=12δcpδdq(T2(1)*)ijac(T2(1))ijadΔωaij12δkpδlq(T2(1)*)ilab(T2(1))ikabΔωabikl  +14δkpδkq(T2(1)*)ijab(T2(1))ijabΔωabijk=12δcpδdqij||acad||ijΔωaijΔijacΔijad12δkpδlqil||abab||ikΔωabiklΔilabΔikab  +14δkpδkqij||abab||ijΔωabijk(Δijab)2.
(A34)

The terms in the (N + 1)-electron sector are expanded as

z̃N+1,μq(0)*G̃N+1,μν(2)z̃N+1,νp(0)=z̃N+1,μq(0)*G̃N+1,μκ(1)ṼN+1,κν(1)G̃N+1,νν(0)z̃N+1,νp(0)  z̃N+1,μq(0)*G̃N+1,μν(1)EN,0(1)G̃N+1,νν(0)z̃N+1,νp(0)  δμνz̃N+1,μq(0)*G̃N+1,μμ(0)EN,0(2)G̃N+1,μμ(0)z̃N+1,μp(0)=12δdqδcpid||abab||icΔωdΔωiabΔωc+14δcqδcpij||abab||ij(Δωc)2Δωijabc  12δdqδcpij||acad||ijΔωdΔωijacdΔωc14δcqδcpij||abab||ij(Δωc)2Δijab,
(A35)
z̃N+1,μq(0)*G̃N+1,μν(1)z̃N+1,νp(1)=12δcqδkpic||ab(T2(1))ikabΔωcΔωiab12δdqδcpij||ac(T2(1))ijadΔωdΔωijacd  +14δcqδcpij||ab(T2(1))ijabΔωcΔωijabc=12δcqδkpic||abab||ikΔωcΔωiabΔikab  12δdqδcpij||acad||ijΔωdΔωijacdΔijad+14δcqδcpij||abab||ijΔωcΔωijabcΔijab,
(A36)
z̃N+1,μq(0)*G̃N+1,μν(0)z̃N+1,νp(2)=δcqδkp(T1(2))kcΔωc=12δcqδkpij||ka(T2(1))ijacΔωcΔkc12δcqδkpic||ab(T2(1))ikabΔωcΔkc=12δcqδkpij||akac||ijΔωcΔkcΔijac12δcqδkpic||abab||ikΔωcΔkcΔikab,
(A37)

and

z̃N+1,μq(1)*G̃N+1,μν(0)z̃N+1,νp(1)=12δlqδkp(T2(1)*)ilab(T2(1))ikabΔωiab12δdqδcp(T2(1)*)ijac(T2(1))ijadΔωijacd  +14δcqδcp(T2(1)*)ijab(T2(1))ijabΔωijabc=12δlqδkpil||abab||ikΔωiabΔilabΔikab  12δdqδcpij||acad||ijΔωijacdΔijacΔijad+14δcqδcpij||abab||ijΔωijabc(Δijab)2.
(A38)

Of the last three renormalization terms in Eq. (A30), only the following is nonzero (since D(1) = 0):

D(2)Gpq(0)=14δpq(T2(1)*)ijab(T2(1))ijabΔωp=14δpqij||abab||ij(Δijab)2Δωp.
(A39)

Let us first show that unlinked terms mutually cancel with one another, leaving only linked and thus size-consistent contributions. The reader is reminded that an unlinked term is diagrammatically defined in Table I; algebraically, it is a simple product of two or more factors (ignoring the denominator), at least one of which is an extensive quantity (i.e., energy in this case).

For instance, the second and fourth terms in Eq. (A31) are unlinked because their numerators are a simple product of ij||abab||ij and δkpδkq with no common summation index, with the former factor being a part of an extensive quantity. Note that the fourth term is a renormalization term containing EN,0(2). Similarly, the last term of Eq. (A32) is unlinked (and so is its complex conjugate), whereas Eq. (A33) is linked. The last term of Eq. (A34) and the renormalization term [Eq. (A39)] are also unlinked. We can then isolate these unlinked contributions in the (N − 1)-electron sector and show, using some trivial identity involving the denominators, that they sum to zero,

{Gkk(2)}UL=14ij||abab||ij(Δkω)2Δijkωab+14ij||abab||ij(Δkω)2Δijab14ij||abab||ijΔkωΔijkωabΔijab14ij||abab||ijΔijabΔijkωabΔkω14ij||abab||ijΔijkωab(Δijab)2+14ij||abab||ij(Δijab)2Δkω
(A40)
=0,
(A41)

where “UL” means unlinked contributions only. In the left-hand side of Eq. (A41), the second and last terms (with the opposite sign from the rest if the denominator factors are all brought to the Δijab form) are the renormalization terms.

We can also show, after some algebra, the unlinked terms in the (N + 1)-electron sector to also sum to zero,

{Gcc(2)}UL=14ij||abab||ij(Δωc)2Δωijabc14ij||abab||ij(Δωc)2Δijab+14ij||abab||ijΔωcΔωijabcΔijab+14ij||abab||ijΔijabΔωijabcΔωc+14ij||abab||ijΔωijabc(Δijab)214ij||abab||ij(Δijab)2Δωc
(A42)
=0.
(A43)

Equation (A41) leaves only linked terms for the occupied-occupied block of the second-order one-particle Green’s function,

Gkl(2)=12il||abab||ikΔωkΔωabiklΔikab+12il||abab||ikΔilabΔωabiklΔωl12il||abab||ikΔilabΔωabiklΔikab12il||abab||ikΔωlΔωabiklΔωk+12il||abab||ikΔilabΔωiabΔikab+12ij||akal||ijΔωkΔωaijΔωl=12il||abab||ikΔωkΔωiabΔωl+12al||ijij||akΔωkΔωaijΔωl.
(A44)

The linkedness of each term can be algebraically verified by noting that there is at least one common summation index shared by every pair of numerator factors.

Likewise, the virtual-virtual block consists of linked terms only thanks to Eq. (A43),

Gcd(2)=12ij||acad||ijΔωcΔωijacdΔijac12ij||acad||ijΔijadΔωijacdΔωd12ij||acad||ijΔijadΔωijacdΔijac12ij||acad||ijΔωdΔωijacdΔωc+12ij||acad||ijΔijadΔωaijΔijac+12id||abab||icΔωcΔωiabΔωd=12ad||ijij||acΔωcΔωaijΔωd+12id||abab||icΔωcΔωiabΔωd.
(A45)

Lastly, the occupied-virtual block, which does not contain unlinked contributions to begin with, is given by

Gkc(2)=12ij||akac||ijΔωkΔkcΔijac+12ij||akac||ijΔωkΔωaijΔijac+12ij||akac||ijΔωcΔkcΔijac12ic||abab||ikΔωcΔkcΔikab12ic||abab||ikΔωcΔωiabΔikab+12ic||abab||ikΔωkΔkcΔikab=12ac||ijij||akΔωkΔωaijΔωc+12ic||abab||ikΔωkΔωiabΔωc.
(A46)

Together, we find

Gpq(2)=12iq||abab||ipΔωpΔωiabΔωq+12aq||ijij||apΔωpΔωaijΔωq,
(A47)

the substitution of which into Eq. (23) yields

Σpq(2)={(G(0))1G(2)(G(0))1}pq{(G(0))1Σ(1)(G(0))1Σ(1)(G(0))1}pq=ΔωpGpq(2)Δωq=12iq||abab||ipΔωiab+12aq||ijij||apΔωaij,
(A48)

where Σpq(1)=0 [Eq. (A28)] is used. The final result agrees with Eq. (24).

2. Diagrammatic commentary

Here, we consider the diagrammatic representations of some of the foregoing algebraic derivation steps to illustrate the algebraic machinery ensuring the linkedness and irreducibleness of the self-energy diagrams, as a subtext for the linked- and irreducible-diagram theorems discussed in  Appendixes B and  C.

Let us first consider Eq. (A41), which gives us a glimpse of the operation of the factorization theorem of Frantz and Mills,92 which consolidates unlinked contributions differing only in denominators into the ones with factored denominators that cancel exactly with the renormalization terms. Its diagrammatic expression is depicted in Fig. 15. The diagrams in this figure are ordered identically as the terms in this equation, and the one-to-one correspondence can be verified by translating each diagram into the corresponding algebraic expression by using the interpretation rules in Table II amended with two additional rules in Table VII.

FIG. 15.

The diagrammatic representation of Eq. (A41). The terms and diagrams are ordered in the same way.

FIG. 15.

The diagrammatic representation of Eq. (A41). The terms and diagrams are ordered in the same way.

Close modal
TABLE VII.

Amendments to Table II for the algebraic interpretation of one-particle Green’s-function diagrams in the Hugenholtz style. The rules do not apply to insertion diagrams, which are, however, simply the product of the algebraic interpretations of constituent subdiagrams with specified resolvent lines.

(1Label each dangling line with a hole or particle index 
 Consider all possible time orderings of their termini 
(3A terminus of each of the dangling lines 
 is viewed as a vertex in this context 
(1Label each dangling line with a hole or particle index 
 Consider all possible time orderings of their termini 
(3A terminus of each of the dangling lines 
 is viewed as a vertex in this context 

We can recognize that the diagrammatic equation in Fig. 15 is a concatenation of two applications of the factorization theorem with insertion,31,34,91,92 as depicted in Figs. 16 and 17. The reader is referred to the work of Manne,91 Harris et al.,31 or Shavitt and Bartlett34 for an overview of the factorization theorem and its utility in the time-independent proof of the linked-diagram theorem of MBPT.

FIG. 16.

An application of the factorization theorem with insertion to diagrams A1 and A3 dictates that their sum is equal to a diagram with insertion (A2), whose value is the product of the values of its disconnected parts (A2).

FIG. 16.

An application of the factorization theorem with insertion to diagrams A1 and A3 dictates that their sum is equal to a diagram with insertion (A2), whose value is the product of the values of its disconnected parts (A2).

Close modal
FIG. 17.

Another application of the factorization theorem with insertion to diagrams A4 and A5.

FIG. 17.

Another application of the factorization theorem with insertion to diagrams A4 and A5.

Close modal

Specifically, diagrams A1 and A3 encompass all time orderings (vertical positions) of the bottom V vertex of the closed part relative to the p and q vertexes (the dangling-line termini) of the open part, while maintaining the position of the top V vertex of the closed part immediately below the top (p) vertex (the higher terminus) of the open part. As per the factorization theorem,31,34,91,92 the sum of these two diagrams is equal to the insertion of the closed part into the open part immediately below the top (p) vertex of the latter (diagram A2). The open part has two resolvent lines in between adjacent vertexes, one above and one below the insertion point. The value of insertion diagram A2 is, in turn, the simple product of the values of its open subdiagram (with two resolvent lines in this case) and closed subdiagram (A2). The validity of this assertion can be easily verified algebraically by noting that the terms differ only in denominators and they satisfy

1ΔkωΔijkωabΔkω+1ΔkωΔijkωabΔijab=1ΔkωΔijabΔkω,
(A49)

which is, however, unnecessary as it is guaranteed by the factorization theorem.

Likewise, diagrams A4 and A5 cover all time orderings of the bottom (q) vertex (the lower terminus) of the open part, while keeping its top (p) vertex immediately below the top V vertex of the closed part. Their sum is, therefore, the insertion of the open part into the closed part (A6) with two resolvent lines above and below the insertion point. The result of this insertion is equal to the product of two disconnected subdiagrams (A6).

Next, we consider Eqs. (A44), (A45), and (A46). The diagrammatic versions of Eqs. (A44) and (A46) are shown in Figs. 18 and 19, respectively.

FIG. 18.

The diagrammatic representation of Eq. (A44). The terms and diagrams are ordered in the same way.

FIG. 18.

The diagrammatic representation of Eq. (A44). The terms and diagrams are ordered in the same way.

Close modal
FIG. 19.

The diagrammatic representation of Eq. (A46). The terms and diagrams are ordered in the same way. That diagrams A15, A16, and A17 sum to A21 can be shown by an application of the original factorization theorem for diagrams with two open parts by first deleting the areas enclosed by dashed boxes and then later restoring them.

FIG. 19.

The diagrammatic representation of Eq. (A46). The terms and diagrams are ordered in the same way. That diagrams A15, A16, and A17 sum to A21 can be shown by an application of the original factorization theorem for diagrams with two open parts by first deleting the areas enclosed by dashed boxes and then later restoring them.

Close modal

The six linked terms in the left-hand side of Eq. (A44), differing only in the denominators, are consolidated into two terms with simpler denominators (the right-hand side). The same is true with Eqs. (A45) and (A46). This simplification, or the denominator factorization, can be easily verified algebraically and is essential during the subsequent step of extracting the self-energy from the one-particle Green’s function [Eq. (A48)], where the reciprocals of the zeroth-order Green’s function (Δωp) cancel with the identical factors in the denominator of the second-order Green’s function. It justifies part of diagrammatic rule 1 of “trimming” of dangling lines in Fig. 1. Note that the said cancellation does not occur for individual terms in the left-hand side of Eq. (A44), (A45), or (A46); only their sums, after the factorization, have the appropriate forms of the denominators that are canceled by the zeroth-order Green’s functions.

Diagrammatically, the left-hand side of each of these equations (before the denominator factorization) consists of the Green’s-function diagrams whose resolvent lines cut across many lines including dangling lines. Their sum, however, is two diagrams in which the resolvent lines are split and each dangling line has its own resolvent line intersecting only it but no other line. Again, this occurs systematically owing to the factorization theorem.

In Fig. 19, corresponding to Eq. (A46), the first three diagrams (A15, A16, and A17) in the left-hand side sum up to diagram A21 in the right-hand side, while the sum of diagrams A18, A19, and A20 is equal to diagram A22. That diagrams A15, A16, and A17 sum to A21 is a straightforward application of the original (i.e., without insertion) factorization theorem31,34,91,92 for diagrams with two open parts. First, the areas enclosed by dashed boxes are deleted temporarily so that each diagram becomes unlinked with two open parts (which does not alter the resolvent lines or denominators). The three diagrams cover all time orderings of vertexes (including the dangling-line termini) while maintaining the open parts aligned at the bottom. They, therefore, add up to become a product of two separate open parts with each resolvent line spanning only the part it belongs to, as per the original factorization theorem. At this point, the deleted areas (which are algebraically a common multiplicative factor) can be restored, forming diagram A21. The same logic proves that diagrams A18, A19, and A20 sum to A22.

In Fig. 18, the first five diagrams (A7–A11) sum to diagram A13, whereas diagram A12 already has a convenient form for trimming of dangling lines and is, therefore, unchanged (A14). That the sum of A7 through A11 is equal to A13 can be shown in three steps, as depicted in Figs. 20–22. First (Fig. 20), applying the factorization theorem to diagrams A7, A9, and A11 with the areas in the dashed boxes deleted and later restored, their sum is shown to be equal to diagram A13. Second (Fig. 21), doing the same to diagrams A8 and A10, we obtain diagram A13 as their sum. Third (Fig. 22), the sum of diagrams A13 and A13 is shown to be equal to diagram A13 by again applying the factorization theorem.

FIG. 20.

An application of the factorization theorem to diagrams A7, A9, and A11 with the areas enclosed by dashed boxes temporarily deleted proves that their sum is equal to diagram A13.

FIG. 20.

An application of the factorization theorem to diagrams A7, A9, and A11 with the areas enclosed by dashed boxes temporarily deleted proves that their sum is equal to diagram A13.

Close modal
FIG. 21.

Same as Fig. 20 but to diagrams A8 and A10.

FIG. 21.

Same as Fig. 20 but to diagrams A8 and A10.

Close modal
FIG. 22.

Same as Fig. 20 but to diagrams A13 and A13.

FIG. 22.

Same as Fig. 20 but to diagrams A13 and A13.

Close modal

In summary, the exact cancellation of the unlinked terms such as occurring in Eq. (A41) and the factorization of zeroth-order Green’s function for dangling lines in Eq. (A44), facilitating the trimming of the latter, are not coincidence but an inevitable, systematic consequence of the factorization theorem, which underlies the linked- and irreducible-diagram theorems of the one-particle Green’s functions and self-energies.

A perturbation correction to the one-particle Green’s function is linked at any order. This is considered well known and certainly widely accepted with a proof in a time-dependent picture provided by Gell-Mann and Low.90 Here, we show the same using only time-independent (frequency-domain) arguments, which may be easier to follow by quantum chemists concerned with solving time-independent Schrödinger equations. It is not a rigorous mathematical proof, however, but an outline of one with concrete examples. We heavily borrow notations, concepts, and proofs from the closely related linked-diagram theorem of MBPT in a time-independent picture as described by Frantz and Mills,92 by Manne,91 by Harris et al.,31 and by Shavitt and Bartlett.34 The original proof of the theorem for MBPT is due to Goldstone,94 on the basis of the work by Brueckner.93 The first time-independent version of the proof is due to Hugenholtz.95 The reader is also referred to Kutzelnigg and Mukherjee,25 who discussed time-independent Rayleigh–Schrödinger and Brillouin–Wigner perturbation expansions of one- and two-particle Green’s functions.

It is sufficient to show the linkedness of either the (N + 1)- or (N − 1)-electron sector of the nth-order correction to the one-particle Green’s function. Therefore, by dividing the latter as

G(n)=GN+1(n)+GN1(n),
(B1)

we will focus only on the (N + 1)-electron sector. Also, by introducing a shifted perturbation operator,

W^(1)=V^(1)EN,0(1),
(B2)

we can rewrite recursion II for the many-particle Green’s function [Eq. (61)] as

G̃N+1,μν(n)=G̃N+1,μκ(n1)W̃N+1,κν(1)G̃N+1,νν(0)i=2nG̃N+1,μν(ni)EN,0(i)G̃N+1,νν(0)
(B3)

with

W̃N+1,μν(1)=ṼN+1,μν(1)δμνEN,0(1).
(B4)

This trick merely decreases the number of mutually canceling unlinked terms involving EN,0(1) and reduces the nonessential clutter in formalisms; see p. 42 of Shavitt and Bartlett.34 Einstein’s convention is used.

Whereas this is not absolutely necessary, we can also exploit the proven fact that the MBPT(n) wave function is linked at any n and is written as

|n|ΨN,0(n)={(R^V^(1))n|ΦN,0(0)}=|{Vn},
(B5)

where R^ is the resolvent operator already defined by Eq. (A3) and a pair of braces “{}” indicates the linkedness. Then, the MBPT energies, whose linkedness is also completely established, are written in Brueckner’s bracket notation34,91,93 as

EN,0(1)=0|V^(1)|0=V,
(B6)
EN,0(2)=0|V^(1)R^V^(1)|0=VV,
(B7)
EN,0(3)=0|V^(1)R^W^(1)R^V^(1)|0=VVV,
(B8)
EN,0(4)=0|V^(1)R^W^(1)R^W^(1)R^V^(1)|0  0|V^(1)R^V^(1)|00|V^(1)R^R^V^(1)|0=VVVVVVVV={VVVV},
(B9)

etc., where braces again mean that vertexes enclosed by them are linked; they are analogous to (but not the same as) “superbrackets” of Paldus and Čížek.11 In the most simplified notation (the rightmost term in each equation), the resolvent operators and the distinction between W^(1) and V^(1) are suppressed. Using this simplified notation of MBPT wave functions, we can rewrite recursion II [Eq. (68)] for the one-particle Green’s function as

{GN+1(n)}pq=i,j,k0i+j+k=ni|q^G^N+1(j)p^|ki,k1,j0i+j+k=ni|k{GN+1(j)}pq,
(B10)

where G^N+1(n) is the operator form of G̃N+1,μν(n), an example of which is given in Eq. (A22) for n = 1, and the summation range in the second term is narrower with the aid of n|0=δn0.

As in MBPT, the objective here is to show that every unlinked contribution in MBGF is exactly canceled by the renormalization terms, which in turn correspond to Brueckner bracket insertions. Unlike in MBPT, unlinked diagrams in MBGF are of one of the following two types: (i) a closed disconnected part is an energy diagram or (ii) it is an overlap diagram [a term in D(n) of Eq. (68) or i|k in Eq. (B10)]. Ultimately, we want to show that every unlinked diagram with a disconnected energy part [type (i)] is canceled by the renormalization terms containing an MBPT energy factor, i.e., the second term of Eq. (B3). We also want to show that an unlinked diagram with a disconnected overlap part [type (ii)] is canceled by the renormalization terms containing a one-particle Green’s function, i.e., the second term of Eq. (B10).

The zeroth-order one-particle Green’s function is written as

{GN+1(0)}pq=0|q^G^N+1(0)p^|0=qp,
(B11)

which is trivially linked. The rightmost is the compressed expression in which the Green’s function in between the adjacent operators is suppressed to avoid clutter.

The first-order one-particle Green’s function which is written as

{GN+1(1)}pq=0|q^G^N+1(0)W^(1)G^N+1(0)p^|0=qVp
(B12)

is also linked by virtue of the fact that W^(1) has been defined as V^(1) minus its internal contraction and is by itself a linked operator (notwithstanding the simplified bracket notation in the rightmost that suppresses the distinction between the two).

In  Appendix A, we have painstakingly derived the second-order Green’s function expression and explicitly showed its linkedness. Here, we illustrate the systematic nature of the cancellation of unlinked terms. In the bracket notation, we have

{GN+1(2)}pq=0|q^G^N+1(2)p^|0+0|q^G^N+1(1)p^|1  +1|q^G^N+1(1)p^|0+0|q^G^N+1(0)p^|2  +2|q^G^N+1(0)p^|0+1|q^G^N+1(0)p^|1  1|1{GN+1(0)}pq
(B13)
=qVVpqVVp+qVpV+VqVp+qp{VV}+{VV}qp+VqpVVqpV
(B14)
=all linked.
(B15)

In the last equality, we have used

qp{VV}UL={VV}qpUL=0,
(B16)
qVVpUL+qVpVUL=qVVp,
(B17)
VqVpUL+VqpVUL=VqpV,
(B18)

where the “UL” subscript extracts unlinked contributions. Equations (B17) and (B18) correspond exactly to the diagrammatic equations in Figs. 16 and 17, respectively. The left-hand side of Eq. (B17) covers all possible time orderings of the lower (i.e., right) V vertex (relative to the p vertex or the lower terminus of the dangling line), which is to be contracted with the upper (i.e., left) V vertex with its position kept immediately below the q vertex (the upper terminus). As per the factorization theorem, their sum is then equal to the insertion of an MBPT(2) energy subdiagram, VV, at the upper V position (the right-hand side). This corresponds to the deletion of the energy-unlinked diagram by the renormalization term in Eq. (B3). Equation (B18) in turn shows that the sum of the overlap-unlinked terms with all possible time orderings of the p vertex while keeping the position of the q vertex immediately below the upper V vertex is equal to the insertion of a zeroth-order Green’s function subdiagram, qp, at the same position. It corresponds to the deletion of the overlap-unlinked diagram by the renormalization term.

Moving on to the third order, we have

{GN+1(3)}pq=0|q^G^N+1(3)p^|0+0|q^G^N+1(2)p^|1+1|q^G^N+1(2)p^|0+0|q^G^N+1(1)p^|2+2|q^G^N+1(1)p^|0+1|q^G^N+1(1)p^|1+0|q^G^N+1(0)p^|3+3|q^G^N+1(0)p^|0+1|q^G^N+1(0)p^|2+2|q^G^N+1(0)p^|11|1{GN+1(1)}pq1|2{GN+1(0)}pq2|1{GN+1(0)}pq
(B19)
=qVVVpqVVVpqVVVpqVVVp+qVVpVqVVpV+VqVVpVqVVp+qVp{VV}+{VV}qVp+VqVpV+qp{VVV}+{VVV}qp+Vqp{VV}+{VV}qpVVqVpVVqp{VV}{VV}qpV
(B20)
=qVVVp+qVVpV+qVp{VV}qVVVp+VqVVp+VqVpV+Vqp{VV}Vqp{VV}+{VV}qVp+{VV}qpV{VV}qpV+qp{VVV}L+{VVV}qpL
(B21)
=all linked,
(B22)

where “L” means linked. The penultimate equality is based on the following identities: Since the first-order Green’s function is zero, we have

VqVpV=0.
(B23)

For the same reason, we deduce

qVVVp=qVVVp=0.
(B24)

Furthermore, n|0=δn0 implies

qVVpV=VqVVp=0,
(B25)
qp{VVV}UL={VVV}qpUL=0.
(B26)

The remaining terms are organized into four groups, each occupying one line in the right-hand side of Eq. (B21). The terms in a group share several upper vertexes in an identical order with the rest of the vertexes permuted in all possible time orders.

The first group has three terms that all start with “qV” with the rest of vertexes scrambled in three distinct orders, which is subtracted by a matching bracket insertion at the V position of the “qV” motif. The sum of the unlinked contributions of the first three terms is, as per the factorization theorem, equal to the bracket insertion,

qVVVpUL+qVVpVUL+qVp{VV}UL=qVVVp.
(B27)

The equation is diagrammatically depicted in Fig. 23, indicating that the three unlinked contributions (B1, B2, and B3) consist of two disconnected parts with identical topology, with the top vertex of the right disconnected part located immediately below the top vertex of the left disconnected part, while the rest of the vertexes of the right part are time ordered in all possible ways. As per the factorization theorem, their sum is the insertion diagram (B4), whose value is equal to that of the bracket insertion, which is, in turn, the product of the values of the disconnected parts with their resolvent lines spanning only one part (B4). Since the algebraic interpretations of the unlinked diagrams differ only in the denominators, this identity can be easily confirmed algebraically also by noting

1ΔmωΔmijωabΔmklωabΔmω+1ΔmωΔmijωabΔmklωabΔklab+1ΔmωΔmijωabΔijabΔklab=1ΔijabΔklab(Δmω)2,
(B28)

for one possible pattern of index labeling. This is unnecessary as it is guaranteed by the factorization theorem.

FIG. 23.

The diagrammatic representation of Eq. (B27). The first, second, and third terms in the left-hand side are diagrams B1, B2, and B3, respectively, the sum of which is insertion diagram B4, which is equal to the renormalization term (diagram B4), representing the right-hand side. Hole/particle distinctions are suppressed.

FIG. 23.

The diagrammatic representation of Eq. (B27). The first, second, and third terms in the left-hand side are diagrams B1, B2, and B3, respectively, the sum of which is insertion diagram B4, which is equal to the renormalization term (diagram B4), representing the right-hand side. Hole/particle distinctions are suppressed.

Close modal

Likewise, the following three terms [the second line of Eq. (B21)] share the “Vq” motif with the rest of vertexes permuted in all possible time orders and form another group, the unlinked contributions of which sum up to the bracket insertion in the same line, again according to the factorization theorem

VqVVpUL+VqVpVUL+Vqp{VV}UL=Vqp{VV},
(B29)

which is depicted diagrammatically in Fig. 24. The terms in the third group in the third line of Eq. (B21) (having the “{VV}q” motif) satisfy the relation

{VV}qVpUL+{VV}qpVUL={VV}qpV,
(B30)

whose diagrammatic version is drawn in Fig. 25. The fourth line of Eq. (B21) is the immediate result of Eq. (B26). Therefore, each line of Eq. (B21) has linked contributions only, proving the size-consistency of the third-order one-particle Green’s function.