Thermalization in a Hartree ensemble approximation to quantum field dynamics

For homogeneous initial conditions, Hartree (gaussian) dynamical approximations are known to have problems with thermalization, because of insufficient scattering. We attempt to improve on this by writing an arbitrary density matrix as a superposition of gaussian pure states and applying the Hartree approximation to each member of such an ensemble. Particles can then scatter via their back-reaction on the typically inhomogeneous mean fields.
Starting from initial states which are far from equilibrium we numerically compute the time evolution of particle distribution functions and observe that they indeed display approximate thermalization on intermediate time scales by approaching a Bose-Einstein form. However, for very large times the distributions drift towards classical-like equipartition.