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.