Finite-temperature many-body perturbation theory for electrons: Algebraic recursive definitions, second-quantized derivation, linked-diagram theorem, general-order algorithms, and grand canonical and canonical ensembles

Hirata, So

doi:10.1063/5.0061384

A comprehensive and detailed account is presented for the finite-temperature many-body perturbation theory for electrons that expands in power series all thermodynamic functions on an equal footing. Algebraic recursions in the style of the Rayleigh–Schrödinger perturbation theory are derived for the grand potential, chemical potential, internal energy, and entropy in the grand canonical ensemble and for the Helmholtz energy, internal energy, and entropy in the canonical ensemble, leading to their sum-over-states analytical formulas at any arbitrary order. For the grand canonical ensemble, these sum-over-states formulas are systematically transformed to sum-over-orbitals reduced analytical formulas by the quantum-field-theoretical techniques of normal-ordered second quantization and Feynman diagrams extended to finite temperature. It is found that the perturbation corrections to energies entering the recursions have to be treated as a nondiagonal matrix, whose off-diagonal elements are generally nonzero within a subspace spanned by degenerate Slater determinants. They give rise to a unique set of linked diagrams—renormalization diagrams—whose resolvent lines are displaced upward, which are distinct from the well-known anomalous diagrams of which one or more resolvent lines are erased. A linked-diagram theorem is introduced that proves the size-consistency of the finite-temperature many-body perturbation theory at any order. General-order algorithms implementing the recursions establish the convergence of the perturbation series toward the finite-temperature full-configuration-interaction limit unless the series diverges. The normal-ordered Hamiltonian at finite temperature sheds light on the relationship between the finite-temperature Hartree–Fock and first-order many-body perturbation theories.

I. INTRODUCTION

Strong electron correlation^1–13 is said to occur in systems with many low-lying excited states, rendering their ground states quasidegenerate. Not only are they difficult to characterize theoretically and therefore a worthy computational challenge,^14,15 but they are also technologically important, serving as a basis of useful materials whose structures and properties can change dramatically upon external stimuli. Electronic excitations, phase transitions, and crossovers underlying these large, abrupt changes may also occur thermally, giving rise to fascinating phenomena such as Mott transitions,^16,17 Peierls distortion,^18–20 and high-T_c superconductivity.^{16,17,20–23}

It can be imagined that in such systems, thermally averaged properties are more accurately computable than are the zero-temperature properties of a minimum-energy state, whose mean-field wave function tends to be unstable.^24–27 It may even be argued that a finite-temperature treatment is more realistic and meaningful for systems whose very attractiveness derives from its large response to perturbations including temperature variations. For these and other reasons, there has been a surge of interest in finite-temperature electron–correlation theories recently.^28–41

Despite being a cornerstone of finite-temperature electron–correlation theories,^42–45 many-body perturbation theory (MBPT) in the grand canonical ensemble^46–52 has not been fully developed, having gaps in its details: The grand potentials and chemical potentials are treated differently⁵³ with no analytical formulas available for the lowest-order corrections to the chemical potential, internal energy, or entropy for a long time.^54,55 The time-dependent, diagrammatic formulation exploiting the isomorphism of the Schrödinger and Bloch equations is elegant at the outset,^46–52 but it becomes quickly inscrutable with the exhaustiveness of diagram enumeration being uncertain. Convergence of the finite-temperature MBPT to the well-established zero-temperature MBPT has been questionable,^48,49,56 compelling Kohn and Luttinger⁴⁸ to conclude at one point that the theory “is in general not correct.”⁴⁸ Convergence toward the exact, i.e., full-configuration-interaction (FCI) limit at finite temperature⁵⁷ is also unverifiable either numerically or analytically because of the lack of a Rayleigh–Schrödinger-type recursion. The relationship between the widely used finite-temperature Hartree–Fock (HF) theory⁵⁸ and finite-temperature MBPT is tenuous, and the physical meaning⁵⁹ of their orbital energies and thus (quasiparticle) energy bands is still unknown.⁶⁰

We recently introduced^54,55 a new finite-temperature MBPT in the grand canonical ensemble, which expands in power series the grand potential, chemical potential, internal energy, and entropy on an equal footing. We obtained sum-over-states and sum-over-orbitals analytical formulas for the first- and second-order perturbation corrections to these thermodynamic functions in a transparent algebraic derivation involving only the combinatorial identities and the energy sum rules of the Hirschfelder–Certain degenerate perturbation theory (HCPT).⁶¹ Sum-over-states formulas in the canonical ensemble were also reported up to the third order.⁶² These formulas were verified by the exact numerical agreement with the benchmark data computed as the λ-derivatives of the finite-temperature FCI (Ref. 57) (where λ is the dimensionless perturbation strength). Using these analytical formulas, we elucidated⁵⁶ the root cause of the nonconvergence at the correct zero-temperature limit,⁴⁸ which has to do with the nonanalytic nature of the Boltzmann factor and thus plagues most any finite-temperature MBPT.

This article is a comprehensive account of this new finite-temperature MBPT that expands all thermodynamic functions in uniform perturbation series. We introduce algebraic recursions in the style of the Rayleigh–Schrödinger perturbation theory,⁶³ generating sum-over-states analytical formulas at any arbitrary order for both grand canonical and canonical ensembles. They are implemented into general-order algorithms that compute the perturbation corrections at any arbitrary order, demonstrating their convergence toward the exact (i.e., finite-temperature FCI) limits.⁵⁷ For the grand canonical ensemble, we introduce systematic methods of converting them into sum-over-orbitals reduced analytical formulas, which will be useful for actual condensed-matter applications. One of the conversion methods is the normal-ordered second quantization,^28,42 which is fully developed in this study, and the other is the Feynman diagrams,^44,47,50,51 the rules of which are stipulated not just for the grand potential but also for the chemical potential and entropy. The second-quantization rules shed light on the relationship between the finite-temperature HF and first-order MBPT.

To derive the correct reduced analytical formulas, it has been found that the perturbation energies entering the recursions have to be treated as a nondiagonal matrix (not as scalars). Off-diagonal elements of this matrix, which are, in general, nonzero within a subspace spanned by degenerate Slater determinants, are shown to give rise to the renormalization terms⁶³ of the Rayleigh–Schrödinger perturbation equations that are not necessarily unlinked. They are responsible for a unique set of linked diagrams whose resolvent lines (i.e., factors of the energy denominators) are shifted but not erased (the existence of “anomalous” diagrams^48,52 with one or more resolvent lines erased is well known). Since such “renormalization” diagrams originate from the degenerate perturbation theory and occur for the first time at the third order, it is unclear whether the time-dependent, diagrammatic derivation^46–51 or density-matrix formulation⁵² correctly takes them into account.

We furthermore prove the linked-diagram theorem of the finite-temperature MBPT, which asserts that all thermodynamic functions in the grand canonical ensemble are diagrammatically linked and therefore size-consistent (size-extensive)⁶⁴ at any perturbation order. The proof is based on the zero-temperature linked-diagram theorem^65–70 (see also Refs. 59, 63, and 71) and the systematic cancellation of unlinked anomalous diagrams. Finally, we document the numerical results of the perturbation corrections to the thermodynamic functions in both ensembles in a wide range of temperatures.

II. ALGEBRAIC RECURSIVE DEFINITIONS

In this section are presented the Rayleigh–Schrödinger-like algebraic recursion relationships of the sum-over-states analytical formulas for the perturbation corrections to all thermodynamic functions for electrons in the grand canonical ensemble. Those for the canonical ensemble are readily inferred from these by restricting the summations to $\bar{N}$ -electron states and setting all chemical potentials to zero. They are relegated to Appendix A. The sum-over-states analytical formulas have limited practical utility except to produce benchmark data. They nevertheless serve as a mathematical basis of second-quantized and diagrammatic derivations of reduced (sum-over-orbitals) analytical formulas expressed in terms of molecular integrals, which are the focus of Secs. III and IV. They also lead to the linked-diagram theorem proven in Sec. V. Translational, rotational, and vibrational degrees of freedom are suppressed.

A. Grand canonical ensemble

The grand partition function for electrons is given by⁵⁷

Ξ = \sum_{I} e^{- β E_{I} + β μ N_{I}},

(1)

where I runs over all states, E_I is the exact (FCI) energy of the Ith state, N_I is the number of electrons of the same state, and $β = {(k_{B} T)}^{- 1}$ (k_B is the Boltzmann constant and T is the temperature). The exact chemical potential, μ, is the one that keeps the average number of electrons at $\bar{N}$ ⁠, the value ensuring the electroneutrality of the system. It is the root of the following equation:

\bar{N} = \frac{1}{β} \frac{\partial}{\partial μ} \ln Ξ .

(2)

The exact thermodynamic functions such as grand potential (Ω), internal energy (U), and entropy (S) are derived from Ξ.

Ω = - \frac{1}{β} \ln Ξ,

(3)

U = - \frac{\partial}{\partial β} \ln Ξ + μ \bar{N},

(4)

S = k_{B} β^{2} \frac{\partial Ω}{\partial β} .

(5)

They bear the following well-known relationship:

Ω = U - μ \bar{N} - T S .

(6)

B. Recursion for Ω⁽ⁿ⁾

The nth-order perturbation correction to any thermodynamic function X (X = Ξ, μ, Ω, U, or S) is given by^53,59

X^{(n)} = \frac{1}{n!} {(\frac{\partial^{n} X (λ)}{\partial λ^{n}}|}_{λ = 0},

(7)

where X(λ) is defined exactly by the finite-temperature FCI (Ref. 57) with a perturbation-scaled Hamiltonian, $\hat{H} = {\hat{H}}_{0} + λ \hat{V}$ ⁠, with λ being the dimensionless perturbation strength.

Using the Taylor expansions of the exponential and logarithm,

e^{a + b} = e^{a} + b e^{a} + \frac{b^{2}}{2!} e^{a} + \frac{b^{3}}{3!} e^{a} + \dots,

(8)

\ln (a + b) = \ln a + \frac{b}{a} - \frac{b^{2}}{2 a^{2}} + \frac{b^{3}}{3 a^{3}} - \dots,

(9)

which are rapidly convergent when a ≫ b, and defining

D_{I}^{(n)} = E_{I}^{(n)} - μ^{(n)} N_{I}

(10)

to minimize clutter, we obtain the recursion for Ξ⁽ⁿ⁾ (n ≥ 1) as

\begin{matrix} \frac{Ξ^{(n)}}{Ξ^{(0)}} & = & (- β) 〈 D_{I}^{(n)} 〉 + \frac{{(- β)}^{2}}{2!} \sum_{i = 1}^{n - 1} 〈 D_{I}^{(i)} D_{I}^{(n - i)} 〉 \\ + \frac{{(- β)}^{3}}{3!} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} 〈 D_{I}^{(i)} D_{I}^{(j)} D_{I}^{(n - i - j)} 〉 \\ + \dots + \frac{{(- β)}^{n}}{n!} 〈 {(D_{I}^{(1)})}^{n} 〉, \end{matrix}

(11)

where the bracket denotes a zeroth-order thermal average,

⟨ X_{I} ⟩ \equiv \frac{\sum_{I} X_{I} e^{- β E_{I}^{(0)} + β μ^{(0)} N_{I}}}{\sum_{I} e^{- β E_{I}^{(0)} + β μ^{(0)} N_{I}}} = \frac{\sum_{I} X_{I} e^{- β D_{I}^{(0)}}}{\sum_{I} e^{- β D_{I}^{(0)}}} .

(12)

State index I is shown explicitly in the bracket to distinguish it from other state indices that can coappear in the summand, as in Eq. (28). A double bracket denotes a thermal average of a Dirac bracket; see, e.g., Eq. (46). In Eq. (119), a Brueckner bracket is represented by the same symbol, but this is appropriate since it too stands for a thermal average.

Postponing until Sec. II C the important discussion of $E_{I}^{(n)}$ ⁠, we can use Eq. (9) to arrive at the following recursion for Ω⁽ⁿ⁾, which, in turn, depends on the above recursion for Ξ⁽ⁿ⁾:

\begin{matrix} - β Ω^{(n)} & = & \frac{Ξ^{(n)}}{Ξ^{(0)}} - \frac{1}{2} \sum_{i = 1}^{n - 1} \frac{Ξ^{(i)}}{Ξ^{(0)}} \frac{Ξ^{(n - i)}}{Ξ^{(0)}} \\ + \frac{1}{3} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} \frac{Ξ^{(i)}}{Ξ^{(0)}} \frac{Ξ^{(j)}}{Ξ^{(0)}} \frac{Ξ^{(n - i - j)}}{Ξ^{(0)}} \\ - \dots - \frac{{(- 1)}^{n}}{n} {(\frac{Ξ^{(1)}}{Ξ^{(0)}})}^{n} \end{matrix}

(13)

for n ≥ 1. Although nontrivial, we can rewrite the above equation as

\begin{matrix} - β Ω^{(n)} & = & \frac{Ξ^{(n)}}{Ξ^{(0)}} - \frac{{(- β)}^{2}}{2!} \sum_{i = 1}^{n - 1} Ω^{(i)} Ω^{(n - i)} \\ - \frac{{(- β)}^{3}}{3!} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} Ω^{(i)} Ω^{(j)} Ω^{(n - i - j)} \\ - \dots - \frac{{(- β)}^{n}}{n!} {(Ω^{(1)})}^{n} . \end{matrix}

(14)

A proof of the equivalence of Eqs. (13) and (14) is given in Appendix B. This is further simplified to a form that does not depend on Ξ⁽ⁿ⁾,

\begin{matrix} Ω^{(n)} & = & 〈 D_{I}^{(n)} 〉 + \frac{(- β)}{2!} \sum_{i = 1}^{n - 1} (〈 D_{I}^{(i)} D_{I}^{(n - i)} 〉 - Ω^{(i)} Ω^{(n - i)}) \\ + \frac{{(- β)}^{2}}{3!} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} (〈 D_{I}^{(i)} D_{I}^{(j)} D_{I}^{(n - i - j)} 〉 - Ω^{(i)} Ω^{(j)} Ω^{(n - i - j)}) \\ + \dots + \frac{{(- β)}^{n - 1}}{n!} \{〈 {(D_{I}^{(1)})}^{n} 〉 - {(Ω^{(1)})}^{n}\} \end{matrix}

(15)

for n ≥ 1.

Next, we describe how $⟨ E_{I}^{(n)} ⟩$ ⁠, $⟨ E_{I}^{(i)} E_{I}^{(n - i)} ⟩$ ⁠, etc., can be evaluated so that the above recursion for Ω⁽ⁿ⁾ can be executed. See Sec. II D for the recursion for μ⁽ⁿ⁾, which is also needed. The above recursion is tantalizing with the linked-diagram theorem of Ω⁽ⁿ⁾, which will be discussed in Sec. V.

C. Recursion for $⟨ E_{I}^{(n)} ⟩$

The perturbation corrections to energy, $E_{I}^{(i)}$ (1 ≤ i ≤ n), and their thermal averages, $⟨ E_{I}^{(i)} ⟩$ ⁠, are crucial quantities entering the recursion for Ω⁽ⁿ⁾. They are identified as those defined by the Hirschfelder–Certain degenerate perturbation theory (HCPT).⁶¹ The more familiar Møller–Plesset perturbation theory (MPPT)⁷² will not suffice here because a finite-temperature MBPT accesses all states, many of whose zeroth-order energies are exactly degenerate. Other degenerate or quasidegenerate perturbation theories may also fall short because HCPT is the proper Rayleigh–Schrödinger perturbation theory for degenerate and nondegenerate reference states, complying with the canonical definition of perturbation theory, i.e., Eq. (7), although there are equivalent degenerate perturbation theories under different names.⁷³

For a nondegenerate reference state,⁶³ HCPT reduces to MPPT. The Schrödinger equation may be written with the perturbatively expanded energy and wave function as

\begin{align} ({\hat{H}}_{0} + λ \hat{V}) (Φ_{I}^{(0)} + λ Φ_{I}^{(1)} + λ^{2} Φ_{I}^{(2)} + \dots) \\ = (E_{I}^{(0)} + λ E_{I}^{(1)} + λ^{2} E_{I}^{(2)} + \dots) (Φ_{I}^{(0)} + λ Φ_{I}^{(1)} + λ^{2} Φ_{I}^{(2)} + \dots) \end{align}

(16)

with intermediate normalization $⟨ Φ_{I}^{(0)} | Φ_{I}^{(n)} ⟩ = δ_{n 0}$ ⁠, where $Φ_{I}^{(0)}$ is the nondegenerate reference (zeroth-order) wave function of the Ith state, which is a single Slater determinant. Expanding and equating terms carrying the nth power of λ, we obtain the recursion for the nth-order correction to the wave function, $Φ_{I}^{(n)}$ ⁠, as

(E_{I}^{(0)} - {\hat{H}}_{0}) Φ_{I}^{(n)} = \hat{V} Φ_{I}^{(n - 1)} - \sum_{i = 1}^{n} E_{I}^{(i)} Φ_{I}^{(n - i)}

(17)

or

Φ_{I}^{(n)} = {\hat{R}}_{I} (\hat{V} Φ_{I}^{(n - 1)} - \sum_{i = 1}^{n - 1} E_{I}^{(i)} Φ_{I}^{(n - i)}),

(18)

where ${\hat{R}}_{I}$ is the resolvent operator,

{\hat{R}}_{I} = \sum_{A}^{denom. \neq 0} \frac{| Φ_{A}^{(0)} 〉 〈 Φ_{A}^{(0)} |}{E_{I}^{(0)} - E_{A}^{(0)}} .

(19)

Here, “denom. ≠ 0” limits the summation to just over the states $Φ_{A}^{(0)}$ that are not degenerate with the Ith state and therefore to the cases where $E_{I}^{(0)} - E_{A}^{(0)} \neq 0$ ⁠. The second term in the right-hand side of Eq. (17) is known as the renormalization term,⁶³ which is entirely diagrammatically unlinked (non-size-consistent) in MPPT. It vanishes exactly by canceling out the unlinked contribution of the same magnitude in the first term, leaving only the linked (size-consistent) contribution as the perturbation correction to the wave function. Multiplying Eq. (17) by $Φ_{I}^{(0) *}$ from left and integrating, we have

E_{I}^{(n)} = ⟨ Φ_{I}^{(0)} | \hat{V} | Φ_{I}^{(n - 1)} ⟩,

(20)

defining the nth-order energy correction according to MPPT.

For M-tuply degenerate reference states (M ≥ 2),⁶¹ HCPT and MPPT differ from each other materially. The perturbation expansion of the Schrödinger equation now becomes an M-by-M matrix equation of the form

\begin{align} ({\hat{H}}_{0} + λ \hat{V}) \{(\begin{matrix} Φ_{1}^{(0)} \\ ⋮ \\ Φ_{M}^{(0)} \end{matrix}) + λ (\begin{matrix} Φ_{1}^{(1)} \\ ⋮ \\ Φ_{M}^{(1)} \end{matrix}) + λ^{2} (\begin{matrix} Φ_{1}^{(2)} \\ ⋮ \\ Φ_{M}^{(2)} \end{matrix}) + \dots\} \\ = \{(\begin{matrix} E_{11}^{(0)} & \dots & E_{M 1}^{(0)} \\ ⋮ & ⋮ \\ E_{1 M}^{(0)} & \dots & E_{M M}^{(0)} \end{matrix}) + λ (\begin{matrix} E_{11}^{(1)} & \dots & E_{M 1}^{(1)} \\ ⋮ & ⋮ \\ E_{1 M}^{(1)} & \dots & E_{M M}^{(1)} \end{matrix})) \\ + (λ^{2} (\begin{matrix} E_{11}^{(2)} & \dots & E_{M 1}^{(2)} \\ ⋮ & ⋮ \\ E_{1 M}^{(2)} & \dots & E_{M M}^{(2)} \end{matrix}) + \dots\} \\ \times \{(\begin{matrix} Φ_{1}^{(0)} \\ ⋮ \\ Φ_{M}^{(0)} \end{matrix}) + λ (\begin{matrix} Φ_{1}^{(1)} \\ ⋮ \\ Φ_{M}^{(1)} \end{matrix}) + λ^{2} (\begin{matrix} Φ_{1}^{(2)} \\ ⋮ \\ Φ_{M}^{(2)} \end{matrix}) + \dots\}, \end{align}

(21)

where all M zeroth-order determinants, ${Φ_{I}^{(0)}}$ (1 ≤ I ≤ M), are degenerate and therefore $E_{I J}^{(0)} = E_{I}^{(0)} δ_{I J}$ ⁠. Collecting terms that are proportional to λⁿ, we obtain

\begin{align} (\begin{matrix} (E_{I}^{(0)} - {\hat{H}}_{0}) Φ_{1}^{(n)} \\ ⋮ \\ (E_{I}^{(0)} - {\hat{H}}_{0}) Φ_{M}^{(n)} \end{matrix}) & = (\begin{matrix} \hat{V} Φ_{1}^{(n - 1)} \\ ⋮ \\ \hat{V} Φ_{M}^{(n - 1)} \end{matrix}) \\ - \sum_{i = 1}^{n} (\begin{matrix} E_{11}^{(i)} & \dots & E_{M 1}^{(i)} \\ ⋮ & ⋮ \\ E_{1 M}^{(i)} & \dots & E_{M M}^{(i)} \end{matrix}) (\begin{matrix} Φ_{1}^{(n - i)} \\ ⋮ \\ Φ_{M}^{(n - i)} \end{matrix}) \end{align}

