Why can arbitrary two qubit density matrices be expressed in this form? In the paper "Violating Bell inequality by mixed spin 1/2 states: necessary and sufficient condition" (http://www.sciencedirect.com/science/article/pii/037596019500214N#) by three Horodecki siblings, it is stated that an arbitrary two qubit density matrix can be written in the following form:
$\rho = \frac{1}{4}(I \otimes I + r \cdot \sigma \otimes I + I \otimes s \cdot \sigma + \Sigma_{i,j=1}^3t_{ij}\sigma_i \otimes \sigma_j)$
where $r$ and $s$ are some vectors.
I understand how $I$ and $\sigma_i$ for $i=1,2,3$ form a basis for operators on a single qubit so taking pairwise cross products of each forms a basis for operators on two qubits. But why does the coefficient of the $I \otimes I$ have to be $\frac{1}{4}$?
As a side question, the paper also says $t_{ij}=tr[\rho\sigma_i \otimes \sigma_j]$. Why is that the case?
Edit: The other page linked doesn't answer either of my questions. They say that $I,\sigma_i$ form a basis, but don't explain why $\frac{1}{4}$ is the coefficient. It also doesn't mention my second question.
 A: The co-efficient $\frac{1}{4}$ simply ensures that correct normalization so that the probabilities of the different pure states composing the classical mixture (the "mixed state") in question sum to unity: recall that:
$$\rho = \sum\limits_j p_j |\psi_j\rangle\langle \psi_j|$$
and when the pure states are normalized properly, we have ${\rm tr}(|\psi_j\rangle\langle \psi_j|)=1$ so that ${\rm tr}(\rho)=\sum_j p_j=1$.
So you can work out ${\rm tr}(I\otimes I)=4$, and all the other terms in your sum are traceless. Use the identity ${\rm tr}(X\otimes Y) = {\rm tr}(X)\, {\rm tr}(Y)$, and then trace's linearity, to show that, for example, ${\rm tr} (r \cdot \sigma \otimes I )= \sum\limits_{j=1}^3 {\rm tr}(r_j \sigma_j\otimes I)=0$, as the nonextended Pauli matrices ($j=1,\,2\,3$) are all traceless. The Wikipedia Trace Page summarizes many other properties of the trace and it has very handy list for working through stuff like this (most of the properties also hold for infinite dimension spaces, but this doesn't worry us here). I'm also using the convention that $\sigma_0 = I$.
As for your last question, the tensor products of extended Pauli matrices (including the identity matrix) (1) form a basis for the Hermitian $4\times 4$ matrices (as you know) and, crucially, (2) the pairs of the form $\sigma_j\otimes \sigma_k$ ($k=0,\,1,\,2,\,3$) are orthogonal with respect to the inner product defined on Hermitian matrices by ${\rm tr}(X\,Y)={\rm tr}(X^\dagger\,Y)$, so they form an orthonormal basis. So the co-efficients are simply falling out of the expression because you're simply finding components in a Euclidean space with respect to any orthonormal basis.
