Converging finite-temperature many-body perturbation theory that conserves the average number of electrons

TL;DR

This paper introduces a finite-temperature perturbation theory for the grand canonical ensemble that maintains charge neutrality by conserving the average number of electrons at each perturbation order, with analytical formulas verified against benchmark data.

Contribution

It develops a new finite-temperature perturbation framework that conserves electron number and provides analytical formulas for first-order corrections, bridging zero-temperature theories.

Findings

Formulas verified against benchmark data

Reduces to Møller-Plesset theory for nondegenerate ground states

Correct zero-temperature limit for degenerate ground states

Abstract

A finite-temperature perturbation theory for the grand canonical ensemble is introduced that expands chemical potential in a perturbation series and conserves the average number of electrons, ensuring charge neutrality of the system at each perturbation order. Two classes of (sum-over-state and reduced) analytical formulas are obtained in a straightforward, algebraic, time-independent derivation for the first-order corrections to chemical potential, grand potential, and internal energy, with the aid of several identities of the Boltzmann sums also introduced in this study. These formulas are numerically verified against benchmark data from thermal full configuration interaction. In the zero-temperature limit, the finite-temperature perturbation theory reduces analytically to and is consistent with the M{\o}ller-Plesset perturbation theory, but only for a nondegenerate ground state. For a degenerate ground state, its correct zero-temperature limit is not the M{\o}ller-Plesset perturbation theory, but the Hirschfelder-Certain degenerate perturbation theory.

Tables

Table 1. Table 1: Comparison of the first-order corrections to chemical potential ( μ ( 1 ) superscript 𝜇 1 \mu^{(1)} ), grand potential ( Ω ( 1 ) superscript Ω 1 \Omega^{(1)} ), and internal energy ( U ( 1 ) superscript 𝑈 1 U^{(1)} ) obtained with the λ 𝜆 \lambda -variation, sum-over-states formulas, and reduced analytical formulas as a function of temperature ( T 𝑇 T ) for the hydrogen fluoride molecule (0.9168 Å) in the STO-3G basis set.

	${\mu }^{\left(1\right)}/E_{\text{h}}$			${\Omega }^{\left(1\right)}/E_{\text{h}}$			$U^{\left(1\right)}/E_{\text{h}}$
	$\lambda$-variation111Obtained using a central seven-point finite difference. See also Jha and Hirata Jha and Hirata (2019a).	SoS222Obtained using the sum-over-states formulas with $\left\{E_I^{\left(1\right)}\right\}$ computed by the $\lambda$-variation method using forward seven- through nine-point finite differences.	Analytical	$\lambda$-variation111Obtained using a central seven-point finite difference. See also Jha and Hirata Jha and Hirata (2019a).	SoS222Obtained using the sum-over-states formulas with $\left\{E_I^{\left(1\right)}\right\}$ computed by the $\lambda$-variation method using forward seven- through nine-point finite differences.	Analytical	$\lambda$-variation111Obtained using a central seven-point finite difference. See also Jha and Hirata Jha and Hirata (2019a).	SoS222Obtained using the sum-over-states formulas with $\left\{E_I^{\left(1\right)}\right\}$ computed by the $\lambda$-variation method using forward seven- through nine-point finite differences.	Analytical
$T/\text{K}$	Eq. (9)	Eq. (29)	Eq. (48)	Eq. (9)	Eq. (30)	Eq. (46)	Eq. (9)	Eq. (32)	Eq. (49)
${10}^4$	$0.0000$	$0.0000$	$-0.0000$	$-45.9959$	$-45.9961$333This value seems to suffer from a slightly greater numerical error.	$-45.9959$	$-45.9959$	$-45.9959$	$-45.9959$
${10}^5$	$-0.0752$	$-0.0752$	$-0.0752$	$-45.2684$	$-45.2684$	$-45.2684$	$-45.9479$	$-45.9479$	$-45.9479$
${10}^6$	$-0.1690$	$-0.1690$	$-0.1690$	$-44.5256$	$-44.5256$	$-44.5256$	$-46.1767$	$-46.1767$	$-46.1767$
${10}^7$	$-0.2981$	$-0.2981$	$-0.2981$	$-43.1991$	$-43.1991$	$-43.1991$	$-46.2355$	$-46.2355$	$-46.2355$
${10}^8$	$-0.4122$	$-0.4122$	$-0.4122$	$-41.9847$	$-41.9847$	$-41.9847$	$-46.1180$	$-46.1180$	$-46.1180$