(22)

or

\begin{align} (\begin{matrix} Φ_{1}^{(n)} \\ ⋮ \\ Φ_{M}^{(n)} \end{matrix}) & = {\hat{R}}_{I} \{(\begin{matrix} \hat{V} Φ_{1}^{(n - 1)} \\ ⋮ \\ \hat{V} Φ_{M}^{(n - 1)} \end{matrix})) \\ (- \sum_{i = 1}^{n - 1} (\begin{matrix} E_{11}^{(i)} & \dots & E_{M 1}^{(i)} \\ ⋮ & ⋮ \\ E_{1 M}^{(i)} & \dots & E_{M M}^{(i)} \end{matrix}) (\begin{matrix} Φ_{1}^{(n - i)} \\ ⋮ \\ Φ_{M}^{(n - i)} \end{matrix})\}, \end{align}

(23)

where the definition of ${\hat{R}}_{I}$ remains unchanged [Eq. (19)]. Multiplying the Jth row of Eq. (22) by $Φ_{I}^{(0) *}$ from the left and integrating for all I and J (1 ≤ I, J ≤ M), we have, for n ≥ 1,

\begin{align} E^{(n)} & \equiv (\begin{matrix} E_{11}^{(n)} & \dots & E_{1 M}^{(n)} \\ ⋮ & ⋮ \\ E_{M 1}^{(n)} & \dots & E_{M M}^{(n)} \end{matrix}) \\ = (\begin{matrix} 〈 Φ_{1}^{(0)} | \hat{V} | Φ_{1}^{(n - 1)} 〉 & \dots & 〈 Φ_{1}^{(0)} | \hat{V} | Φ_{M}^{(n - 1)} 〉 \\ ⋮ & ⋮ \\ 〈 Φ_{M}^{(0)} | \hat{V} | Φ_{1}^{(n - 1)} 〉 & \dots & 〈 Φ_{M}^{(0)} | \hat{V} | Φ_{M}^{(n - 1)} 〉 \end{matrix}), \end{align}

(24)

which is Hermitian.⁶¹ Owing to the degeneracy, $Φ_{1}^{(0)}$ through $Φ_{M}^{(0)}$ can mix with one another and remain as valid zeroth-order references insofar as they are orthonormal. We seek a unique (up to a phase) set of the zeroth-order wave functions that brings E⁽¹⁾ through E⁽ⁿ⁾ into a diagonal form so that the above M-fold coupled Schrödinger equation [Eq. (21)] is separated into M independent Schrödinger equations of the form of Eq. (16). The eigenvalues of E⁽ⁿ⁾ are then identified as the nth-order HCPT energy corrections,⁶¹ although the foregoing formulation is considerably simpler than that of Ref. 61. Given above is an abridged version that has just enough details to make the recursion for Ω⁽ⁿ⁾ self-contained.

At this point, it may appear hopeless to try to derive compact analytical formulas of $⟨ E_{I}^{(n)} ⟩$ written in terms of molecular integrals because eigenvalues are only procedurally defined and cannot be written in a closed form (for M ≥ 5 according to the Abel–Ruffini theorem). Fortunately, however, we can still derive analytical formulas without knowing the eigenvalues:^54,55 Since these energy corrections are thermally averaged with the same Boltzmann weight within each degenerate subspace, we only need the sum of the eigenvalues (not the individual eigenvalues) in order to calculate the thermal average correctly. This sum is equal to the sum of diagonal elements owing to the trace invariance for a cyclic permutation of a matrix product,

Tr U^{- 1} E^{(i)} U = Tr E^{(i)},

(25)

Tr U^{- 1} E^{(i)} U U^{- 1} E^{(j)} U = Tr E^{(i)} E^{(j)},

(26)

etc., where U is the unitary matrix that brings all of E⁽ⁱ⁾ (1 ≤ i ≤ n) into a diagonal form. Each matrix element (before diagonalization) can be written in a closed form, lending the trace and thus thermal average to analytical expressions.

In our previous studies,^54,55 we used the Slater–Condon rules⁷⁴ to evaluate these matrix elements and derived compact analytical formulas for the thermal averages with custom-made combinatorial identities. Such an order-by-order approach reaches an impasse at higher orders, and hence, in this study, we switch to second-quantized and diagrammatic formulations expounded on in the subsequent sections. For the sum-over-states recursion and general-order algorithm based on it, we can still diagonalize E⁽ⁱ⁾ and use its eigenvalues in principle, but we elect to leave E⁽ⁱ⁾ undiagonalized because the diagonalization step is simply superfluous.

The foregoing argument must not be misconstrued to mean that only diagonal elements of E⁽ⁱ⁾ matter. As we demonstrate in Sec. III, off-diagonal elements of E⁽ⁱ⁾ enter thermal averages via the renormalization term [the second term of Eqs. (22) or (23)], giving rise to unique diagrams—the renormalization diagrams with displaced resolvent lines (distinguished from the anomalous diagrams with missing resolvents) appearing for the first time at the third order.

To summarize, ⟨E⁽ⁿ⁾⟩, ⟨E⁽ⁱ⁾E⁽ⁿ⁻ⁱ⁾⟩, etc., are evaluated by thermally averaging the traces of E⁽ⁱ⁾, E⁽ⁱ⁾E⁽ⁿ⁻ⁱ⁾, etc., using the energy recursion [Eq. (24)], which, in turn, depends on the recursion for wave functions [Eq. (23)]. In other words, for the purpose of executing the recursions, we may use the following substitutions, where $E_{I}^{(i)}$ is the Ith eigenvalue, whereas $E_{I J}^{(i)}$ is the IJth element of the matrix E⁽ⁱ⁾:

〈 E_{I}^{(i)} 〉 = 〈 E_{I I}^{(i)} 〉,

(27)

〈 E_{I}^{(i)} E_{I}^{(j)} 〉 = 〈 E_{I J}^{(i)} {\hat{P}}_{I} E_{J I}^{(j)} 〉,

(28)

etc., where the inner projector, ${\hat{P}}_{I}$ ⁠, restricts the summation over J (whose summation symbol is suppressed according to Einstein’s convention) to determinants within the degenerate subspace of I. It is not only advantageous but also necessary to treat the energy corrections that enter the recursions as a nondiagonal matrix rather than scalars.

D. Recursion for μ⁽ⁿ⁾

The recursion for μ⁽ⁿ⁾ comes from⁵⁵

{(\frac{\partial Ω^{(n)}}{\partial μ^{(0)}})}_{μ^{(1)}, \dots, μ^{(n)}} = 0 .

(29)

Substituting this to Eq. (14), we have

\frac{\partial}{\partial μ^{(0)}} \frac{Ξ^{(n)}}{Ξ^{(0)}} = 0 .

(30)

Using the recursion for Ξ⁽ⁿ⁾/Ξ⁽⁰⁾ [Eq. (11)] in the above equation, we arrive at the following recursion for μ⁽ⁿ⁾ (n ≥ 1):

\begin{matrix} μ^{(n)} ⟨ N_{I}^{2} - {\bar{N}}^{2} ⟩ = ⟨ E_{I}^{(n)} (N_{I} - \bar{N}) ⟩ + \frac{(- β)}{2!} \sum_{i = 1}^{n - 1} ⟨ D_{I}^{(i)} D_{I}^{(n - i)} (N_{I} - \bar{N}) ⟩ \\ + \frac{{(- β)}^{2}}{3!} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} ⟨ D_{I}^{(i)} D_{I}^{(j)} D_{I}^{(n - i - j)} (N_{I} - \bar{N}) ⟩ \\ + \dots + \frac{{(- β)}^{n - 1}}{n!} ⟨ {(D_{I}^{(1)})}^{n} (N_{I} - \bar{N}) ⟩, \end{matrix}

(31)

where we used the identity

\frac{\partial}{\partial μ^{(0)}} 〈 X_{I} 〉 = 〈 X_{I} N_{I} 〉 - 〈 X_{I} 〉 〈 N_{I} 〉 .

(32)

Its application to the first term of Eq. (11) leads to

\frac{\partial}{\partial μ^{(0)}} ⟨ D_{I}^{(n)} ⟩ = ⟨ E_{I}^{(n)} (N_{I} - \bar{N}) ⟩ - μ^{(n)} ⟨ N_{I}^{2} - {\bar{N}}^{2} ⟩,

(33)

which was also used in Eq. (31).

E. Recursion for U⁽ⁿ⁾

Differentiating Eq. (4) n times with λ, we obtain

U^{(n)} - μ^{(n)} \bar{N} = {(\frac{\partial (- β Ω^{(n)})}{\partial (- β)})}_{μ^{(0)}, \dots, μ^{(n)}} .

(34)

Substituting Eq. (14) into this, we arrive at

\begin{matrix} U^{(n)} - μ^{(n)} \bar{N} = \frac{\partial}{\partial (- β)} \frac{Ξ^{(n)}}{Ξ^{(0)}} - (- β) \sum_{i = 1}^{n - 1} (U^{(i)} - μ^{(i)} \bar{N}) Ω^{(n - i)} \\ - \frac{{(- β)}^{2}}{2!} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} (U^{(i)} - μ^{(i)} \bar{N}) Ω^{(j)} Ω^{(n - i - j)} \\ - \dots - \frac{{(- β)}^{n - 1}}{(n - 1)!} (U^{(1)} - μ^{(1)} \bar{N}) {(Ω^{(1)})}^{n - 1} . \end{matrix}

(35)

The first term can be expanded as

\begin{matrix} \frac{\partial}{\partial (- β)} \frac{Ξ^{(n)}}{Ξ^{(0)}} & = & ⟨ D_{I}^{(n)} ⟩ - (- β) ⟨ D_{I}^{(n)} ⟩ ⟨ D_{I}^{(0)} ⟩ + (- β) \sum_{i = 1}^{n} ⟨ D_{I}^{(i)} D_{I}^{(n - i)} ⟩ - \frac{{(- β)}^{2}}{2!} \sum_{i = 1}^{n - 1} ⟨ D_{I}^{(i)} D_{I}^{(n - i)} ⟩ ⟨ D_{I}^{(0)} ⟩ \\ + \frac{{(- β)}^{2}}{2!} \sum_{i = 1}^{n - 1} \sum_{j = 1}^{n - i} ⟨ D_{I}^{(i)} D_{I}^{(j)} D_{I}^{(n - i - j)} ⟩ - \frac{{(- β)}^{3}}{3!} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} ⟨ D_{I}^{(i)} D_{I}^{(j)} D_{I}^{(n - i - j)} ⟩ ⟨ D_{I}^{(0)} ⟩ \\ + \frac{{(- β)}^{3}}{3!} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} \sum_{k = 1}^{n - i - j} ⟨ D_{I}^{(i)} D_{I}^{(j)} D_{I}^{(k)} D_{I}^{(n - i - j - k)} ⟩ - \dots - \frac{{(- β)}^{n}}{n!} ⟨ {(D_{I}^{(1)})}^{n} ⟩ ⟨ D_{I}^{(0)} ⟩ + \frac{{(- β)}^{n}}{n!} ⟨ {(D_{I}^{(1)})}^{n} D_{I}^{(0)} ⟩, \end{matrix}

(36)

where we used Eq. (11) and

\frac{\partial}{\partial (- β)} 〈 X_{I} 〉 = 〈 X_{I} D_{I}^{(0)} 〉 - 〈 X_{I} 〉 〈 D_{I}^{(0)} 〉 .

(37)

Therefore, the recursion for U⁽ⁿ⁾ (n ≥ 1) is given by

\begin{matrix} U^{(n)} & = & 〈 E_{I}^{(n)} 〉 + (- β) \sum_{i = 1}^{n} 〈 D_{I}^{(i)} D_{I}^{(n - i)} 〉 \\ - (- β) \sum_{i = 1}^{n} Ω^{(i)} (U^{(n - i)} - μ^{(n - i)} \bar{N}) \\ + \frac{{(- β)}^{2}}{2!} \sum_{i = 1}^{n - 1} \sum_{j = 1}^{n - i} 〈 D_{I}^{(i)} D_{I}^{(j)} D_{I}^{(n - i - j)} 〉 \\ - \frac{{(- β)}^{2}}{2!} \sum_{i = 1}^{n - 1} \sum_{j = 1}^{n - i} Ω^{(i)} Ω^{(j)} (U^{(n - i - j)} - μ^{(n - i - j)} \bar{N}) \\ + \dots + \frac{{(- β)}^{n}}{n!} 〈 {(D_{I}^{(1)})}^{n} D_{I}^{(0)} 〉 \\ - \frac{{(- β)}^{n}}{n!} {(Ω^{(1)})}^{n} (U^{(0)} - μ^{(0)} \bar{N}), \end{matrix}

(38)

where we used $⟨ D_{I}^{(0)} ⟩ = U^{(0)} - μ^{(0)} \bar{N}$ ⁠.

F. Recursion for S⁽ⁿ⁾

The recursion for S⁽ⁿ⁾ is a concatenation of the recursions for Ω⁽ⁿ⁾, μ⁽ⁿ⁾, and U⁽ⁿ⁾ given above because

S^{(n)} = k_{B} β (U^{(n)} - μ^{(n)} \bar{N} - Ω^{(n)}) .

(39)

III. SECOND-QUANTIZED DERIVATION

In this section, we derive the reduced analytical formulas of Ω⁽ⁿ⁾ and μ⁽ⁿ⁾ for the few lowest n’s,^54,55 starting from the sum-over-states analytical formulas obtained from the recursions. We rely on the rules of normal-ordered second quantization at finite temperature,^28,42 expounded on in Appendix C. Second-quantized derivations of U⁽ⁿ⁾ and S⁽ⁿ⁾ follow essentially the same procedure and will not be repeated.

A. Zeroth order

The zeroth-order formulas initiate the recursions. They correspond to the Fermi–Dirac theory discussed in many textbooks.^42–45 We therefore only document the results,^54,55,57

Ω^{(0)} = E_{nuc.} + \frac{1}{β} \sum_{p} \ln f_{p}^{+},

(40)

U^{(0)} = E_{nuc.} + \sum_{p} ϵ_{p} f_{p}^{-},

(41)

S^{(0)} = - k_{B} \sum_{p} (f_{p}^{-} \ln f_{p}^{-} + f_{p}^{+} \ln f_{p}^{+}),

(42)

where E_nuc. is the nucleus–nucleus repulsion energy and $f_{p}^{-}$ (⁠ $f_{p}^{+}$ ⁠) is the Fermi–Dirac occupancy (vacancy) function for the pth spin orbital given by

f_{p}^{-} = \frac{1}{1 + e^{β (ϵ_{p} - μ^{(0)})}},

(43)

f_{p}^{+} = \frac{e^{β (ϵ_{p} - μ^{(0)})}}{1 + e^{β (ϵ_{p} - μ^{(0)})}} .

(44)

Finally, μ⁽⁰⁾ is the root of the following equation embodying the zeroth-order electroneutrality condition:

\bar{N} = \sum_{p} f_{p}^{-} .

(45)

B. First order

According to the recursion for Ω⁽ⁿ⁾ [Eq. (15)], using Einstein’s convention inside thermal averages ⟨⋯⟩, we can write Ω⁽¹⁾ as

Ω^{(1)} = ⟨ D_{I}^{(1)} ⟩ = ⟨ ⟨ I | \hat{V} | I ⟩ ⟩ - μ^{(1)} ⟨ ⟨ I | \hat{N} | I ⟩ ⟩,

(46)

where $⟨ ⟨ I | \hat{V} | I ⟩ ⟩$ is the thermal average of $⟨ I | \hat{V} | I ⟩$ and I stands for $Φ_{I}^{(0)}$ ⁠. Consulting with Appendix C, we can evaluate this using the normal-ordered second quantization as

\begin{matrix} Ω^{(1)} & = & ⟨ ⟨ I | ⟨ E_{I}^{(1)} ⟩ + F_{p q} {{\hat{p}}^{†} \hat{q}} + (1 / 4) ⟨ p q ‖ r s ⟩ {{\hat{p}}^{†} {\hat{q}}^{†} \hat{s} \hat{r}} | I ⟩ ⟩ \\ - μ^{(1)} ⟨ ⟨ I | \bar{N} + {{\hat{p}}^{†} \hat{p}} | I ⟩ ⟩ \end{matrix}

(47)

= ⟨ E_{I}^{(1)} ⟩ - μ^{(1)} \bar{N}

(48)

= \sum_{p} F_{p p} f_{p}^{-} - \frac{1}{2} \sum_{p, q} ⟨ p q ‖ p q ⟩ f_{p}^{-} f_{q}^{-} - μ^{(1)} \bar{N}

(49)

\equiv {⟨ D_{I}^{(1)} ⟩}_{L} .

(50)

In Eq. (47), we substituted the normal-ordered perturbation operator [Eq. (C41)] and number operator [Eq. (C29)] derived in Appendix C 3. Going from Eq. (47) to (48), we invoked the rule that the thermal average of a normal-ordered product of operators is always zero, allowing only the constant terms in $\hat{V}$ and $\hat{N}$ to survive. In Eq. (49), $⟨ E_{I}^{(1)} ⟩$ is defined by Eq. (C36) and F_pq is given by Eq. (C37), whereas ⟨pq‖rs⟩ is an antisymmetrized two-electron integral and μ⁽¹⁾ will be considered below. Therefore, as in the zero-temperature normal ordering,⁶³ an effort has been prepaid to derive the analytical formula for $⟨ E_{I}^{(1)} ⟩$ during the derivation of the normal-ordered form of $\hat{V}$ ⁠. The end result [Eq. (49)] is recognized as diagrammatically linked and thus size-consistent^63,64 and can be symbolically written as Eq. (50) with subscript “L” standing for “linked.” Algebraically, a linked term is the one that is not a simple product of two or more extensive quantities (μ⁽ⁿ⁾, $f_{p}^{\pm}$ ⁠, and F_pq are intensive, while $⟨ E_{I}^{(n)} ⟩$ and $\bar{N}$ are extensive).

We start with the sum-over-states formula of μ⁽¹⁾ [Eq. (31)], which reads

μ^{(1)} 〈 N_{I}^{2} - {\bar{N}}^{2} 〉 = 〈 E_{I}^{(1)} (N_{I} - \bar{N}) 〉 .

(51)

The thermal average in the left-hand side is evaluated as

(52)

which is linked owing to the cancellation of the unlinked terms, ${\bar{N}}^{2}$ ⁠. The right-hand side of Eq. (51) is evaluated as

(53)

which is also linked. Taken together, Eq. (51) can be symbolically rewritten as

μ^{(1)} 〈 N_{I}^{2} - {\bar{N}}^{2} 〉_{L} = 〈 E_{I}^{(1)} (N_{I} - \bar{N}) 〉_{L},

(54)

which reduces to

μ^{(1)} \sum_{p} f_{p}^{-} f_{p}^{+} = \sum_{p} F_{p p} f_{p}^{-} f_{p}^{+} .

(55)

For both Ω⁽¹⁾ and μ⁽¹⁾, we recover the reduced analytical formulas reported earlier.^54,55

C. Second order

From the Ω⁽ⁿ⁾ recursion [Eq. (15)], we have

Ω^{(2)} = 〈 D_{I}^{(2)} 〉 + \frac{(- β)}{2!} (〈 D_{I}^{(1)} D_{I}^{(1)} 〉 - Ω^{(1)} Ω^{(1)}),

(56)

whose components are expanded by using the $⟨ E_{I}^{(n)} ⟩$ recursions [Eqs. (23) and (24)] as

\begin{matrix} 〈 D_{I}^{(2)} 〉 & = & 〈 E_{I}^{(2)} 〉 - μ^{(2)} 〈 N_{I} 〉 \\ = & 〈 E_{I I}^{(2)} 〉 - μ^{(2)} 〈 〈 I | \hat{N} | I 〉 〉 \\ = & 〈 〈 I | \hat{V} {\hat{R}}_{I} \hat{V} | I 〉 〉 - μ^{(2)} \bar{N} \end{matrix}

(57)

and

\begin{matrix} ⟨ D_{I}^{(1)} D_{I}^{(1)} ⟩ & = & ⟨ E_{I}^{(1)} E_{I}^{(1)} ⟩ - 2 μ^{(1)} ⟨ E_{I}^{(1)} N_{I} ⟩ + {(μ^{(1)})}^{2} ⟨ N_{I}^{2} ⟩ \\ = & ⟨ E_{I J}^{(1)} {\hat{P}}_{I} E_{J I}^{(1)} ⟩ - 2 μ^{(1)} ⟨ E_{I}^{(1)} (N_{I} - \bar{N}) ⟩ \\ - 2 μ^{(1)} ⟨ E_{I}^{(1)} ⟩ \bar{N} + {(μ^{(1)})}^{2} ⟨ N_{I}^{2} - {\bar{N}}^{2} ⟩ + {(μ^{(1)})}^{2} {\bar{N}}^{2} \\ = & ⟨ ⟨ I | \hat{V} {\hat{P}}_{I} \hat{V} | I ⟩ ⟩ - 2 μ^{(1)} \sum_{p} F_{p p} f_{p}^{-} f_{p}^{+} + {(μ^{(1)})}^{2} \\ \times \sum_{p} f_{p}^{-} f_{p}^{+} - 2 μ^{(1)} ⟨ E_{I}^{(1)} ⟩ \bar{N} + {(μ^{(1)})}^{2} {\bar{N}}^{2}, \end{matrix}

