Variational Approach for Many-Body Systems at Finite Temperature

We introduce an equation for density matrices that ensures a monotonic decrease of the free energy and reaches a fixed point at the Gibbs thermal. We build a variational approach for many-body systems that can be applied to a broad class of states, including all bosonic and fermionic Gaussian, as well as their generalizations obtained by unitary transformations, such as polaron transformations in electron-phonon systems. We apply it to the Holstein model on 20 x 20 and 50 x 50 square lattices, and predict phase separation between the superconducting and charge-density wave phases in the strong interaction regime.