A finite-temperature many-body perturbation theory is presented, which expands in power series the electronic grand potential, chemical potential, internal energy, and entropy on an equal footing. Sum-over-states and sum-over-orbitals analytical formulas for the second-order perturbation corrections to these thermodynamic properties are obtained in a time-independent, nondiagrammatic, algebraic derivation, relying on the sum rules of the Hirschfelder–Certain degenerate perturbation energies in a degenerate subspace as well as nine algebraic identities for the zeroth-order thermal averages of one- through four-indexed quantities and products thereof. They reproduce numerically exactly the benchmark data obtained as the numerical derivatives of the thermal-full-configuration-interaction results for a wide range of temperatures.
It is shown analytically and numerically that the finite-temperature many-body perturbation theory in the grand canonical ensemble has zero radius of convergence at zero temperature when the energy ordering or degree of degeneracy for the ground state changes with the perturbation strength. When the degeneracy of the reference state is partially or fully lifted at the first-order Hirschfelder-Certain degenerate perturbation theory, the second-order grand potential and internal energy diverge as $T\ensuremath{\rightarrow}0$. Contrary to earlier suggestions of renormalizability by the chemical potential $\ensuremath{\mu}$, this nonconvergence, first suspected by Kohn and Luttinger, is caused by the nonanalytic nature of the Boltzmann factor ${e}^{\ensuremath{-}E/{k}_{\text{B}}T}$ at
T=0, also plaguing the canonical ensemble, which does not involve $\ensuremath{\mu}$. The finding reveals a fundamental flaw in perturbation theory, which is deeply rooted in the mathematical limitation of power-series expansions and is unlikely to be removed within its framework.