(58)

where Eqs. (52) and (53) were used. The last two terms of Eq. (58) are unlinked, growing quadratically with size, which are, however, easily seen to be canceled by the corresponding contributions in the Ω⁽¹⁾Ω⁽¹⁾ term of Eq. (56).

Let us evaluate $⟨ E_{I}^{(2)} ⟩$ using the normal-ordered second-quantization rules for resolvent operators described in Appendix C 4. Using Einstein’s convention of implied summations, we write

(59)

where i, j, a, b, as well as p through w run over all spin orbitals and “× 16” means that there are 16 distinct, equal-valued full contraction patterns. The constant term $⟨ E_{I}^{(1)} ⟩$ in $\hat{V}$ [see Eq. (C41)] does not contribute because the resolvent ${\hat{R}}_{I}$ erases all determinants in the degenerate subspace of I (see Appendix C 4). The latter is also responsible for the restriction (“denom. ≠ 0”) of the summations to nonzero denominators. Although we did not examine the linkedness of the second-order HCPT energy corrections (be they eigenvalues or the trace of E⁽²⁾), we could nonetheless establish the linkedness of their thermal average, $⟨ E_{I}^{(2)} ⟩$ ⁠.

Next, we evaluate $⟨ E_{I}^{(1)} E_{I}^{(1)} ⟩$ ⁠, recalling the rules for the inner projector ${\hat{P}}_{I}$ (see Appendix C 4),

(60)

Unlike in $⟨ E_{I}^{(2)} ⟩$ ⁠, the inner projector ${\hat{P}}_{I}$ allows the constant part $⟨ E_{I}^{(1)} ⟩$ of $\hat{V}$ to survive, resulting in the unlinked contribution ${⟨ E_{I}^{(1)} ⟩}_{L} {⟨ E_{I}^{(1)} ⟩}_{L}$ ⁠, which is canceled by the corresponding contribution in the Ω⁽¹⁾Ω⁽¹⁾ term of Eq. (56). It is also responsible for the restrictions (“denom. = 0”) to the summations to the cases whose fictitious denominator (ϵ_p − ϵ_q in the first sum or ϵ_p + ϵ_q −ϵ_r − ϵ_s in the second sum) is zero. See the discussion regarding Eq. (28) and Appendix C 4 for the justification of this rule.

Substituting all these results into Eq. (56), we observe exact mutual cancellations of all unlinked terms, reproducing the reduced analytical formula reported earlier,⁵⁵

\begin{matrix} Ω^{(2)} & = & \sum_{p, q}^{denom. \neq 0} \frac{| F_{p q} |^{2} f_{p}^{-} f_{q}^{+}}{ϵ_{p} - ϵ_{q}} + \frac{1}{4} \sum_{p, q, r, s}^{denom. \neq 0} \frac{| 〈 p q ‖ r s 〉 |^{2} f_{p}^{-} f_{q}^{-} f_{r}^{+} f_{s}^{+}}{ϵ_{p} + ϵ_{q} - ϵ_{r} - ϵ_{s}} \\ - \frac{β}{2} \sum_{p, q}^{denom. = 0} | F_{p q} |^{2} f_{p}^{-} f_{q}^{+} - \frac{β}{8} \sum_{p, q, r, s}^{denom. = 0} | 〈 p q ‖ r s 〉 |^{2} f_{p}^{-} f_{q}^{-} f_{r}^{+} f_{s}^{+} \\ + β μ^{(1)} \sum_{p} F_{p p} f_{p}^{-} f_{p}^{+} - \frac{β}{2} {(μ^{(1)})}^{2} \sum_{p} f_{p}^{-} f_{p}^{+} - μ^{(2)} \bar{N} \end{matrix}

(61)

\equiv 〈 D_{I}^{(2)} 〉_{L} + \frac{(- β)}{2!} 〈 D_{I}^{(1)} D_{I}^{(1)} 〉_{L},

(62)

which is linked.

As per the μ⁽ⁿ⁾ recursion [Eq. (31)], we have

μ^{(2)} ⟨ N_{I}^{2} - {\bar{N}}^{2} ⟩ = ⟨ E_{I}^{(2)} (N_{I} - \bar{N}) ⟩ + \frac{(- β)}{2!} ⟨ D_{I}^{(1)} D_{I}^{(1)} (N_{I} - \bar{N}) ⟩,

(63)

of which the left-hand side has already been evaluated in Eq. (52). The second term in the right-hand side is further expanded as

\begin{matrix} 〈 D_{I}^{(1)} D_{I}^{(1)} (N_{I} - \bar{N}) 〉 & = & 〈 E_{I}^{(1)} E_{I}^{(1)} (N_{I} - \bar{N}) 〉 \\ - 2 μ^{(1)} 〈 E_{I}^{(1)} N_{I} (N_{I} - \bar{N}) 〉 \\ + {(μ^{(1)})}^{2} 〈 N_{I} N_{I} (N_{I} - \bar{N}) 〉 . \end{matrix}

(64)

The first term of Eq. (63) is evaluated as

(65)

\begin{matrix} = & \sum_{p, q}^{denom. \neq 0} \frac{| F_{p q} |^{2} f_{p}^{-} f_{q}^{+} f_{p}^{+}}{ϵ_{p} - ϵ_{q}} - \sum_{p, q}^{denom. \neq 0} \frac{| F_{p q} |^{2} f_{p}^{-} f_{q}^{+} f_{q}^{-}}{ϵ_{p} - ϵ_{q}} + \sum_{p, q, s}^{denom. \neq 0} \frac{F_{p q} 〈 q s ‖ p s 〉 f_{p}^{-} f_{q}^{+} f_{s}^{-} f_{s}^{+}}{ϵ_{p} - ϵ_{q}} \\ + \sum_{p, q, r}^{denom. \neq 0} \frac{〈 p q ‖ r q 〉 F_{r p} f_{p}^{-} f_{q}^{-} f_{q}^{+} f_{r}^{+}}{ϵ_{p} - ϵ_{r}} + \frac{1}{2} \sum_{p, q, r, s}^{denom. \neq 0} \frac{| 〈 p q ‖ r s 〉 |^{2} f_{p}^{-} f_{q}^{-} f_{r}^{+} f_{s}^{+} f_{p}^{+}}{ϵ_{p} + ϵ_{q} - ϵ_{r} - ϵ_{s}} - \frac{1}{2} \sum_{p, q, r, s}^{denom. \neq 0} \frac{| 〈 p q ‖ r s 〉 |^{2} f_{p}^{-} f_{q}^{-} f_{r}^{+} f_{s}^{+} f_{r}^{-}}{ϵ_{p} + ϵ_{q} - ϵ_{r} - ϵ_{s}} \end{matrix}

(66)

\equiv {⟨ E_{I}^{(2)} (N_{I} - \bar{N}) ⟩}_{L} .

(67)

The constant part [Eq. (C41)] in the first $\hat{V}$ vanishes because the resolvent ${\hat{R}}_{I}$ annihilates it. The constant part in the second $\hat{V}$ does not contribute, either, because it cannot form a valid contraction with $\hat{N} - \bar{N} = {{\hat{p}}^{†} \hat{p}}$ ⁠. As a result, the above thermal average is linked. The first term of Eq. (64) is evaluated similarly as

(68)

\begin{matrix} = & 2 ⟨ E_{I}^{(1)} ⟩ \sum_{p} F_{p p} f_{p}^{-} f_{p}^{+} + \sum_{p, q}^{denom. = 0} | F_{p q} |^{2} f_{p}^{-} f_{q}^{+} f_{p}^{+} - \sum_{p, q}^{denom. = 0} | F_{p q} |^{2} f_{p}^{-} f_{q}^{+} f_{q}^{-} + \sum_{p, q, s}^{denom. = 0} F_{p q} ⟨ q s ‖ p s ⟩ f_{p}^{-} f_{q}^{+} f_{s}^{-} f_{s}^{+} \\ + \sum_{p, q, r}^{denom. = 0} ⟨ p q ‖ r q ⟩ F_{r p} f_{p}^{-} f_{q}^{-} f_{q}^{+} f_{r}^{+} + \frac{1}{2} \sum_{p, q, r, s}^{denom. = 0} | ⟨ p q ‖ r s ⟩ |^{2} f_{p}^{-} f_{q}^{-} f_{r}^{+} f_{s}^{+} f_{p}^{+} - \frac{1}{2} \sum_{p, q, r, s}^{denom. = 0} | ⟨ p q ‖ r s ⟩ |^{2} f_{p}^{-} f_{q}^{-} f_{r}^{+} f_{s}^{+} f_{r}^{-} \end{matrix}

(69)

\equiv 2 {⟨ E_{I}^{(1)} ⟩}_{L} {⟨ E_{I}^{(1)} (N_{I} - \bar{N}) ⟩}_{L} + {⟨ E_{I}^{(1)} E_{I}^{(1)} (N_{I} - \bar{N}) ⟩}_{L} .

(70)

The first term is unlinked, as each factor is extensive, making the product grow quadratically with size. This unlinked term will be canceled exactly (see below).

Likewise, the thermal average in the second term of Eq. (64) is evaluated as

(71)

\begin{matrix} \equiv 〈 E_{I}^{(1)} 〉_{L} 〈 N_{I}^{2} - {\bar{N}}^{2} 〉_{L} + \bar{N} 〈 E_{I}^{(1)} (N_{I} - \bar{N}) 〉_{L} \\ + 〈 E_{I}^{(1)} N_{I} (N_{I} - \bar{N}) 〉_{L} . \end{matrix}

(72)

The first two terms are unlinked.

The thermal average in the third term of Eq. (64) becomes

(73)

\equiv 2 \bar{N} 〈 N_{I}^{2} - {\bar{N}}^{2} 〉_{L} + 〈 N_{I} N_{I} (N_{I} - \bar{N}) 〉_{L} .

(74)

The first term is unlinked.

Substituting these into Eq. (63), we again observe systematic cancellations of unlinked terms with the aid of Eq. (54), leaving

\begin{matrix} μ^{(2)} {⟨ N_{I}^{2} - {\bar{N}}^{2} ⟩}_{L} = {⟨ E_{I}^{(2)} (N_{I} - \bar{N}) ⟩}_{L} + \frac{(- β)}{2!} {⟨ E_{I}^{(1)} E_{I}^{(1)} (N_{I} - \bar{N}) ⟩}_{L} \\ + β μ^{(1)} {⟨ E_{I}^{(1)} N_{I} (N_{I} - \bar{N}) ⟩}_{L} + \frac{(- β)}{2!} {(μ^{(1)})}^{2} \\ \times {⟨ N_{I} N_{I} (N_{I} - \bar{N}) ⟩}_{L} \\ = {⟨ E_{I}^{(2)} (N_{I} - \bar{N}) ⟩}_{L} + \frac{(- β)}{2!} {⟨ D_{I}^{(1)} D_{I}^{(1)} (N_{I} - \bar{N}) ⟩}_{L}, \end{matrix}

(75)

which translates to

\begin{matrix} μ^{(2)} \sum_{p} f_{p}^{-} f_{p}^{+} & = & \sum_{p, q}^{denom. \neq 0} \frac{| F_{p q} |^{2} f_{p}^{-} f_{q}^{+} (f_{p}^{+} - f_{q}^{-})}{ϵ_{p} - ϵ_{q}} + \sum_{p, q, r}^{denom. \neq 0} \frac{(F_{p q} ⟨ q r ‖ p r ⟩ + ⟨ p r ‖ q r ⟩ F_{q p}) f_{p}^{-} f_{q}^{+} f_{r}^{-} f_{r}^{+}}{ϵ_{p} - ϵ_{q}} \\ + \frac{1}{4} \sum_{p, q, r, s}^{denom. \neq 0} \frac{| ⟨ p q ‖ r s ⟩ |^{2} f_{p}^{-} f_{q}^{-} f_{r}^{+} f_{s}^{+} (f_{p}^{+} + f_{q}^{+} - f_{r}^{-} - f_{s}^{-})}{ϵ_{p} + ϵ_{q} - ϵ_{r} - ϵ_{s}} \\ - \frac{β}{2} \sum_{p, q}^{denom. = 0} | F_{p q} |^{2} f_{p}^{-} f_{q}^{+} (f_{p}^{+} - f_{q}^{-}) - \frac{β}{2} \sum_{p, q, r}^{denom. = 0} (F_{p q} ⟨ q r ‖ p r ⟩ + ⟨ p r ‖ q r ⟩ F_{q p}) f_{p}^{-} f_{q}^{+} f_{r}^{-} f_{r}^{+} \\ - \frac{β}{8} \sum_{p, q, r, s}^{denom. = 0} | ⟨ p q ‖ r s ⟩ |^{2} f_{p}^{-} f_{q}^{-} f_{r}^{+} f_{s}^{+} (f_{p}^{+} + f_{q}^{+} - f_{r}^{-} - f_{s}^{-}) \\ + β μ^{(1)} \sum_{p} F_{p p} f_{p}^{-} f_{p}^{+} (f_{p}^{+} - f_{p}^{-}) + β μ^{(1)} \sum_{p, q} ⟨ p q ‖ p q ⟩ f_{p}^{-} f_{p}^{+} f_{q}^{-} f_{q}^{+} \\ - \frac{β}{2} {(μ^{(1)})}^{2} \sum_{p} f_{p}^{-} f_{p}^{+} (f_{p}^{+} - f_{p}^{-}) . \end{matrix}

(76)

This is identified as Eq. (69) of Ref. 55.

D. Third order

We shall derive the reduced analytical formulas for Ω⁽³⁾ using diagrams in Sec. IV. Here, we focus on its most important contribution, ⟨E⁽³⁾⟩, since this is where linked renormalization terms appear for the first time, giving rise to a new, unique class of renormalization diagrams.

According to Eq. (24), the third-order correction to the energy matrix of HCPT is written as

E_{I J}^{(3)} = 〈 I | \hat{V} {\hat{R}}_{I} \hat{V} {\hat{R}}_{I} \hat{V} | J 〉 - \sum_{K} 〈 I | \hat{V} {\hat{R}}_{I} {\hat{R}}_{I} \hat{V} | K 〉 E_{K J}^{(1)},

(77)

where K runs over all states that are degenerate with I and J. The second term is the renormalization term of the Rayleigh–Schrödinger perturbation theory.⁶³

If I is nondegenerate, HCPT reduces to MPPT, and therefore, I = J = K. The first term then consists of the linked and unlinked contributions and is written as

⟨ I | \hat{V} {\hat{R}}_{I} \hat{V} {\hat{R}}_{I} \hat{V} | I ⟩ = {⟨ I | \hat{V} {\hat{R}}_{I} \hat{V} {\hat{R}}_{I} \hat{V} | I ⟩}_{L} + {⟨ I | \hat{V} {\hat{R}}_{I} {\hat{R}}_{I} \hat{V} | I ⟩}_{L} E_{I}^{(1)},

(78)

with the unlinked (second) term canceling exactly the renormalization term [the second term of Eq. (77)], leaving

E_{I}^{(3)} = {⟨ I | \hat{V} {\hat{R}}_{I} \hat{V} {\hat{R}}_{I} \hat{V} | I ⟩}_{L} .

(79)

When I is degenerate, the cancellation of unlinked terms may or may not be as straightforward in individual matrix elements, $E_{I J}^{(3)}$ ⁠. However, the cancellation occurs completely analogously and straightforwardly in the thermal average, $⟨ E_{I}^{(3)} ⟩$ ⁠. Using the trace invariance of E⁽³⁾, we have

\begin{matrix} ⟨ E_{I}^{(3)} ⟩ & = & ⟨ E_{I I}^{(3)} ⟩ \\ = & ⟨ ⟨ I | \hat{V} {\hat{R}}_{I} \hat{V} {\hat{R}}_{I} \hat{V} | I ⟩ ⟩ - ⟨ ⟨ I | \hat{V} {\hat{R}}_{I} {\hat{R}}_{I} \hat{V} | K ⟩ E_{K I}^{(1)} ⟩ \\ = & ⟨ ⟨ I | \hat{V} {\hat{R}}_{I} \hat{V} {\hat{R}}_{I} \hat{V} | I ⟩ ⟩ - ⟨ ⟨ I | \hat{V} {\hat{R}}_{I} {\hat{R}}_{I} \hat{V} | K ⟩ ⟨ K | \hat{V} | I ⟩ ⟩ \\ = & ⟨ ⟨ I | \hat{V} {\hat{R}}_{I} \hat{V} {\hat{R}}_{I} \hat{V} | I ⟩ ⟩ - ⟨ ⟨ I | \hat{V} {\hat{R}}_{I} {\hat{R}}_{I} \hat{V} {\hat{P}}_{I} \hat{V} | I ⟩ ⟩, \end{matrix}

(80)

where summation symbols are again suppressed. Each term of the last line is then decomposed into linked and unlinked contributions. Recalling the normal-ordered form of $\hat{V}$ [Eq. (C41)], we can decompose each term as

(81)

(82)

where the unlinked contraction patterns are indicated by the staple symbols in the second terms, which can be seen to have an equal value. They therefore cancel each other in Eq. (80), leaving linked ⟨E⁽³⁾⟩,

⟨ E_{I}^{(3)} ⟩ = {⟨ ⟨ I | \hat{V} {\hat{R}}_{I} \hat{V} {\hat{R}}_{I} \hat{V} | I ⟩ ⟩}_{L} - {⟨ ⟨ I | \hat{V} {\hat{R}}_{I} {\hat{R}}_{I} \hat{V} {\hat{P}}_{I} \hat{V} | I ⟩ ⟩}_{L} .

(83)

Unlike in MPPT, the renormalization term in HCPT has both unlinked and linked contributions with the former canceling exactly the unlinked contribution in the parent term, while the latter is nonzero and therefore should not be overlooked or neglected. In the above example, this linked renormalized contribution has one resolvent (⁠ ${\hat{R}}_{I}$ ⁠) operator shifted from in between the last two $\hat{V}$ operators to in between the first two $\hat{V}$ operators. Diagrammatically, it will correspond to a closed, connected diagram with no resolvent line between a pair of adjacent vertexes but with two resolvent lines crammed in between the other pair of adjacent vertexes (see Sec. IV C). This unusual diagram with a displaced resolvent occurs at T > 0 and should be distinguished from the anomalous diagrams^48,52 with missing resolvents. In Sec. IV, we will document the reduced analytical formula of $⟨ E_{I}^{(3)} ⟩$ obtained with diagrams.

IV. FEYNMAN DIAGRAMS

The diagrammatic rules are obtained by a systematic modification of the zero-temperature counterpart^42,44,63 and justified by the one-to-many correspondence with the foregoing normal-ordered second-quantization logic at finite temperature. The rules to generate diagrams and to interpret them algebraically are given in Tables I and II, respectively.

TABLE I.

Rules to generate all Ω⁽ⁿ⁾ diagrams in the Hugenholtz style (n ≥ 2). See the text for generating diagrammatic equations for μ⁽ⁿ⁾ and S⁽ⁿ⁾ from the Ω⁽ⁿ⁾ equation.

(1)

Place n vertexes in an unambiguous vertical order using three types of vertexes: filled-circle two-line (F_pq) vertex, numbered two-line (μ) vertex, and filled-circle four-line (⟨pq‖rs⟩) vertex. Count an μ^(m) vertex as m vertexes.

(2)

Connect all vertexes with lines to form a closed connected diagram in all topologically distinct manner. Each two-line (four-line) vertex should be connected with one outgoing and one incoming (two outgoing and two incoming) lines.

(3)

Insert zero through n − 1 resolvent (wiggly) lines into the diagram in all 2ⁿ⁻¹ ways with either zero or one resolvent in between each pair of adjacent vertexes. Call the diagram with n − 1 resolvents a normal diagram and the one with n − 2 or less resolvents an anomalous diagram.

(4)

Starting with each diagram, shift upward one or more resolvents onto existing resolvents. Call the diagram with at least one shifted resolvent a renormalization diagram.

TABLE II.

Rules to transform an Ω⁽ⁿ⁾, μ⁽ⁿ⁾, or S⁽ⁿ⁾ diagram in the Hugenholtz style to an algebraic formula (n ≥ 1).

(1)

Label lines with indices p, q, r, s, etc.

(2)

Associate each filled-circle two-line vertex with F_pq, where p (q) is the outgoing (incoming) line label.

(3)

Associate each vertex with number n with − μ⁽ⁿ⁾δ_pq, where p (q) is the outgoing (incoming) line label.

(4)