Equations

$$\Xi =\sum _{I=1}^{2^m}{e}^{-\beta E_I+\beta \mu N_I},$$

$$\bar{N}=\frac{{\sum }_{I=1}^{2^m}{N}_I{e}^{-\beta E_I+\beta \mu N_I}}{{\sum }_{I=1}^{2^m}{e}^{-\beta E_I+\beta \mu N_I}},$$

$$\Omega$$

$$U$$

$$\Omega$$

$$\Xi$$

$$\mu$$

$$E_I$$

$$X^{(n)}$$

$$e^{a+b}$$

$$\ln (a+b)$$

$${\Xi }^{(0)}$$

$${\xi }_I^{(0)}=e^{-\beta E_I^{(0)}+\beta {\mu }^{(0)}{N}_I}$$

$$E_I^{(0)}=E_{\text{nuc.}}+\sum _i^{\text{occ.}}{ϵ}_i^{(0)};N_I=\sum _i^{\text{occ.}}1,$$

$${\Xi }^{(0)}$$

$$=$$

$$\bar{N}$$

$$=$$

$$f_i^{-}=\frac{1}{1+e^{\beta ({ϵ}_i^{(0)}-{\mu }^{(0)})}}.$$

$${\Omega }^{(0)}$$

$$U^{(0)}$$

$$f_i^{+}=1-f_i^{-}.$$

$$\frac{{\sum }_{I=1}^{2^m}{a}_I{e}^{-\beta E_I^{(0)}+\beta {\mu }^{(0)}{N}_I}}{{\sum }_{I=1}^{2^m}{e}^{-\beta E_I^{(0)}+\beta {\mu }^{(0)}{N}_I}}$$

$${\Xi }^{(1)}$$

$${\xi }_I^{(1)}={\xi }_I^{(0)}(-\beta E_I^{(1)}+\beta {\mu }^{(1)}{N}_I),$$

$$\bar{N}$$

$${\mu }^{(1)}$$

$${\Omega }^{(1)}$$

$$U^{(1)}$$

$$\frac{{\sum }_{I=1}^{2^m}{\xi }_I^{(0)}(E_I^{(0)}-{\mu }^{(0)}{N}_I)}{{\Xi }^{(0)}}=U^{(0)}-{\mu }^{(0)}\bar{N},$$

$$\frac{{\sum }_{i=1}^m{a}_i{b}_i{e}^{{\nu }_i}+{\sum }_{i<j}^m(a_i+a_j)(b_i+b_j)e^{{\nu }_i}{e}^{{\nu }_j}+\cdots +{\sum }_{i_1<\cdots <i_m}^m(a_{i_1}+\cdots +a_{i_m})(b_{i_1}+\cdots +b_{i_m})e^{{\nu }_{i_1}}\cdots e^{{\nu }_{i_m}}}{(1+e^{{\nu }_1})(1+e^{{\nu }_2})\cdots (1+e^{{\nu }_m})}$$

$$=\sum _{i=1}^m{a}_i{b}_i{f}_i^{-}{f}_i^{+}+\sum _{i=1}^m{a}_i{f}_i^{-}\sum _{i=1}^m{b}_i{f}_i^{-};(\text{Identity II}).$$

$$\frac{{\sum }_{i<j}^m{c}_{ij}{e}^{{\nu }_i}{e}^{{\nu }_j}+{\sum }_{i<j<k}^m(c_{ij}+c_{ik}+c_{jk})e^{{\nu }_i}{e}^{{\nu }_j}{e}^{{\nu }_k}+\cdots +{\sum }_{i_1<\cdots <i_m}^m(c_{i_1{i}_2}+\cdots +c_{i_{m-1}{i}_m})e^{{\nu }_{i_1}}\cdots e^{{\nu }_{i_m}}}{(1+e^{{\nu }_1})(1+e^{{\nu }_2})\cdots (1+e^{{\nu }_m})}$$

Full text

Keywords: thermodynamics; many-body perturbation theory; temperature; chemical potential; grand canonical ensemble; grand potential; internal energy

I Introduction

