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 $\Delta $MP*n* 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 $\Delta $MP*n*, which agrees with MBGF in the diagonal and frequency-independent approximations at $1\u2264n\u22643$ but converges at the full-configuration-interaction (FCI) limit at $n=\u221e$ (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 $\Delta $MP*n* 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 $\Delta $MP*n* and MBGF in the diagonal and frequency-independent approximations at $1\u2264n\u22643$ is algebraically proven, also ascribing the differences at *n* = 4 to the so-called semi-reducible and linked-disconnected diagrams.

## I. INTRODUCTION

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 mapping^{41} 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 dependence^{46–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 three^{15,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 expansion^{9,17,18} of the self-energy. Whereas the Gell-Mann–Low theorem^{87,89,90} offers a time-dependent perturbation-theoretical justification of this diagrammatic derivation, it is an oblique^{91} 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 out^{55} that there were six additional third-order self-energy diagrams, the so-called $\Sigma (\u221e)$ 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 formalism^{53} 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 implemented^{97} what we call the $\Delta $MP*n* method, originated by Pickup and Goscinski^{7} and extended by Chong *et al.*,^{98–100} in a determinant-based, general-order algorithm. In this method, the *n*th-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\xb11)$-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 confirmed^{97} that $\Delta $MP*n* reduces to MBGF(*n*) with the self-energy in the diagonal, frequency-independent approximation at $1\u2009\u2264\u2009n\u2009\u2264\u20093$ but converges at the nonapproximated MBGF(*n*) at $n=\u221e$, 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 definition^{88} of a perturbation theory, which equates the *n*th-order perturbation correction of a given quantity to the *n*th derivative with respect to $\lambda $ of the same quantity computed exactly with a perturbation-scaled Hamiltonian, $\u0124=H^0\u2009+\u2009\lambda 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 $\lambda $-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-energy^{9} 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 theorems^{11,31,34,91–95} of MBPT in the time-independent form of Manne^{91} hold for the one-particle Green’s function, justifying its linked-diagrammatic perturbation expansion and proving its size-consistency.^{101,102} The factorization theorem^{31,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 $\lambda $-variation method, using only slightly modified FCI code and finite-difference approximations for the $\lambda $-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 $\lambda $-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 $\Delta $MP*n* method^{97} 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 $\Delta $MP*n* and MBGF at several low perturbation orders.

## II. THEORY

### A. Definitions

The exact one-particle Green’s function^{9,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 ($\omega $), namely,

where $\Psi N,\mu $ and $EN,\mu $ are the exact (FCI) *N*-electron (unnormalized) wave function and energy for the $\mu $th state ($\mu =0$ being the ground state) and $p^\u2020$ ($q^$) is a creation (annihilation) operator of an electron in the *p*th (*q*th) 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

where $\mu $ runs over all exact $(N\xb11)$-electron states. This indicates that the exact Green’s function diverges whenever $\omega $ agrees with an exact electron-detachment ($EN,0\u2212EN\u22121,\mu $) or electron-attachment ($EN+1,\mu \u2212EN,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\xb11)$-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 partitioning^{104} of the exact Hamiltonian $\u0124$,

where the zeroth-order Hamiltonian $H^(0)$ is the Fock operator. At this stage, $\lambda $, 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 $\lambda $-variation method). The zeroth-order or HF Green’s function is defined as

where $\Delta \omega p=\omega \u2212\mathit{\epsilon}p$ and $\mathit{\epsilon}p$ is the energy of the *p*th canonical HF molecular spinorbital (MO). The zeroth-order wave function, $\Phi 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 “$\Psi $” for a (potentially) multi-Slater-determinant wave function and the letter “$\Phi $” 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 $\omega $ 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, $\Sigma (\omega )$, which together satisfy the so-called Dyson equation

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 $\Sigma (\omega )$ in terms of $G(\omega )$ and $G(0)(\omega )$. In other words, $\Sigma (\omega )$ does not contain any new information which is not already encoded in $G(\omega )$, and we need to define $G(\omega )$ or its perturbative approximation before we can know $\Sigma (\omega )$ (not the other way around).

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

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

Furthermore, since

where $1$ is an *m*-by-*m* unit matrix and $\mathit{\epsilon}$ is the diagonal Fock matrix ($\mathit{\epsilon}pq=\delta pq\mathit{\epsilon}p$) in the canonical HF MO basis, Eq. (8) is equivalent to an *m*-by-*m* matrix eigenvalue equation of the form

or

where $U$ is a unitary matrix and $\omega $ 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, $\Sigma (\omega )$, 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 $\Sigma (\omega )$, 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 $\Sigma (\omega )$ 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 $\omega $ appears both in the left- and right-hand sides, which is solved typically by repeated diagonalizations in an iterative root search for $\omega $. 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,

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 $\omega =\mathit{\epsilon}p$ in the left-hand side,

which can be solved by one diagonalization for each state for which $\omega \u2248\mathit{\epsilon}p$ (valid only for Koopmans states). In the diagonal, frequency-independent approximation, $\omega $ is obtained by one-shot evaluation (without diagonalization or root search) of the left-hand side of

This is also appropriate only for Koopmans states.

However, the most important approximation is the perturbation expansion of $\Sigma (\omega )$. Since $\Sigma (\omega )$ is defined by $G(\omega )$ through the Dyson equation, the perturbation series of $\Sigma (\omega )$ is derived from the same of $G(\omega )$. 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

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

in which the term proportional to $\lambda 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 $\lambda $, we obtain

etc., where the frequency argument is omitted for brevity. Generally, for $n\u22651$, we have

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,

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}

### B. Diagrams

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 *n*th-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.

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^\u2020$ | |

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 $\omega $ | |

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^\u2020$ | |

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 $\omega $ | |

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 *n*th-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 theorem^{11,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))\u22121$ 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 $\Sigma (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.

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.

If we used only rules 1 and 2, we could be overlooking a whole class of diagrams. Rule 3 (“vertex insertion”) dictates that additional *n*th-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*) diagrams^{34} and cut open a bubble implicit in each of the Fock operator vertexes.^{105}

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,

with

(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 $\u27e8left out, right out||left in, right in\u27e9$ |

(4) | Between each adjacent pair of vertexes, draw a resolvent line |

(5) | Associate each resolvent line with $1\u2215\Delta ijk\u2026abc\u2026$, if the number of |

intersections is even; $1\u2215\Delta ijk\u2026\omega ab\u2026$ if the number of hole | |

intersections exceeds the number of particle intersections; $1\u2215\Delta \omega ij\u2026abc\u2026$ | |

otherwise. Here, $i,j,k,\u2026$ ($a,b,c,\u2026$) 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 $\u27e8left out, right out||left in, right in\u27e9$ |

(4) | Between each adjacent pair of vertexes, draw a resolvent line |

(5) | Associate each resolvent line with $1\u2215\Delta ijk\u2026abc\u2026$, if the number of |

intersections is even; $1\u2215\Delta ijk\u2026\omega ab\u2026$ if the number of hole | |

intersections exceeds the number of particle intersections; $1\u2215\Delta \omega ij\u2026abc\u2026$ | |

otherwise. Here, $i,j,k,\u2026$ ($a,b,c,\u2026$) 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 |

### C. $\Delta $MP*n* for Koopmans and non-Koopmans states

Two of the present authors with two coauthors generalized^{97} the $\Delta $MP*n* method of Pickup and Goscinski^{7} and of Chong *et al.*^{98–100} to higher orders. It defines the *n*th-order perturbation correction to the $\mu $th electron binding energy, denoted by $\Sigma \xaf\mu (n)$, as the MBPT(*n*) correction difference between the *N*- and $(N\xb11)$-electron states using the same *N*-electron HF reference wave function,

for $n\u22651$, where $\mu $ labels an $(N\xb11)$-electron state with the MBPT(*n*) correction of $EN\xb11,\mu (n)$. It is well known^{7} that this difference reduces (at the formalism level) to the *n*th-order self-energy in the diagonal, frequency-independent approximation at *n* = 2, the fact used for pedagogical purposes by Szabo and Ostlund^{30} 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\u2009\xb1\u20091)$-electron systems,

Here, the *n*th-order correction to the wave function of the $\mu $th *M*-electron state, $\Psi M,\mu (n)$, is given by the recursion equation,

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

where $\Phi M,\mu (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 $\Delta $MP*n* at any arbitrary order (*n*) for Koopmans states only. Using this method, we numerically confirmed that the $\Delta $MP*n* correction ($\Sigma \xaf\mu (n)$) agrees identically with the *n*th-order self-energy in the diagonal, frequency-independent approximation for $1\u2009\u2264\u2009n\u22643$,

where *p* labels the spinorbital that is vacated or filled in the $\mu $th Koopmans state. As $n\u2192\u221e$, $\Delta $MP*n* approaches MBGF(*n*) with the nondiagonal, frequency-dependent self-energy. We speculated^{97} that this switch occurs because two unexpected classes of diagrams (the semi-reducible diagrams and linked-disconnected diagrams) are included in $\Delta $MP*n*, which, respectively, correct the diagonal and frequency-independent approximations systematically. Examples of such diagrams are given in Fig. 5.

One may intuitively understand that diagram 13 recuperates the effect of off-diagonal $\Sigma pq(2)(\mathit{\epsilon}p)$ at the fourth order, which is already included in MBGF(2) but not in $\Delta $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 $\omega $ because

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 $\Delta $MP*n* and approximate MBGF(*n*) at $1\u2009\u2264\u2009n\u2009\u2264\u20093$ 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 $\Delta $MP*n* 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\xb11)$-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,

where $\nu $ runs over all states having the same zeroth-order energy with the $\mu $th state and $UM,\mu \nu (0)$ is a unitary matrix element. These zeroth-order wave functions are determined by the recursion plus some additional “consistency” conditions^{103} on wave functions, which satisfy

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 $\Delta $MP*n* 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 $\lambda $,

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}