The F_pq and μ⁽¹⁾ vertexes can be consolidated into a shaded two-line vertex to reduce the number of diagrams. If this simplification is used, associate the vertex with F_pq − μ⁽¹⁾δ_pq.

(5)

Associate each filled-circle four-line vertex with ⟨pq‖rs⟩, where p, q, r, and s are the left-outgoing, right-outgoing, left-incoming, and right-incoming line labels, respectively.

(6)

Associate a downgoing (upgoing) line with

f_{p}^{-}

(⁠

f_{p}^{+}

⁠). The directionality of a bubble is judged at the farthest point from the vertex or from its point of contact.

(7)

Associate each dashed bubble with β(ϵ_p − μ⁽⁰⁾) additionally.

(8)

Associate each resolvent line with 1/(ϵ_p + ϵ_q ⋯ − ϵ_r −ϵ_s…), where p, q, … (r, s, …) are the labels of the downgoing (upgoing) lines intersecting the resolvent line.

(9)

Sum over all line indices. Restrict the summations to cases that existing denominator factors are nonzero and fictitious (due to missing resolvents) denominator factors are zero.

(10)

Multiply 1/n! for each set of n equivalent lines. Two lines are equivalent when they start from a same vertex and end at a same vertex as well as have the identical line directions. Ignore attached bubbles when judging the equivalence.

(11)

Multiply (−1)^h+l to the diagram with h downgoing lines and l loops. A bubble counts as a loop.

(12)

Multiply (−β)ⁿ/(n + 1)! to the diagram with n missing resolvents. Further multiply with n to an entropy diagram with no dashed bubble.

(13)

Multiply (−1)ⁿ to the diagram with n shifted resolvents.

Each vertex graphically represents a matrix element of the normal-ordered Hamiltonian or number operator. Each line (edge) connecting two vertexes corresponds to a Wick contraction. That only full contractions yield a nonzero thermal average translates to the diagrammatic rule that all lines must be terminated by vertexes at both ends in order to form a closed diagram (no dangling line). An internal contraction is depicted as a “bubble,” i.e., a contraction of a vertex with itself forming a short loop. In a Ω⁽ⁿ⁾ diagram, the bubbles are seen only as a part of the normal-ordered Hamiltonian or number operator such as $⟨ E_{I}^{(1)} ⟩$ ⁠, F_pq, and $\bar{N}$ ⁠. This is because outside of these operator definitions, internal contractions are zero in a thermal average (see Appendix C). An unlinked diagram (“unlinked” is synonymous with “disconnected” for a closed diagram) corresponds to a full contraction pattern that results in a product of two or more extensive scalars with no shared summation index [e.g., the second terms of Eqs. (81) and (82)]. Section V proves that all thermodynamic functions of the finite-temperature MBPT are represented by linked diagrams only.

A. First order

The diagrammatic rules are not so helpful at the zeroth and first orders, but for the sake of completeness, we show in Fig. 1 the diagrammatic formulation of Ω⁽¹⁾. Figure 2 defines the relevant vertexes.

FIG. 1.

View large Download slide

Diagrammatic equations of Ω⁽¹⁾.

FIG. 2.

View large Download slide

Definitions of vertexes (see also Appendix C 3).

Using the interpretation rules in Table II, we can transform each diagram in Fig. 1 as

Ω 1 a = {(- 1)}^{1 + 1} \sum_{p} F_{p p} f_{p}^{-},

(84)

Ω 1 b = {(- 1)}^{2 + 2} {(\frac{1}{2})}^{1} \sum_{p, q} ⟨ p q ‖ p q ⟩ f_{p}^{-} f_{q}^{-},

(85)

Ω 1 c = {(- 1)}^{1 + 1} \sum_{p} (- μ^{(1)}) f_{p}^{-} = - μ^{(1)} \bar{N} .

(86)

Summing them with the signs indicated in the figure, we recover the reduced analytical formula of Ω⁽¹⁾ [Eq. (49)].

We can rewrite the diagrammatic equation in the first line of Fig. 1 into the second line, which will turn out to be more convenient for the diagrammatic derivation of μ⁽¹⁾ and S⁽¹⁾ (see below). An open-circle, two-line vertex stands for h_pq −ϵ_pδ_pq. The algebraic interpretation of Ω1d is then

Ω 1 d = {(- 1)}^{1 + 1} \sum_{p} (h_{p p} - ϵ_{p}) f_{p}^{-},

(87)

leading to the same reduced analytical formula of Ω⁽¹⁾.

The diagrammatic rules for μ⁽ⁿ⁾ are based on Eq. (29) with

\frac{\partial f_{p}^{\pm}}{\partial μ^{(0)}} = \mp β f_{p}^{-} f_{p}^{+} .

(88)

The diagrammatic equation for μ⁽ⁿ⁾ is therefore obtained by “differentiating” each line in the diagrammatic equation of Ω⁽ⁿ⁾ with respect to μ⁽⁰⁾. In view of Eq. (88), this process can be represented as attaching to each line a bubble with the opposite directionality (which is judged at the farthest point from its contact with the line). Note that the point of contact of the bubble is not a vertex, and therefore the bubble and line share the same spin orbital index.

Figure 3 is the diagrammatic equation for μ⁽¹⁾. It is based on the second line of Fig. 1 rather than on the first line. This is because the latter may hide the fact that F_pq contains $f_{r}^{-}$ that also needs to be differentiated with μ⁽⁰⁾. Each diagram is algebraically interpreted by the same rules of Table II as

μ 1 a = {(- 1)}^{1 + 2} \sum_{p} (h_{p p} - ϵ_{p}) f_{p}^{-} f_{p}^{+},

(89)

μ 1 b = {(- 1)}^{2 + 3} {(\frac{1}{2})}^{1} \sum_{p, q} 〈 p q ‖ p q 〉 f_{p}^{-} f_{q}^{-} f_{q}^{+},

(90)

μ 1 c = {(- 1)}^{2 + 3} {(\frac{1}{2})}^{1} \sum_{p, q} 〈 p q ‖ p q 〉 f_{p}^{-} f_{q}^{-} f_{p}^{+},

(91)

μ 1 d = (- μ^{(1)}) {(- 1)}^{1 + 2} \sum_{p} f_{p}^{-} f_{p}^{+} .

(92)

They lead to the same reduced analytical equation for μ⁽¹⁾ as Eq. (55).

FIG. 3.

View large Download slide

Diagrammatic equation of μ⁽¹⁾.

The diagrammatic equation for S⁽ⁿ⁾ utilizes the following identity obtained by differentiating Eq. (5) with λ n times:

T S^{(n)} = β {(\frac{\partial Ω^{(n)}}{\partial β})}_{μ^{(0)}, \dots, μ^{(n)}},

(93)

with

β \frac{\partial f_{p}^{\pm}}{\partial β} = \pm β (ϵ_{p} - μ^{(0)}) f_{p}^{-} f_{p}^{+} .

(94)

Therefore, the diagrammatic rules for S⁽ⁿ⁾ can be obtained by “differentiating” the Ω⁽ⁿ⁾ diagram equation with β, keeping in mind that anomalous diagrams also carry a power of β. A β-differentiation of a line is diagrammatically depicted as attaching to the line a dashed bubble with the opposite directionality.

At the first order, Ω⁽¹⁾ does not have anomalous diagrams, and so only lines (or more precisely their $f_{p}^{-}$ factors) need to be differentiated with β. The diagrammatic equation for S⁽¹⁾ is drawn in Fig. 4. Each diagram is interpreted according to Table II as

S 1 a = {(- 1)}^{1 + 2} \sum_{p} (h_{p p} - ϵ_{p}) β (ϵ_{p} - μ^{(0)}) f_{p}^{-} f_{p}^{+},

(95)

S 1 b = {(- 1)}^{2 + 3} {(\frac{1}{2})}^{1} \sum_{p, q} 〈 p q ‖ p q 〉 f_{p}^{-} f_{q}^{-} β (ϵ_{q} - μ^{(0)}) f_{q}^{+},

(96)

S 1 c = {(- 1)}^{2 + 3} {(\frac{1}{2})}^{1} \sum_{p, q} 〈 p q ‖ p q 〉 f_{p}^{-} f_{q}^{-} β (ϵ_{p} - μ^{(0)}) f_{p}^{+},

(97)

S 1 d = (- μ^{(1)}) {(- 1)}^{1 + 2} \sum_{p} f_{p}^{-} β (ϵ_{p} - μ^{(0)}) f_{p}^{+},

(98)

leading to

T S^{(1)} = - β \sum_{p} (F_{p p} - μ^{(1)}) (ϵ_{p} - μ^{(0)}) f_{p}^{-} f_{p}^{+}

(99)

= - β \sum_{p} (F_{p p} - μ^{(1)}) ϵ_{p} f_{p}^{-} f_{p}^{+} .

(100)

The second equality follows from Eq. (55). In general, for the purpose of deriving the correct S⁽ⁿ⁾ formulas, μ⁽⁰⁾ in Eq. (94) can be replaced by any constant (such as zero) by virtue of Eq. (29).

FIG. 4.

View large Download slide

Diagrammatic equation of TS⁽¹⁾.

The diagrammatic equation and rules for U⁽ⁿ⁾ are a concatenation of those for Ω⁽ⁿ⁾, μ⁽ⁿ⁾, and S⁽ⁿ⁾ given above because

U^{(n)} = Ω^{(n)} + μ^{(n)} \bar{N} + T S^{(n)} .

(101)

The U⁽¹⁾ formula obtained as a combination of Eqs. (49) and (100) is the same as Eq. (61) of Ref. 55.

B. Second order

Following the rules in Table I, we obtain eight diagrams that consist in Ω⁽²⁾, which are shown in Fig. 5. They are interpreted algebraically according to Table II, leading to the same result as Eq. (61).

FIG. 5.

View large Download slide

Diagrammatic equation of Ω⁽²⁾.

The first two (Ω2a and Ω2b) are isomorphic to the usual MBPT(2) energy diagrams (Ω2b being the non-HF term⁶³) with exactly one resolvent (wiggly) line in between the pair of vertexes. The third and fourth (Ω2c and Ω2d) are the corresponding anomalous diagrams,^48,52 whose resolvent line is missing. In these diagrams, the summations over spin orbital indices are limited to those cases whose fictitious denominators are zero. Instead of being divided by the energy denominators, they are divided by 2k_BT to restore the dimension of energy.

Diagrams Ω2e through Ω2g are also unique to the finite-temperature MBPT. The absence of a resolvent line implies that they are technically counted as anomalous diagrams. The summations in these diagrams are therefore also limited to zero fictitious denominators, but this restriction is automatically fulfilled because the μ⁽ⁿ⁾ vertex is zero whenever p ≠ q.

We can compress the diagrammatic equations of Ω⁽ⁿ⁾ by introducing skeleton diagrams⁶³ and a modified Fock operator defined in Fig. 6. A skeleton diagram is the one that is stripped of line directions, line labels, and resolvents. Single skeleton diagram Ω2A in Fig. 7 stands for Ω2a plus Ω2c, which differ from each other by the presence and absence of the resolvent. Skeleton diagram Ω2B in Fig. 8 is the sum of five diagrams. A replacement of the F_pq vertex by the F_pq − μ⁽¹⁾δ_pq vertex has no effect on the normal diagram (Ω2b) because its summation is limited to p ≠ q, but it consolidates the remaining four anomalous diagrams into one. This process is closely related to the Luttinger–Ward prescription⁴⁹ for the Kohn–Luttinger conundrum.^48,56 Together, we can simplify the diagrammatic equation of Ω⁽²⁾ into Fig. 9.

FIG. 6.

View large Download slide

Definitions of vertexes.

FIG. 7.

View large Download slide

A skeleton diagram of Ω⁽²⁾.

FIG. 8.

View large Download slide

Another skeleton diagram of Ω⁽²⁾.

FIG. 9.

View large Download slide

Diagrammatic equation of Ω⁽²⁾ in terms of skeleton diagrams.

The diagrammatic equation of μ⁽²⁾ is drawn in Fig. 10. Translating them into algebraic formulas using the rules in Table II, we arrive at the same reduced analytical formula of μ⁽²⁾ as Eq. (69) of Ref. 55.

FIG. 10.

View large Download slide

Diagrammatic equation of μ⁽²⁾.

The diagrams defining S⁽²⁾ are shown in Fig. 11. Evaluating them according to the rules in Table II, we obtain the reduced analytical formula that is consistent with the U⁽²⁾ formula, i.e., Eq. (74) of Ref. 55. Note that the evaluation rules are slightly different between diagrams S2a and S2f and diagrams Ω2c and Ω2d–Ω2g (cf. rule 12 of Table II) despite their identical appearances. Care must be exercised when combining diagrammatic equations of Ω⁽ⁿ⁾ and S⁽ⁿ⁾ to form a concatenated equation of U⁽ⁿ⁾.

FIG. 11.

View large Download slide

Diagrammatic equation of S⁽²⁾.

C. Third order

The diagrammatic equation of Ω⁽³⁾ is given in Fig. 12. It is the sum of eight skeleton diagrams [cf. the zero-temperature MBPT(3) diagrams in page 133 of Shavitt and Bartlett⁶³] plus two anomalous diagrams involving μ⁽²⁾ and a bubble diagram of μ⁽³⁾. Each skeleton diagram, in turn, consists of five to 15 diagrams differing from one another in line direction and/or placement of the resolvents, shown in Figs. 13–18. In each row of the right-hand sides of these diagram equations, the first diagram (Ω3x1) is the normal (parent) diagram. The second diagram (Ω3x2) is the renormalization diagram with a resolvent shifted upward. The third through fifth diagrams are an anomalous diagram with one or two missing resolvents.

FIG. 12.

View large Download slide

Diagrammatic equation of Ω⁽³⁾.

FIG. 13.

View large Download slide

Definition of the Ω3A skeleton diagram.

FIG. 14.

View large Download slide

Definition of the Ω3B skeleton diagram.

FIG. 15.

View large Download slide

Definition of the Ω3D skeleton diagram.

FIG. 16.

View large Download slide

Definition of the Ω3E skeleton diagram.

FIG. 17.

View large Download slide

Definition of the Ω3G skeleton diagram.

FIG. 18.

View large Download slide

Definition of the Ω3H skeleton diagram.

The diagrams in the first row of the right-hand side of Fig. 13 are interpreted algebraically as

\begin{matrix} Ω 3 a 1 & = & {(- 1)}^{2 + 2} {(\frac{1}{2})}^{3} \\ \times \sum_{p, q, r, s, t, u}^{\begin{array}{c} denom. 1 \neq 0 \\ denom. 2 \neq 0 \end{array}} \frac{⟨ p s ‖ q r ⟩ ⟨ q r ‖ t u ⟩ ⟨ t u ‖ p s ⟩ f_{p}^{-} f_{q}^{+} f_{r}^{+} f_{s}^{-} f_{t}^{+} f_{u}^{+}}{(ϵ_{p} + ϵ_{s} - ϵ_{q} - ϵ_{r}) (ϵ_{p} + ϵ_{s} - ϵ_{t} - ϵ_{u})}, \end{matrix}

(102)

\begin{matrix} Ω 3 a 2 & = & {(- 1)}^{1} {(- 1)}^{2 + 2} {(\frac{1}{2})}^{3} \\ \times \sum_{p, q, r, s, t, u}^{\begin{array}{c} denom. 1 \neq 0 \\ denom. 2 = 0 \end{array}} \frac{⟨ p s ‖ q r ⟩ ⟨ q r ‖ t u ⟩ ⟨ t u ‖ p s ⟩ f_{p}^{-} f_{q}^{+} f_{r}^{+} f_{s}^{-} f_{t}^{+} f_{u}^{+}}{{(ϵ_{p} + ϵ_{s} - ϵ_{q} - ϵ_{r})}^{2}}, \end{matrix}

(103)

\begin{matrix} Ω 3 a 3 & = & \frac{(- β)}{2!} {(- 1)}^{2 + 2} {(\frac{1}{2})}^{3} \\ \times \sum_{p, q, r, s, t, u}^{\begin{array}{c} denom. 1 = 0 \\ denom. 2 \neq 0 \end{array}} \frac{⟨ p s ‖ q r ⟩ ⟨ q r ‖ t u ⟩ ⟨ t u ‖ p s ⟩ f_{p}^{-} f_{q}^{+} f_{r}^{+} f_{s}^{-} f_{t}^{+} f_{u}^{+}}{(ϵ_{p} + ϵ_{s} - ϵ_{t} - ϵ_{u})}, \end{matrix}

(104)

\begin{matrix} Ω 3 a 4 & = & \frac{(- β)}{2!} {(- 1)}^{2 + 2} {(\frac{1}{2})}^{3} \\ \times \sum_{p, q, r, s, t, u}^{\begin{array}{c} denom. 1 \neq 0 \\ denom. 2 = 0 \end{array}} \frac{⟨ p s ‖ q r ⟩ ⟨ q r ‖ t u ⟩ ⟨ t u ‖ p s ⟩ f_{p}^{-} f_{q}^{+} f_{r}^{+} f_{s}^{-} f_{t}^{+} f_{u}^{+}}{(ϵ_{p} + ϵ_{s} - ϵ_{q} - ϵ_{r})}, \end{matrix}

(105)

and

\begin{matrix} Ω 3 a 5 & = & \frac{{(- β)}^{2}}{3!} {(- 1)}^{2 + 2} {(\frac{1}{2})}^{3} \\ \times \sum_{p, q, r, s, t, u}^{\begin{array}{c} denom. 1 = 0 \\ denom. 2 = 0 \end{array}} ⟨ p s ‖ q r ⟩ ⟨ q r ‖ t u ⟩ ⟨ t u ‖ p s ⟩ f_{p}^{-} f_{q}^{+} f_{r}^{+} f_{s}^{-} f_{t}^{+} f_{u}^{+}, \end{matrix}

(106)

where “denom.1” and “denom.2” denote the first and second denominator factors, respectively, of the normal diagram [Eq. (102)]. Owing to this provision, none of the above summations encounters a division by zero.

Unlike in the zero-temperature MBPT, where the renormalization term vanishes upon canceling unlinked contribution in the parent diagram, in the finite-temperature MBPT, it has both linked and unlinked contributions, the latter canceling exactly the unlinked contribution in the parent diagram, while the former persists as the renormalization diagram Ω3a2 [Eq. (103)] with its resolvent shifted upward. Unlike in anomalous diagrams, no resolvent is deleted in a renormalization diagram, and therefore, it does not have to be divided by k_BT to restore the correct dimension of energy. Instead, it is multiplied by the parity of (−1) raised to the power of the number of shifted resolvents (which is equal to the number/nestedness of diagram insertions).⁶³ It is not clear if the time-dependent, diagrammatic derivation^46–51 or the density-matrix formulation⁵² of the conventional finite-temperature MBPT properly accounts for such diagrams, the neglect of which makes it erroneous at the third and higher orders and nonconvergent to the finite-temperature FCI.⁵⁷ The renormalization diagrams may be related to the linked–disconnected diagrams in ΔMPn.⁷⁵

Diagrams Ω3a3, Ω3a4, and Ω3a5 are anomalous diagrams of which one or two resolvents are missing. The summation is correspondingly restricted to the cases where the missing denominator factors are zero. For each missing denominator factor, the diagram is divided by k_BT times some integer factor. The presence of these anomalous diagrams is considered well known.^48,52

The remaining normal diagrams (Ω3x1) are interpreted algebraically as

\begin{matrix} Ω 3 b 1 & = & {(- 1)}^{4 + 2} {(\frac{1}{2})}^{3} \\ \times \sum_{p, q, r, s, t, u}^{\begin{array}{c} denom. 1 \neq 0 \\ denom. 2 \neq 0 \end{array}} \frac{⟨ q r ‖ p s ⟩ ⟨ t u ‖ q r ⟩ ⟨ p s ‖ t u ⟩ f_{p}^{+} f_{q}^{-} f_{r}^{-} f_{s}^{+} f_{t}^{-} f_{u}^{-}}{(ϵ_{q} + ϵ_{r} - ϵ_{p} - ϵ_{s}) (ϵ_{t} + ϵ_{u} - ϵ_{p} - ϵ_{s})}, \end{matrix}

(107)

\begin{matrix} Ω 3 c 1 & = & {(- 1)}^{3 + 2} \\ \times \sum_{p, q, r, s, t, u}^{\begin{array}{c} denom. 1 \neq 0 \\ denom. 2 \neq 0 \end{array}} \frac{⟨ q s ‖ p r ⟩ ⟨ t r ‖ q u ⟩ ⟨ p u ‖ t s ⟩ f_{p}^{+} f_{q}^{-} f_{r}^{+} f_{s}^{-} f_{t}^{-} f_{u}^{+}}{(ϵ_{q} + ϵ_{s} - ϵ_{p} - ϵ_{r}) (ϵ_{t} + ϵ_{s} - ϵ_{p} - ϵ_{u})}, \end{matrix}

(108)

