Let $\gamma_i\ (i=1,2,\ldots N)$ be the Dirac gamma matrices satisfying the Clifford algebra $$\gamma_i\gamma_j+\gamma_j\gamma_i=2\delta_{ij} I\ \ (i,j=1,2,\ldots,N).$$ Then the tensor products $\gamma_i\otimes\gamma_i$ commute with each other: $$[\gamma_i\otimes\gamma_i, \gamma_j\otimes\gamma_j]=0.$$ This means all of $\gamma_i\otimes\gamma_i$ can be simultaneously diagonalized by some similarity transformation $S$. I wonder what the characterizations of such diagonalization are, such as the arrangement of eigenvalues of $\gamma_i\otimes\gamma_i$ after diagonalization, the (fast) computational algorithm for $S$ etc. I find surprisingly little about this elementary construction online, any reference will also be appreciated.

Update: I have basically solved problem, I will write an answer to it in the future if time permits.