In the previous Chapter Jha and Hirata (2019a), we showed that the finite-temperature perturbation theory for electrons described in a number of modern textbooks Thouless (1990); Mattuck (1992); March et al. (1995); Fetter and Walecka (2003) and excellent reviews Bloch and De Dominicis (1958); Balian et al. (1961); Bloch (1965); Santra and Schirmer (2017) does not reproduce the benchmark data obtained as the $\lambda$-derivatives of the exact finite-temperature theory Kou and Hirata (2014) with the perturbation-scaled Hamiltonian, $\hat{H}=\hat{H}_{0}+\lambda\hat{V}$. The discrepancy was traced Jha and Hirata (2019a); White and Chan (2018a) not so much to mathematical issues but to a misuse of the ensemble. The theory given in textbooks (beyond the zeroth-order Fermi–Dirac theory) adopts a grand canonical ensemble that holds chemical potential fixed at some arbitrary value and allows the average number of electrons to fluctuate, violating charge neutrality as a basic tenet of equilibrium thermodynamics Levin et al. (2014). Note that the energy of a system with long-range unscreened interactions such as a charged plasma is not thermodynamically extensive Fisher and Ruelle (1966); Dyson and Lenard (1967); Hirata et al. (2012). As a consequence, the perturbation theory described in textbooks does not converge at the exact limit of the neutral system. Worse still, beyond the zeroth order, it describes a massively charged system to which equilibrium thermodynamics is no longer applicable.

We need a finite-temperature perturbation theory that also expands the chemical potential in a perturbation series so as to ensure charge neutrality at any perturbation order and thereby converges at the exact, thermal full-configuration-interaction (FCI) Kou and Hirata (2014) limits. In this Chapter, we present such a perturbation theory and derive analytical formulas of its first-order corrections to the chemical potential, grand potential, and internal energy. Two classes of analytical formulas were derived in an easy-to-follow, time-independent, algebraic (nondiagrammatic) manner: The sum-over-states formulas are correct in all cases, while the more compact, reduced formulas are valid only for a nondegenerate ground state at any temperature or for a degenerate or nondegenerate ground state at a nonzero temperature. We show that both formulas reproduce the benchmark data obtained as the $\lambda$-derivatives of the thermal FCI results in the case of a nondegenerate ground state.

There have been a persistent question as to whether the zero-temperature limits of a finite-temperature perturbation theory for internal energy should agree with the Møller–Plesset (MP) perturbation theory Shavitt and Bartlett (2009) for energy at respective orders. Kohn and Luttinger Kohn and Luttinger (1960) were the first to pose a related (if not the same) question, taking the finite-temperature second-order perturbation theory for a homogeneous electron gas (having a degenerate ground state) as an example. We show that, for a nondegenerate ground state, the finite-temperature perturbation theory introduced here is consistent with the MP perturbation theory in the sense that the first-order internal-energy correction reduces to the first-order MP energy correction both analytically and numerically. For a degenerate ground state, however, the finite-temperature perturbation theory is not consistent with the MP perturbation theory; it instead reduces to and is consistent with the Hirschfelder–Certain (HC) degenerate perturbation theory Hirschfelder and Certain (1974).

II Finite-temperature perturbation theory in grand canonical ensemble

The electronic grand partition function Kou and Hirata (2014) per molecule of an ideal gas of identical molecules at given temperature $T$ is

[TABLE]

where $\beta=(k_{\text{B}}T)^{-1}$, $E_{I}$ is the FCI energy of the $I$th state and $N_{I}$ is the number of electrons in the same state. The summation over a capital letter index ($I$) runs over all (exponentially many) $2^{m}$ states spanned by $m$ spinorbital basis functions. Chemical potential $\mu$ is determined by solving

[TABLE]

where $\bar{N}$ is a given average number of electrons that keeps the system electrically neutral. These two are the equations of state to be solved simultaneously for $\Xi$ and $\mu$. A valid grand canonical ensemble for electrons holds $\bar{N}$ fixed at the value that maintains neutrality of the system by adjusting $\mu$ (not vice versa). This is not a choice but a nonnegotiable requirement, without which the system is massively charged, causing the energy per particle to be divergent Fisher and Ruelle (1966); Dyson and Lenard (1967); Hirata et al. (2012) and equilibrium thermodynamics to break down Levin et al. (2014).