\begin{matrix} Ω 3 d 1 & = & {(- 1)}^{2 + 2} {(\frac{1}{2})}^{1} \\ \times \sum_{p, q, r, s, t}^{\begin{array}{c} denom. 1 \neq 0 \\ denom. 2 \neq 0 \end{array}} \frac{(F_{p q} - μ^{(1)} δ_{p q}) 〈 q s ‖ r t 〉 〈 r t ‖ p s 〉 f_{p}^{-} f_{q}^{+} f_{r}^{+} f_{s}^{-} f_{t}^{+}}{(ϵ_{p} - ϵ_{q}) (ϵ_{p} + ϵ_{s} - ϵ_{r} - ϵ_{t})}, \end{matrix}

(109)

\begin{matrix} Ω 3 e 1 & = & {(- 1)}^{3 + 2} {(\frac{1}{2})}^{1} \\ \times \sum_{p, q, r, s, t}^{\begin{array}{c} denom. 1 \neq 0 \\ denom. 2 \neq 0 \end{array}} \frac{(F_{q p} - μ^{(1)} δ_{q p}) 〈 r t ‖ q s 〉 〈 p s ‖ r t 〉 f_{p}^{+} f_{q}^{-} f_{r}^{-} f_{s}^{+} f_{t}^{-}}{(ϵ_{q} - ϵ_{p}) (ϵ_{r} + ϵ_{t} - ϵ_{p} - ϵ_{s})}, \end{matrix}

(110)

\begin{matrix} Ω 3 f 1 & = & {(- 1)}^{2 + 2} {(\frac{1}{2})}^{1} \\ \times \sum_{p, q, r, s, t}^{\begin{array}{c} denom. 1 \neq 0 \\ denom. 2 \neq 0 \end{array}} \frac{〈 p r ‖ q s 〉 (F_{s t} - μ^{(1)} δ_{s t}) 〈 q t ‖ p r 〉 f_{p}^{-} f_{q}^{+} f_{r}^{-} f_{s}^{+} f_{t}^{+}}{(ϵ_{p} + ϵ_{r} - ϵ_{q} - ϵ_{s}) (ϵ_{p} + ϵ_{r} - ϵ_{q} - ϵ_{t})}, \end{matrix}

(111)

\begin{matrix} Ω 3 g 1 & = & {(- 1)}^{3 + 2} {(\frac{1}{2})}^{1} \\ \times \sum_{p, q, r, s, t}^{\begin{array}{c} denom. 1 \neq 0 \\ denom. 2 \neq 0 \end{array}} \frac{〈 q s ‖ p r 〉 (F_{t s} - μ^{(1)} δ_{t s}) 〈 p r ‖ q t 〉 f_{p}^{+} f_{q}^{-} f_{r}^{+} f_{s}^{-} f_{t}^{-}}{(ϵ_{q} + ϵ_{s} - ϵ_{p} - ϵ_{r}) (ϵ_{q} + ϵ_{t} - ϵ_{p} - ϵ_{r})}, \end{matrix}

(112)

\begin{matrix} Ω 3 h 1 & = & {(- 1)}^{2 + 1} \\ \times \sum_{p, q, r, s}^{\begin{array}{c} denom. 1 \neq 0 \\ denom. 2 \neq 0 \end{array}} \frac{(F_{q p} - μ^{(1)} δ_{q p}) (F_{r s} - μ^{(1)} δ_{r s}) ⟨ p s ‖ r q ⟩ f_{p}^{+} f_{q}^{-} f_{r}^{-} f_{s}^{+}}{(ϵ_{q} - ϵ_{p}) (ϵ_{q} + ϵ_{r} - ϵ_{p} - ϵ_{s})}, \end{matrix}

(113)

\begin{matrix} Ω 3 i 1 & = & {(- 1)}^{2 + 1} \\ \times \sum_{p, q, r, s}^{\begin{array}{c} denom. 1 \neq 0 \\ denom. 2 \neq 0 \end{array}} \frac{(F_{p q} - μ^{(1)} δ_{p q}) ⟨ r q ‖ p s ⟩ (F_{s r} - μ^{(1)} δ_{s r}) f_{p}^{-} f_{q}^{+} f_{r}^{-} f_{s}^{+}}{(ϵ_{p} - ϵ_{q}) (ϵ_{r} - ϵ_{s})}, \end{matrix}

(114)

\begin{matrix} Ω 3 j 1 & = & {(- 1)}^{1 + 1} \\ \times \sum_{p, q, r}^{\begin{array}{c} denom. 1 \neq 0 \\ denom. 2 \neq 0 \end{array}} \frac{(F_{p q} - μ^{(1)} δ_{p q}) (F_{q r} - μ^{(1)} δ_{q r}) (F_{r p} - μ^{(1)} δ_{r p}) f_{p}^{-} f_{q}^{+} f_{r}^{+}}{(ϵ_{p} - ϵ_{q}) (ϵ_{p} - ϵ_{r})}, \end{matrix}

(115)

and

\begin{matrix} Ω 3 k 1 & = & {(- 1)}^{2 + 1} \sum_{p, q, r}^{\begin{array}{c} denom. 1 \neq 0 \\ denom. 2 \neq 0 \end{array}} \\ \times \frac{(F_{q p} - μ^{(1)} δ_{q p}) (F_{r q} - μ^{(1)} δ_{r q}) (F_{p r} - μ^{(1)} δ_{p r}) f_{p}^{+} f_{q}^{-} f_{r}^{-}}{(ϵ_{q} - ϵ_{p}) (ϵ_{r} - ϵ_{p})} . \end{matrix}

(116)

The algebraic formulas for the accompanying renormalization and anomalous diagrams can be inferred from these formulas by a systematic modification of the denominators, summation restrictions, and prefactors. Therefore, despite the numerousness of these third-order diagrams, their algebraic evaluations and computer implementation are not as horrendous as they seem.

In the zero-temperature MBPT, Ω3B and Ω3C (Ω3E and Ω3F also) are a pair of conjugated diagrams,⁶³ whose values are the complex conjugate of each other. In the finite-temperature formalism, renormalization diagrams destroy such symmetry at a diagram level (not at the whole theory level), and we need to evaluate Ω3B and Ω3C separately, for they have values that are not simply related to each other.

The last three diagrams of Fig. 12 are interpreted as

Ω 3 m = {(- 1)}^{1 + 1} \frac{(- β)}{2!} \sum_{p} (F_{p p} - μ^{(1)}) (- μ^{(2)}) f_{p}^{-} f_{p}^{+},

(117)

Ω 3 o = {(- 1)}^{1 + 1} \sum_{p} (- μ^{(3)}) f_{p}^{-} = - μ^{(3)} \bar{N} .

(118)

Diagrams Ω3m and Ω3n are a pair of conjugate diagrams (since they are not renormalization diagrams), and the algebraic formula of the latter is the same as Eq. (117).

D. Fourth order

At the fourth order, enumerating and interpreting diagrams manually is no longer practicable, and a computerized scheme⁷⁶ may be necessary (not pursued in this study). Here, we focus on how to generate renormalization diagrams, which may be the least trivial step of the diagram enumeration. The nth-order anomalous diagrams are straightforwardly spawned by simply eliminating resolvents in all 2ⁿ⁻¹ − 1 ways.

Figure 19 illustrates how all renormalization diagrams are generated systematically from a normal diagram. The latter has exactly one resolvent line in between every pair of adjacent vertexes. The corresponding renormalization diagrams are obtained by shifting one or more (up to n − 2) resolvent lines upward in all distinct ways. However, it may be safer to use the method of the Brueckner brackets,^66,70 which enumerates all renormalization terms by indicating how various operators are Wick contracted in such a way that the result is unlinked (in the case of MPPT). In our case, many reference states are degenerate and these brackets are not necessarily unlinked, as emphasized earlier.

FIG. 19.

View large Download slide

A normal diagram (the leftmost) and its renormalization diagrams in $⟨ E_{I}^{(4)} ⟩$ and their corresponding Brueckner brackets. Each wiggly line represents a resolvent. Line directions and line labels are suppressed. Signs in front of the Brueckner brackets are absorbed in the algebraic interpretations of the respective diagrams.

For instance, the first renormalization diagram (the second diagram) in Fig. 19 corresponds to the Brueckner bracket,

(119)

where V stands for $\hat{V}$ ⁠, P stands for ${\hat{P}}_{I}$ ⁠, R stands for ${\hat{R}}_{I}$ ⁠, and ⟨⋯⟩ is a thermal average, i.e., ⟨⟨I|⋯|I⟩⟩. Associating V’s from left to right with the vertexes from top to bottom, Eq. (119) is the diagram in which there is one resolvent in between the first (top) and second vertexes, two in between the second and third, and none in between the third and fourth (bottom). In other words, the bottom resolvent has been shifted up to the middle. The parity associated with this bracket is (−1) because there is one insertion (one inner bracket),^63,70 and this factor is properly accounted for by rule 13 of Table II. The fourth diagram in Fig. 19 is an example of double insertions,^63,70 whose parity is (−1)². Generally, an n-tuple insertion or n-fold nested insertion carries a factor of (−1)ⁿ.

There should also be anomalous–renormalization diagrams.

V. LINKED-DIAGRAM THEOREM

We have established in Eqs. (50) and (62) that Ω⁽¹⁾ and Ω⁽²⁾ are diagrammatically linked and thus size-consistent. We may summarize this as

Ω^{(1)} = 〈 D_{I}^{(1)} 〉 = 〈 D_{I}^{(1)} 〉_{L},

(120)

\begin{matrix} Ω^{(2)} & = & 〈 D_{I}^{(2)} 〉 + \frac{(- β)}{2!} (〈 D_{I}^{(1)} D_{I}^{(1)} 〉 - Ω^{(1)} Ω^{(1)}) \\ = & 〈 D_{I}^{(2)} 〉_{L} + \frac{(- β)}{2!} (〈 D_{I}^{(1)} D_{I}^{(1)} 〉 - 〈 D_{I}^{(1)} 〉_{L} 〈 D_{I}^{(1)} 〉_{L}) \\ = & 〈 D_{I}^{(2)} 〉_{L} + \frac{(- β)}{2!} 〈 D_{I}^{(1)} D_{I}^{(1)} 〉_{L}, \end{matrix}

(121)

where subscript L means linked. These can be used to show that Ω⁽³⁾ is also linked,

\begin{matrix} Ω^{(3)} & = & ⟨ D_{I}^{(3)} ⟩ + \frac{(- β)}{2!} (⟨ D_{I}^{(1)} D_{I}^{(2)} ⟩ - Ω^{(1)} Ω^{(2)}) \\ + \frac{(- β)}{2!} (⟨ D_{I}^{(2)} D_{I}^{(1)} ⟩ - Ω^{(2)} Ω^{(1)}) \\ + \frac{{(- β)}^{2}}{3!} (⟨ D_{I}^{(1)} D_{I}^{(1)} D_{I}^{(1)} ⟩ - Ω^{(1)} Ω^{(1)} Ω^{(1)}) \end{matrix}

(122)

\begin{matrix} = & {⟨ D_{I}^{(3)} ⟩}_{L} + \frac{(- β)}{2!} (⟨ D_{I}^{(1)} D_{I}^{(2)} ⟩ - {⟨ D_{I}^{(1)} ⟩}_{L} {⟨ D_{I}^{(2)} ⟩}_{L}) \\ + \frac{(- β)}{2!} (⟨ D_{I}^{(2)} D_{I}^{(1)} ⟩ - {⟨ D_{I}^{(2)} ⟩}_{L} {⟨ D_{I}^{(1)} ⟩}_{L}) \\ + \frac{{(- β)}^{2}}{3!} (⟨ D_{I}^{(1)} D_{I}^{(1)} D_{I}^{(1)} ⟩ - 3 {⟨ D_{I}^{(1)} ⟩}_{L} {⟨ D_{I}^{(1)} D_{I}^{(1)} ⟩}_{L}) \\ (- {⟨ D_{I}^{(1)} ⟩}_{L} {⟨ D_{I}^{(1)} ⟩}_{L} {⟨ D_{I}^{(1)} ⟩}_{L}) \end{matrix}

(123)

\begin{matrix} = & {⟨ D_{I}^{(3)} ⟩}_{L} + \frac{(- β)}{2!} {⟨ D_{I}^{(1)} D_{I}^{(2)} ⟩}_{L} + \frac{(- β)}{2!} {⟨ D_{I}^{(2)} D_{I}^{(1)} ⟩}_{L} \\ + \frac{{(- β)}^{2}}{3!} {⟨ D_{I}^{(1)} D_{I}^{(1)} D_{I}^{(1)} ⟩}_{L} . \end{matrix}

(124)

They lead us to speculate that the recursion of Ω⁽ⁿ⁾ [Eq. (15)] is always reducible to a linked form,

\begin{matrix} Ω^{(n)} & = & 〈 D_{I}^{(n)} 〉 + \frac{(- β)}{2!} \sum_{i = 1}^{n - 1} (〈 D_{I}^{(i)} D_{I}^{(n - i)} 〉 - Ω^{(i)} Ω^{(n - i)}) \\ + \frac{{(- β)}^{2}}{3!} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} (〈 D_{I}^{(i)} D_{I}^{(j)} D_{I}^{(n - i - j)} 〉) - (Ω^{(i)} Ω^{(j)} Ω^{(n - i - j)}) \\ + \dots + \frac{{(- β)}^{n - 1}}{n!} (〈 {(D_{I}^{(1)})}^{n} 〉 - {(Ω^{(1)})}^{n}) \end{matrix}

(125)

\begin{matrix} = & 〈 D_{I}^{(n)} 〉_{L} + \frac{(- β)}{2!} \sum_{i = 1}^{n - 1} 〈 D_{I}^{(i)} D_{I}^{(n - i)} 〉_{L} \\ + \frac{{(- β)}^{2}}{3!} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} 〈 D_{I}^{(i)} D_{I}^{(j)} D_{I}^{(n - i - j)} 〉_{L} \\ + \dots + \frac{{(- β)}^{n - 1}}{n!} 〈 {(D_{I}^{(1)})}^{n} 〉_{L} . \end{matrix}

(126)

This is indeed the case and has already been implicit in the diagrammatic rules presented in Sec. IV. A proof of this assertion is given here as the linked-diagram theorem at finite temperature.

The proof consists of two parts: First, we show that ⟨D⁽ⁿ⁾⟩ is always linked by relying heavily on the linked-diagram theorem at zero temperature.^63,65–71 Second, we prove that the unlinked terms in the subsequent sums in Eq. (125) are canceled out exactly. Note that the cancellation does not complete within each parenthesis, and we have to consider the recursion holistically.

Using the ⟨E⁽ⁿ⁾⟩ recursion in Sec. II C, we can rewrite ⟨D⁽ⁿ⁾⟩ as

〈 D^{(n)} 〉 = 〈 E^{(n)} 〉 - μ^{(n)} \bar{N}

(127)

= 〈 〈 Φ_{I}^{(0)} | \hat{V} | Φ_{I}^{(n - 1)} 〉 〉 - {(μ^{(n)} \bar{N})}_{L}

(128)

\begin{matrix} = & 〈 〈 Φ_{I}^{(0)} | \hat{V} {\hat{R}}_{I} \hat{V} | Φ_{I}^{(n - 2)} 〉 〉 \\ - \sum_{i = 1}^{n - 2} 〈 \sum_{J} 〈 Φ_{I}^{(0)} | \hat{V} {\hat{R}}_{I} | Φ_{J}^{(n - 1 - i)} 〉 E_{J I}^{(i)} 〉 - {(μ^{(n)} \bar{N})}_{L}, \end{matrix}

(129)

where $μ^{(n)} \bar{N}$ is linked despite the fact that it is a product of two scalars. This is because μ⁽ⁿ⁾ is intensive.⁶⁴

As emphasized, at finite temperature, the second term (the renormalization term) consists of both unlinked and linked contributions. The unlinked contribution comes from the diagonal elements $E_{I I}^{(i)}$ ⁠, making the corresponding summand a product of two extensive scalars. Its off-diagonal elements $E_{J I}^{(i)}$ (J ≠ I) are linked to the $⟨ Φ_{I}^{(0)} | \hat{V} {\hat{R}}_{I} | Φ_{J}^{(n - 1 - i)} ⟩$ part via J, which is degenerate with but distinct from I.

Therefore, we can further rewrite the salient portion of the right-hand side of Eq. (129) as

\begin{align} 〈 〈 Φ_{I}^{(0)} | \hat{V} {\hat{R}}_{I} \hat{V} | Φ_{I}^{(n - 2)} 〉 〉 - \sum_{i = 1}^{n - 2} 〈 \sum_{J} 〈 Φ_{I}^{(0)} | \hat{V} {\hat{R}}_{I} | Φ_{J}^{(n - 1 - i)} 〉 E_{J I}^{(i)} 〉 \\ = 〈 〈 Φ_{I}^{(0)} | \hat{V} {\hat{R}}_{I} \hat{V} | Φ_{I}^{(n - 2)} 〉 〉_{L} - \sum_{i = 1}^{n - 2} 〈 \sum_{J} 〈 Φ_{I}^{(0)} | \hat{V} {\hat{R}}_{I} | Φ_{J}^{(n - 1 - i)} 〉 E_{J I}^{(i)} 〉_{L} \\ + 〈 〈 Φ_{I}^{(0)} | \hat{V} {\hat{R}}_{I} \hat{V} | Φ_{I}^{(n - 2)} 〉 〉_{U} - \sum_{i = 1}^{n - 2} 〈 \sum_{J} 〈 Φ_{I}^{(0)} | \hat{V} {\hat{R}}_{I} | Φ_{J}^{(n - 1 - i)} 〉 E_{J I}^{(i)} 〉_{U} \\ = 〈 〈 Φ_{I}^{(0)} | \hat{V} {\hat{R}}_{I} \hat{V} | Φ_{I}^{(n - 2)} 〉 〉_{L} - \sum_{i = 1}^{n - 2} 〈 \sum_{J \neq I} 〈 Φ_{I}^{(0)} | \hat{V} {\hat{R}}_{I} | Φ_{J}^{(n - 1 - i)} 〉 E_{J I}^{(i)} 〉_{L} \\ + 〈 〈 Φ_{I}^{(0)} | \hat{V} {\hat{R}}_{I} \hat{V} | Φ_{I}^{(n - 2)} 〉 〉_{U} - \sum_{i = 1}^{n - 2} 〈 〈 Φ_{I}^{(0)} | \hat{V} {\hat{R}}_{I} | Φ_{I}^{(n - 1 - i)} 〉 E_{I I}^{(i)} 〉_{U}, \end{align}

(130)

where subscript U means unlinked. The identical line of logic proving the linked-diagram theorem at zero temperature can be used to show that the last two terms cancel each other out, leaving only the linked contributions, i.e.,

\begin{matrix} 〈 〈 Φ_{I}^{(0)} | \hat{V} {\hat{R}}_{I} \hat{V} | Φ_{I}^{(n - 2)} 〉 〉 - \sum_{i = 1}^{n - 2} 〈 \sum_{J} 〈 Φ_{I}^{(0)} | \hat{V} {\hat{R}}_{I} | Φ_{J}^{(n - 1 - i)} 〉 E_{J I}^{(i)} 〉 \\ = 〈 〈 Φ_{I}^{(0)} | \hat{V} {\hat{R}}_{I} \hat{V} | Φ_{I}^{(n - 2)} 〉 〉_{L} \\ - \sum_{i = 1}^{n - 2} 〈 \sum_{J \neq I} 〈 Φ_{I}^{(0)} | \hat{V} {\hat{R}}_{I} | Φ_{J}^{(n - 1 - i)} 〉 E_{J I}^{(i)} 〉_{L} . \end{matrix}

(131)

The reader is referred to Chap. 6 of Ref. 63 or Ref. 70 for the time-independent proof of the linked-diagram theorem of the zero-temperature MBPT, which can be easily seen to be applicable to the thermal averages with no modification.

Next, we prove that the second and subsequent terms in Eq. (125) are linked by virtue of a systematic cancellation of all unlinked terms between $⟨ D_{I}^{(i_{1})} \dots D_{I}^{(i_{k})} ⟩$ (2 ≤ k ≤ m) and the products of Ω⁽ⁱ⁾’s. We borrow equations in Appendix B only for notational clarity.

Let us define A⁽ⁿ⁾ and B⁽ⁿ⁾ by

A^{(n)} = (- β) Ω^{(n)},

(132)

\begin{matrix} B^{(n)} & = & (- β) 〈 D_{I}^{(n)} 〉 + \frac{{(- β)}^{2}}{2!} \sum_{i = 1}^{n - 1} 〈 D_{I}^{(i)} D_{I}^{(n - i)} 〉 \\ + \frac{{(- β)}^{3}}{3!} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} 〈 D_{I}^{(i)} D_{I}^{(j)} D_{I}^{(n - i - j)} 〉 \\ + \dots + \frac{{(- β)}^{n}}{n!} 〈 {(D_{I}^{(1)})}^{n} 〉 . \end{matrix}

(133)

Equation (125) can then be identified as Eq. (B2) with these definitions of A⁽ⁿ⁾ and B⁽ⁿ⁾. Furthermore, if we divide B⁽ⁿ⁾ into the linked and unlinked parts,

