  A normal ordered exponential parametrization is used to obtain equations for
thermal one-and two-particle reduced density matrices, as well as free
energies, partition functions and entropy for both Fermionic (electronic) and
Bosonic (vibrational) Hamiltonians. A first principles derivation of the
equations, relying only on a simple Wick's theorem and starting from the
differential equation $\frac{d \hat{D}}{d \beta}= - (\hat{H}-\mu
\hat{N})\hat{D}$, is presented that yields a differential equation for the
amplitudes representing density cumulants, as well as the grand potential. In
contrast to other approaches reported in the literature the theory does not use
perturbation theory in the interaction picture and an integral formulation as a
starting point, but rather requires a propagation of the resulting differential
equation for the amplitudes. While the theory is applicable to general classes
of many-body problems in principle, here, the theory is illustrated using
simple model systems. For one-body Fermionic Hamiltonians, Fermi-Dirac one-body
reduced density matrices are recovered for the grand-canonical formulation. For
multidimensional harmonic oscillators numerically exact results are obtained
using the thermal normal ordered exponential (TNOE) approach. As an application
of the related time-dependent formulation numerically exact
time-autocorrelation functions and absorption spectra are obtained for harmonic
Franck Condon problems. These examples illustrate the basic soundness of the
scheme and are used for pedagogical purposes. Other approaches in the
literature are only discussed briefly and no detailed comparative discussion is
attempted.
