  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.