B^{(n)} = {(B^{(n)})}_{L} + {(B^{(n)})}_{U},

(134)

the unlinked part is a sum of products of the linked part, having the following exponential structure:⁷⁷

\begin{matrix} {(B^{(n)})}_{U} & = & \frac{1}{2!} \sum_{i = 1}^{n - 1} {(B^{(i)})}_{L} {(B^{(n - i)})}_{L} \\ + \frac{1}{3!} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} {(B^{(i)})}_{L} {(B^{(j)})}_{L} {(B^{(n - i - j)})}_{L} \\ + \dots + \frac{1}{n!} {(B^{(1)})}_{L}^{n} . \end{matrix}

(135)

This means

\begin{matrix} B^{(n)} & = & {(B^{(n)})}_{L} + \frac{1}{2!} \sum_{i = 1}^{n - 1} {(B^{(i)})}_{L} {(B^{(n - i)})}_{L} \\ + \frac{1}{3!} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} {(B^{(i)})}_{L} {(B^{(j)})}_{L} {(B^{(n - i - j)})}_{L} + \dots + \frac{1}{n!} {(B^{(1)})}_{L}^{n} . \end{matrix}

(136)

Comparing this with Eq. (B3), we immediately surmise $A^{(n)} = {(B^{(n)})}_{L}$ ⁠, proving the linkedness of Ω⁽ⁿ⁾.

VI. NUMERICAL BENCHMARKS

A. General-order algorithms

The recursion relationships for the grand canonical ensemble in Sec. II and those for the canonical ensemble in Appendix A have been implemented in a determinant-based FCI program. The algorithms are literal translations of those equations. They enable calculations of μ⁽ⁿ⁾, Ω⁽ⁿ⁾, and U⁽ⁿ⁾ in the grand canonical ensemble and F⁽ⁿ⁾ and U⁽ⁿ⁾ in the canonical ensemble for an ideal gas of the smallest molecules at any arbitrary perturbation order, while S⁽ⁿ⁾ can be inferred from them. These calculations are more expensive than FCI itself and are not intended for realistic applications. Instead, they are meant to furnish benchmark data^{53,59,78–83} to assess the validity of the theories and the correctness of more efficient algorithms^84–88 based on reduced analytical formulas to be developed in the future. The latter algorithms should have the same cost scaling as their respective zero-temperature counterparts and be applicable to solids in the grand canonical ensemble.^89,90 Practical utility of the canonical ensemble, which does not lend itself to reduced analytical formulas, is doubtful at this point.

We have also implemented the reduced analytical formulas of Ω⁽ⁿ⁾, μ⁽ⁿ⁾, U⁽ⁿ⁾, and S⁽ⁿ⁾ in the grand canonical ensemble in the range of 0 ≤ n ≤ 3. We already implemented the analytical formulas for n ≤ 2 in our previous studies.^54,55 That of Ω⁽³⁾ derived in Sec. IV C has been manually implemented. The analytical formulas for μ⁽³⁾, U⁽³⁾, and S⁽³⁾ have then been directly computer-programmed by systematic modifications of the Ω⁽³⁾ code. The derivations at higher n would be extremely tedious and need to be computerized.⁷⁶

We have repeated the λ-variation calculations in both ensembles up to a higher order than previously reported.^53,62 A λ-variation calculation^53,59 determines the several lowest-order perturbation corrections by numerical differentiation with respect to λ of the thermodynamic functions obtained by the finite-temperature FCI (Ref. 57) with a perturbation-scaled Hamiltonian, $\hat{H} = {\hat{H}}_{0} + λ \hat{V}$ ⁠. The lower-order data already reported in our previous studies are reproduced by the present calculations. The new data at higher orders are not necessarily precise because of the round-off and finite-difference errors, which have been detected by the disagreements among different finite-difference formulas (not shown).

B. Grand canonical ensemble

Tables III–V compile the thermodynamic functions in the grand canonical ensemble at T = 10⁵, 10⁶, and 10⁷ K, respectively, of an ideal gas of the identical hydrogen fluoride molecules (0.9168 Å) in the STO-3G basis set.^53,57,62 The gas consists of infinitely many nonrotating rigid hydrogen fluoride molecules that are not interacting with one another, except to exchange energies and electrons. The system and temperatures considered are unrealistic, only serving as a convenient test case for formalisms and computer codes, although a similar atomic calculation might be relevant to the Saha ionization equation.^91,92 The zero-temperature canonical HF wave function for the singlet ground state has been used as the reference.

TABLE III.

Perturbation corrections to the thermodynamic functions of an ideal gas of the hydrogen fluoride molecules (0.9168 Å in STO-3G) in the grand canonical ensemble at T = 10⁵ K.

	Ω⁽ⁿ⁾/E_h			μ⁽ⁿ⁾/E_h			U⁽ⁿ⁾/E_h
n	Recursion^a	Analytical^b	λ-variation^c	Recursion^a	Analytical^b	λ-variation^c	Recursion^a	Analytical^b	λ-variation^c
0	−55.636 56	−55.636 56	−55.636 56	0.272 24	0.272 24	0.272 24	−52.016 59	−52.016 59	−52.016 59
1	−45.268 43	−45.268 43	−45.268 43	−0.075 19	−0.075 19	−0.075 19	−45.947 86	−45.947 86	−45.947 86
2	−2.581 48	−2.581 48	−2.581 47	0.231 98	0.231 98	0.231 98	0.098 41	0.098 41	0.098 41
3	4.413 31	4.413 31	4.413 56	−0.421 77	−0.421 77	−0.421 80	−0.146 04	−0.146 04	−0.146 02
4	−9.729 34		−9.738 64	0.927 40		0.928 33	−0.171 27		−0.171 96
5	22.276 04		22.180 99	−2.130 89		−2.120 63	0.723 67		0.719 65
6	−53.385 26			5.116 03			−2.949 46
7	129.978 20			−12.461 29			10.603 52
8	−311.271 00			29.812 98			−38.
9	713.882 81			−67.969 33			121.
10	−1661.6			147.877 02			3888.
$\sum_{0}^{10}$	−1268.9			101.179 18			3881.
FCI^d	−102.106 59			0.295 68			−98.049 38

	Ω⁽ⁿ⁾/E_h			μ⁽ⁿ⁾/E_h			U⁽ⁿ⁾/E_h
n	Recursion^a	Analytical^b	λ-variation^c	Recursion^a	Analytical^b	λ-variation^c	Recursion^a	Analytical^b	λ-variation^c
0	−55.636 56	−55.636 56	−55.636 56	0.272 24	0.272 24	0.272 24	−52.016 59	−52.016 59	−52.016 59
1	−45.268 43	−45.268 43	−45.268 43	−0.075 19	−0.075 19	−0.075 19	−45.947 86	−45.947 86	−45.947 86
2	−2.581 48	−2.581 48	−2.581 47	0.231 98	0.231 98	0.231 98	0.098 41	0.098 41	0.098 41
3	4.413 31	4.413 31	4.413 56	−0.421 77	−0.421 77	−0.421 80	−0.146 04	−0.146 04	−0.146 02
4	−9.729 34		−9.738 64	0.927 40		0.928 33	−0.171 27		−0.171 96
5	22.276 04		22.180 99	−2.130 89		−2.120 63	0.723 67		0.719 65
6	−53.385 26			5.116 03			−2.949 46
7	129.978 20			−12.461 29			10.603 52
8	−311.271 00			29.812 98			−38.
9	713.882 81			−67.969 33			121.
10	−1661.6			147.877 02			3888.
$\sum_{0}^{10}$	−1268.9			101.179 18			3881.
FCI^d	−102.106 59			0.295 68			−98.049 38

^a

MBPT(n) recursions. The entries with fewer significant figures suffer from large round-off errors.

^b

MBPT(n) reduced analytical formulas.

^c

Central seven-point formula with Δλ = 0.01. See also Ref. 53.

^d

Finite-temperature full configuration interaction. See also Ref. 57.

TABLE IV.

Same as Table III but at T = 10⁶ K.

	Ω⁽ⁿ⁾/E_h			μ⁽ⁿ⁾/E_h			U⁽ⁿ⁾/E_h
n	Recursion^a	Analytical^b	λ-variation^c	Recursion^a	Analytical^b	λ-variation^c	Recursion^a	Analytical^b	λ-variation^c
0	−105.947 53	−105.947 53	−105.947 53	3.961 30	3.961 30	3.961 30	−50.596 35	−50.596 35	−50.596 35
1	−44.525 64	−44.525 64	−44.525 64	−0.168 96	−0.168 96	−0.168 96	−46.176 65	−46.176 65	−46.176 65
2	−0.964 31	−0.964 31	−0.964 25	0.085 09	0.085 09	0.085 09	−0.219 84	−0.219 84	−0.219 83
3	0.249 39	0.249 39	0.257 74	−0.022 70	−0.022 70	−0.023 53	0.064 64	0.064 64	0.066 48
4	−0.073 81		−0.452 86	0.006 76		0.044 67	−0.023 89		−0.107 20
5	0.022 96		−15.641 35	−0.002 10		1.564 33	0.009 45		−3.442 42
6	−0.006 99			0.000 63			−0.003 73
7	0.001 87			−0.000 17			0.001 40
8	−0.000 32			0.000 03			−0.000 48
9	−0.000 05			0.000 01			0.000 14
10	0.000 09			−0.000 01			−0.000 02
$\sum_{0}^{10}$	−151.244 36			3.859 89			−96.945 33
FCI^d	−151.244 40			3.859 90			−96.945 34

	Ω⁽ⁿ⁾/E_h			μ⁽ⁿ⁾/E_h			U⁽ⁿ⁾/E_h
n	Recursion^a	Analytical^b	λ-variation^c	Recursion^a	Analytical^b	λ-variation^c	Recursion^a	Analytical^b	λ-variation^c
0	−105.947 53	−105.947 53	−105.947 53	3.961 30	3.961 30	3.961 30	−50.596 35	−50.596 35	−50.596 35
1	−44.525 64	−44.525 64	−44.525 64	−0.168 96	−0.168 96	−0.168 96	−46.176 65	−46.176 65	−46.176 65
2	−0.964 31	−0.964 31	−0.964 25	0.085 09	0.085 09	0.085 09	−0.219 84	−0.219 84	−0.219 83
3	0.249 39	0.249 39	0.257 74	−0.022 70	−0.022 70	−0.023 53	0.064 64	0.064 64	0.066 48
4	−0.073 81		−0.452 86	0.006 76		0.044 67	−0.023 89		−0.107 20
5	0.022 96		−15.641 35	−0.002 10		1.564 33	0.009 45		−3.442 42
6	−0.006 99			0.000 63			−0.003 73
7	0.001 87			−0.000 17			0.001 40
8	−0.000 32			0.000 03			−0.000 48
9	−0.000 05			0.000 01			0.000 14
10	0.000 09			−0.000 01			−0.000 02
$\sum_{0}^{10}$	−151.244 36			3.859 89			−96.945 33
FCI^d	−151.244 40			3.859 90			−96.945 34

^a

MBPT(n) recursions.

^b

MBPT(n) reduced analytical formulas.

^c

Central seven-point formula with Δλ = 0.01. See also Ref. 53.

^d

Finite-temperature full configuration interaction. See also Ref. 57.

TABLE V.

Same as Table III but at T = 10⁷ K.

	Ω⁽ⁿ⁾/E_h			μ⁽ⁿ⁾/E_h			U⁽ⁿ⁾/E_h
n	Recursion^a	Analytical^b	λ-variation^c	Recursion^a	Analytical^b	λ-variation^c	Recursion^a	Analytical^b	λ-variation^c
0	−686.708 14	−686.708 14	−686.708 14	47.150 12	47.150 12	47.150 12	−45.789 11	−45.789 11	−45.789 11
1	−43.199 11	−43.199 11	−43.199 11	−0.298 11	−0.298 11	−0.298 11	−46.235 54	−46.235 54	−46.235 54
2	−0.196 96	−0.196 96	−0.196 91	0.017 74	0.017 74	0.017 74	−0.032 60	−0.032 60	−0.032 58
3	0.009 51	0.009 51	0.053 65	−0.000 88	−0.000 88	−0.005 29	0.001 79	0.001 79	0.013 33
4	−0.000 53		1.457 17	0.000 05		−0.145 72	−0.000 13		0.380 65
5	0.000 03		−44.483 27	−0.000 00		4.448 33	0.000 01		−11.631 01
6	−0.000 00			0.000 00			−0.000 00
7	0.000 00			−0.000 00			0.000 00
8	−0.000 00			0.000 00			−0.000 00
9	0.000 00			−0.000 00			0.000 00
10	−0.000 00			0.000 00			−0.000 00
$\sum_{0}^{10}$	−730.095 19			46.868 92			−92.055 57
FCI^d	−730.095 19			46.868 92			−92.055 57

	Ω⁽ⁿ⁾/E_h			μ⁽ⁿ⁾/E_h			U⁽ⁿ⁾/E_h
n	Recursion^a	Analytical^b	λ-variation^c	Recursion^a	Analytical^b	λ-variation^c	Recursion^a	Analytical^b	λ-variation^c
0	−686.708 14	−686.708 14	−686.708 14	47.150 12	47.150 12	47.150 12	−45.789 11	−45.789 11	−45.789 11
1	−43.199 11	−43.199 11	−43.199 11	−0.298 11	−0.298 11	−0.298 11	−46.235 54	−46.235 54	−46.235 54
2	−0.196 96	−0.196 96	−0.196 91	0.017 74	0.017 74	0.017 74	−0.032 60	−0.032 60	−0.032 58
3	0.009 51	0.009 51	0.053 65	−0.000 88	−0.000 88	−0.005 29	0.001 79	0.001 79	0.013 33
4	−0.000 53		1.457 17	0.000 05		−0.145 72	−0.000 13		0.380 65
5	0.000 03		−44.483 27	−0.000 00		4.448 33	0.000 01		−11.631 01
6	−0.000 00			0.000 00			−0.000 00
7	0.000 00			−0.000 00			0.000 00
8	−0.000 00			0.000 00			−0.000 00
9	0.000 00			−0.000 00			0.000 00
10	−0.000 00			0.000 00			−0.000 00
$\sum_{0}^{10}$	−730.095 19			46.868 92			−92.055 57
FCI^d	−730.095 19			46.868 92			−92.055 57

^a

MBPT(n) recursions.

^b

MBPT(n) reduced analytical formulas.

^c

Central seven-point formula with Δλ = 0.01. See also Ref. 53.

^d

Finite-temperature full configuration interaction. See also Ref. 57.

In all cases, Ω⁽ⁿ⁾, μ⁽ⁿ⁾, and U⁽ⁿ⁾ computed by the sum-over-states (recursion) and sum-over-orbitals (reduced analytical) formulas agree with each other for at least ten decimal places, verifying the formulas and computer programs mutually.

At T = 10⁵ K, they display a clear sign of divergence. This may not be surprising in light of the fact that HCPT for many states are already divergent.⁵⁹ This is a different type of divergence than discovered by Kohn and Luttinger⁴⁸ and analyzed by us.⁵⁶ The latter is concerned with the zero-temperature limit, at which the radius of convergence is zero under some circumstances,⁵⁶ whereas the former should have a finite radius of convergence, the rate of which may be accelerated by, e.g., the Padé approximant.^79,93–95 At lower temperatures (not shown), U⁽ⁿ⁾ reproduces E⁽ⁿ⁾ for the reference state computed by the zero-temperature MBPT and is much less prone to divergence (see, however, Ref. 96). At T = 10⁶ and 10⁷ K, the thermodynamic functions are rapidly convergent at the respective finite-temperature FCI values. At 10⁷ K, the tenth-order approximations are within 10⁻⁷ E_h of FCI, testifying that the finite-temperature MBPT presented here is a converging series of approximation toward exactness. This would not be the case if the renormalization diagrams were overlooked or neglected.

At T = 10⁵ K, the λ-variation method turns out to be surprisingly accurate. It remains somewhat useful up to the fifth order. At T = 10⁶ and 10⁷ K, it breaks down at lower orders, hardly serving as a benchmark even for the second-order corrections. The stability of the λ-variation method may be anticorrelated with the convergence of the perturbation series. At lower temperatures (not shown), the finite-temperature FCI (and hence the λ-variation also) becomes unstable because of the inherent difficulty of precisely determining μ. At T = 0, any value of μ falling in between the highest-occupied (HOMO) and lowest-unoccupied molecular orbital (LUMO) energies satisfies the electroneutrality condition, making it increasingly difficult for μ to reach the correct zero-temperature limit,^56,57 which is the midpoint of HOMO and LUMO. Nevertheless, for the temperatures used in these tables, analytical results and the λ-variation benchmarks are in reasonable agreement.

Figures 20 and 21 plot the perturbation approximations to U and S and their FCI values as a function of temperature. It can be seen that the perturbation approximations are so accurate that they are indiscernible from the FCI curves at both low and high temperatures. Only at mid-temperatures (10⁴ ≤ T/K ≤ 10⁶) are the differences noticeable. This is the same temperature domain where these thermodynamic functions exhibit strong temperature dependence and the perturbation series tend to diverge. Still, low-order perturbation theories seem to work well, correcting the vast majority of the discrepancies between the zeroth-order (Fermi–Dirac) theory and FCI in U (U⁽⁰⁾ is not shown in Fig. 20 because it is far outside the graph). For S, the first-order approximation overcorrects the errors in S⁽⁰⁾ in this temperature domain, while the second- and third-order approximations roughly halve the errors. Given the divergence of the underlying perturbation series at these temperatures, the lower-order perturbation theories are promising.

FIG. 20.

View large Download slide

U⁽⁰⁾ + ⋯ + U⁽ⁿ⁾ (1 ≤ n ≤ 3) as a function of temperature for the same system as Table III.

FIG. 21.

View large Download slide

S⁽⁰⁾ + ⋯ + S⁽ⁿ⁾ (0 ≤ n ≤ 3) as a function of temperature for the same system as Table III.

C. Canonical ensemble

Thermodynamic functions in the canonical ensemble have been calculated for the same system as Sec. VI B. The results at T = 10⁵, 10⁶, and 10⁷ K are summarized in Tables VI and VII, and VIII, respectively. Only the sum-over-states (recursion) formulas and λ-variation method are available for this ensemble.

TABLE VI.

Perturbation corrections to the thermodynamic functions of an ideal gas of the hydrogen fluoride molecules (0.9168 Å in STO-3G) in the canonical ensemble at T = 10⁵ K.

	F⁽ⁿ⁾/E_h		U⁽ⁿ⁾/E_h
n	Recursion^a	λ-variation^b	Recursion^a	λ-variation^b
0	−52.671 66	−52.671 66	−52.264 52	−52.264 52
1	−46.163 15	−46.163 15	−45.694 42	−45.694 42
2	−0.146 57	−0.146 57	−0.021 51	−0.021 51
3	−0.052 39	−0.052 39	−0.166 48	−0.166 48
4	0.001 62	0.001 62	−0.088 26	−0.088 26
5	0.011 17	0.011 21	0.019 89	0.020 29
6	0.003 59		0.056 64
7	−0.003 46		−0.5
8	0.045 90		16.
9	−0.390 63		−562.
10	−36.6		22 464.
$\sum_{0}^{10}$	−136.0		21 820.
FCI^c	−99.020 43		−98.178 36

	F⁽ⁿ⁾/E_h		U⁽ⁿ⁾/E_h
n	Recursion^a	λ-variation^b	Recursion^a	λ-variation^b
0	−52.671 66	−52.671 66	−52.264 52	−52.264 52
1	−46.163 15	−46.163 15	−45.694 42	−45.694 42
2	−0.146 57	−0.146 57	−0.021 51	−0.021 51
3	−0.052 39	−0.052 39	−0.166 48	−0.166 48
4	0.001 62	0.001 62	−0.088 26	−0.088 26
5	0.011 17	0.011 21	0.019 89	0.020 29
6	0.003 59		0.056 64
7	−0.003 46		−0.5
8	0.045 90		16.
9	−0.390 63		−562.
10	−36.6		22 464.
$\sum_{0}^{10}$	−136.0		21 820.
FCI^c	−99.020 43		−98.178 36

^a

MBPT(n) recursions. The entries with fewer significant figures suffer from large round-off errors.

^b

Central seven-point formula with Δλ = 0.01. See also Ref. 62.

^c

Finite-temperature full configuration interaction. See also Ref. 57.

TABLE VII.

Same as Table VI but at T = 10⁶ K.

	F⁽ⁿ⁾/E_h		U⁽ⁿ⁾/E_h
n	Recursion^a	λ-variation^b	Recursion^a	λ-variation^b
0	−62.555 50	−62.555 50	−50.622 84	−50.622 84
1	−46.778 59	−46.778 59	−46.716 63	−46.716 63
2	−0.016 47	−0.016 47	−0.034 20	−0.034 20
3	0.000 30	0.000 30	0.000 90	0.000 90
4	0.000 00	0.000 00	0.000 00	0.000 00
5	−0.000 00	0.000 01	−0.000 00	0.000 02
6	0.000 00		0.000 00
7	0.000 00		−0.000 00
8	−0.000 00		0.000 00
9	−0.000 00		−0.000 00
10	0.000 00		−0.000 00
$\sum_{0}^{10}$	−109.350 26		−97.372 78
FCI^c	−109.350 26		−97.372 78

	F⁽ⁿ⁾/E_h		U⁽ⁿ⁾/E_h
