  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.
