If I have a quantum mechanical system (QMS) whose density matrix is at temperature $T$ and I put it in contact with another quantum mechanical system (QMS2) at temperature $T$ again.
However, if I try to model this dynamically I'd be tempted to use the sudden approximation. Let the probability of QMS being in energy state $| E_i \rangle$ is $e^{- \beta E_i}$ and QMS2 being in energy $| E_j \rangle$ is $e^{- \beta E_j}$. Hence, the probability of the below outcome
$$ \Big (H_{QMS} \otimes 1+ 1 \otimes H_{QMS2} \Big ) |E_i ,E_j\rangle \to \Big (H_{QMS} \otimes 1 + 1 \otimes H_{QMS2} + H_{int}\Big )|\lambda_l ,\lambda_k \rangle$$
(where $H_{int}$ is the interaction Hamiltonian and $|\lambda_l ,\lambda_k \rangle$ are arbitrary states of the combined system) is given by:
$$ e^{-\beta(E_i +E_j)} |\langle E_i | \lambda_l \rangle \langle E_j | \lambda_k \rangle|^2 $$
However, by the zeroth law of thermodynamics I know that these systems should remain in thermal equilibrium as well. I'm uncertain as how to show that is the case here?