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.