How to linearly expand density matrices with "physical" basis elements? If we have a two qubit system, can we find a convenient set $\{\hat e_i\}$ for the $4 \times 4$ density matrix $\hat \rho$, such that $\rho$ can be writen as a linear combination of the elements in the set, and each element of the set a physical state, i.e. the $\mathrm{Tr} \{ \hat e_i\} = 1$, $\hat e^{\dagger} = \hat e$ and $\hat e_i \ge 0$?
For example, two qubits can be expanded as:
$$\hat \rho = \frac{1}{4} \sum_{i,j = 0}^3 S_{ij} \sigma_i \otimes \sigma_j $$
Where $\sigma_i$ belongs to the vector $\vec \sigma = \{1_2, \sigma_1, \sigma_s, \sigma_3 \}$, i.e. the $2 \times 2$ identity matrix followed by the Pauli matrices.  However, not all of the 16 terms in this expansion are physical states. Most of them have trace zero.
We can also expand $\hat \rho$ as:
$$ \hat \rho = \sum_{i, j, k, l = 0}^{1} A_{ijkl} | i \rangle \langle j | \otimes | k \rangle \langle l | $$
But again, the basis states in this expansion are not physical, most of them have trace 0.
I'm looking for an expansion $\hat \rho = \sum_{i=0}^{15} \hat e_i$ where all of $\hat e_i$'s are physical. Is it possible, or is there any resource out there which talks about this?
EDIT
A single qubit can be written as a linear combination of $\frac{1_2 + \sigma_i}{2}$ (of course, there'll be conditions on the coefficients) where $i$ runs from 1 to 3. All three of these are physical.I'm looking for something similar for two qubits.
 A: There is several options, as long as you don't require positive weights:

*

*Choose a SIC-POVM.


*Take your favorite hermitian basis -- for instance, the tensor products of the Pauli matrices (including the identity), and add the identity of them until the resulting operators are positive semidefinite. (In case you chose the Pauli products $P_i$, this means your basis has e.g. entries $I+P_i$, together with the identity matrix.) If you want that they have trace 1: Divide by the trace.
A: An interesting choice is to use so-called measurement frames. In brief, the idea is that given any informationally-complete POVM with elements $\Pi_i$, you can find another set of operators, call them $F_i$, such that any state can be decomposed as $\rho=\sum_i \operatorname{Tr}(\Pi_i \rho) F_i$. This makes the coefficients "physical" because they just become the output probabilities. This is useful in tomography applications and many other things. A good reference is (A. J. Scott 2006).
This probably also works for not informationally-complete POVMs, but some modifications in the formalism might be required.
