Why do we need mixed states in quantum mechanics? I am trying to understand the necessity of density matrices and the notion of "mixed states" in quantum mechanics (I read all the other posts about this, I promise).
As far as I understand, one could motivate these notions as follows:
Let $H_1$ and $H_2$ be Hilbert spaces of two quantum systems $S_1$ and $S_2$. Moreover, let $H_1 \otimes H_2$ be the Hilbert space of the composite system $S$. Now, let $x \in H_1 \otimes H_2$ be a pure state, by which I mean just some arbitrary (unit) vector. If $x$ is entangled (i.e. cannot be written as $x=x_1\otimes x_2$) then we cannot in a reasonable way reduce $x$ to pure a state (again, read (unit) vectors) in $H_1$ or $H_2$. A way out of this is to introduce density matrices in order to relax the "state"-notion just enough that the "reducing to subsystem" is well-defined by taking the partial trace.
My question is now: Why would we want to "reduce" to $S_1$ or $S_2$ in the first place? If we want to measure some $S_1$-observable $A$ in $S$, we could just apply $A \otimes I$ to $x$ itself which should give the same expectation values etc., right? Of course, I'm hiding here that $x$ itself could be a density matrix but I feel like my argument works inductively. Are there other benefits of introducing "mixed states" which I am not seeing?
 A: In addition to what has already been said/written, you need to consider mixed states, as otherwise the entropy of your system is always either zero (like von Neumann entropy) or ill-defined, as pure states are by definition the solutions of a perfectly time-reversible, deterministic, unitary Schroedinger equation.
