Gaussian resolutions for equilibrium density matrices

A Gaussian resolution method for the computation of equilibrium density matrices rho(T) for a general multidimensional quantum problem is presented. The variational principle applied to the ``imaginary time'' Schroedinger equation provides the equations of motion for Gaussians in a resolution of rho(T) described by their width matrix, center and scale factor, all treated as dynamical variables. The method is computationally very inexpensive, has favorable scaling with the system size and is surprisingly accurate in a wide temperature range, even for cases involving quantum tunneling. Incorporation of symmetry constraints, such as reflection or particle statistics, is also discussed.