A finite-temperature perturbation theory for the grand canonical ensemble is introduced that expands the chemical potential in a perturbation series and conserves the average number of electrons, ensuring charge neutrality of the system at each perturbation order. Two (sum-over-state and reduced) classes of analytical formulas are obtained in a straightforward, algebraic, time-independent derivation for the first-order corrections to the 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. For a nondegenerate ground state, the finite-temperature perturbation theory reduces analytically to and is consistent with the Møller–Plesset perturbation theory as temperature (T) tends to zero. For a degenerate ground state, it should instead reduce to the Hirschfelder–Certain degenerate perturbation theory as T → 0.