The density operator
To introduce the density operator, let us first consider a pure state, i.e. not a statistical mixture of states.  The initial conditions are completely known.  Let {|u_n>} be an orthonormal eigenbasis of some operator.  Most often this operator is chosen to be H.  Let the initial state of the system at time t be
|ψ_k(t)> =  ∑_nc_n^k(t)|u_n>,  c_n^k(t) = <u_n|ψ_k(t)>.
The system is in a well defined state, but if more than one of the c_n^k(t) are non-zero, it is not in an eigenstate of H, but in a coherent state. 
Let <ψ_k|ψ_k> = 1, then ∑_n|c_n^k(t)|² = 1.

The mean value of an observable A at time t is
<A>(t) = <ψ_k(t)|A|ψ_k(t)> =  ∑_npc_n^k*(t)c_p^k(t)A_np, where A_np = <u_n|A|u_p>.

We now introduce the pure state density operator ρ^k(t) = |ψ_k(t)><ψ_k(t)|,
i.e. the projector onto the state |ψ_k(t)>, and the density matrix with matrix elements
ρ_np^k(t) = <u_n|ρ^k(t)|u_p> = <u_n|ψ_k(t)><ψ_k(t)|u_p> = c_n^k(t)c_p^k*(t).

We will now show that knowing the density operator (or the elements of the density matrix) completely suffices to characterize the state of the system.  All of the predictions of Quantum Mechanics that can be made knowing |ψ_k(t)> can also be made knowing ρ^k(t).

In terms of the density operator ρ^k(t) we express:

• (a) conservation of probability,
  <ψ_k|ψ_k> = ∑_n<ψ_k|u_n><u_n|ψ_k> = ∑_n<u_n|ψ_k><ψ_k|u_n> = ∑_nρ_nn^k(t) = Tr{ρ^k(t)} = 1.
  (d/dt)Tr{ρ^k(t)} = 0.
• (b) the mean value of the observable A,
  <A> =  ∑_npc_n^k*(t)c_p^k(t)A_np = ∑_npρ_pn^k(t)A_np =  ∑_np<u_p|ρ^k(t)|u_n><u_n|A|u_p>
  = ∑_p<u_p|ρ^k(t)A|u_p> = Tr{ρ^k(t)A}.
• (c) the time evolution of the system,
  (d/dt)ρ^k(t) = (d|ψ_k>/dt) <ψ_k| +  |ψ_k>(d<ψ_k|/dt)
  = (1/(iħ)H|ψ_k><ψ_k| - (1/(iħ)|ψ_k><ψ_k|H
  = (iħ)^-1[H,ρ^k(t)].
• (d) the probability that a measurement of A will yield the eigenvalue a_n.
  P_k(a_n) = <ψ_k(t)|P_n|ψ_k(t)> = Tr{ρ^k(t)P_n},
  where P_n is the projector into the eigensubspace of A associated with the eigenvalue a_n.  The probability is the mean value of P_n.
  We have therefore shown that knowing ρ^k(t) suffices to characterize the quantum state of the system.

Why is the density matrix formulation useful?
The important equations,
Tr{ρ^k(t)} = 1,
<A> = Tr{ρ^k(t)A},
P_k(a_n) = <ψ_k(t)|P_n|ψ_k(t) = Tr{ρ^k(t)P_n},
(d/dt)ρ^k(t) = (iħ)^-1[H,ρ^k(t)],
are all linear, while the corresponding equations for the first three of these equations in terms of |ψ_k(t)>,
d|ψ(x)|²/dt = 0,
<A> = ∫ψ*(x) A ψ(x)dx
P_k(a_n) = |<ψ_a|A||ψ>|² = |∫ψ_a*(x) A ψ(x) dx|²
are quadratic.  The equation corresponding to the fourth equation is the Schroedinger equation.

The linear equations can easily be generalized to describe a statistical mixture of states and express all physical predictions about such a system in terms of a density operator.

Note:

• If {|u_n>} is an eigenbasis of H and |ψ_k(t)> is an eigenstate of H, then the matrix ρ_np^k(t) of the pure state density operator is diagonal.  All off-diagonal elements are zero.
  |ψ_k> = |u_k>, ρ_k = |u_k><u_k|.
• If |ψ_k(t)> is a coherent state then ρ_np^k(t) has off-diagonal elements. 
  | ψ_k> = ∑_nc_n|u_n>, ρ^k = ∑_nn'c_n*c_n'|u_n'><u_n|.
• We can always choose {|u_n>} to be an eigenbasis of the density operator itself.
  Then ρ_np^k(t) is diagonal.

In experiments, we are often dealing with systems about which we have incomplete information.

Example:

• A C⁴⁺ beam emerges from an accelerator.  The C ions have been stripped of most of their electrons by having been past through a thin foil.  Most ions are in their ground state, but some are in an excited, metastable state.  For each ion we only know that it may be in a state |ψ_k> with probability p_k,  ∑_kp_k = 1.  We are dealing with a statistical mixture of states.

For such a system we introduce the density operator ρ(t) =  ∑_kp_kρ^k(t), and the density matrix with matrix elements
ρ_pn(t) = ∑_kp_k<u_p|ρ^k(t)|u_n> = ∑_kp_kc_p^k(t)c_n^k*(t).
In terms of the density operator we express:

• (a) conservation of probability,
  Tr{ρ(t)} = Tr{∑_kp_kρ^k(t))} = ∑_kp_kTr{ρ^k(t)) = ∑_kp_k = 1.
• (b) the probability that a measurement of the observable A will yield the eigenvalue a_n,
  P(a_n) = ∑_kp_kP_k(a_n) = ∑_kp_kTr{ρ^k(t)P_n} = Tr{∑_kp_kρ^k(t)P_n} = Tr{ρ(t)P_n},
  where P_n is the projector into the eigensubspace of A associated with the eigenvalue a_n.
• (c) the mean value of the observable A,
  <A>(t)  = ∑_na_nP(a_n) = ∑_na_nTr{ρ(t)P_n} = Tr{ρ(t)∑_na_nP_n} = Tr{ρ(t)A},
  since ∑_na_nP_n = A,  A|ψ> =  ∑_nA|a_n><a_n|ψ> = ∑_na_n|a_n><a_n|ψ> = ∑_na_nP_n|ψ>,
  where {|a_n>} is an eigenbasis of A.
• (d) the time evolution of the system.
  The linearity of (d/dt)ρ^k(t) = (iħ)^-1[H,ρ^k(t)] implies that (d/dt)ρ(t) = (iħ)^-1[H,ρ(t)].

Properties of ρ(t)
ρ(t) is not a projector in general; ρ(t)² ≠ ρ(t) in general.

Example:

• Let |ψ_k> = |u_k> with probability p_k, ρ =  ∑_k'p_k'|u_k'><u_k'|, 
  ρ_kk =  <u_k|∑_k'p_k'|u_k'><u_k'|u_k> = p_k,  ρ_kn(k ≠ n) = 0, Tr{ρ} = ∑_kp_k = 1, Tr{ρ²} = ∑_kp_k²  ≤ 1.

ρ(t) is a Hermitian operator (ρ^† = ρ from the definition of ρ) and
<ψ|ρ||ψ> = ∑_kp_k<ψ|ψ_k><ψ_k|ψ> = ∑_kp_k|<ψ|ψ_k>|² ≥ 0.
ρ(t) is a positive operator.