  By combining conventional finite-temperature many-body perturbation theory
with cluster expansions, we develop a systematic method to carry out high order
arbitrary temperature perturbative calculations on the computer. The method is
well suited to studying the thermodynamic properties of quantum disordered and
quantum critical phases at finite temperatures. As an application, we calculate
the magnetic susceptibility, internal energy and specific heat of the bilayer
Heisenberg model. It is shown that for a wide range of coupling constants these
expansions show excellent convergence at all temperatures. Comparing the direct
series (without extrapolations) for the bulk susceptibility to Quantum Monte
Carlo simulations we find an almost perfect agreement between the two methods
even at the quantum critical coupling separating the dimerized and
antiferromagnetic phases. The convergence fails only at very low temperatures,
which are also difficult to reach by Quantum Monte Carlo simulations.
