  A unified theory is presented for finite-temperature many-body perturbation
expansions of the anharmonic vibrational contributions to thermodynamic
functions: the free energy, internal energy, and entropy. The theory is
diagrammatically size-consistent at any order, as ensured by the linked-diagram
theorem proved here, and thus applicable to molecular gases and solids on an
equal footing. It is also a basis-set-free formalism, just like its underlying
Bose-Einstein theory, capable of summing anharmonic effects over an infinite
number of states analytically. It is formulated by the
Rayleigh-Schrodinger-style recursions, generating sum-over-states formulas for
the perturbation series, which unambiguously converges at the
finite-temperature vibrational full-configuration-interaction limits. Two
strategies are introduced to reducing these sum-over-states formulas into
compact sum-over-modes analytical formulas. One is a purely algebraic method
that factorizes each many-mode thermal average into a product of one-mode
thermal averages, which are then evaluated by the thermal Born-Huang rules.
Canonical forms of these rules are proposed, dramatically expediting the
reduction process. The other is finite-temperature normal-ordered second
quantization, which is fully developed in this study, including a proof of
thermal Wick's theorem and the derivation of a normal-ordered vibrational
Hamiltonian at finite temperature. The latter naturally defines a
finite-temperature extension of size-extensive vibrational self-consistent
field theory. These reduced formulas can be represented graphically as Feynman
diagrams with resolvent lines, which include anomalous and renormalization
diagrams. Two order-by-order and one general-order algorithms of computing
these perturbation corrections are implemented and applied up to the eighth
order. The results show no signs of Kohn-Luttinger-type nonconvergence.