n	Recursion^a	λ-variation^b	Recursion^a	λ-variation^b
0	−62.555 50	−62.555 50	−50.622 84	−50.622 84
1	−46.778 59	−46.778 59	−46.716 63	−46.716 63
2	−0.016 47	−0.016 47	−0.034 20	−0.034 20
3	0.000 30	0.000 30	0.000 90	0.000 90
4	0.000 00	0.000 00	0.000 00	0.000 00
5	−0.000 00	0.000 01	−0.000 00	0.000 02
6	0.000 00		0.000 00
7	0.000 00		−0.000 00
8	−0.000 00		0.000 00
9	−0.000 00		−0.000 00
10	0.000 00		−0.000 00
$\sum_{0}^{10}$	−109.350 26		−97.372 78
FCI^c	−109.350 26		−97.372 78

^a

MBPT(n) recursions.

^b

Central seven-point formula with Δλ = 0.01. See also Ref. 62.

^c

Finite-temperature full configuration interaction. See also Ref. 57.

TABLE VIII.

Same as Table VI but at T = 10⁷ K.

	F⁽ⁿ⁾/E_h		U⁽ⁿ⁾/E_h
n	Recursion^a	λ-variation^b	Recursion^a	λ-variation^b
0	−176.803 50	−176.803 50	−46.002 78	−46.002 78
1	−46.857 43	−46.857 43	−46.845 16	−46.845 16
2	−0.002 43	−0.002 43	−0.003 71	−0.003 71
3	0.000 03	0.000 03	0.000 05	0.000 05
4	−0.000 00	−0.000 00	−0.000 00	−0.000 00
5	0.000 00	0.000 01	0.000 00	0.000 01
6	−0.000 00		−0.000 00
7	0.000 00		0.000 00
8	0.000 00		0.000 00
9	0.000 00		0.000 00
10	0.000 00		0.000 00
$\sum_{0}^{10}$	−223.663 34		−92.851 59
FCI^c	−223.663 34		−92.851 59

	F⁽ⁿ⁾/E_h		U⁽ⁿ⁾/E_h
n	Recursion^a	λ-variation^b	Recursion^a	λ-variation^b
0	−176.803 50	−176.803 50	−46.002 78	−46.002 78
1	−46.857 43	−46.857 43	−46.845 16	−46.845 16
2	−0.002 43	−0.002 43	−0.003 71	−0.003 71
3	0.000 03	0.000 03	0.000 05	0.000 05
4	−0.000 00	−0.000 00	−0.000 00	−0.000 00
5	0.000 00	0.000 01	0.000 00	0.000 01
6	−0.000 00		−0.000 00
7	0.000 00		0.000 00
8	0.000 00		0.000 00
9	0.000 00		0.000 00
10	0.000 00		0.000 00
$\sum_{0}^{10}$	−223.663 34		−92.851 59
FCI^c	−223.663 34		−92.851 59

^a

MBPT(n) recursions.

^b

Central seven-point formula with Δλ = 0.01. See also Ref. 62.

^c

Finite-temperature full configuration interaction. See also Ref. 57.

Owing to the much fewer number of states involved in the canonical ensemble, the thermodynamic functions are less prone to divergence and the λ-variation method tends to be more accurate. At T = 10⁶ and 10⁷ K, the tenth-order perturbation theory agrees with FCI for five and ten decimal places, respectively, which attests to the convergence of the perturbation series to exactness at any temperature unless they diverge.

VII. CONCLUSIONS

We have fully developed the finite-temperature MBPT that expands all thermodynamic functions in uniform perturbation series. Both grand canonical and canonical ensembles are considered, although only the former lends itself to more drastic mathematical simplifications, making it useful for condensed-matter applications. The Rayleigh–Schrödinger-like algebraic recursions have been obtained for both, yielding their sum-over-states analytical formulas of all thermodynamic functions at any arbitrary order.

For the grand canonical ensemble, the sum-over-states formulas have been transformed to sum-over-orbitals analytical formulas by the method of normal-ordered second quantization at finite temperature. The rules of Wick contractions and various operators in the normal-ordered form have been derived (as opposed to postulated) and how they are practically applied to integral and thermal-average evaluations has been illustrated. Feynman diagrammatic rules are then introduced as a straightforward graphical representation of Wick contractions as lines and of molecular integrals as vertexes. Nowhere in this formulation do we resort to time-dependent logic^46–51 or human intuitions to exhaustively enumerate all valid diagrams, which include renormalization, anomalous, and anomalous–renormalization diagrams.

Normal ordering regroups the Hamiltonian operator into the finite-temperature HF energy, finite-temperature Fock matrix elements times one-electron operators, and two-electron integrals times two-electron operators. This natural partitioning underscores the significance of the finite-temperature HF theory as the foundation of all converging finite-temperature electron–correlation theories. Its energy expression is identified with $⟨ E_{I}^{(0)} ⟩ + ⟨ E_{I}^{(1)} ⟩$ ⁠, which differs from either U⁽⁰⁾ + U⁽¹⁾ or Ω⁽⁰⁾ + Ω⁽¹⁾ at T > 0.

In the nondegenerate zero-temperature MBPT, the renormalization terms cancel exactly the unlinked contributions of the equal magnitudes in the parent terms, ensuring the size-consistency of MBPT at any order, as proven by the linked-diagram theorem.^65–70 In the finite-temperature MBPT, which sums over all possible states many of which are degenerate, the renormalization terms consist of both linked and unlinked contributions. The unlinked contributions cancel the same in the parent terms, restoring the linkedness of the perturbation corrections at any order, while nonzero linked renormalization terms persist. They correspond to diagrams whose resolvent lines are displaced upward in all possible ways, in which there are two or more resolvents in between some pairs of adjacent vertexes and none in between others. These renormalization diagrams are distinguished from the well-known anomalous diagrams whose resolvent lines are deleted. They both are essential numerically, and had it not for them, the perturbation series would not converge at the finite-temperature FCI limit. It is unclear whether the conventional time-dependent, diagrammatic formulation^46–51 or density-matrix derivation⁵² is aware of such diagrams.

We have proven the linked-diagram theorem of the finite-temperature MBPT. The theorem is based on the linked-diagram theorem at zero temperature^65–70 and the systematic cancellation of unlinked anomalous diagrams.

We have derived reduced analytical formulas for the perturbation corrections to the thermodynamic functions in the grand canonical ensemble up to the third order. We have also implemented general-order algorithms in both ensembles and calculated the benchmark perturbation corrections up to the tenth order. At intermediate temperatures, where thermodynamic functions vary with temperature considerably, the perturbation series tend to diverge, but elsewhere, they converge rapidly toward the finite-temperature FCI limits, establishing the correctness of the theory.

In a separate study,⁵⁶ we showed that the finite-temperature MBPT at the first and second orders does not converge at the correct zero-temperature limit when the reference wave function is qualitatively different from the exact one, confirming Kohn and Luttinger.⁴⁸ This originates from the nonanalyticity of the Boltzmann factor and is therefore expected to plague most any finite-temperature MBPT, e.g., the one in the canonical ensemble also.⁶² Nonconvergence should persist at the third and higher orders under the same condition.

This study and another reported in Ref. 59 have been conducted to lay a firm mathematical foundation of perturbation theories,^97,98 which continue to be a workhorse for efficient, size-consistent,^64,99 converging ab initio electronic structure calculations^100,111 for large complex molecules and solids.¹¹² Unlike time-dependent, diagrammatic expositions,^42–45 which are often incomplete and/or intractable, our derivations follow the transparent, linear, time-independent, algebraic logic, involving only elementary calculus and combinatorial identities. They justify the second-quantized and diagrammatic rules just as tools for expediency rather than relying on them as a graphical method of derivations that are often based on human intuitions detached from mathematics. Our derivation strategy can therefore be universally applied to any perturbation theory insofar as the exact limit is known, reliably leading to recursions, order-by-order analytical formulas, and general-order algorithms.

ACKNOWLEDGMENTS

The author thanks Professor Rodney J. Bartlett for teaching him the quantum-field-theoretical techniques used in this study. This work was supported by the Center for Scalable, Predictive methods for Excitation and Correlated phenomena (SPEC), which is funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, as a part of the Computational Chemical Sciences Program, and also by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Grant No. DE-SC0006028.

DATA AVAILABILITY

The data that support the findings of this study are available within the article.

APPENDIX A: RECURSIONS IN THE CANONICAL ENSEMBLE

Thermodynamic functions in the canonical ensemble considered are the canonical partition function (Z), Helmholtz energy (F), internal energy (U), and entropy (S), which are related to one another by

Z = \sum_{I} e^{- β E_{I}},

(A1)

F = - \frac{1}{β} \ln Z,

(A2)

U = - \frac{\partial}{\partial β} \ln Z,

(A3)

and

F = U - T S,

(A4)

where the I-summation goes over all spin states with a fixed total number (⁠ $\bar{N}$ ⁠) of electrons. The recursion for Z⁽ⁿ⁾ is given by

\begin{matrix} \frac{Z^{(n)}}{Z^{(0)}} & = & (- β) 〈 E^{(n)} 〉 + \frac{{(- β)}^{2}}{2!} \sum_{i = 1}^{n - 1} 〈 E^{(i)} E^{(n - i)} 〉 \\ + \frac{{(- β)}^{3}}{3!} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} 〈 E^{(i)} E^{(j)} E^{(n - i - j)} 〉 \\ + \dots + \frac{{(- β)}^{n}}{n!} 〈 {(E^{(1)})}^{n} 〉, \end{matrix}

(A5)

where ⟨⋯⟩ is defined differently from Eq. (12) only in this appendix as

⟨ X_{I} ⟩ \equiv \frac{\sum_{I} X_{I} e^{- β E_{I}^{(0)}}}{\sum_{I} e^{- β E_{I}^{(0)}}} .

(A6)

The recursion for F⁽ⁿ⁾ is inferred from Eq. (15) as

\begin{matrix} F^{(n)} & = & 〈 E^{(n)} 〉 + \frac{(- β)}{2!} \sum_{i = 1}^{n - 1} (〈 E^{(i)} E^{(n - i)} 〉 - F^{(i)} F^{(n - i)}) \\ + \frac{{(- β)}^{2}}{3!} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} (〈 E^{(i)} E^{(j)} E^{(n - i - j)} 〉) - (F^{(i)} F^{(j)} F^{(n - i - j)}) \\ + \dots + \frac{{(- β)}^{n - 1}}{n!} \{〈 {(E^{(1)})}^{n} 〉 - {(F^{(1)})}^{n}\} . \end{matrix}

(A7)

Likewise, U⁽ⁿ⁾ is defined recursively by

\begin{matrix} U^{(n)} & = & 〈 E_{I}^{(n)} 〉 + (- β) \sum_{i = 1}^{n} 〈 E_{I}^{(i)} E_{I}^{(n - i)} 〉 - (- β) \sum_{i = 1}^{n} F^{(i)} U^{(n - i)} \\ + \frac{{(- β)}^{2}}{2!} \sum_{i = 1}^{n - 1} \sum_{j = 1}^{n - i} 〈 E_{I}^{(i)} E_{I}^{(j)} E_{I}^{(n - i - j)} 〉 \\ - \frac{{(- β)}^{2}}{2!} \sum_{i = 1}^{n - 1} \sum_{j = 1}^{n - i} F^{(i)} F^{(j)} U^{(n - i - j)} \\ + \dots + \frac{{(- β)}^{n}}{n!} 〈 {(E_{I}^{(1)})}^{n} E_{I}^{(0)} 〉 - \frac{{(- β)}^{n}}{n!} {(F^{(1)})}^{n} U^{(0)}, \end{matrix}

(A8)

which implies the recursion for S⁽ⁿ⁾ = k_Bβ(U⁽ⁿ⁾ − F⁽ⁿ⁾).

These have been obtained by restricting the summations to $\bar{N}$ -electron states and setting μ⁽ⁿ⁾ = 0 for all n in the corresponding recursions for the grand canonical ensemble.

APPENDIX B: A PROOF OF EQUIVALENCE OF EQS. (13) AND (14)

Generally, we prove the equivalence of two recursions,

\begin{matrix} A^{(n)} & = & B^{(n)} - \frac{1}{2} \sum_{i = 1}^{n - 1} B^{(i)} B^{(n - i)} + \frac{1}{3} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} B^{(i)} B^{(j)} B^{(n - i - j)} \\ - \dots - \frac{{(- 1)}^{n}}{n} {(B^{(1)})}^{n} \end{matrix}

(B1)

and

\begin{matrix} A^{(n)} & = & B^{(n)} - \frac{1}{2!} \sum_{i = 1}^{n - 1} A^{(i)} A^{(n - i)} - \frac{1}{3!} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} A^{(i)} A^{(j)} A^{(n - i - j)} \\ - \dots - \frac{1}{n!} {(A^{(1)})}^{n} . \end{matrix}

(B2)

Equations (13) and (14) are a special case of these in which A⁽ⁿ⁾ = −βΩ⁽ⁿ⁾ and B⁽ⁿ⁾ = Ξ⁽ⁿ⁾/Ξ⁽⁰⁾. We can rewrite the second equation as

\begin{matrix} B^{(n)} & = & A^{(n)} + \frac{1}{2!} \sum_{i = 1}^{n - 1} A^{(i)} A^{(n - i)} + \frac{1}{3!} \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} A^{(i)} A^{(j)} A^{(n - i - j)} \\ + \dots + \frac{1}{n!} {(A^{(1)})}^{n} . \end{matrix}

(B3)

Substituting this into the right-hand side of Eq. (B1), we can convert it into a form expressed entirely in terms of A⁽ⁱ⁾ as

\begin{matrix} A^{(n)} & = & C_{A^{(n)}} A^{(n)} + \sum_{i = 1}^{n - 1} C_{A^{(i)} A^{(n - i)}} A^{(i)} A^{(n - i)} \\ + \sum_{i = 1}^{n - 2} \sum_{j = 1}^{n - i - 1} C_{A^{(i)} A^{(j)} A^{(n - i - j)}} A^{(i)} A^{(j)} A^{(n - i - j)} \\ + \dots + C_{{(A^{(1)})}^{n}} {(A^{(1)})}^{n} . \end{matrix}

(B4)

We shall show that $C_{A^{(n)}} = 1$ and $C_{A^{(i_{1})} A^{(i_{2})} \dots A^{(i_{m})}} = 0$ (m ≥ 2).

By inspection, we find $C_{A^{(n)}} = 1$ ⁠, $C_{A^{(i)} A^{(n - i)}} = 1 / 2! - 1 / 2 = 0$ ⁠, and $C_{A^{(i)} A^{(j)} A^{(n - i - j)}} = 1 / 3! - (1 / 2) (1 / 2!) - (1 / 2) (1 / 2!) + 1 / 3 = 0$ ⁠. In general, the coefficient multiplying $A^{(i_{1})} A^{(i_{2})} \dots A^{(i_{m})}$ (m ≥ 2) is given as

\begin{matrix} C_{A^{(i_{1})} A^{(i_{2})} \dots A^{(i_{m})}} \\ = \frac{1}{m!} - \frac{1}{2} \sum_{i = 1}^{m - 1} \frac{1}{i! (m - i)!} + \frac{1}{3} \sum_{i = 1}^{m - 2} \sum_{j = 1}^{m - i - 1} \frac{1}{i! j! (m - i - j)!} \\ + \dots - \frac{{(- 1)}^{m}}{m} {(\frac{1}{1!})}^{m} \\ = \frac{1}{m!} (\frac{1!}{1} \{\begin{matrix} m \\ 1 \end{matrix}\} - \frac{2!}{2} \{\begin{matrix} m \\ 2 \end{matrix}\} + \frac{3!}{3} \{\begin{matrix} m \\ 3 \end{matrix}\}) \\ + (\dots - \frac{{(- 1)}^{m} m!}{m} \{\begin{matrix} m \\ m \end{matrix}\}), \end{matrix}

(B5)

where $\{\begin{matrix} m \\ k \end{matrix}\}$ is the Stirling number of the second kind (the number of ways m objects are partitioned into k nonempty sets).¹¹³ It is defined as

\{\begin{matrix} m \\ k \end{matrix}\} = \frac{1}{k!} \sum_{m_{1}, \dots, m_{k}} \frac{m!}{\prod_{i = 1}^{k} m_{i}!},

(B6)

where the summation goes over all natural numbers m₁, …, m_k that satisfy m₁ + ⋯ + m_k = m. There is an identity involving the Stirling number,

\sum_{k = 1}^{m} {(- 1)}^{k} (k - 1)! \{\begin{matrix} m \\ k \end{matrix}\} = 0,

(B7)

which implies $C_{A^{(i_{1})} A^{(i_{2})} \dots A^{(i_{m})}} = 0$ (m ≥ 2), proving the equivalence.

Identity (B7) can, in turn, be proven by induction using the recursion

\{\begin{matrix} n + 1 \\ k \end{matrix}\} = k \{\begin{matrix} n \\ k \end{matrix}\} + \{\begin{matrix} n \\ k - 1 \end{matrix}\} .

(B8)

Equation (B7) is true for m = 2. If it is true for m = 2, …, n, it also holds true for m = n + 1 and therefore for any m ≥ 2 because

\begin{matrix} \sum_{k = 1}^{n + 1} {(- 1)}^{k} (k - 1)! \{\begin{matrix} n + 1 \\ k \end{matrix}\} \\ = \sum_{k = 1}^{n + 1} {(- 1)}^{k} k! \{\begin{matrix} n \\ k \end{matrix}\} + \sum_{k = 1}^{n + 1} {(- 1)}^{k} (k - 1)! \{\begin{matrix} n \\ k - 1 \end{matrix}\} \\ = \sum_{k = 1}^{n + 1} {(- 1)}^{k} k! \{\begin{matrix} n \\ k \end{matrix}\} + \sum_{k = 0}^{n} {(- 1)}^{k + 1} k! \{\begin{matrix} n \\ k \end{matrix}\} \\ = \sum_{k = 1}^{n} {(- 1)}^{k} k! \{\begin{matrix} n \\ k \end{matrix}\} + \sum_{k = 1}^{n} {(- 1)}^{k + 1} k! \{\begin{matrix} n \\ k \end{matrix}\} = 0, \end{matrix}

(B9)

where we used $\{\begin{matrix} n \\ 0 \end{matrix}\} = \{\begin{matrix} n \\ n + 1 \end{matrix}\} = 0$ ⁠.

Finally, the recursion for the Stirling number [Eq. (B8)] can be justified as follows: The left-hand side counts the number of ways n + 1 objects are partitioned into k nonempty sets. The first term in the right-hand side is the number of ways in which the (n + 1)th object is absorbed into one of the existing k nonempty sets of the first n objects. The second term is the number of ways in which the (n + 1)th object forms a new (kth) set containing only it and the first n objects are partitioned into k − 1 nonempty sets.

One may view Maclaurin series of the right-hand sides of the following equations as a mnemonic for Eqs. (B1) and (B3), respectively,

A^{(n)} = \ln (B^{(n)} + 1),

(B10)

B^{(n)} = \exp (A^{(n)}) - 1,

(B11)

which, however, may not be taken literally.

APPENDIX C: NORMAL-ORDERED SECOND QUANTIZATION AT FINITE TEMPERATURE

In this appendix, we fully develop the normal-ordered second quantization at finite temperature^28,42 within the time-independent framework in a transparent and pedagogical manner. It serves as a basis of the diagrammatic derivation and the linked-diagram theorem. The reader is referred to the work of Shavitt and Bartlett⁶³ for normal-ordered second quantization at zero temperature and to the work of March, Young, and Sampanthar⁴² for an abridged exposition of its finite-temperature counterpart. The following expounds on the latter with enough details to be actually used in deriving the second-order formulas in Sec. III.

1. Normal ordering at finite temperature

A product of creation (⁠ ${\hat{p}}^{†}$ ⁠) and annihilation (⁠ $\hat{p}$ ⁠) operators in the normal order (denoted by {⋯}) is defined [cf. Eq. (9.3.2) of Ref. 42] as

{{\hat{p}}^{†} \hat{p}} \equiv f_{p}^{+} {\hat{p}}^{†} \hat{p} - f_{p}^{-} \hat{p} {\hat{p}}^{†},

(C1)

{\hat{p} {\hat{p}}^{†}} \equiv f_{p}^{-} \hat{p} {\hat{p}}^{†} - f_{p}^{+} {\hat{p}}^{†} \hat{p},

(C2)

where $f_{p}^{-}$ and $f_{p}^{+}$ are the Fermi–Dirac distribution functions given by Eqs. (43) and (44). When p ≠ q,

{{\hat{p}}^{†} \hat{q}} \equiv {\hat{p}}^{†} \hat{q},

(C3)

{\hat{p} {\hat{q}}^{†}} \equiv \hat{p} {\hat{q}}^{†} .

