Strong thermalization of the two-component Bose-Hubbard model at finite temperatures

We study thermalization of a two-component Bose-Hubbard model by exact diagonalization. Initially the two components do not interact and are each at equilibrium but with different temperatures. As the on-site inter-component interaction is turned on, perfect thermalization occurs. Remarkably, not merely those simple "realistic" physical observables thermalize but even the density matrix of the \textit{whole} system---the time-averaged density matrix of the system can be well approximated by that of a canonical ensemble. A conjecture about this fact is put forward.