What is the intuition behind density matrix? What gives density matrix the expressive power to be able to represent mixture of pure states?
For example if $|\Omega\rangle$ is 50-50 mixture of $|\Psi\rangle=\frac{1}{\sqrt{2}}(|u\rangle+|d\rangle)$ and $|\Phi\rangle=\frac{1}{\sqrt{2}}(|u\rangle-|d\rangle)$ we can not represent such as $|\Omega\rangle=\frac{1}{2}(|\Psi\rangle+|\Phi\rangle)$, but we can represent it as $\rho=\frac{1}{2}(|\Psi\rangle\langle\Psi|+|\Phi\rangle\langle\Phi|)$.
Intuitively what makes $|\Psi\rangle\langle\Psi|$ able to represent such as opposed to $|\Psi\rangle$?
 A: Here's one way you can think of it - we will specialize to a finite dimensional Hilbert space $\mathcal H$ for simplicity.  Abstractly, the state of a quantum system should be identified with a probability distribution for every possible measurement one could make.
A given self-adjoint operator $A$ can be expanded as $A=\sum_i \lambda_i \Pi_i$, where $\lambda_i$ is the $i^{th}$ eigenvalue and $\Pi_i$ the orthogonal projection operator onto the $i^{th}$ eigenspace.  With this in mind, a ray $\boldsymbol \Psi$ in the underlying Hilbert space defines a probability distribution as follows. Pick any vector $\psi$ from that ray; then the probability of measuring $A$ to have value $\lambda_i$ is simply
$$\frac{\langle \psi, \Pi_i \psi\rangle}{\langle\psi,\psi\rangle}$$
With this as motivation, we define an event as an orthogonal projector $\Pi$, and the probability of that event to be
$$P_\Psi(\Pi):= \frac{\langle \psi,\Pi \psi\rangle}{\langle \psi,\psi\rangle} \in [0,1]$$
In that sense, $P_\Psi$ constitutes a probability measure on the set $\mathcal P(\mathcal H)$ of orthogonal projection operators. We now ask whether it is the most general form of probability measure; this question is answered in the negative by Gleason's theorem, which says that the set of all probability measures is in one-to-one correspondence with the set of positive-semidefinite, self-adjoint operators $\rho$ with $\mathrm{Tr}(\rho)=1$ (i.e. density operators), with the corresponding probability measure being given by $\Pi \mapsto \mathrm{Tr}(\Pi \rho)$. The probability measures which arise from rays in $\mathcal H$ - which correspond to so-called pure states - form a strict subset of all probability measures.

The proof of Gleason's theorem is technical and beyond the scope of a PSE answer in my opinion. However, it is useful to see how the set of pure states is insufficient.  Consider $\mathcal H= \mathbb C^2$ - does there exist a state of the system whose probability measure is uniform, such that for any observable $A$ with two distinct outcomes (eigenvalues) the probability of each is 1/2?
If we restrict ourselves to pure states, then the answer is no which can be seen easily.  Consider $P_\Psi$ for some arbitrary ray $\Psi$, let $\psi$ be a normalized vector in that ray, and let $\phi$ be a normalized vector in the orthogonal complement of $\Psi$.  Then $\{|\psi\rangle\langle\psi|,|\phi\rangle\langle\phi|\}$ is a set of orthogonal projectors which sum to the identity, yet $P_\Psi(|\psi\rangle\langle\psi|)=1\neq 1/2$ and $P_\Psi(|\phi\rangle\langle\phi|) = 0\neq 1/2$.
As a result, if we want to generalize to allow for such statistical ensembles, we must enlarge the state space from the set of rays in $\mathcal H$ to the set of density operators on $\mathcal H$. It's not hard to show that the latter objects can be written $\rho = \sum_i p_i |\psi_i\rangle\langle \psi_i|$ where $\sum_i p_i = 1$ and $\{\psi_i\}$ is some collection of unit vectors in $\mathcal H$.
A: $|\Omega\rangle=\frac{1}{\sqrt{2}}(|\Psi\rangle+|\Phi\rangle)$ is a superposition state, but it's a state on its own and will evolve accordingly. If it described a system undergoing two-slit diffraction, with $|\Psi\rangle$ being the wavefunction of one slit and $|\Phi\rangle$ that of the other one, the evolution of $|\Omega\rangle$ would give you the two-slit interference pattern, because both slits are contributing. The particle goes through both slits, since it's a superposition.
$\rho=\frac{1}{2}(|\Psi\rangle\langle\Psi|+|\Phi\rangle\langle\Phi|)$ is a statistical mixture, meaning that the system is either in $|\Psi\rangle$ or in $|\Phi\rangle$. That is you know that the particle has gone through one of the two slits, you just don't know which one. The density matrix allows to evolve the two states at the same time, but each of them is evolving as a single state, not as a superposition.
