A path-integral based non-perturbative thermal cluster cumulant (PI-TCC) method is presented for thermal averages of quantum many-body hamiltonians. This provides a PI-expression (y) of the quantum fluctuation component of the trace of Boltzmann operator with a gaussian measure, involving the fluctuation variable y around each point x(0) of the classical centroid path. The T-ordered products of y-variables in U, are evaluated efficiently via a 'Wick's theorem' for the y-variables, which enables us to express (y) as the y-average of a y-normal ordered exponential, where only the number component of cumulant survives. The full gpf is obtained as an integral over x(0). The gpf of the quartic an harmonic oscillator is computed as an illustration, using an optimal reference frequency Ohm(x(0)) at each x(0).