(C4)

A normal-ordered product of two creation operators or of two annihilation operators is also the same as the original order.

In analogy to the zero-temperature case, the normal ordering at finite temperature is designed so that the thermal average of a normal-ordered product is always zero. Using Einstein’s convention of implied summations of repeated spin orbital indices, we confirm

\begin{matrix} 〈 〈 I | x_{p q} {{\hat{p}}^{†} \hat{q}} | I 〉 〉 & \equiv & \frac{\sum_{I} 〈 I | x_{p q} {{\hat{p}}^{†} \hat{q}} | I 〉 e^{- β D_{I}^{(0)}}}{\sum_{I} e^{- β D_{I}^{(0)}}} \\ = & \frac{\sum_{I} 〈 I | x_{p p} {{\hat{p}}^{†} \hat{p}} | I 〉 e^{- β D_{I}^{(0)}}}{\sum_{I} e^{- β D_{I}^{(0)}}} \\ = & \frac{\sum_{I} 〈 I | x_{p p} f_{p}^{+} {\hat{p}}^{†} \hat{p} | I 〉 e^{- β D_{I}^{(0)}}}{\sum_{I} e^{- β D_{I}^{(0)}}} \\ - \frac{\sum_{I} 〈 I | x_{p p} f_{p}^{-} \hat{p} {\hat{p}}^{†} | I 〉 e^{- β D_{I}^{(0)}}}{\sum_{I} e^{- β D_{I}^{(0)}}} \\ = & \sum_{p} x_{p p} f_{p}^{+} f_{p}^{-} - \sum_{p} x_{p p} f_{p}^{-} f_{p}^{+} \\ = & 0, \end{matrix}

(C5)

where we used the Boltzmann sum identities I and V of Ref. 55 in the penultimate equality.

Next, we define a Wick contraction (denoted by a staple symbol) as the difference between the operator product and its normal-ordered counterpart,

(C6)

(C7)

where we used the anticommutation rules of fermion creation and annihilation operators. When p ≠ q,

(C8)

(C9)

A Wick contraction of two creation or of two annihilation operators is also zero.

The concept of the normal ordering can be generalized to any even number of operators. For instance, a product of four operators in the normal order is defined as (p ≠ q),

{{\hat{p}}^{†} {\hat{q}}^{†} \hat{q} \hat{p}} = - {{\hat{p}}^{†} {\hat{q}}^{†} \hat{p} \hat{q}} = {{\hat{p}}^{†} \hat{p} {\hat{q}}^{†} \hat{q}}

(C10)

\equiv (f_{p}^{+} {\hat{p}}^{†} \hat{p} - f_{p}^{-} \hat{p} {\hat{p}}^{†}) (f_{q}^{+} {\hat{q}}^{†} \hat{q} - f_{q}^{-} \hat{q} {\hat{q}}^{†})

(C11)

\begin{matrix} = & f_{p}^{+} f_{q}^{+} {\hat{p}}^{†} {\hat{q}}^{†} \hat{q} \hat{p} - f_{p}^{+} f_{q}^{-} {\hat{p}}^{†} \hat{q} {\hat{q}}^{†} \hat{p} - f_{p}^{-} f_{q}^{+} \hat{p} {\hat{q}}^{†} \hat{q} {\hat{p}}^{†} \\ + f_{p}^{-} f_{q}^{-} \hat{p} \hat{q} {\hat{q}}^{†} {\hat{p}}^{†} . \end{matrix}

(C12)

When pq ≠ rs or pq ≠ sr,

{{\hat{p}}^{†} {\hat{q}}^{†} \hat{r} \hat{s}} \equiv {\hat{p}}^{†} {\hat{q}}^{†} \hat{r} \hat{s} .

(C13)

The thermal average of a normal-ordered product of any number of operators is also always zero by construction. This can be verified for the above product as follows:

\begin{align} 〈 〈 I | x_{p q p q} {{\hat{p}}^{†} {\hat{q}}^{†} \hat{q} \hat{p}} | I 〉 〉 \\ = 〈 〈 I | x_{p q p q} f_{p}^{+} f_{q}^{+} {\hat{p}}^{†} {\hat{q}}^{†} \hat{q} \hat{p} | I 〉 〉 - 〈 〈 I | x_{p q p q} f_{p}^{+} f_{q}^{-} {\hat{p}}^{†} \hat{q} {\hat{q}}^{†} \hat{p} | I 〉 〉 \\ - 〈 〈 I | x_{p q p q} f_{p}^{-} f_{q}^{+} \hat{p} {\hat{q}}^{†} \hat{q} {\hat{p}}^{†} | I 〉 〉 + 〈 〈 I | x_{p q p q} f_{p}^{-} f_{q}^{-} \hat{p} \hat{q} {\hat{q}}^{†} {\hat{p}}^{†} | I 〉 〉 \\ = \sum_{p, q} x_{p q p q} f_{p}^{+} f_{q}^{+} f_{p}^{-} f_{q}^{-} - \sum_{p, q} x_{p q p q} f_{p}^{+} f_{q}^{-} f_{p}^{-} f_{q}^{+} \\ - \sum_{p, q} x_{p q p q} f_{p}^{-} f_{q}^{+} f_{p}^{+} f_{q}^{-} + \sum_{p, q} x_{p q p q} f_{p}^{-} f_{q}^{-} f_{p}^{+} f_{q}^{+} \\ = 0, \end{align}

(C14)

where Boltzmann sum identities III, V, and IX of Ref. 55 were used.

For a product of more than two operators, there are two types of contractions: partial and full contractions. Together, they are defined as

(C15)

where {⋯} in the right-hand side gathers the uncontracted operators (kept in the original order) and P is the number of permutations necessary to reorder the whole operator sequence from left to right. If there are no operators left in {⋯}, it is a full contraction; otherwise, it is a partial contraction. A full contraction is a real number, while a partial contraction is a normal-ordered operator multiplied by a real number. The parity of a full contraction is +1 (−1) for an even (odd) number of intersections of the staple symbols of Wick contractions.⁶³

2. Wick’s theorem at finite temperature

Wick’s theorem states that a product of operators is the sum of its normal-ordered product and all of its partial and full contractions.⁶³ This holds true at zero and nonzero temperatures. For a product of two operators, the theorem is just a restatement of the definitions of Wick contractions,

(C16)

(C17)

For a product of four operators, the theorem asserts

(C18)

We can verify this for the following simpler but nontrivial case by inspection:

(C19)

(C20)

\begin{matrix} = & f_{p}^{+} f_{q}^{+} {\hat{p}}^{†} {\hat{q}}^{†} \hat{q} \hat{p} - f_{p}^{+} f_{q}^{-} {\hat{p}}^{†} \hat{q} {\hat{q}}^{†} \hat{p} - f_{p}^{-} f_{q}^{+} \hat{p} {\hat{q}}^{†} \hat{q} {\hat{p}}^{†} \\ + f_{p}^{-} f_{q}^{-} \hat{p} \hat{q} {\hat{q}}^{†} {\hat{p}}^{†} + f_{p}^{-} (f_{q}^{+} {\hat{q}}^{†} \hat{q} - f_{q}^{-} \hat{q} {\hat{q}}^{†}) \\ + f_{q}^{-} (f_{p}^{+} {\hat{p}}^{†} \hat{p} - f_{p}^{-} \hat{p} {\hat{p}}^{†}) + f_{p}^{-} f_{q}^{-} \end{matrix}

(C21)

\begin{matrix} = & f_{p}^{+} f_{q}^{+} {\hat{p}}^{†} {\hat{q}}^{†} \hat{q} \hat{p} - f_{p}^{+} f_{q}^{-} {\hat{p}}^{†} (1 - {\hat{q}}^{†} \hat{q}) \hat{p} \\ - f_{p}^{-} f_{q}^{+} {\hat{q}}^{†} (1 - {\hat{p}}^{†} \hat{p}) \hat{q} + f_{p}^{-} f_{q}^{-} (1 - {\hat{p}}^{†} \hat{p}) (1 - {\hat{q}}^{†} \hat{q}) \\ + f_{p}^{-} (f_{q}^{+} {\hat{q}}^{†} \hat{q} - f_{q}^{-} \hat{q} {\hat{q}}^{†}) + f_{q}^{-} (f_{p}^{+} {\hat{p}}^{†} \hat{p} - f_{p}^{-} \hat{p} {\hat{p}}^{†}) + f_{p}^{-} f_{q}^{-} \end{matrix}

(C22)

\begin{matrix} = & (f_{p}^{+} f_{q}^{+} + f_{p}^{+} f_{q}^{-} + f_{p}^{-} f_{q}^{+} + f_{p}^{-} f_{q}^{-}) {\hat{p}}^{†} {\hat{q}}^{†} \hat{q} \hat{p} \\ + f_{p}^{-} f_{q}^{-} - f_{p}^{-} f_{q}^{-} ({\hat{p}}^{†} \hat{p} + \hat{p} {\hat{p}}^{†}) \\ - f_{p}^{-} f_{q}^{-} ({\hat{q}}^{†} \hat{q} + \hat{q} {\hat{q}}^{†}) + f_{p}^{-} f_{q}^{-} \end{matrix}

(C23)

= {\hat{p}}^{†} {\hat{q}}^{†} \hat{q} \hat{p} .

(C24)

When pq ≠ rs or pq ≠ sr, Eq. (C13) embodies the theorem. In general situations, the theorem can be proven by induction. See Ref. 63 for an outline of the proof for the zero-temperature case.

Since the thermal average of any normal-ordered product of operators vanishes, a nonzero thermal average arises solely from full contractions, e.g.,

(C25)

However, more frequently, we need to evaluate the thermal average of a product of two or more normal-ordered products. We show that such a product is the sum of concatenated normal-ordered product plus all of its partial and full contractions excluding internal contractions. An internal contraction is the one that involves at least one contraction within an original normal-ordered product. The following is an illustration of this rule:

(C26)

Internal contractions in this example are the ones that contract ${\hat{p}}^{†}$ with $\hat{q}$ and/or ${\hat{r}}^{†}$ with $\hat{s}$ ⁠, which are therefore excluded from the right-hand side. This rule can be justified by comparing Eq. (C18) with

(C27)

The last three terms are internal contractions, which cancel out the same in Eq. (C18) when equated with the above. Therefore, the thermal average of a product of normal-ordered products comes solely from its full contractions excluding internal contractions, e.g.,

(C28)

3. Hamiltonian and number operators

Let us express the number operator in the normal-ordered second quantization at finite temperature,

(C29)

where we used the zeroth-order electroneutrality condition [Eq. (45)] in the last equality. Taking the thermal average of the above equation, we immediately find $⟨ \hat{N} ⟩ = \bar{N}$ ⁠.

The Hamiltonian operator in the second quantization is written as^63,74

\hat{H} = E_{nuc.} + \sum_{p, q} h_{p q} {\hat{p}}^{†} \hat{q} + \frac{1}{4} \sum_{p, q, r, s} ⟨ p q ‖ r s ⟩ {\hat{p}}^{†} {\hat{q}}^{†} \hat{s} \hat{r},

(C30)

where E_nuc. is the nucleus–nucleus repulsion energy, h_pq is the integral of the one-electron part of the Hamiltonian operator, and ⟨pq‖rs⟩ is the antisymmetrized two-electron integral. This can be finite-temperature normal ordered as

(C31)

where E_HF(T) and f_pq(T) are the finite-temperature Hartree–Fock energy and Fock matrix element,⁵⁸ defined by

E_{HF} (T) = E_{nuc.} + \sum_{p} h_{p p} f_{p}^{-} + \frac{1}{2} \sum_{p, q} 〈 p q ‖ p q 〉 f_{p}^{-} f_{q}^{-},

(C32)

f_{p q} (T) = h_{p q} + \sum_{r} 〈 p r ‖ q r 〉 f_{r}^{-} .

(C33)

Therefore, E_HF(T) is identified as the following sum:

E_{HF} (T) = 〈 E_{I}^{(0)} 〉 + 〈 E_{I}^{(1)} 〉,

(C34)

where the thermal averages of the zeroth- and first-order energy corrections^54,55 are given by

⟨ E_{I}^{(0)} ⟩ = E_{nuc.} + \sum_{p} ϵ_{p} f_{p}^{-},

(C35)

⟨ E_{I}^{(1)} ⟩ = \sum_{p} F_{p p} f_{p}^{-} - \frac{1}{2} \sum_{p, q} ⟨ p q ‖ p q ⟩ f_{p}^{-} f_{q}^{-},

(C36)

with F_pq being the difference between the finite-temperature Fock and diagonal zero-temperature Fock matrix elements,

F_{p q} = f_{p q} (T) - ϵ_{p} δ_{p q} .

(C37)

Equation (C34) mirrors the well-known identity⁷⁴ in the zero-temperature case, E_HF = E⁽⁰⁾ + E⁽¹⁾, but it does not agree with either U⁽⁰⁾ + U⁽¹⁾ or Ω⁽⁰⁾ + Ω⁽¹⁾ at T > 0. They are instead related to one another by^54,55

\begin{align} U^{(0)} + U^{(1)} & = E_{nuc.} + \sum_{p} ϵ_{p} f_{p}^{-} + \sum_{p} F_{p p} f_{p}^{-} \\ - \frac{1}{2} \sum_{p, q} 〈 p q ‖ p q 〉 f_{p}^{-} f_{q}^{-} - β \sum_{p} (F_{p p} - μ^{(1)}) ϵ_{p} f_{p}^{-} f_{p}^{+} \\ = E_{HF} (T) + T S^{(1)}, \end{align}

(C38)

\begin{align} Ω^{(0)} + Ω^{(1)} & = E_{nuc.} + \frac{1}{β} \sum_{p} \ln f_{p}^{+} + \sum_{p} F_{p p} f_{p}^{-} \\ - \frac{1}{2} \sum_{p, q} 〈 p q ‖ p q 〉 f_{p}^{-} f_{q}^{-} - μ^{(1)} \bar{N} \\ = E_{HF} (T) - (μ^{(0)} + μ^{(1)}) \bar{N} - T S^{(0)} . \end{align}

(C39)

We employ the Møller–Plesset partitioning⁷² of the Hamiltonian. The zeroth-order Hamiltonian and perturbation operators are then written in the normal-ordered form as

{\hat{H}}_{0} \equiv E_{nuc.} + \sum_{p} ϵ_{p} {\hat{p}}^{†} \hat{p} = ⟨ E_{I}^{(0)} ⟩ + \sum_{p} ϵ_{p} {{\hat{p}}^{†} \hat{p}},

(C40)

\hat{V} \equiv ⟨ E_{I}^{(1)} ⟩ + \sum_{p, q} F_{p q} {{\hat{p}}^{†} \hat{q}} + \frac{1}{4} \sum_{p, q, r, s} ⟨ p q ‖ r s ⟩ {{\hat{p}}^{†} {\hat{q}}^{†} \hat{s} \hat{r}} .

(C41)

4. Projector and resolvent operator

In order to understand the second-quantization rule for the resolvent operator, let us first examine the effect of inserting the resolution-of-the-identity in Eq. (C28). Using Einstein’s implied summation of repeated spin orbital and state indices, we can rewrite Eq. (C28) as

(C42)

where Q runs over all determinants, and hence, ∑|Q⟩⟨Q| = 1. In the second equality, we wrote $| Q 〉 = {{\hat{t}}^{†} \hat{u}} | I 〉$ ⁠. For the last equality to hold, we should regard a chain of contractions, p-u-s or q-t-r, as a single contraction evaluated as $f_{p}^{-}$ or $f_{q}^{+}$ ⁠, respectively, rather than view each as three consecutive contractions evaluated as $f_{p}^{-} f_{p}^{+} f_{p}^{-}$ or $f_{q}^{+} f_{q}^{-} f_{q}^{+}$ ⁠, which is erroneous.

Then, the rule for evaluating thermal averages involving a resolvent operator is essentially the same. To explain the rule by a similar example, we may consider

(C43)

where A stands for a determinant outside the degenerate subspace of I, which explains why the summations over p and q must be restricted (“denom. ≠ 0”) to nonzero denominators. In other words, ${\hat{R}}_{I}$ acts as an outer projector.

Let us define the inner projector ${\hat{P}}_{I} = \sum | J 〉〈 J |$ ⁠, where J runs over only those determinants that are degenerate with I,

(C44)

where the summations over p and q are now restricted (“denom. = 0”) to the cases whose fictitious denominator (ϵ_p − ϵ_q) is zero.

REFERENCES

1.

Y.

Kurashige

and

T.

Yanai

,

J. Chem. Phys.

130

,

234114

(

2009

).

https://doi.org/10.1063/1.3152576

Crossref

PubMed

2.

D.

Zgid

and

G. K.-L.

Chan

,

J. Chem. Phys.

134

,

094115

(

2011

).

https://doi.org/10.1063/1.3556707

Crossref

PubMed

3.

G. H.

Booth

,

D.

Cleland

,

A. J. W.

Thom

, and

A.

Alavi

,

J. Chem. Phys.

135

,

084104

(

2011

).

https://doi.org/10.1063/1.3624383

Crossref

PubMed

4.

J. K.

Ellis

,

R. L.

Martin

, and

G. E.

Scuseria

,

J. Chem. Theory Comput.

9

,

2857

(

2013

).

https://doi.org/10.1021/ct400307d

Crossref

PubMed

5.

L.

Bytautas

,

G. E.

Scuseria

, and

K.

Ruedenberg

,

J. Chem. Phys.

143

,

094105

(

2015

).

https://doi.org/10.1063/1.4929904

Crossref

PubMed

6.

Y.

Ohtsuka

and

S.

Ten-no

,

J. Chem. Phys.

143

,

144116

(

2015

).

https://doi.org/10.1063/1.4936412

Crossref

7.

G. K.-L.

Chan

,

A.

Keselman

,

N.

Nakatani

,

Z.

Li

, and

S. R.

White

,

J. Chem. Phys.

145

,

014102

(

2016

).

https://doi.org/10.1063/1.4955108

Crossref

PubMed

8.

A. A.

Holmes

,

N. M.

Tubman

, and

C. J.

Umrigar

,

J. Chem. Theory Comput.

12

,

3674

(

2016

).

https://doi.org/10.1021/acs.jctc.6b00407

Crossref

PubMed

9.

Y. Y.

He

,

M. P.

Qin

,

H.

Shi

,

Z. Y.

Lu

, and

S. W.

Zhang

,

Phys. Rev. B

99

,

045108

(

2019

).

https://doi.org/10.1103/physrevb.99.045108

Crossref

10.

A.

Baiardi

and

M.

Reiher

,

J. Chem. Phys.

152

,

040903

(

2020

).

https://doi.org/10.1063/1.5129672

Crossref

PubMed

11.

N. H.

Stair

and

F. A.

Evangelista

,

J. Chem. Phys.

153

,

104108

(

2020

).

https://doi.org/10.1063/5.0014928

Crossref

PubMed

12.

Z.-H.

Cui

,

T.

Zhu

, and

G. K.-L.

Chan

,

J. Chem. Theory Comput.

16

,

119

(

2020

).

https://doi.org/10.1021/acs.jctc.9b00933

Crossref

PubMed

13.

J. E.

Deustua

,

J.

Shen

, and

P.

Piecuch

,

J. Chem. Phys.

154

,

124103

(

2021

).

https://doi.org/10.1063/5.0045468

Crossref

PubMed

14.

M.

Motta

,

D. M.

Ceperley

,

G.

Finite-temperature many-body perturbation theory for electrons: Algebraic recursive definitions, second-quantized derivation, linked-diagram theorem, general-order algorithms, and grand canonical and canonical ensembles

I. INTRODUCTION

II. ALGEBRAIC RECURSIVE DEFINITIONS

A. Grand canonical ensemble

B. Recursion for Ω(n)

C. Recursion for ⟨EI(n)⟩

D. Recursion for μ(n)

E. Recursion for U(n)

F. Recursion for S(n)

III. SECOND-QUANTIZED DERIVATION

A. Zeroth order

B. First order

C. Second order

D. Third order

IV. FEYNMAN DIAGRAMS

A. First order

B. Second order

C. Third order

D. Fourth order

V. LINKED-DIAGRAM THEOREM

VI. NUMERICAL BENCHMARKS

A. General-order algorithms

B. Grand canonical ensemble

C. Canonical ensemble

VII. CONCLUSIONS

ACKNOWLEDGMENTS

DATA AVAILABILITY

APPENDIX A: RECURSIONS IN THE CANONICAL ENSEMBLE

APPENDIX B: A PROOF OF EQUIVALENCE OF EQS. (13) AND (14)

APPENDIX C: NORMAL-ORDERED SECOND QUANTIZATION AT FINITE TEMPERATURE

1. Normal ordering at finite temperature

2. Wick’s theorem at finite temperature

3. Hamiltonian and number operators

4. Projector and resolvent operator

REFERENCES

B. Recursion for Ω⁽ⁿ⁾

C. Recursion for $⟨ E_{I}^{(n)} ⟩$

D. Recursion for μ⁽ⁿ⁾

E. Recursion for U⁽ⁿ⁾

F. Recursion for S⁽ⁿ⁾