Also another way to think about is the following (here I follow the 5th Volume of Landau&Lifshitz). Let $\rho$ be a density matrix of a system, and the system can be subdivided into 2 parts, described by denisty matrices $\rho_1, \rho_2$ with weak interaction, such that they can be considered independent, but at the same time there is thermodynamical equilibrium between them. Therefore the density matrix is product of density matrices of each subsystem:
$$
\rho = \rho_1 \rho_2 \Rightarrow \ln \rho = \ln \rho_1 + \ln \rho_2
$$
Therefore $\ln \rho$ is an additive integral of motion. The system has additive following integrals of motion $E, P, J$ - total energy, momentum, angular momentum, respectively. So, $\ln \rho $ has to be some linear combination of them. Going to the frame of reference, which is comoving and corotating with the system concerned, we may set $P, J = 0$. $\beta$, the inverse temperature $1/T$, has the meaning of proportionality coefficient.
$$
\ln \rho = -\beta E
$$