Grand potential $\Omega$ and internal energy $U$ per molecule are related to $\Xi$ by

[TABLE]

and also to each other by

[TABLE]

where $S$ is entropy.

We expand $\Xi$ and $\mu$ as well as $\Omega$ and $U$ in converging perturbation series,

[TABLE]

and so forth, and seek to write each correction in terms of the perturbation corrections of the energy,

[TABLE]

where $E_{I}^{(n)}$ is the $n$th-order correction to the $I$th-state energy according to the HC degenerate perturbation theory Hirschfelder and Certain (1974), of which the MP perturbation theory Shavitt and Bartlett (2009) for a nondegenerate ground state is a special case.

They conform to the canonical definition of a perturbation theory, in which the $n$th-order correction of quantity $X$ is

[TABLE]

where $X(\lambda)$ is the value of $X$ (such as $E_{I}$, $\Xi$, $\mu$, $\Omega$, and $U$) calculated by FCI with the Hamiltonian $\hat{H}_{0}+\lambda\hat{V}$. We adopt the Møller–Plesset partitioning and, therefore, the zeroth-order Hamiltonian, $\hat{H}_{0}$, is the sum of the (zero-temperature) Fock operators plus the nuclear repulsion energy. We emphasize again that the perturbation corrections to energies that exactly match this canonical definition are furnished, in the most general case, by the HC degenerate perturbation theory and not by the MP perturbation theory (because many excited, ionized, and electron-attached states are exactly degenerate at the zeroth order even if the ground state is not).

We have obtained numerically these corrections at lower orders and shown that the finite-temperature perturbation theory described in textbooks does not reproduce such benchmark data Jha and Hirata (2019a) because it neglects to expand $\mu$ in a perturbation series, thereby allowing $\bar{N}$ to fluctuate and causing the system to be massively charged. We call the method that evaluates these derivatives by finite-difference formulas the $\lambda$-variation method Hirata et al. (2017); Jha and Hirata (2019a). In this study, we will perform the differentiation analytically to arrive at the corresponding analytical formulas for $\mu^{(1)}$, $\Omega^{(1)}$, and $U^{(1)}$.

In the following, we use the Taylor series,

[TABLE]

which are rapidly convergent when $a\gg b$. We will also use four identities of the Boltzmann sums, which will be introduced as needed.

III Zeroth order

The zeroth-order theory Fetter and Walecka (2003) in the grand canonical ensemble is correct in the sense that chemical potential $\mu^{(0)}$ is adjusted to keep the average number of electrons constant at $\bar{N}$. We review a time-independent and algebraic (as opposed to time-dependent or diagrammatic) derivation of this theory as a preparation for its application to the first-order theory.

III.1 Fermi–Dirac theory

The zeroth-order grand partition function is given by

[TABLE]

with

[TABLE]

and

[TABLE]

where $E_{\text{nuc.}}$ is the nuclear repulsion energy and $\epsilon_{i}^{(0)}$ is the $i$th canonical Hartree–Fock (HF) spinorbital energy. The summation over small letter index ($i$) runs over spinorbitals occupied in the $I$th Slater determinant, which is a solution of the zeroth-order Schrödinger equation. The occupied orbitals should not be confused with those occupied in the HF ground state; they vary depending on state $I$.

This can be simplified to

[TABLE]

where products over $i$ run over all $m$ spinorbitals spanned by the basis set. Equation (2) also reduces to

[TABLE]

where $f_{i}^{-}$ is the Fermi–Dirac occupancy function of the $i$th spinorbital,

[TABLE]

Using $\Xi^{(0)}$ of equation (16) with $\mu^{(0)}$ determined by solving equation (20), we obtain the zeroth-order grand potential and internal energy as

[TABLE]

where $f_{i}^{+}$ is the Fermi–Dirac vacancy function of the $i$th spinorbital,

[TABLE]

III.2 Boltzmann-sum identity I

It is instructive to rederive equation (20) in a more general manner. Let $a_{I}=\sum_{i}^{\text{occ.}}a_{i}$. Using a notation, $\nu_{i}=\beta(\mu^{(0)}-\epsilon_{i}^{(0)})$, and canceling the common factor of $e^{-\beta E_{\text{nuc.}}}$ in the numerator and denominator, we can prove an identity,

