Why is the partial differentiation of density operator with respect to time zero, for an ensemble in thermal equilibrium? Sakurai initially says that density operator evolves with time because state kets evolve with time. But for an ensemble in thermal equilibrium, its partial differentiation is zero.
As far as I know, thermal equilibrium means the temperature is same everywhere and no heat flows between the particles. How is it anyway related to time? Why does the density operator stop evolving?
 A: The equality of temperatures is a condition ensuring thermal equilibrium. But equilibrium primarily refers to a stable time independent condition of the system. In particular, all the average values of every observable must be stationary. This can be achieved only by stationary density matrices.
A: $\newcommand{\ket}[1]{\left|#1\right>}$
Here is a not very rigorous explanation. In equilibrium, any macroscopic quantity A should have a constant value. This means the corresponding operator A should have a constant mean. For convenience, I use the same symbol here for both the quantity and the operator.
Suppose that the Hamiltonian operator H is independent of time. From the Schrodinger equation
\begin{equation}
i \hbar \frac{\partial \ket{\psi}}{\partial t} = H\ket{\psi},
\end{equation}
we have
$$\ket{\psi}=\exp \left( -\frac{iHt}{\hbar} \right) \ket{\psi(0)}.$$
Let $\{\ket{u_k}\}$ be a basis which is also a set of eigenkets of H, then we have
$$\langle A\rangle=Tr(\rho A)$$
$$=\sum_{j,k} \langle u_j|\rho|u_k\rangle\langle u_k|A|u_j\rangle.$$
Since $\{\ket{u_k}\}$ are eigenkets of H, we have
$$\exp\left( \frac{iHt}{\hbar} \right) |u_j\rangle = \exp\left( \frac{i E_j t}{\hbar} \right) |u_j\rangle.$$
This gives us
$$\langle u_j|\rho|u_k\rangle = \exp( \frac{i(E_k-E_j)t}{\hbar} ) 
\langle u_j|\rho(0)|u_k\rangle,$$
which yields
$$\langle A\rangle=\sum_{j,k}\exp( \frac{i(E_k-E_j)t}{\hbar} ) 
\langle u_j|\rho(0)|u_k\rangle A_{kj}.$$
Then we can see that the "safest" way to make such an average independent of time for any A is to let $\langle u_j|\rho(0)|u_k\rangle=0$ for all $j \neq k$. This gives
$$\rho=\rho(0)=\sum_j \rho_j |u_j\rangle\langle u_j|$$.
From here, we can see two results at the same time:
(1) $\rho$ is independent of time.
(2) $\rho$ is diagonalized in $\{\ket{u_k}\}.$
From (2), we can see that $\rho$ and $H$ can be simulateously diagonalized, so they commute. This conclusion can also be derived by using (1) and the von Neumann equation
$$i\hbar \frac{\partial \rho}{\partial t} = [H, \rho]$$
A: At thermal equilibrium, the density matrix can be written as
\begin{equation}
\hat{\rho}_{\text{eq}} = \frac{e^{-\beta \hat{H}}}{\mathcal{Z}},
\end{equation}
where $\mathcal{Z} = \text{Tr}\ e^{-\beta \hat{H}}$ is the partition function.
The states are distributed thermally, and the density operator can be written diagonally in the basis of energy eigenstates, so it commutes with the Hamiltonian. Since the time evolution of the density operator is proportional to its commutator with the Hamiltonian, evolution stops if the commutator vanishes.
Eigenstates of the Hamiltonian are stationary states and do not evolve in time. Since the density matrix is diagonal, we simply have different probabilities for different stationary states with no interference between them.
I should also connect how the commutator with $\hat{H}$ relates to what you stated. As you say, the density operator evolves with time because state kets evolve with time. To see this, write the density operator in terms of an outer product of bras and kets; take the time derivative, applying the product rule; and apply the Schrödinger equation (and its Hermitian conjugate for the bras). From this, the commutator naturally pops out. If these bras and kets are stationary states, then all of the time derivatives vanish.
