Is the density operator a mathematical convenience or a 'fundamental' aspect of quantum mechanics? In quantum mechanics, one makes the distinction between mixed states and pure states.  A classic example of a mixed state is a beam of photons in which 50% have spin in the positive $z$-direction and 50% have spin in the positive $x$-direction.  Note that this is not the same as a beam of photons, 100% of which are in the state
$$
\frac{1}{\sqrt{2}}[\lvert z,+\rangle + \lvert x,+\rangle].
$$
It seems, however, that at least in principle, we could describe this beam of particles as a single pure state in a 'very large' Hilbert space, namely the Hilbert space that is the tensor product of all the $\sim 10^{23}$ particles (I think that's the proper order of magnitude at least).
So then, is the density operator a mathematical convenience, or are there other aspects of quantum mechanics that truly require the density operator to be a 'fundamental' object of the theory?
(If what I mean by this is at all unclear, please let me know in the comments, and I will do my best to clarify.)
 A: My impression from the literature is that physicists are still divided on this question.  The Quantum Information Theory camp says the latter, but the Quantum Optics people say the former.  
A related, but distinct, issue is whether one regards the concept of «open system» as a mere mathematical convenience, or as a fundamental concept.  This issue has received less attention.  It is striking that in Dirac's basic textbook, «system» means «closed system», just as it always did in Hamiltonian mechanics: there is no other kind of system.  Obviously in order to do some practical calculations in physical situations where dissipation or decoherence plays a role, the device of an open system is convenient.  The issue of whether it is fundamental is still not settled.
Both issues are related to one's sympathies or lack thereof with «decoherence» as the answer to the puzzle of quantum measurement.  In spite of some claims to the contrary, decoherence is not yet accepted by consensus.  See Stephen Adler's well argued paper pointing out the serious problems with decoherence when taken as the solution to the quantum measurement problem (everyone agrees that decoherence is sometimes a convenient mathematical device) 

Why Decoherence has not Solved the Measurement Problem: A Response to P. W. Anderson. Stephen L. Adler. Stud. Hist. Philos. Sci. Part B 34 no. 1, pp. 135-142 (2003). doi:10.1016/S1355-2198(02)00086-2, arXiv:quant-ph/0112095.

See also my own (which sometimes, on a bad day, doesn't even convince me) 

Thermodynamic Limits, Non-commutative Probability, and Quantum Entanglement. Joseph F. Johnson. In Proceedings of the 3rd International Symposium, Quantum Theory and Symmetries (Cincinnati, USA, 10-14 September 2003). arXiv:quant-ph/0507017

and the excellent paper by the great Roger Balian and two assistants, 

The quantum measurement process in an exactly solvable model. Armen E. Allahverdyan, Roger Balian and Theo M. Nieuwenhuizen. AIP Conf. Proc. 750, pp. 26-34 (Foundations of probability and physics - 3, Vaxjo, Sweden, 7-12 June 2004). doi:10.1063/1.1874554 arXiv:cond-mat/0408316

all of which should be studied in order to make up your mind on the subject.  Balian, as all statmech people, relies heavily on the use of the density matrix without enquiring too closely as to its foundational status.  But since there is real important physics in their paper, it (and its precursor zipped files posted at http://www.chicuadro.es/ , H. Green, Nuovo Cimento 9 (1958), 880, written under Schroedinger's influence while at Dublin---but note that decades later, having long escaped Schroedinger's influence, Green adopted an opposing viewpoint) is of key importance to answering your question.
A: Your question is ambiguous. If you are asking if the density operator $\hat{\rho}$ formalism is more fundamental than the state vector $|\Psi\rangle$ formalism, the response is yes, because the density operator formalism applies to open quantum systems for the which no state vector exists. Moreover, it has been shown in last decades that there exists a kind of isolated quantum systems (LPSs) for the which the equivalence between density operators and state vectors is broken and a new kind of solutions (outside the Hilbert space) appear which only can be described in terms of density operators. Some recent Solvay conferences were devoted to such systems.
If you are asking if the density operator formalism is the only way to study such systems (LPS, open...), the response is no, because such systems can be studied using alternatives to the density operator formalism: e.g. in the Wigner-Moyal formulation an open quantum system is described by a Wigner distribution $\rho_W$. However, notice that in both cases the Wigner-Moyal formalism can describe systems beyond the scope of the state vector formalism.
A: No macroscopic quantum system is described by a pure state. For example, notions like temperature or pressure, which apply to macroscopic systems do not even exist for systems described by a pure state. The description of macroscopic objects (discussed in statistical mechanics) is always in terms of a density operator (or the essentially equivalent notion of a positive linear functional on a suitable $C^*$-algebra of observables.
However, what is considered as fundamental depends on the fundamentals upon which one erects QM. In the standard textbook foundations, the density operator is therefore not fundamental. However, the axioms for QM become much nicer if one takes the density operator as fundamental; see the topic ''Postulates for the formal core of quantum mechanics'' in Chapter A1: Fundamental concepts in quantum mechanics of my theoretical physics FAQ at http://arnold-neumaier.at/physfaq/physics-faq.html
Moreover, when analyzing actual experiments on the foundations of quantum mechanics, one must represent all states as mixed states to get correct models. These models are in the Markov approximation of the form of a Lindblad equation for the state of the system, and in the general case of a similar nature but with non-Markovian memory terms. The latter gives the exact dynamics. In both cases, the dynamics is that of a density operator, never that of a wave function.
The usual wave function formalism appears just to be an idealization in which one ignores the effect of the environment (rigorously allowed only at the temperature 0K) and thus arrives at marginally simpler equations.
[added Nov.20,2016:] We cannot really decide whether the universe is described by a wave function or a density operator, since we have far too little information to pin it down. What we observe points however to a density operator, since a wave function has no temperature, while the universe has one. In any case, any proper subsystem of the universe must be obtained by tracing out the deleted part and hence is given by a density operator, except if the subsystem is so tiny that we can force it into an eigenstate of a complete set of observables.
A: Let me try to convince you that the density operator is a mathematical convenience and not a fundamental aspect of quantum mechanics by describing a very general setup for states and observables in both classical and quantum mechanics. This may not directly answer your question, but hopefully it will settle whatever motivated this question. 
Briefly,


*

*Observables $A$ are a Poisson *-algebra. This means a complex algebra $A$ equipped with an antilinear involution $a \mapsto a^{\ast}$ such that $(ab)^{\ast} = b^{\ast} a^{\ast}$ as well as a Lie bracket $\{ -, - \} : A \otimes A \to A$ which is also a derivation with respect to multiplication (and maybe some compatibility between these two structures). Classical examples occur when $A$ is the algebra of complex-valued smooth functions on a symplectic manifold $M$, the involution is pointwise conjugation, and $\{ -, - \}$ is the usual Poisson bracket, and quantum examples occur when $A$ is the algebra of linear operators on a Hilbert space $H$, the involution is adjoint, and $\{ -, - \}$ is the commutator. 

*States are *-linear functionals $\mathbb{E} : A \to \mathbb{C}$ on $A$ such that $\mathbb{E}(1) = 1$ and such that $\mathbb{E}(a^{\ast} a) \ge 0$ for all $a$. Classical examples occur when $\mathbb{E}$ is integration against a probability measure on a symplectic manifold $M$. Pure quantum examples occur when $\mathbb{E}(a) = \langle \psi, a \psi \rangle$ for some unit vector $\psi$ in some Hilbert $^{\ast}$-representation of $A$ and you can get more examples by taking linear combinations or more generally integrals of these.


I wrote down what amounts to a justification of this formalism in a somewhat long series of blog posts:
http://qchu.wordpress.com/2011/07/16/the-heisenberg-picture-of-quantum-mechanics/
http://qchu.wordpress.com/2011/08/14/poisson-algebras-and-the-classical-limit/
http://qchu.wordpress.com/2012/08/18/noncommutative-probability/
http://qchu.wordpress.com/2012/09/09/finite-noncommutative-probability-the-born-rule-and-wave-function-collapse/
The Poisson structure is only necessary to understand time evolution; if you just want to understand states, you can safely ignore the first two posts. 
Density operators come into the picture as follows. $A$ is sometimes equipped with a canonical linear functional (which may not be everywhere defined); a normalized version of the corresponding state (when it exists) should be thought of as the "uniform distribution." Classical examples occur when this linear functional is given by integration against Liouville measure on a symplectic manifold $M$, and quantum examples occur when this linear functional is given by trace. If $\text{tr}$ denotes this functional, then you can use it to write down a distinguished class of states of the form
$$\mathbb{E}(a) = \text{tr}(\rho a)$$
for some $\rho \in A$. This is the density operator of the state. If $\text{tr}$ satisfies a suitable nondegeneracy condition, $\rho$ can be uniquely recovered from the state it defines. The axioms defining a state require $\rho$ to be self-adjoint and have trace $1$, and it must also have suitable positivity properties (for example in the case of a finite-dimensional matrix algebra it must have non-negative eigenvalues). 
You are absolutely correct that every state can be described using a vector in a suitably large Hilbert space (here I am using "Hilbert space" in the physicist's sense, which I understand to really be "inner product space"). This is a corollary of a version of the Gelfand-Naimark-Segal construction which is explained in the third post above. 
However, it is worth mentioning that there is an intrinsic definition of pure state in the theory of operator algebras: namely, the space of states has a natural convex structure and you can define a pure state to be an extreme point of the space of states. 
So, summary: what is fundamental is a linear functional on the algebra of observables, and everything else comes down to finding convenient ways to write down and analyze such linear functionals. 
A: This is an interesting as still open question. 
@Quiaochu Yan The point is that we are using a mathematical model to understand how the world is working. We are using vectors and matrices to make prevision of some feature of the nature. You cannot answer to the question "is the mathematical model that we are using describing completely the nature or not?" staying in the mathematical model! In this case experimental observations are needed.
What physicist are trying to do is to extend quantum mechanics, removing for instance one axioms, and to look if it is still a good mathematical theory, that is if it still in same way  make right prevision. This is analog to how we have proven that classical mechanics is not complete. These trials are in the same trend to understand if a density matrix is a fundamental aspect of quantum mechanics or not.
A: I'll throw in two "classical" citations from  Landau-Lifshitz, book 5, chapter 5:

The averaging by means of the statisitcal matrix ... has a twofold nature. 
  It comprises both the averaging due to the probalistic nature of the quantum description (even when as complete as possible) 
  and the statistical averaging necessiated by the incompleteness of our information concerning the object considered....
  It must be borne in mind, however, that these constituents cannot be separated; the whole averaging procedure is carried 
  out as a single operation, and cannot be represented as the result of succesive averagings, 
  one purely quantum-mechanical and the other purely statistical.

And then this:

It must be emphasised that the averaging over various $\psi$ states, which we have used in order to 
  illustrate the transition from a complete to an incomplete quantum-mechanical description has only 
  a very formal significance. In particular, it would be quite incorrect to suppose that the description 
  by means of the density matrix signifies that the subsystem can be in various $\psi$ states with various 
  probabilities and that the average is over these probabilities. Such a treatment would be in conflict 
  with the basic principles of quantum mechanics.

Taking these into account, I'd say that density matrix is actually more "fundamental" than just pure states description. 