[TABLE]

by mathematical induction. Equation (20) is nothing but Identity I with $a_{i}=1$ for any $i$.

IV First order

IV.1 Sum-over-states analytical formulas

An application of equation (10) to equation (1) leads to the following expression for the first-order correction to grand partition function $\Xi$:

[TABLE]

with

[TABLE]

where $E_{I}^{(1)}$ is the first-order correction to the energy of the $I$th state according to the HC degenerate perturbation theory Hirschfelder and Certain (1974).

The same procedure to equation (2) leads to the equation for $\mu^{(1)}$ as

[TABLE]

which can be solved for $\mu^{(1)}$ as

[TABLE]

Using equation (11), we obtain the first-order corrections to the grand potential and internal energy as

[TABLE]

In the last equality, we used

[TABLE]

which can be readily proven with Identity I [equation (25)].

The third and fourth terms of equation (32) contain products of extensive quantities such as $\Omega^{(1)}U^{(0)}$, which are individually non-size-consistent Hirata (2011). However, it can be shown (see below) that they mutually cancel exactly across the two terms, leaving only rigorously extensive contributions; the first-order finite-temperature perturbation theory is size-consistent in spite of these terms.

The aforementioned equations are the correct and complete (if not compact) analytical formulas of the first-order finite-temperature perturbation theory in the grand canonical ensemble applicable to either a degenerate or nondegenerate ground state at any temperature down to zero. Equation (29) is first solved for $\mu^{(1)}$ using $\{E_{I}^{(1)}\}$ unambiguously defined and calculable by the HC degenerate perturbation theory Hirschfelder and Certain (1974). Once $\mu^{(1)}$ is determined, $\Xi^{(1)}$ is evaluated with equation (26). Next, $\Omega^{(1)}$ and $U^{(1)}$ are computed with equations (30) and (32), respectively, using $\mu^{(1)}$ and $\{E_{I}^{(1)}\}$ as well as the zeroth-order quantities. We have verified these formulas numerically (see the Appendix).

IV.2 Boltzmann-sum identities II, III, and IV

Here, we introduce three more identities for the Boltzmann sums, which are useful for deriving compact reduced formulas. The second identity reads

[TABLE]

The third identity is

[TABLE]

where $c_{ij}=c_{ji}$ for any $i$ and $j$. Note that $c_{ii}$ is absent in any summation. Identities II and III can be proven by mathematical induction with Identity I. Finally, the fourth reads

[TABLE]

which can also be proven by mathematical induction with the aid of Identities I, II, and III. These proofs are elementary and not reproduced here.

IV.3 Reduced analytical formulas

The foregoing analytical formulas for the first-order perturbation theory involve sum over exponentially many states, and are inconvenient to use. For a system with a nondegenerate ground state or at a nonzero temperature, we can reduce them into much more compact formulas involving only molecular integrals in the spirit of the Fermi–Dirac theory or (incorrect) diagrammatic finite-temperature perturbation theory given in textbooks.

First, we rewrite equation (30) as

[TABLE]

Here, $E_{I}^{(1)}$ is the first-order HC perturbation correction to the energy Hirschfelder and Certain (1974) and is written as

[TABLE]

where $\Psi_{I}^{(0)}$ is the zeroth-order wave function and $\hat{V}$ is the perturbation operator in the Møller–Plesset partitioning.

For nondegenerate $I$th state, $\Psi_{I}^{(0)}$ is a single Slater determinant $\Phi_{I}$, making the evaluation of the above integral trivial:

[TABLE]

with

[TABLE]

where ‘occ.’ again refers to the spinorbitals occupied in the $I$th Slater determinant, whereas ‘HOMO’ stands for the highest-occupied spinorbital of the $\bar{N}$-electron ground state (hereafter simply the ground state); the last summation runs over spinorbitals occupied in the ground state. Each four-index integral is a usual antisymmetrized two-electron integral Shavitt and Bartlett (2009).

When the $I$th state is degenerate (in the sense that there are $d$ states sharing the same zeroth-order energy), $\Psi_{I}^{(0)}$ is no longer a Slater determinant, but a linear combination of $d$ degenerate Slater determinants, $\{\Phi_{J}\}$ Hirschfelder and Certain (1974).

[TABLE]

