In quantum computation and quantum information it is very common to use e.g. the effect of a Hadamard matrix $H$ over $2n$ spins. Using (I think it is called the Kroenecker product in mathematical literature) the tensor product $\otimes$, one can write the Hadamard matrix for example for two spins as
$$ H\otimes H. $$
I want to ask precisely the opposite question. Is there a theorem(s) that ensures that, e.g. $M\in U(4)$ could be written as
$$ M = M_1\otimes M_2, $$
where (perhaps) $M_i\in U(2)$?
There is a rich culture of 4-way arrays, that is square matrices whose elements are square matrices, but I am only a physicist, so I can mention what we mugs in the trenches do faced with these problems, routinely.
Assume, first, as you did for your 4×4 matrix, that it is of the form $M\otimes m$ for M and m diagonalizable 2×2 matrices--no preteses of generality here. That is, it diagonalizes to the 4×4 diagonal matrix $$\begin{pmatrix} A&0\\ 0&B \end{pmatrix} \otimes \begin{pmatrix} a&0\\ 0&b \end{pmatrix}= \begin{pmatrix} Aa&0&0&0\\ 0&Ab&0&0 \\ 0&0&Ba&0 \\ 0&0&0&Bb \end{pmatrix}. $$ Its eigenvalues are not independent, but satisfy the condition $\lambda_1 \lambda_4=\lambda_2 \lambda_3 $. A similarity transformation of this 4×4 matrix, however, may permute these eigenvalues, and hence sour this relation.
So, given a 4×4 matrix you can diagonalize, it suffices if the above condition on eignevalues is satisfied, but is not really necessary: if you found such a pairing among any two pairs, you might go backwards, rearrange them in such a desirable ordering as above, incorporate the permuting matrices into the original diagonalizing transformation, and be done. You might pick a simple example of Pauli matrices. In physics, Dirac's $\gamma^\mu$ matrices are routinely represented and handled this way.
I'd be shocked if multilinear many-array people did not have elegant schemes for this sort of thing, but, here, this is what the grunts do.