### D. $\lambda $-variation

The $\lambda $-variation method^{88} can evaluate virtually any perturbation series *numerically* up to high orders. It computes the *n*th-order perturbation correction of any given target quantity as the *n*th $\lambda $-derivative of the same quantity calculated by the FCI method with a scaled Hamiltonian, $H^(0)+\lambda V^(1)$.

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

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

### E. Algebraic recursive definition I

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

Using the resolution-of-the-identity with exact $(N\xb11)$-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

Every factor except $p^\u2020$, $q^$, and $\omega $ is dependent on $\lambda $. The $\lambda $-derivatives of the wave functions and energies of the *N*- and $(N\xb11)$-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 $\lambda $-derivative of the second factor in the numerator of each term,

where $GN\xb11,\mu \nu (n)$ (the frequency argument is omitted for brevity) is the *n*th-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\xb11)$-electron states. Hereafter, a “Green’s function” without a qualifier refers to the *one-particle* Green’s function, $G$ or *G*_{pq}.

Since the bra and ket wave functions are the eigenfunctions of $\u0124$, 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, $\u27e8\Psi N,0|\Psi N,0\u27e9$, 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\xb11)$-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}

with

we obtain the recursion for $GN\xb11,\mu \nu (n)$ as

with

which is initiated with the zeroth-order quantities given by