where $U_{JI}$ is an element of a unitary matrix, which can only be determined by a highly involved procedure of the HC degenerate perturbation theory Hirschfelder and Certain (1974). It is more a rule than an exception that the zeroth-order energies of Slater determinants display degeneracy.

However, thanks to the orthonormality of $\{\Psi_{I}^{(0)}\}$ within each degenerate group, the following equality holds:

[TABLE]

where the summation over $J$ is taken for all $d$ degenerate states. In other words, although $E_{J}^{(1)}$ for the $J$th individual state cannot be written simply as the right-hand side of the above equation, its average over all $d$ degenerate states can be Santra and Schirmer (2017).

Therefore, if all $E_{J}^{(1)}$ from a degenerate group are summed with an equal weight, one obtains the correct final formula by pretending as if equation (39) were true for each state. This is indeed the case with equation (37) if the ground state is nondegenerate or temperature is nonzero. For a degenerate ground state at zero temperature, the states that are degenerate with the ground state are summed with an unequal weight and the use of equation (42) becomes inappropriate. We will consider this case in Sec. VI, but in this section we confine ourselves to a nondegenerate ground state or nonzero temperature.

We can then write equation (39) as

[TABLE]

with

[TABLE]

and $c_{ii}=0$. Using Identities I and III [equation (25) and (35)], we can reduce equation (37) into

[TABLE]

where the summations run over all spinorbitals.

With the same choice of $\{a_{i}\}$ and $\{c_{ij}\}$ in equation (44) as well as $b_{i}=1$ for any $i$, equation (28) for $\mu^{(1)}$ is simplified with Identities II and IV [equation (34) and (36)] as follows:

[TABLE]

which can be solved for $\mu^{(1)}$ as

[TABLE]

The last equation can also be obtained by a direct application of Identities I–IV to equation (29).

Likewise, applying Identities I–IV to the sum-over-states formula for $U^{(1)}$ [equation (32)], we first observe exact mutual cancellation of non-size-consistent terms and then arrive at a compact analytical formula,

[TABLE]

which no longer contains a non-size-consistent term that is a product of two or more extensive quantities. The same expression can be obtained alternatively Santra and Schirmer (2017) by returning to equation (32), namely,

[TABLE]

and using

[TABLE]

Care must be exercised that partial differentiation with $\beta$ acts only on $f_{i}^{-}$, $f_{i}^{+}$, and $\bar{N}$ via equation (45), but not on $\mu^{(1)}$.

It should be noted that the foregoing reduced formulas are valid only for a nondegenerate ground state or at a nonzero temperature, whereas the sum-over-states formulas and underlying HC degenerate perturbation theory should be correct in all cases.

The corresponding reduced formulas Santra and Schirmer (2017) for $\Omega^{(1)}$ and $U^{(1)}$ according to the textbook ansatz wherein $\mu$ is held fixed at $\mu^{(0)}$ are special cases of the more general expressions reported here and are readily obtained by substituting $\mu^{(1)}=0$ in equations (46) and (49), respectively.

V Comments on related methods

The finite-temperature HF Thouless (1990); Mattuck (1992); March et al. (1995); Fetter and Walecka (2003), finite-temperature Green’s function Phillips and Zgid (2014); Rusakov and Zgid (2016), and thermal coupled-cluster Sanyal et al. (1993); Mandal et al. (2001, 2002); White and Chan (2018b) theories as well as thermal FCI theory Kou and Hirata (2014) all have a provision to adjust $\mu^{(0)}$ so as to restore the correct average number of electrons, making the calculation physically sound and their results more in line with the correct values Jha and Hirata (2019a).

It is then also possible to adjust $\mu^{(0)}$ by some root-finding procedure in the framework of the textbook finite-temperature perturbation theory to restore the net charge neutrality Santra and Schirmer (2017). Then, the theory forms a converging series of approximations towards thermal FCI. However, this series may not be considered a canonical perturbation series defined by equation (9) because adjusting $\mu^{(0)}$ amounts to altering the partitioning of the Hamiltonian at each perturbation order, obscuring the order of perturbation corrections (in the same sense that HF theory involving rotation of orbitals cannot be equated to the first-order perturbation theory). In other words, this series differs from the series presented in this work that gives closed expressions for the first-order corrections to all thermodynamic quantities of interest, which furthermore agree exactly with the $\lambda$-variation results.

