What is the reduced density matrix that is $|0\rangle$ or $|1\rangle$ with probabilities $1/2$? $
\newcommand{\bra}[1]{\langle{#1}\vert}
\newcommand{\ket}[1]{\vert{#1}\rangle}
$I was trying to solve the following:

Suppose a system is produced in state |0⟩ with probability p0=1/2 and
in state |−⟩ with probability p1=1/2. What is the resulting density
matrix?

I am given several options but my answer is none of ones suggested here
Here is what I tried:
\begin{array}{l}
\ket{\psi }  = \frac{1}{2}\left(\ket{0}  + \ket{-}\right)\\
\rho  = \ket{\psi }\bra{\psi }\\
\frac{1}{2} \times \frac{1}{2}\left(  \left(\ket{0}  + \ket{-}\right)  \otimes  \left(\bra{0}  + \bra{-}\right)\right)\\
\frac{1}{4}\left(  \ket{0}\bra{0}  + \ket{0}\bra{-}  + \ket{-}\bra{0}  + \ket{-}\bra{-}  \right)\\
\ket{0}\bra{0}  = \left[\begin{matrix}
1 & 0\\
0 & 0
\end{matrix}\right]\\
\ket{0}\bra{-}  = 1/\sqrt{2}\left[\begin{matrix}
1 & 1\\
0 & 0
\end{matrix}\right]\\
\ket{-}\bra{0}  = 1/\sqrt{2}\left[\begin{matrix}
1 & 0\\
1 & 0
\end{matrix}\right]\\
\ket{-}\bra{-}  = 1/2\left[\begin{matrix}
1 & 1\\
1 & 1
\end{matrix}\right]
\end{array}
when I sum these 4 matrices none of them is the answer given in the online exercice.
What exactly am I doing wrong?
 A: The mistake is considering probability amplitudes and superposition states instead of probabilities and mixed states. In general, it is incorrect to say that a state of the form
\begin{equation}
|\psi\rangle=\sqrt{p_0}|\psi_0\rangle+\sqrt{p_1}|\psi_1\rangle
\end{equation} has a probability $p_0$ of being found to be in state $|\psi_0\rangle$ and probability $p_0$ of being found to be in state $|\psi_1\rangle$. The only time that such a definitive statement can be made is when the two states on the right hand side are orthogonal; i.e., $\langle \psi_1| \psi_0\rangle=0$. Indeed, when the orthogonality condition is not met, the state $|\psi\rangle$ written here is not normalized for $p_0+p_1=1$.
The correct thing to do is to express our density matrix as an ensemble of pure states. This means that we write down the density matrices that we would have gotten from each of the pure states, and we sum them together with the weights corresponding to their probabilities. So, in this case, we would write
$$\rho=p_0 |0\rangle\langle 0|+p_1 |-\rangle\langle -|=\frac{1}{2} |0\rangle\langle 0|+\frac{1}{2} |-\rangle\langle -|.$$
This type of "convex combination" could be extended to any number of states and any set of probabilities for the states, including a continuous probability distribution for a continuous set of pure states (but that is irrelevant here). If you do this with the states you have calculated, you will get one of the answers from your picture! But I warn you to look into the minus sign in the definition of the state $|-\rangle=(|0\rangle-|1\rangle)/\sqrt{2}$ because it will change the result to another answer from your picture.
Finally, if we take a general density matrix $\rho$, the probability that it is measured to be in state $|\phi\rangle$ is given by \begin{equation}\mathrm{Tr}\left(|\phi\rangle\langle\phi|\rho\right)=\langle\phi|\rho|\phi\rangle.
\end{equation} You cannot directly use this method to verify that you have gotten the correct solution because, even if the state is prepared in $|\psi_0\rangle$, there is still some nonzero probability that it will be meausured to be in any other state $|\phi\rangle$ with $\langle\phi|\psi_0\rangle\neq 0$.
A: One fundamental reason to introduce the density matrix is to describe systems that are classical "mixtures" of states, where each state is associated with a "classical" probability. Such a "mixture" is fundamentally different from a superposition state. Superposition states do not correspond to a classical system where you have a probabilty to find the system in state $|n\rangle$ with probabilty $p_n$.
This is easily illustrated when we look at the time dependent expectation value of an operator of a superposition of energy eigenstates. A superposition state is for example
$$
|\psi\rangle = c_a|a\rangle +c_b|b\rangle 
$$
The time dependent state would be
$$
|\psi(t)\rangle  = c_ae^{-iE_at/\hbar }|a\rangle + c_b e^{-iE_bt/\hbar }|b\rangle
$$
Now lets take a look at the expectation value of an arbitrary hermitian operator $\hat O$
$$\begin{aligned}
\langle \psi(t)|\hat O|\psi(t) \rangle &=\left( c_a^*e^{iE_at/\hbar }\langle a| + c_b^* e^{iE_bt/\hbar }\langle b|\right)|\hat O| \left( c_ae^{-iE_at/\hbar }|a\rangle + c_b e^{-iE_bt/\hbar }|b\rangle \right)\\
\langle \psi(t)|\hat O|\psi(t) \rangle  &= |c_a|^2\langle a|\hat O|a \rangle  + |b|^2\langle b|\hat O|b\rangle + 2\Re\left(c_a^*c_be^{i(E_a-E_b)t\hbar}\langle a|\hat O|b\rangle  \right)
\end{aligned}$$
If we assume for simplicity that $c_a,c_b,\langle a|\hat O|b\rangle \in \Re $ and define $\omega_{ab}=(E_a- E_b)/\hbar $  we can simplify to
$$
\langle \psi(t)|\hat O|\psi(t) \rangle = |c_a|^2\langle a|\hat O|a \rangle  + |b|^2\langle b|\hat O|b\rangle + 2c_ac_b \cos(\omega_{ab} t)  \langle a|\hat O|b\rangle 
$$
This result for a time dependent expectation value of a superposition looks a bit like the sum of the expectation value associated with state a $\langle a| \hat O|a\rangle\equiv O_a $ weighted with probabilty $p_a \equiv |c_a|^2$ plus the the expectation value of state b $O_b$ weighted with $p_b$ if it wasn't for the time dependent oscillating factor
$2c_ac_b \cos(\omega_{ab} t)  \langle a|\hat O|b\rangle $ that causes the expectation value to oscillate and vary with time. This is contrary to a classical system where you have fixed probabilities of $p_a$ and $p_b$. With "classical" probabilities you wouldn't expect any oscillation with time of the expectation value and rather assume that the expectation value is simply the weighted sum of expectation values $\text{"classical sum"} = |c_a|^2O_a + |c_b|^2O_b$.
Which brings us to the important point, namely that it is impossible to describe such a system with a single state, also called pure state,
$$\begin{aligned}
\langle \psi(t)|\hat O|\psi(t)\rangle  &\not = |c_a|^2O_a + |c_b|^2O_b\\
 & \not= p_a O_a + p_b O_b\\
\langle \psi(t)|\hat O|\psi(t)\rangle  &\not= |c_a|^2\langle a|\hat O|a\rangle + |c_b|^2\langle b|\hat O|b\rangle
\end{aligned}$$
which yields an expectation value that looks like it was formed by a "classical" sum of probability weighted states.
To describes systems that are simple sums of classical probabilities we need the formalism of density matrices/operators. Only within this formalism we can model quantum systems that are "classical" probability weighted sums.
The problem with your ansatz is that you have constructed only superposition states and you haven't made use of the property of density matrices to describe "classical" probabilities.
The density matrix formalism allows us to construct probability weighted sums like this
$$
\hat \rho = \sum p_n |n\rangle \langle n|
$$
This allows us to construct a density matrix in a basis of two states like this,
$$
\rho_{\text{classical}} = \left( \begin{matrix}|c_a|^2 & 0\\ 0 & |c_b|^2 \end{matrix}\right)
$$
while a single superposition state will always lead to matrix of the form
$$
\rho_{\text{pure state}} = \left( \begin{matrix}|c_a|^2 & c_ac_b^*\\ c_a^*c_b & |c_b|^2 \end{matrix}\right)
$$
I hope this helps you to understand the question, which is trivial to answer if you get the basics and the "why" of density matrix/operator formalism.
