Majorana-like representation for mixed symmetric states? Is there a generalization of the Majorana representation of pure symmetric $n$-qubit states to mixed states (made of pure symmetric $n$-qubit)?
By Majorana representation I mean the decomposition of a state 
$$|\psi\rangle = \text{normalization} \times \sum_{perm} |\eta_1\rangle |\eta_2\rangle \cdots |\eta_n\rangle,$$
where $|\eta_k\rangle$ are uniquely determined (up to a global phase in each and the permutation) qubit states. Typically, they are presented as points on the Bloch sphere.
When it comes to a desired generalization - by it I mean a unique set of invariants and covariants of $\text{SU}(2)$ assigned to each density matrix. The uniqueness (up to permutations) is crucial, otherwise one can just spectrally decompose the density matrix 
$$\rho = \sum_i p_i |\psi_i\rangle\langle\psi_i|$$
and apply the standard Majorana representation to each of its eigenvectors (obtaining invariants $\{ p_i \}_{i\in\{1,\ldots,n\}}  $ and covariants $\{ |\eta_{i,j}\rangle \}_{i,j\in\{1,\ldots,n\},}$ ). However, for the case of the eigenvalue degeneracy it is no longer unique.
A naive dimension counting $(n+1)^2-1=2n + n^2$ gives a hint that there may be $n$ covariant points plus an invariant $n\times n$ matrix. 
An auxiliary question is if there is a unique representation of $k$-dimensional subspaces of pure symmetric $n$-qubit states? (If there is one, then it is possible to 'fix' the mentioned approach with the spectral decomposition.) 
 A: A generalization of the Majorana representation to $N+1$ dimensional density matrices can be performed as follows:
(By a generalization I mean a representation of the density matrix by means of a certain number of points on the Bloch sphere (not necessarily independent) + a polarization vector belonging to the $N+1$ dimensional probability simplex)
I'll consider first the generic case (multiplicity free eigenvalues) described in the question for completeness. 
Every $N+1$ dimensional matrix can be written as:
$\rho = \sum_{i=0}^{N}p_i\theta_i$
where $p$ is the eigenvalue vector $p\in \Delta^n$ (the probability simplex in $\mathbb{R}^{n+1}$)
and $\theta_i$ are the one dimensional projectors on the the eigenvectors:
$\theta_i^2=\theta_i$ 
satisfying the orthonormality constraints.
$tr(\theta_i\theta_j)=\delta_{ij}$
The Majorana representation 
$\mathrm{Sym}^N(\mathbb{C}P^1)  \approx \mathbb{C}P^N$
allows to express the first projector in terms of $N$ points on the Bloch sphere, and the second one in terms of $N-1$ points because the orthonormality constraints and the third in terms of $N-2$ points etc.
Dimension count
$2 \times (0+1+ . . .+N) + N =(N+1)^2-1$
The case of an eigenvalue of multiplicity $1<M<N$
In this case the projector on the degeneracy eigenspace is of dimension greater than 1. The orbit of these projectors is the Grassmannian $Gr(M,N)$.
In order to perform a Majorana representation of this projector, we first embed the Grassmannian into a complex projective space by means of  a Plucker embedding:
$Gr(M,N) \rightarrow \mathbb{C}P^{{N \choose M}-1}$
and then perform the Majorana map. 
Now, both the Majorana map and the Plucker embedding are known explicitely, which makes the the whole construction possible. In the degenerate cases, the dimension count is less than that of the space of all density matrices.
Update:
This update is intended to provide a partial answer for Piotr's comment
Remark: I should have remarked that material on the classification of the density matrix orbits according to their eigenvalue degeneracy in this answer is based on
Geometry of quantum states by Bengtsson and Życzkowski (Mainly on chapter 8)
The case of the Grassmannian $Gr(M,N)$ (which is the unitary orbit of $N$ dimensional density matrices with two distinct eigenvalues, one of the of multiplicity $M

The rigid $SU(2)$ which acts on the qubits of the Majorana representation of $\mathbb{C}P^{{N \choose M}-1}$ is a subgroup of the isometry group $SU({N \choose M})$ of the complex projective space. One should not expect apriori that it is a subgroup of the isometry group $SU(N)$ of the Grassmannian.
However, I checked the the simplest case $Gr(2,4) \rightarrow \mathbb{C}P^5$. This Grassmannian is given by the Plucker relation given by the quadric $w \wedge w = 0$ in the homogeneous coordinates $w = \sum_{i=j=1, i< j}^4w_{ij} e_i\wedge e_j$ of  $\mathbb{C}P^5$, ($e_i$ form a basis of $\mathbb{C}^4$).
After the identification of the homogeneous coordinates $w$ with the homogeneous coordinates of the Majorana representation such that on both the action of the $SU(2)$ generator $\sigma_3$ is diagonal, it turned out that the Plucker relation is invariant under the whole $SU(2)$ action. In other words there exists an $SU(4)$ transformation $U$ such that this action can be implemented as: $U^{-1} \theta U$, where $\theta$ is the projector onto the two dimensional degeneracy space of the density matrix. This result is new (and very interesting) to me and I don't have a deep understanding of its reason. I hope it generalizes to all cases.
The generic case of distinct eigenvalues:
In this case, the orbit of the quantum states is the flag manifold $Fl(N)$. In this case there are $N-1$ distinct  Bloch spheres corresponding to the hierarchy of eigenvectors as given in the main answer.
A rigid $SU(2)$ action on the qubits of the Majorana representation of each of the $n$-th eigenvector is realized as a spin $\frac{N-n+1}{2}$ representation on the subspace orthogonal to the higher in hierarchy eigenvalues. 
In addition there will be a corresponding action on all lower in hierarchy eigenvectors. I think that in order to feel comfortable with the construction, one can work out explicitely the three dimensional case with the flag manifold $Fl(3)$  expressed as a $\mathbb{C}P^1$ bundle over $\mathbb{C}P^2$