VI Zero-temperature limit and degeneracy

VI.1 Zeroth order

In the zero-temperature ($\beta\to\infty$) limit, we can analytically show Kou and Hirata (2014) and numerically confirm Jha and Hirata (2019a)

[TABLE]

where HOMO again stands for the highest-occupied spinorbital of the $\bar{N}$-electron ground state and LUMO the lowest-unoccupied spinorbital.

In this limit, the Fermi–Dirac occupancy function is

[TABLE]

assuming that the spinorbitals are in an ascending order of energy. However, a special consideration Pederson and Jackson (1991) may become necessary when the ground state is degenerate. In this case,

[TABLE]

leaving $f_{\text{HOMO}}^{-}$ and $f_{\text{LUMO}}^{-}$ (as well as the Fermi–Dirac functions of all other degenerate spinorbitalss) undefined, although their sum must take a certain value that maintains the average number of electrons at $\bar{N}$, say,

[TABLE]

in the case of doubly degenerate HOMOs sharing one electron on average. However, a partitioning between $f_{\text{HOMO}}^{-}$ and $f_{\text{LUMO}}^{-}$ is indeterminate.

Nevertheless, we can unambiguously determine the zero-temperature limits of $\Omega^{(0)}$ and $U^{(0)}$ as

[TABLE]

because they do not depend on the partitioning of the Fermi–Dirac occupancy functions among the degenerate HOMOs. These are consistent with the zeroth-order MP perturbation theory and also with each other according to equation (5).

VI.2 First order

Next, we consider the zero-temperature limits of the first-order corrections. We divide the problem into two cases: a nondegenerate ground state versus a degenerate ground state.

In a nondegenerate ground state, the reduced analytical formulas such as equations (47), (46), and (49) are valid at any temperature down to zero. Using equation (54) and the identity Kohn and Luttinger (1960); Hirata and He (2013); Santra and Schirmer (2017),

[TABLE]

we find

[TABLE]

the last of which is the familiar first-order MP perturbation energy correction at zero temperature, $E_{\text{MP1}}$. These are numerically verified in the Appendix. For a nondegenerate ground state, therefore, there is no inconsistency between the finite-temperature perturbation theory and MP perturbation theory.

When the ground state is $d$-fold degenerate ($d>1$), at zero temperature, they are no longer summed at an equal Boltzmann weight, and the reduced analytical formulas derived with equation (42) become invalid. We, therefore, have to revert to the sum-over-states analytical formulas. However, at first glance, even the sum-over-states formulas such as equation (26) appear to suffer from the same issue of ambiguous Boltzmann weights; at zero temperature, $\beta=\infty$ and $E_{1}^{(0)}-\mu^{(0)}N_{1}=0$, rendering the exponent $-\beta E_{1}^{(0)}+\beta\mu^{(0)}N_{1}$ indeterminate.

The correct zero-temperature behavior of the exponents can be inferred by considering the $\lambda$-dependence of the thermal FCI energies or the closely related HC degenerate perturbation corrections to the energies of degenerate states. States that are degenerate at the zeroth order or $\lambda=0$ become nondegenerate or have a lesser degree of degeneracy as $\lambda$ is raised from zero in thermal FCI. The perturbation corrections to the internal energy at $T=0$ are, by definition, the $\lambda$-derivatives of the lowest energy (not the average energy of the degenerate ground states at $\lambda=0$). Correspondingly, the HC degenerate perturbation corrections to energies are generally different from one another for degenerate states, gradually lifting degeneracy as the perturbation order is raised Hirschfelder and Certain (1974). Then, for the first-order correction to be consistent with the $\lambda$-derivative of the corresponding thermal FCI result, the sum over states should be interpreted as being dominated by the zeroth-order state whose first-order energy correction is the most negative. Applying this logic to $U^{(1)}$ [equation (32)], we predict

[TABLE]

where $E_{\text{HC1}}$ is the most negative among the first-order HC degenerate perturbation corrections to the degenerate zeroth-order ground states.

