Two qubits system in polar co-ordinates I know that I can write a single qubit state in terms of polar co-ordinates $(r,\theta,\phi)$ on a Bloch sphere.
\begin{equation} 
\rho = 
\begin{pmatrix}
\frac{1+r \cos\theta}{2} &\frac{r \exp(-i\phi)\sin\theta}{2} \\
\frac{r \exp(i\phi)\sin\theta}{2} &\frac{1-r \cos\theta}{2} \nonumber
\end{pmatrix}
\end{equation}
Are this kind of polar decomposition exist for two qubits system?
 A: A valid density operator is any Hermitian, trace 1, matrix (with complex entries) and all eigenvalues between 0 and 1. Any two qubit system may be represented therefore by a Hermitian, trace 1 4x4 matrix.
Your qubit representation could be rewritten, more suggestively as:
\begin{align}
   \rho &= \frac{1}{2}\left(\operatorname{I} + a_1 \sigma_x + a_2 \sigma_y + a_3\sigma_z \right)\\ 
&= \frac{1}{2}\left(\operatorname{I} + \vec{a} \cdot\vec{\sigma}\right)
\end{align}
Where $\operatorname{I}$ is the 2x2 identity matrix and the $\sigma_i$ are the Pauli matrices which are a basis for the space of Hermitian, trace 0 2x2 complex matrices.
Similarly to qubits we can decompose any n-dimensional quantum system using the identity and any basis for the space of Hermitian, traceless nxn complex matrices, a good choice of basis is the Gell-Mann matrices, and a construction in dimension n is given here.
Since they are Hilbert-Schmidt orthogonal you can find the coefficients $a_i$ for your state $\rho$ by taking the inner product:
$$a_i = \operatorname{tr}\sigma_i \rho$$
If you like you can then rewrite the vector $\vec{a}$ in any coordinate system you like (including polars).
There is a problem however because although every valid density operator may be rewritten in this way (with $\vec{a}$ being a n-dimensional sub-normalised vector) it is not the case that every operator of that form will be a valid density operator, the problem is that some of your eigenvalues will become negative and others will become greater than 1 to balance this out. This doesn't happen for qubits (any operator of the form $\rho$ I wrote above will be a valid density operator) and it makes working in higher dimensions fairly annoying sometimes.
A: Yes.  We can do this for any number of qubits  using N dimensional spherical coordinates.  For two qubits  we can write a general density matrix as a linear combination of direct products of Pauli matrices and identity , $$  \rho = \sum_{ij=0}^3 a_{ij}~ \sigma_{i} \otimes \sigma_{j}$$.
Here $\sigma_{0} = I$, and the rest are the usual Pauli matrices. For a valid density matrix $a_{00} = 1/4$. It can be shown that all valid density matrices lie inside a higher dimensional sphere centered at the coordinate corresponding to $a_{00}$. And these points can be represented using higher dimensional spherical coordinates. But unlike the Bloch sphere for a single qubit, all points do not represent valid density matrices as there can be matrices with negative eigenvalues. In all higher dimensions it is a convex body inside this sphere. A detailed explanation can be found here and here.
A: A one qubit state may be written in general as 
\begin{equation}
|\psi_1\rangle=\alpha|0\rangle+\beta|1\rangle
\end{equation}
where $\alpha,\beta \in \mathbb C$ and there is the further restriction that $\langle \psi_1|\psi_1\rangle=|\alpha|^2+|\beta|^2=1$. However, only the relative phase between $|0\rangle$ and $|1\rangle$ is physically meaningful (quantum states are rays in Hilbert space), so we may factor out a phase $e^{i\theta}$.
Let's choose that phase so that $a=e^{-i\vartheta}\alpha$ is real while $b=e^{-i\vartheta} \beta$ can remain complex. Our state is then
\begin{equation}
|\psi_1\rangle=e^{i\vartheta}(a|0\rangle+b|1\rangle).
\end{equation}
The previous condition is now $\langle \psi_1|\psi_1\rangle=|a|^2+|b|^2=1$. We have three real parameters and one constraint and so we may say that our state $|\psi_1\rangle$ lives on the unit 2-sphere.
Analagously, a two qubit state may be written in general as
\begin{equation}
|\psi_2\rangle=a|00\rangle+b|01\rangle+c|10\rangle+d|11\rangle
\end{equation}
where we can choose $a\in \mathbb R$ and $b,c,d \in \mathbb C$ with the further restriction that $|a|^2 + |b|^2 + |c|^2 + |d|^2 =1$. In this case there are seven free parameters and one constraint so a two-qubit system would live on the unit 6-sphere and have a corresponding higher dimensional polar decomposition.