A: Philosophically, the reasoning is exactly the same as why it's often very useful to use classical statistical mechanics instead of just treating every system as a single point in a gigantic many-body classical phase space. As a practical matter, we can rarely exactly determine the value of all of the relevant degrees of freedom in a system - and moreover, many of those degrees of freedom are often unimportant anyway for the kinds of questions that we want to ask.
(I'll always remember a great comment from my professor in a stat mech class: "Everyone always says that the reason to use quantum stat mech is that the Schrodinger is too hard to solve exactly for a many-body system. But that's far from the only reason. Suppose that someone magically could give you a 500-exabyte hard drive storing an exact many-body wavefunction for a macroscopic object. What would you do with it?")
Due the locality of the laws of physics, it's very hard to measure highly nonlocal observables. So it's very difficult, even in principle, to fully reconstruct a wavefunction for a large system. But we can often measure local quantities, like the values of single-qubit observables. It's extremely convenient to have a mimimal description of a subsystem that only describes the outcomes of those local (physically realistic) observables with a minimal amount of information (i.e. the minimal number of degrees of freedom). The density-matrix formalism does that beautifully. In a completely nonlocal universe where nonlocal obervables are just as natural to consider as local ones, you might well be right that density matrices wouldn't be that useful - but that is not the universe that we live in.
More practically, density matrices are just really, really useful in practice. But like many (most?) useful physics concepts, it's very hard to explain why they're useful in the abstract. You need to get your hands dirty and work out a lot of problems or do a lot of experiments. That's the only way to truly appreciate why it's often extremely useful to be able to ignore (or, in fancy probability language, "marginalize over") degrees of freedom that are difficult or impossible to determine in practice.
A: Imagine that Alice has a qubit in the state $$|\psi\rangle = \alpha |0\rangle + \beta|1\rangle.$$ She measures her qubit in the computational basis $\{|0\rangle, |1\rangle\}$, obtaining the outcome $|0\rangle$ with probability $p_0 = |\alpha|^2$ and the outcome $|1\rangle$ with probability $ p_1 = |\beta|^2= 1-|\alpha|^2$. She now hands the qubit to Bob, without telling him the measurement outcome. What state does Bob assign to the qubit?
There is no pure state that Bob can write down that can describe the system in his hands.* But, assuming Bob knows that the initial state was $|\psi\rangle$ and that the measurement basis was $\{|0\rangle, |1\rangle\}$, he can write down the mixed state
$$ \rho = p_0 |0\rangle\langle 0| + p_1 |1\rangle \langle 1|,$$
which, unlike a pure state, faithfully predicts the outcome of any measurement Bob cares to do.
Any workable theory of physics must incorporate ignorance: one does not need to invoke mischievous human agents like Alice to see this. Real measuring devices are imperfect, and only give us results with some uncertainty. The only feasible way of describing this uncertainty in the quantum state of the measured system is with a mixed state. For exactly the same reason, we use probability theory when analysing experiments in classical physics.

*That is, if this procedure is repeated many times – Alice repeatedly hands Bob identically prepared qubits that he can measure in any way he likes – there is no pure state that will accurately describe the outcome of all measurements Bob could choose to do. To see this, imagine that Bob tries to assign the pure state $$|\phi\rangle = a|0\rangle + b|1\rangle.$$ Bob now measures in the basis $\{|\phi\rangle, |\phi^\perp\rangle\}$, where  $|\phi^\perp\rangle = b|0\rangle - a|1\rangle$ is orthogonal to $|\phi\rangle$. He expects to obtain outcome $|\phi\rangle$ with probability one and never obtain outcome $|\phi^\perp\rangle$, due to quantum interference. Alas, feckless Bob will not see interference: he obtains outcome $|\phi^\perp\rangle$ with probability
$$ p_{\rm failure} = p_0|b|^2 + p_1|a|^2.$$
Clearly, $p_{\rm failure}\neq 0$ unless $p_0=0,1$ and $a = 0,1$.
A: There are a lot of reasons.  Here's one.  We sometimes want to couple a small quantum system to a very large one, like a two-level atom interacting with the (quantized) electromagnetic field (perhaps in its vacuum state so that we are modeling spontaneous emission).  The large system has many degrees of freedom, and so writing down the coupled-system Hamiltonian and using that to solve for anything like the dynamics of the coupled system is impossible.
However, under suitable conditions (weak coupling, etc.), we can make approximations that lead to an effective equation of motion for just the system degrees of freedom.  Through the interaction of the system with the large environment, there is necessarily decoherence, and so we have to treat the system's state via a density operator.  Keyword: quantum master equation.
A: The counterpart of density matrix in classical physics is probability distribution.
In classical physics one uses the description in terms of probability distributions if one describes open system dynamics, where the initial conditions for the environmental degrees of freedom are not exactly known and one can only give some statistical information about it (eg. by Brownian motion one could specify the temperature of the solution i.e. give information about average velocity of the solution's molecules).
This is also true in quantum physics, i.e. if you do not have exact knowledge about some degrees of freedom, you can describe it in terms of density matrix, where you build a mixture of all possible states with a corresponding probabilities
$\rho = \sum\limits_i p_i | \Phi_i\rangle \langle \Phi_i |$.
However, in quantum physics the density operator is even more fundamental than probability distributions in classical physics (at least in this context), as even by starting from the pure system and environmental state
$|\Phi_{tot}(0) \rangle =|\Phi_{S}(0) \rangle|\Phi_{E}(0) \rangle $
(which would in classical physics corresponding to exact knowledge of initial conditions for open system and environment), the total state can be entangled at some time
$|\Phi_{tot}(t) \rangle \neq|\Phi_{S}(t) \rangle|\Phi_{E}(t) \rangle $.
By tracing out the environmental degrees of freedom (partial trace), one obtains the reduced density operator for open quantum system, which is mixed state. Note this is due entanglement between system and environment, i.e. is purely quantum phenomenon.
A: There is a very simple reason for introducing the density matrix: for some systems the Hamiltonian doesn't exist, so there is no Schrödinger equation and no wavefunction or state vector.
The common occurence of this is thermodynamics, where the system $S$ is interacting with some environment $E$, which we obviously can't describe microscopically but, at beast, by macroscopic quantities like temperature, pressure, etc. The Hamiltonian of $S$ can't be defined because its energy depends on the state of $E$. Classically the coordinates of $S$ don't pin down the interaction terms, so the function $H(p,q)$ can't be defined.
