  We introduce a non-linear differential flow equation for density matrices
that provides a monotonic decrease of the free energy and reaches a fixed point
at the Gibbs thermal state. We use this equation to build a variational
approach for analyzing equilibrium states of many-body systems and demonstrate
that it can be applied to a broad class of states, including all bosonic and
fermionic Gaussian states, as well as their generalizations obtained by unitary
transformations, such as polaron transformations, in electron-phonon systems.
We benchmark this method with a BCS lattice Hamiltonian and apply it to the
Holstein model in two dimensions. For the latter, our approach reproduces the
transition between the BCS pairing regime at weak interactions and the
polaronic regime at stronger interactions, displaying phase separation between
superconducting and charge-density wave phases.
