  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.
