If quantum Liouville's equation describes nonequilibrium evolution of the density matrix, why its derivation assumes time-independent probabilities? For an ensemble characterized by the states $\{|n\rangle\}$ with probabilities $p_n$ at time $t$, the density operator is defined as $$\rho\equiv \sum\limits_n p_n|n\rangle\langle n|.$$ Assuming that the states $\{|n\rangle\}$ evolve according to the Schrodinger's equation but the probabilities $p_n$ are time-independent, we obtain the quantum Liouville's equation $$i\hslash\frac{\partial\rho}{\partial t}+[\rho,H]_-=0.$$ This equation, in general, describes the nonequilibrium situation and only when, $\partial\rho/\partial t=0$, it describes equilibrium.


*

*But if the density matrix above, in general, is a descriptor of nonequilibrium situation, why $p_n$'s are not taken to be time-dependent?  

 A: 
But if the density matrix above, in general, is a descriptor of nonequilibrium situation, why $p_n$'s are not taken to be time-dependent? 

Quantities $p_n$ are not probabilities of some physical quantity measurement results (which generally change in time), they are weights of representation of the pure states in the ensemble (members of the ensemble that can be described using single ket). For example, if the Hamiltonian operator does not depend on time and $\psi_n$'s are not its eigenfunctions (they could be eigenfunctions of some spin operator), the projectors $|\psi_n\rangle\langle\psi_n|$ will change in time, but $p_n$'s will not.
In other words, changing $p_n$'s in time would mean changing the ensemble in time (weights of pure states $|\psi_n\rangle$ in the mix). Usually one density operator symbol, as a function of time, is supposed to describe only one ensemble.
The definition
$$
\hat{\rho}(t) = \sum_n p_n |\psi_n(t)\rangle\langle \psi_n(t) |
$$
is therefore mathematically a definition of an operator at any time $t$, based on probabilities $p_n$ that are assumed valid for all of the time.
The purpose of this definition is to have a quantity (matrix) that can give expected average of some quantity for this ensemble at any later time $t$.
Note.
The above density operator depends on time purely via dependence of $\psi_n$'s on time and so the usual way to write down the von Neumann equation:
$$
\frac{\partial \hat{\rho}}{\partial t} = \frac{1}{i\hbar}[\hat{H},\hat{\rho}]
$$
is misleading, because there is only one way the time derivative of $\hat{\rho}$ can be calculated - it is the full time derivative. So the preferable way to write the equation is
$$
\frac{d \hat{\rho}}{d t} = \frac{1}{i\hbar}[\hat{H},\hat{\rho}].
$$
So far the above was valid for ensemble with fixed weights of known pure states. This is the case for various spin/light polarization experiments.
In case we want to apply this density matrix idea to a system of many possible states (such as gas in a box), in principle we can:
$$
\hat{\rho}(t) =  \sum_n p_n |\psi_n(t)\rangle\langle \psi_n(t) |
$$
For a real macroscopic system such as gas in a box it is close to impossible to know the states $|\psi_n\rangle$ and their probabilities $p_n$. But if we assume that the gas is initially in thermodynamic equilibrium, we may start with the plausible assumption that the above $\hat{\rho}$ (only at time $t_0$)  may be expressed as
$$
\hat{\rho}(t_0) = \sum_n p_n^{B} |\phi_n\rangle \langle \phi_n|
$$
that corresponds to ensemble where Hamiltonian eigenstates $|\phi_n\rangle$ are represented with Boltzmannian weights $p_n^B$. Then if the Hamiltonian changes in time $t>t_0$, we can use the von Neumann equation with the new Hamiltonian to derive what happens to $\hat{\rho}$. In this time, $\rho$ in general changes in such a way that we no longer express it as
$$
\hat{\rho} =  \sum_n p_n |\psi_n\rangle\langle \psi_n |.
$$
In theory, you could take $\hat{\rho}(t)$ and try to find out which states $|\psi_n\rangle$ are represented and what are their weights $p_n$, but this is not common - there is no simple ensemble of Hamiltonian eigenstates anymore! Instead, we just work with matrix $\rho_{ik}$ in the original Hamiltonian basis and assume the von Neumann equation is valid for it. The density operator is then less of characteric of some ensemble with fixed weights, and more of a theoretical device to describe evolution of expected averages for a single system.
A: There seem to be a few misconceptions underlying this question.
First of all,
$$\hat{\rho} = \sum_n p_n|n\rangle\langle n|$$
is not a general form of the density matrix. It has this form in equilibrium (where $p_n$ are given by the Boltzmann weights) or when it is a density matrix of a pure state in the diagonal basis.
Secondly, one need not assume that $p_n$ are constant in order to obtain the von Neumann equation. In fact, for a pure state (i.e. a state describable by a wave function) this equation is simply an equivalent of the Schrödinger equation:
$$\hat{\rho}(t) = |\psi(t)\rangle\langle \psi(t)|,$$
where
$$i\hbar\partial_t |\psi(t)\rangle = \hat{H}|\psi(t)\rangle,\\
-i\hbar\partial_t \langle\psi(t)| = \langle\psi(t)|\hat{H}
$$
and therefore
$$i\hbar\partial_t \hat{\rho}(t) =  \hat{H}|\psi(t)\rangle\langle \psi(t)| - |\psi(t)\rangle\langle\psi(t)|\hat{H} = [\hat{H},\hat{\rho}(t)].$$
Now, in a special case where the system at time $t=0$ is described by the density matrix $\hat{\rho} = \sum_n p_n|n\rangle\langle n|$ where $|n\rangle$ are eigenstates of the Hamiltonian, the commutator is zero and the density matrix (and therefore $p_n$) is time-independent. Note that this will work also for non-pure states (where von Neumann equation is postulated rather than derived).