where $EN,0(0)$ or $EN\xb11,\mu (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

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

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

with

and so forth, or more generally,

This can be further simplified with $\u27e8\Psi N,0(n)|\Psi N,0(0)\u27e9=\delta 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 form^{25} 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 technique^{112} to derive these recursions. As Öhrn and Born^{17} implied and Kutzelnigg and Mukherjee^{25} 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 $\lambda $-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 $\Delta $MP*n* and CC theories for electron detachment and attachment is elucidated.

### F. Algebraic recursive definition II

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

Differentiation with respect to $\lambda $ 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 $\lambda $ and their $\lambda $-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

which are shown to obey the following recursion:

with

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

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

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

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

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

for $n\u22651$, 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 $\Sigma (1)$ and $\Sigma (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.

### G. Algebraic recursive definition III

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\xb11)$-electron states,

where $\Psi N\xb11,\mu (0)$ is the zeroth-order HCPT wave function for the $\mu $th state, as given by Eq. (35). For a Koopmans state, it is a single Slater determinant and may be denoted by $\Phi N\xb11,\mu (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 $\Delta $MP*n* expressions, thus clarifying the relationship between $\Delta $MP*n* 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

where the many-particle Green’s function, $G\xaf(n)$ (which differs from either $G(n)$ or $G\u0303(n)$), is defined recursively as

with

where $\Sigma \xaf\nu (n)$ is the $\Delta $MP*n* correction for the $\nu $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, $\Psi N\xb11,\nu (n)$. It is initiated with the zeroth-order quantities

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

for $n\u22651$ (with the last term omitted when *n* = 0) with

## III. GENERAL-ORDER CALCULATIONS

### A. Diagrams

We have developed a symbolic computing program^{124} 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 $\alpha \alpha $ spin block of the self-energy matrix is shown. The $\alpha \beta $ and $\beta \alpha $ blocks are zero, and the $\beta \beta $ block is the same as the $\alpha \alpha $ block.

$\Sigma pq(2)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.010 239 | −0.000 984 | 0.000 000 | −0.000 000 | 0.011 268 |

3 | −0.000 984 | 0.001 304 | 0.000 000 | −0.000 000 | 0.007 050 |

4 | 0.000 000 | 0.000 000 | 0.003 328 | −0.000 000 | 0.000 000 |

5 | −0.000 000 | −0.000 000 | −0.000 000 | 0.003 328 | −0.000 000 |

6 | 0.011 268 | 0.007 050 | 0.000 000 | −0.000 000 | 0.015 802 |

$\Sigma pq(2)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.010 239 | −0.000 984 | 0.000 000 | −0.000 000 | 0.011 268 |

3 | −0.000 984 | 0.001 304 | 0.000 000 | −0.000 000 | 0.007 050 |

4 | 0.000 000 | 0.000 000 | 0.003 328 | −0.000 000 | 0.000 000 |

5 | −0.000 000 | −0.000 000 | −0.000 000 | 0.003 328 | −0.000 000 |

6 | 0.011 268 | 0.007 050 | 0.000 000 | −0.000 000 | 0.015 802 |

$\Sigma pq(3)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.003 632 | −0.004 702 | 0.000 000 | −0.000 000 | 0.008 025 |

3 | −0.004 702 | −0.004 586 | 0.000 000 | −0.000 000 | 0.005 288 |

4 | 0.000 000 | 0.000 000 | 0.012 711 | −0.000 000 | 0.000 000 |

5 | −0.000 000 | −0.000 000 | −0.000 000 | 0.012 711 | −0.000 000 |

6 | 0.008 025 | 0.005 288 | 0.000 000 | −0.000 000 | 0.009 849 |

$\Sigma pq(3)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.003 632 | −0.004 702 | 0.000 000 | −0.000 000 | 0.008 025 |

3 | −0.004 702 | −0.004 586 | 0.000 000 | −0.000 000 | 0.005 288 |

4 | 0.000 000 | 0.000 000 | 0.012 711 | −0.000 000 | 0.000 000 |

5 | −0.000 000 | −0.000 000 | −0.000 000 | 0.012 711 | −0.000 000 |

6 | 0.008 025 | 0.005 288 | 0.000 000 | −0.000 000 | 0.009 849 |

$\Sigma pq(4)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.000 858 | −0.004 281 | 0.000 000 | −0.000 000 | 0.004 964 |

3 | −0.004 281 | −0.004 399 | 0.000 000 | −0.000 000 | 0.003 296 |

4 | 0.000 000 | 0.000 000 | 0.011 417 | −0.000 000 | 0.000 000 |

5 | −0.000 000 | −0.000 000 | −0.000 000 | 0.011 417 | −0.000 000 |

6 | 0.004 964 | 0.003 296 | 0.000 000 | −0.000 000 | 0.005 991 |

$\Sigma pq(4)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.000 858 | −0.004 281 | 0.000 000 | −0.000 000 | 0.004 964 |

3 | −0.004 281 | −0.004 399 | 0.000 000 | −0.000 000 | 0.003 296 |

4 | 0.000 000 | 0.000 000 | 0.011 417 | −0.000 000 | 0.000 000 |

5 | −0.000 000 | −0.000 000 | −0.000 000 | 0.011 417 | −0.000 000 |

6 | 0.004 964 | 0.003 296 | 0.000 000 | −0.000 000 | 0.005 991 |

$\Sigma pq(5)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | 0.000 316 | −0.003 118 | 0.000 000 | −0.000 000 | 0.003 026 |

3 | −0.003 118 | −0.003 129 | 0.000 000 | −0.000 000 | 0.001 957 |

4 | 0.000 000 | 0.000 000 | 0.009 122 | −0.000 000 | 0.000 000 |

5 | −0.000 000 | −0.000 000 | −0.000 000 | 0.009 122 | −0.000 000 |

6 | 0.003 026 | 0.001 957 | 0.000 000 | −0.000 000 | 0.003 596 |

$\Sigma pq(5)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | 0.000 316 | −0.003 118 | 0.000 000 | −0.000 000 | 0.003 026 |

3 | −0.003 118 | −0.003 129 | 0.000 000 | −0.000 000 | 0.001 957 |

4 | 0.000 000 | 0.000 000 | 0.009 122 | −0.000 000 | 0.000 000 |

5 | −0.000 000 | −0.000 000 | −0.000 000 | 0.009 122 | −0.000 000 |

6 | 0.003 026 | 0.001 957 | 0.000 000 | −0.000 000 | 0.003 596 |

### B. $\lambda $-variation

Table IV lists perturbation corrections to the self-energy of the BH molecule, the same system as Table III but obtained with the $\lambda $-variation method. The symmetric seven-point finite-difference formulas^{109,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 $\lambda \u2009=\u20090$ with the grid spacing of $\Delta \lambda =0.01$. Very tight convergence of the FCI roots is essential for more precise results.

$\Sigma pq(2)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.010 239 | 0.000 984 | 0.000 000 | 0.000 000 | −0.011 268 |

3 | 0.000 984 | 0.001 304 | 0.000 000 | 0.000 000 | 0.007 050 |

4 | 0.000 000 | 0.000 000 | 0.003 328 | −0.000 000 | 0.000 000 |

5 | 0.000 000 | 0.000 000 | −0.000 000 | 0.003 328 | −0.000 000 |

6 | −0.011 268 | 0.007 050 | 0.000 000 | −0.000 000 | 0.015 802 |

$\Sigma pq(2)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.010 239 | 0.000 984 | 0.000 000 | 0.000 000 | −0.011 268 |

3 | 0.000 984 | 0.001 304 | 0.000 000 | 0.000 000 | 0.007 050 |

4 | 0.000 000 | 0.000 000 | 0.003 328 | −0.000 000 | 0.000 000 |

5 | 0.000 000 | 0.000 000 | −0.000 000 | 0.003 328 | −0.000 000 |

6 | −0.011 268 | 0.007 050 | 0.000 000 | −0.000 000 | 0.015 802 |

$\Sigma pq(3)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.003 632 | 0.004 702 | −0.000 000 | −0.000 000 | −0.008 025 |

3 | 0.004 702 | −0.004 586 | −0.000 000 | −0.000 000 | 0.005 288 |

4 | −0.000 000 | −0.000 000 | 0.012 711 | 0.000 000 | −0.000 000 |

5 | −0.000 000 | −0.000 000 | 0.000 000 | 0.012 711 | −0.000 000 |

6 | −0.008 025 | 0.005 288 | −0.000 000 | −0.000 000 | 0.009 849 |

$\Sigma pq(3)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.003 632 | 0.004 702 | −0.000 000 | −0.000 000 | −0.008 025 |

3 | 0.004 702 | −0.004 586 | −0.000 000 | −0.000 000 | 0.005 288 |

4 | −0.000 000 | −0.000 000 | 0.012 711 | 0.000 000 | −0.000 000 |

5 | −0.000 000 | −0.000 000 | 0.000 000 | 0.012 711 | −0.000 000 |

6 | −0.008 025 | 0.005 288 | −0.000 000 | −0.000 000 | 0.009 849 |

$\Sigma pq(4)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.000 858 | 0.004 281 | −0.000 000 | −0.000 000 | −0.004 964 |

3 | 0.004 281 | −0.004 398 | 0.000 000 | −0.000 000 | 0.003 296 |

4 | −0.000 000 | 0.000 000 | 0.011 416 | 0.000 000 | 0.000 000 |

5 | −0.000 000 | −0.000 000 | 0.000 000 | 0.011 417 | 0.000 000 |

6 | −0.004 964 | 0.003 296 | 0.000 000 | 0.000 000 | 0.005 994 |

$\Sigma pq(4)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.000 858 | 0.004 281 | −0.000 000 | −0.000 000 | −0.004 964 |

3 | 0.004 281 | −0.004 398 | 0.000 000 | −0.000 000 | 0.003 296 |

4 | −0.000 000 | 0.000 000 | 0.011 416 | 0.000 000 | 0.000 000 |

5 | −0.000 000 | −0.000 000 | 0.000 000 | 0.011 417 | 0.000 000 |

6 | −0.004 964 | 0.003 296 | 0.000 000 | 0.000 000 | 0.005 994 |

$\Sigma pq(5)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | 0.000 312 | 0.003 120 | 0.000 000 | 0.000 000 | −0.003 028 |

3 | 0.003 120 | −0.003 131 | 0.000 000 | 0.000 000 | 0.001 958 |

4 | 0.000 000 | 0.000 000 | 0.009 137 | 0.000 000 | 0.000 000 |

5 | 0.000 000 | 0.000 000 | 0.000 000 | 0.009 126 | 0.000 004 |

6 | −0.003 028 | 0.001 958 | 0.000 000 | 0.000 004 | 0.003 593 |

$\Sigma pq(5)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | 0.000 312 | 0.003 120 | 0.000 000 | 0.000 000 | −0.003 028 |

3 | 0.003 120 | −0.003 131 | 0.000 000 | 0.000 000 | 0.001 958 |

4 | 0.000 000 | 0.000 000 | 0.009 137 | 0.000 000 | 0.000 000 |

5 | 0.000 000 | 0.000 000 | 0.000 000 | 0.009 126 | 0.000 004 |

6 | −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^{−5} *E*_{h} 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).

### C. $\Delta $MP*n*

The $\Delta $MP*n* method^{97} was extended to non-Koopmans states by implementing a general-order HCPT method^{103} using the determinant-based algorithm.^{66} Our implementation was verified by the $\lambda $-variation method, which computed up to the fourth-order HCPT corrections to the energies of non-Koopmans states. They were evaluated as the $\lambda $-derivatives of the FCI energies with the forward seven- or nine-point finite difference formulas with $\lambda =0$ as the lower end (not the center) of an evenly spaced grid with $\Delta \lambda =0.01$. This asymmetic grid choice was to safeguard against the states that are degenerate at $\lambda =0$ but become non-degenerate and can be ordered differently between $\lambda <0$ and $\lambda >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 $\Delta $MP*n* 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.

We defer the comparison of the $\Delta $MP*n* data with those from various approximations of MBGF until Sec. III D.

### D. Algebraic recursive definitions

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\xb11)$-electron system are then determined by the determinant-based, general-order HCPT algorithm^{103} described in Sec. III C; (iii) with these, *z*’s and *D*’s are generated^{68,72} up to sufficiently high orders; (iv) for a given $\omega $, recursion equations (49) and (50) are evaluated^{68,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 ($\Sigma $) 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 $\lambda $-variation calculations within the numerical precision of the latter. This observation underscores the full mathematical equivalence of the three algebraic recursions, $\lambda $-variation, and diagrammatics.

$\Sigma pq(2)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.010 239 | 0.000 984 | 0.000 000 | 0.000 000 | −0.011 268 |

3 | 0.000 984 | 0.001 304 | 0.000 000 | −0.000 000 | 0.007 050 |

4 | 0.000 000 | 0.000 000 | 0.003 328 | 0.000 000 | −0.000 000 |

5 | 0.000 000 | −0.000 000 | −0.000 000 | 0.003 328 | −0.000 000 |

6 | −0.011 268 | 0.007 050 | −0.000 000 | −0.000 000 | 0.015 802 |

$\Sigma pq(2)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.010 239 | 0.000 984 | 0.000 000 | 0.000 000 | −0.011 268 |

3 | 0.000 984 | 0.001 304 | 0.000 000 | −0.000 000 | 0.007 050 |

4 | 0.000 000 | 0.000 000 | 0.003 328 | 0.000 000 | −0.000 000 |

5 | 0.000 000 | −0.000 000 | −0.000 000 | 0.003 328 | −0.000 000 |

6 | −0.011 268 | 0.007 050 | −0.000 000 | −0.000 000 | 0.015 802 |

$\Sigma pq(3)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.003 632 | 0.004 702 | 0.000 000 | 0.000 000 | −0.008 025 |

3 | 0.004 702 | −0.004 586 | 0.000 000 | −0.000 000 | 0.005 288 |

4 | 0.000 000 | 0.000 000 | 0.012 711 | 0.000 000 | −0.000 000 |

5 | 0.000 000 | −0.000 000 | 0.000 000 | 0.012 711 | −0.000 000 |

6 | −0.008 025 | 0.005 288 | −0.000 000 | −0.000 000 | 0.009 849 |

$\Sigma pq(3)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.003 632 | 0.004 702 | 0.000 000 | 0.000 000 | −0.008 025 |

3 | 0.004 702 | −0.004 586 | 0.000 000 | −0.000 000 | 0.005 288 |

4 | 0.000 000 | 0.000 000 | 0.012 711 | 0.000 000 | −0.000 000 |

5 | 0.000 000 | −0.000 000 | 0.000 000 | 0.012 711 | −0.000 000 |

6 | −0.008 025 | 0.005 288 | −0.000 000 | −0.000 000 | 0.009 849 |

$\Sigma pq(4)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.000 858 | 0.004 281 | 0.000 000 | 0.000 000 | −0.004 964 |

3 | 0.004 281 | −0.004 399 | 0.000 000 | −0.000 000 | 0.003 296 |

4 | 0.000 000 | 0.000 000 | 0.011 417 | 0.000 000 | −0.000 000 |

5 | 0.000 000 | −0.000 000 | 0.000 000 | 0.011 417 | −0.000 000 |

6 | −0.004 964 | 0.003 296 | −0.000 000 | −0.000 000 | 0.005 991 |

$\Sigma pq(4)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | −0.000 858 | 0.004 281 | 0.000 000 | 0.000 000 | −0.004 964 |

3 | 0.004 281 | −0.004 399 | 0.000 000 | −0.000 000 | 0.003 296 |

4 | 0.000 000 | 0.000 000 | 0.011 417 | 0.000 000 | −0.000 000 |

5 | 0.000 000 | −0.000 000 | 0.000 000 | 0.011 417 | −0.000 000 |

6 | −0.004 964 | 0.003 296 | −0.000 000 | −0.000 000 | 0.005 991 |

$\Sigma pq(5)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | 0.000 316 | 0.003 118 | 0.000 000 | 0.000 000 | −0.003 026 |

3 | 0.003 118 | −0.003 129 | 0.000 000 | −0.000 000 | 0.001 957 |

4 | 0.000 000 | 0.000 000 | 0.009 122 | 0.000 000 | −0.000 000 |

5 | 0.000 000 | −0.000 000 | 0.000 000 | 0.009 122 | −0.000 000 |

6 | −0.003 026 | 0.001 957 | −0.000 000 | −0.000 000 | 0.003 596 |

$\Sigma pq(5)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | 0.000 316 | 0.003 118 | 0.000 000 | 0.000 000 | −0.003 026 |

3 | 0.003 118 | −0.003 129 | 0.000 000 | −0.000 000 | 0.001 957 |

4 | 0.000 000 | 0.000 000 | 0.009 122 | 0.000 000 | −0.000 000 |

5 | 0.000 000 | −0.000 000 | 0.000 000 | 0.009 122 | −0.000 000 |

6 | −0.003 026 | 0.001 957 | −0.000 000 | −0.000 000 | 0.003 596 |

$\Sigma pq(6)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | 0.000 709 | 0.002 080 | 0.000 000 | 0.000 000 | −0.001 836 |

3 | 0.002 080 | −0.002 031 | 0.000 000 | −0.000 000 | 0.001 124 |

4 | 0.000 000 | 0.000 000 | 0.007 036 | 0.000 000 | −0.000 000 |

5 | 0.000 000 | −0.000 000 | 0.000 000 | 0.007 036 | −0.000 000 |

6 | −0.001 836 | 0.001 124 | −0.000 000 | −0.000 000 | 0.002 117 |

$\Sigma pq(6)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | 0.000 709 | 0.002 080 | 0.000 000 | 0.000 000 | −0.001 836 |

3 | 0.002 080 | −0.002 031 | 0.000 000 | −0.000 000 | 0.001 124 |

4 | 0.000 000 | 0.000 000 | 0.007 036 | 0.000 000 | −0.000 000 |

5 | 0.000 000 | −0.000 000 | 0.000 000 | 0.007 036 | −0.000 000 |

6 | −0.001 836 | 0.001 124 | −0.000 000 | −0.000 000 | 0.002 117 |

$\Sigma pq(7)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | 0.000 725 | 0.001 318 | 0.000 000 | 0.000 000 | −0.001 089 |

3 | 0.001 318 | −0.001 282 | 0.000 000 | −0.000 000 | 0.000 628 |

4 | 0.000 000 | 0.000 000 | 0.005 286 | 0.000 000 | −0.000 000 |

5 | 0.000 000 | −0.000 000 | 0.000 000 | 0.005 286 | −0.000 000 |

6 | −0.001 089 | 0.000 628 | −0.000 000 | −0.000 000 | 0.001 207 |

$\Sigma pq(7)$ . | q = 2
. | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|

p = 2 | 0.000 725 | 0.001 318 | 0.000 000 | 0.000 000 | −0.001 089 |

3 | 0.001 318 | −0.001 282 | 0.000 000 | −0.000 000 | 0.000 628 |

4 | 0.000 000 | 0.000 000 | 0.005 286 | 0.000 000 | −0.000 000 |

5 | 0.000 000 | −0.000 000 | 0.000 000 | 0.005 286 | −0.000 000 |

6 | −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 $\Delta $MP*n* electron binding energy of the same orbital is also included.

n
. | Full . | Diagonal . | $\omega $-indep. . | Diagonal, $\omega $-indep. . | $\Delta $MPn
. |
---|---|---|---|---|---|

0 | −0.246 54 | −0.246 54^{a} | −0.246 54 | −0.246 54^{a} | −0.246 54^{a} |

1 | −0.246 54 | −0.246 54^{a} | −0.246 54 | −0.246 54^{a} | −0.246 54^{a} |

2 | −0.244 11 | −0.244 07^{a} | −0.244 05 | −0.244 00^{a} | −0.244 00^{a} |

3 | −0.247 69 | −0.247 61^{a} | −0.247 74 | −0.247 66^{a} | −0.247 66^{a} |

4 | −0.251 13 | −0.251 13 | −0.251 40 | −0.251 40 | −0.251 52 |

5 | −0.253 45 | −0.253 56 | −0.253 92 | −0.254 04 | −0.254 18 |

6 | −0.254 86 | −0.255 06 | −0.255 47 | −0.255 71 | −0.255 71 |

7 | −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 |

$\u221e$ | −0.257 00^{a} | −0.257 37^{a} | −0.257 92 | −0.258 37^{a} | −0.257 00^{a} |

n
. | Full . | Diagonal . | $\omega $-indep. . | Diagonal, $\omega $-indep. . | $\Delta $MPn
. |
---|---|---|---|---|---|

0 | −0.246 54 | −0.246 54^{a} | −0.246 54 | −0.246 54^{a} | −0.246 54^{a} |

1 | −0.246 54 | −0.246 54^{a} | −0.246 54 | −0.246 54^{a} | −0.246 54^{a} |

2 | −0.244 11 | −0.244 07^{a} | −0.244 05 | −0.244 00^{a} | −0.244 00^{a} |

3 | −0.247 69 | −0.247 61^{a} | −0.247 74 | −0.247 66^{a} | −0.247 66^{a} |

4 | −0.251 13 | −0.251 13 | −0.251 40 | −0.251 40 | −0.251 52 |

5 | −0.253 45 | −0.253 56 | −0.253 92 | −0.254 04 | −0.254 18 |

6 | −0.254 86 | −0.255 06 | −0.255 47 | −0.255 71 | −0.255 71 |

7 | −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 |

$\u221e$ | −0.257 00^{a} | −0.257 37^{a} | −0.257 92 | −0.258 37^{a} | −0.257 00^{a} |

^{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 $\omega =\mathit{\epsilon}p$ because they contain dangling lines (and/or articulation lines) with a divergent resolvent factor of $1\u2215\Delta p\omega $ (see Appendix C). Therefore, the self-energies in the frequency-independent approximation at $\omega =\mathit{\epsilon}p$ were determined by a linear interpolation from those on an evenly spaced grid of $\omega $’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 $n\u22652$. These differences persist in the $n=\u221e$ 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 $\omega $. 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, $\omega $, 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 $\omega $-line are the MBGF(20) electron binding energies in the diagonal approximation. On the other hand, the values of these curves at $\omega =\mathit{\epsilon}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.

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 $\Delta $MP*n* self-energy agrees with the one in the diagonal, frequency-independent approximation at $1\u2264n\u22643$ but converges at the full self-energy in the FCI limit. These data substantiate our speculation^{97} that $\Delta $MP*n* does not agree with MBGF(*n*) or any of its approximation in the range of $4\u2264n<\u221e$ 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 $\Delta $MP*n*. For this tiny system, the diagonal approximation has no effect on the LUMO electron binding energy.

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, $\Delta $MP*n* 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, $\Delta $MP*n* 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 $\omega $ 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 $\omega $ itself. The intersections between them seen in $\omega <\u22121.2\u2009Eh$ are the roots corresponding to non-Koopmans states, whereas those around $\omega \u2248\u22120.25\u2009Eh$ are the Koopmans state (the area enlarged in Fig. 7).

Figure 11 plots the same but for the LUMO. In this figure, non-Koopmans solutions are visible both in $\omega <\u22120.7\u2009Eh$ and $\omega >1.2\u2009Eh$, while the Koopmans roots are concentrated in a small area around $\omega \u22480.27\u2009Eh$. Although the perturbative self-energies beyond the first order are invariably divergent at many $\omega $’s, the intersections (roots) can never be divergent for the simple reason that they must also be on the diagonal $\omega $-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=\u221e$, it cannot be divergent even at $n<\u221e$.

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 $\omega \u2248\u22120.8\u2009Eh$ in Fig. 11 as a function of the perturbation order, together with the $\Delta $MP*n* 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). $\Delta $MP*n*, 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 $\Delta $MP*n* over MBGF(*n*), but it should be recalled that the former is frequently divergent, while the latter (without the frequency-independent approximation) is not.

Figures 13 and 14 plot the electron binding energies of the HOMO−1 and LUMO, respectively, of the H_{2}O 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, $\Delta $MP*n* converges less rapidly than full MBGF, reaffirming that the opposite case found in Fig. 8 is accidental.

## IV. CONCLUDING REMARKS

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 Manne^{91} 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 $\Delta $MP*n* and MBGF(*n*) in the diagonal, frequency-independent approximation at $1\u2264n\u22643$ 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, $\lambda $-variation, and $\Delta $MP*n*, all but the last yielding the identical, nondivergent perturbation series, mutually verifying the formulations and implementations. In fact, the $\lambda $-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 $n\u22653$ 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); $\Delta $MP*n* 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*.

## ACKNOWLEDGMENTS

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.

### APPENDIX A: DERIVATION OF $\Sigma (1)$ AND $\Sigma (2)$

#### 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 theorem^{34} facilitated by a symbolic computing program^{126}).

Let us first consider the expressions of $z\u0303$’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

where

and

where $\Phi 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

where

and

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\u0303$’s contracted with appropriate determinants as follows:

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

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

and

Next, we construct the order-by-order algebraic formula of $G\u0303$ (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 $\Sigma $ (the self-energy) via Eq. (23).

The zeroth-order many-particle Green’s functions, $G\u0303N\xb11,\mu \nu (0)$, are already given by Eq. (65) and are diagonal. This immediately leads to

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

which are further expanded into an operator form as

and

These sums do not terminate until all terms involving up to the *n*h(*n* − 1)p and *n*h(*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

In the last equality, we have used *D*^{(1)} = 0,

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

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

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

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

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

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

and

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

and

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

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 $\delta kp\delta 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,

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 $\Delta ij\u2026ab\u2026$ 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,

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

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),

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

Together, we find

the substitution of which into Eq. (23) yields

#### 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.

($1\u2032$) | Label each dangling line with a hole or particle index |

Consider all possible time orderings of their termini | |

($3\u2032$) | A terminus of each of the dangling lines |

is viewed as a vertex in this context |

($1\u2032$) | Label each dangling line with a hole or particle index |

Consider all possible time orderings of their termini | |

($3\u2032$) | A 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 Bartlett^{34} for an overview of the factorization theorem and its utility in the time-independent proof of the linked-diagram theorem of MBPT.

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\u2032$). 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\u2032$ 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

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\u2032$) 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.

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 ($\Delta \omega 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 theorem^{31,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\u2032$. Second (Fig. 21), doing the same to diagrams A8 and A10, we obtain diagram $A13\u2032\u2032$ as their sum. Third (Fig. 22), the sum of diagrams $A13\u2032$ and $A13\u2032\u2032$ is shown to be equal to diagram A13 by again applying the factorization theorem.

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.

### APPENDIX B: LINKED-DIAGRAM THEOREM

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 *n*th-order correction to the one-particle Green’s function. Therefore, by dividing the latter as

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

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

with

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

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

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

where $G^N+1(n)$ is the operator form of $G\u0303N+1,\mu \nu (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=\delta 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

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

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

In the last equality, we have used

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

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

For the same reason, we deduce

Furthermore, $n|0=\delta n0$ implies

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\u2009$” 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\u2009$” motif. The sum of the unlinked contributions of the first three terms is, as per the factorization theorem, equal to the bracket insertion,

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\u2032$), 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

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

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

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