This equation indicates that there is no inconsistency between the finite-temperature perturbation theory and HC degenerate perturbation theory, but there is an inconsistency (reminiscent of the Kohn–Luttinger inconsistency Kohn and Luttinger (1960); Hirata and He (2013); Santra and Schirmer (2017)) between the finite-temperature perturbation theory and MP perturbation theory because $E_{\text{HC1}}$ is not equal to $E_{\text{MP1}}$. The former, $E_{\text{HC1}}$, is determined by the complex, but unambiguous procedure Hirschfelder and Certain (1974) outlined by Hirschfelder and Certain and cannot be expressed simply as the right-hand side of equation (62).

The compact formulas for $U^{(1)}$ [equation (49)] and even for $\Omega^{(1)}$ [equation (46)] appear to reduce analytically to the MP1 expression at zero temperature in all cases including a degenerate ground state. We consider such an interpretation incorrect for three reasons: (1) these formulas are derived using an assumption [equation (42)] that is not valid for a degenerate ground state at zero temperature, (2) $\mu^{(1)}$ is not necessarily zero in this case, and (3) the correct zero-temperature limit is not the MP perturbation theory but the HC degenerate perturbation theory for a degenerate ground state.

We should note that the foregoing argument in the case of degenerate ground states has not been numerically verified unlike the conclusion for nondegenerate ground states, which is numerically confirmed in the Appendix. It is based on an assumption of the availability of strictly degenerate reference wave functions, which are not uniquely determined by the zeroth-order theory. However, only the rate of convergence of the HC degenerate perturbation series and thus of the finite-temperature many-body perturbation theory varies depending on the reference, while the converged limits should be the same correct thermal FCI limits.

VII Summary

Original contributions of this work are summarized as follows:

(1) We have introduced a converging finite-temperature perturbation theory in the grand canonical ensemble that expands chemical potential in a perturbation series and conserves the average number of electrons, thereby ensuring charge neutrality of the system at each perturbation order.

(2) We have obtained the sum-over-states analytical formulas for the first-order corrections of the chemical potential, grand potential, and internal energy, which are valid for a degenerate or nondegenerate ground state and at any temperature down to zero. They are given as equations (29), (30), and (32), respectively. These formulas depend crucially on the HC degenerate perturbation theory, whose details are described elsewhere Hirschfelder and Certain (1974).

(3) We have reported the compact, reduced analytical formulas for the same quantities as equations (48), (46), and (49), respectively, which are valid for a nondegenerate ground state or at a nonzero temperature. We have established a straightforward, time-independent, nondiagrammatic derivation strategy of these formulas by introducing several identities for Boltzmann sums in the grand canonical ensemble.

(4) We shed new light on the issue of the zero-temperature limit of the converging finite-temperature perturbation theory. The finite-temperature perturbation theory introduced here is consistent with the MP perturbation theory for a nondegenerate ground state. For a degenerate ground state, the finite-temperature theory is consistent with the HC degenerate perturbation theory, but not with the MP perturbation theory. For a degenerate ground state, the compact analytical formulas such as equation (49) are correct at $T>0$, but suddenly become incorrect at $T=0$; the sum-over-states formulas should be used instead.

VIII Future outlook

This theory should be extended to the second and higher orders in the grand canonical ensemble as well as to the canonical ensemble Jha and Hirata (2019b). These results should be generalized to a recursive algebraic definition in the style of the Rayleigh–Schödinger perturbation theory, if possible, which should permit a general-order implementation and establish and justify diagrammatic rules for the grand potential and internal energy as well as their size-consistency through a time-indepdendent linked-diagram theorem Shavitt and Bartlett (2009); Hirata et al. (2017).

Acknowledgements.

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

Appendix A Numerical verification

Table 1 documents the first-order corrections to $\mu$, $\Omega$, and $U$ in a converging perturbation series that expands $\mu$ perturbatively and maintains the average number of electrons at $\bar{N}$ at any perturbation order or at any value of perturbation strength ($\lambda$). The $\lambda$-variation data were already reported in the previous Chapter Jha and Hirata (2019a).

At a lower temperature ($10^{4}$ K), a determination of $\mu^{(0)}$ by a bisection method becomes technically difficult, causing some numerical noise. Otherwise, the three methods are mutually in exact agreement, verifying the correctness of formulas and computational methods. The results numerically confirm the zero-temperature limits at $T=10^{4}$ K. The perturbation theory described in textbooks that holds $\mu$ fixed (e.g., at $\mu^{(0)}$) gives completely different values except at the low-temperature limit.