Eigenvalues of tensor products of Pauli matrices $\sigma_\mu\otimes\sigma_\nu (\mu,\nu\in\{x,y,z\})$ always come in pairs, as a consequence of the fact that eigenvalues of Pauli matrices are $\pm 1$. As a reminder, Pauli matrices are defined to be $$\sigma_x= \begin{pmatrix} 0 & 1 \\ 1 & 0 \\ \end{pmatrix}, \sigma_y= \begin{pmatrix} 0& -i \\ i& 0 \\ \end{pmatrix}, \sigma_z= \begin{pmatrix} 1 & 0\\ 0 & -1\\ \end{pmatrix} $$
Adding in the identity matrix $\sigma_0=\begin{pmatrix} 1 & 0\\ 0 & 1\\ \end{pmatrix}$ and change letter subscipts $\{x,y,z\}$ into numbers $\{1,2,3\}$ for notational convenience, we use the tensor products of the four matrices to generate $4\times 4$ matrices of the form $\sigma_i\otimes\sigma_j$. There are $16$ of them in total.
Now arbitrary linear combinations of them form a $16$-dimensional linear space, and any vector in it must have real eigenvalues since it must be a Hermitian matrix. However, it does not necessarily need to have paired eigenvalues $\pm\lambda$.
My questions is, what are the subspaces that produce paired eigenvalues? I feel it has something to do with group theory since Pauli matrices are generator of $SU(2)$, but I am not very good at representation theory and not sure if there is any theoretical approach/ method to handle problems like this?
I have some computational results. For instance, if my algorithm is not wrong, the dimension of the largest subspace is $10$. If we use a $4\times4$ indicator matrix to indicate which $\sigma_i\otimes\sigma_j$ is in the subspace. They are the following, a total of $6$ of them. Here I show one of them
$$\begin{pmatrix} 0 & 0 & 0 & 1 \\ 1 & 1 & 1 & 0 \\ 1 & 1 & 1 & 0 \\ 1 & 1 & 1 & 0 \end{pmatrix}$$ which means that any linear combinations of the form $$a_{03}\sigma_0\otimes\sigma_3+\sum_{i\in\{1,2,3\},j\in\{0,1,2\}} a_{ij}\sigma_{i}\otimes\sigma_j $$ have paired eigenvalues.
And it follows immediately that any linear subspaces of those $6$ subspaces again have the property of producing paired eigenvalues. A quick sanity check that the top left element cannot be $1$ in any of these indicator matrices, since $\sigma_0\otimes\sigma_0=\mathbb{I}_4$ does not produce paired eigenvalues $\pm\lambda$.
And that is not the end of the story. It seems that there are other small subspaces that are not subspaces of these $6$ largest ones. For subspaces with dimension $9$, there are $75$ in total, and it means that apart from the $10\times 6=60$ ones that we get immediately from the largest ones, there are $15$ that are new. (It reminds me of the number of generators of $SO(6)$, though they may be totally unrelated.)
We can continue doing this: there are $60$ maximal "paired-eigenvalues" subspaces of dimension $8$. By "maximal" I mean they are not subspaces of any other high-dimensional "paired-eigenvalues" subspaces (i.e., they cannot be extended). For dimension $7$, there is no maximal one. For $6$, there are $17$...
One possible problem in the results I get above is that I am only considering linear subspaces of the form $a_{ij}\sigma_i\otimes\sigma_j$ (summed over $i$ and $j$ with arbitrary $a_{ij}$). There are $2^{16}$ of them and correspond 1-to-1 to indicator matrices. There might be other subspaces that have this "paired eigenvalues" property with fine-tuned coefficients $a_{ij}$ that cannot be represented by an indicator matrix.
For instance, the vector $(\sigma_x\otimes\sigma_x+\sigma_y\otimes\sigma_y)$ also forms a $1$-dimensional subspace, but it cannot be represented by our indicator matrix.