Is a mixture of Bell states entangled? Setup: 
I have a two qubit system $\mathbb{C}^2 \otimes \mathbb{C}^2$ an we consider the Bell states
$\newcommand{\ket}[1]{\left|#1\right>}$
$\newcommand{\bra}[1]{\left<#1\right|}$
$\DeclareMathOperator{\spn}{span}$
$$\begin{aligned}
\ket{\Phi^+} &= \frac{1}{\sqrt{2}} ( \ket{00} + \ket{11} ) & \ket{\Psi^+} &= \frac{1}{\sqrt{2}} ( \ket{01} + \ket{10} ) \\
\ket{\Phi^-} &= \frac{1}{\sqrt{2}} ( \ket{00} - \ket{11} ) & \ket{\Psi^-} &= \frac{1}{\sqrt{2}} ( \ket{01} - \ket{10} )
\end{aligned} $$
I decompose the Hilbert space into
$$\mathbb{C}^2 \otimes \mathbb{C}^2 = \mathcal{H}_+ \oplus \mathcal{H}_- $$
where $\mathcal{H}_+ = \spn \{ \ket{\Phi^+}, \ket{\Phi^-}, \ket{\Psi^+} \}$ and $\mathcal{H}_- = \spn \{ \ket{\Psi^-} \}$, the spaces of symmetric and antisymmetric vectors respectively. By $P_+$ and $P_-$ we denote the projections onto these subspaces and $\rho_+ = \frac{1}{3} P_+$ and $ \rho_- = P_- = \ket{\Psi^-}\bra{\Psi^-}$ are two states. Clearly, $\rho_-$ is entangled.
Questions:
Is the projection $P_+ = \ket{\Phi^+}\bra{\Phi^+} + \ket{\Phi^-}\bra{\Phi^-} + \ket{\Psi^+}\bra{\Psi^+}$? 
Is $\rho_+$ separable? According to what I found (for example, this wikipedia page), it should be, but what is the decomposition into product states? Also, how do I interpret this? It starts as a mixture of entangled states, but it turns out to be separable? Intuitively, that feels off to me.
 A: Notice that a projector onto the symmetric states is, up to normalization, a Casimir of the symmetric representation of $SU(2)$. The quadratic Casimir is $\sum_a T_{\rm{product}}^{a} T_{\rm{product}}^a = J^2$. This gives $j(j+1)$ - it will annihilate the anti-symmetric ('spin $0$') states. Hence, up to normalization this is just the projector we're looking for!
So can this be 'separated'?
Note that in terms of fundamental generators  $T^{a}_{\rm{product}} = T^a \otimes I + I \otimes T^a $ where $T^a = \sigma^a /2$ (a Pauli matrix). Hence: $$ J^2 = \sum_{a=1}^{3} \left(T^a T^a \otimes I+ I \otimes T^a T^a +2 T^a \otimes T^a \right)\\= \frac{3}{2} I\otimes I  +2 \sum_{a=1}^{3} T^a \otimes T^a \\=  \sum_{a=1}^{3}\left((\frac{I}{2}-T^a)\otimes(\frac{I}{2}-T^a)+(\frac{I}{2}+T^a)\otimes(\frac{I}{2}+T^a) \right)  $$
Notice that each of the factors in each term is a Hermitian matrix with non-negative Eigenvalues and trace 1, making them all legal density matrices. I arrived at this precisely by trying to construct such density matrices. In fact they are pure-states - projections unto 'up' or 'down' along each of the 3 axes.

"Also, how do I interpret this?"

Less sure about this one. I suspect the issue is that while you can prepare this state by 'preparing 3 entangled states with equal probability', measuring one spin can still give you no more information about the other than is achievable without entanglement, and so there is no 'observable' entanglement.
By this I don't mean that there is entanglement. There isn't. I just mean that what we understand as entanglement must be a physical observable that can be experimentally verifiable, as opposed to an artifact of writing the density matrix in a certain way.
