  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\"{o}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 downwards, 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 towards the finite-temperature
full-configuration-interaction limit unless the series diverges. Normal-ordered
Hamiltonian at finite temperature sheds light on the relationship between the
finite-temperature Hartree--Fock and first-order many-body perturbation
theories.
