Can permutation similarities change order of Kronecker products? Consider the Kronecker products
$${\bf R_1 \otimes R_2} \text{ and } {\bf R_2 \otimes R_1}$$
We can show in the special case they both reside in $\mathbb R^{2\times 2}$ that we can find permutation matrix $\bf P$ so that:
$$(\bf R_1 \otimes R_2) - P(\bf R_2 \otimes R_1)P = 0$$
In fact this $\bf P$ is rather simple:
$${\bf P} = \begin{bmatrix}1&&&\\&&1&\\&1&&\\&&&1\end{bmatrix}$$
and also happens to be quite restricted, as ${\bf P}^{-1} = {\bf P}^{T} = \bf P$
So this could potentially in fact be a similarity transformation $\bf A = P^{-1}BP$
Can we prove that for any sizes 
$${\bf R_1} \in \mathbb R^{N\times N}, {\bf R_2} \in \mathbb R^{M\times M}$$
we can always find such a $\bf P$?
I also found one by trial-and-error for $(N,M) = (2,3)$ and transpose.

(As an engineer, of course I am more interested of construction than existence, but I would be satisfied with existence proof).
 A: The matrices you are looking for are called commutation matrices, defined as the unique permutation matrix that satisfies $$K_{m,n} \cdot {\rm vec}\{A\} = {\rm vec}\{A^{\rm T}\}$$ for an arbitrary $m \times n$ matrix $A$. 
It can then be shown that $K_{m,n} \cdot (A\otimes B) \cdot K_{p,q} = B \otimes A$ where $A$ is $m \times p$ and $B$ is $n \times q$. In your case, since $R_1$ and $R_2$ are square, you could have $P=K_{N,M}$, as you already observed in the special case.
The commutation matrices are introduced in [MN79] and further studied in the book [MN88], which is a pretty nice read.
[MN79] Magnus, Jan R.; Neudecker, H., The commutation matrix: Some properties and applications, Ann. Stat. 7, 381-394 (1979). ZBL0414.62040.
[MN88]  Magnus, Jan R.; Neudecker, Heinz, Matrix differential calculus with applications in statistics and econometrics, Wiley Series in Probability and Mathematical Statistics. Applied Probability and Statistics. Chichester etc.: John Wiley & Sons. XVII, 393 p.; \textsterling 24.50 (1988). ZBL0651.15001.
A: The answer I have found to be the most likely (by brute force calculations on random matrices of a large span of sizes $N,M$) is to choose $P$ so that it performs transpose on the vectorization on the vector space spanned by the vectorized matrices of size $N \times M$. I still have no proof why it works, so any proof would be welcome.
