Background
Somebody asked a question on the Chinese Q&A website ZhiHu(知乎), which roughly translates to the following: given an $n\times n$ matrix $M$ consisting of $n^2$ distinct real numbers, consider the two ways of shuffling its elements:
1). First, apply an intra-row shuffle $\mathcal R$ (which means shuffling the elements within each row for all rows), and then apply an intra-column shuffle $\mathcal C$ (defined similarly) . We get a random matrix $\mathcal C(\mathcal R(M))$.
2). Directly apply an elementwise full shuffle $\mathcal E$ (which means embedding $M$ as a vector in $\Bbb R^{n^2}$ and then shuffling its entries).
The original question is: Does $\mathcal C(\mathcal R(M))$ equal $\mathcal E(M)$ in distribution?
The answer is no. And there are more than one ways to prove it, the quickest perhaps being direct counting: $\mathcal C(\mathcal R(M))$ yields at most $(n!)^{2n}$ outcomes, but $\mathcal E(M)$ yields exactly $(n^2)!$ outcomes. A not-so-difficult analysis using Stirling approximation shows $(n^2)! > (n!)^{2n}$ for all $n\ge 2$. (See this anwer (in Chinese) for example).
The Question
What if we continue to apply additional layers of $\mathcal C$'s and $\mathcal R$'s on $M$? Would the resulting matrix be sufficiently "stochastic" and hence possibly equal $\mathcal E(M)$ in distribution? For example, Does $\mathcal R(\mathcal C(\mathcal R(M)))$ (or $\mathcal C(\mathcal R(\mathcal C(\mathcal R(M))))$ and so on, suppose we can apply more layers if necessary) equal $\mathcal E(M)$ in distribution? Or, at least, can these additional layers of operation span the entirety of outcomes of $\mathcal E(M)$?
My Thought
Counting might still help but is perceivably much more difficult. Because $(n!)^{kn} \gg (n^2)! $ for $k \ge 3$ already, we cannot simply chain-multiply the number of outcomes of each operation. Surely we will have to eliminate a lot of duplicated outcomes to achieve a reasonable estimate of the number of outcomes after multiple $\mathcal C,\mathcal R$. But removing multiplicities itself doesn't seem that easy.
Another possible way to disprove this guess is to use some other algebraic or analytical features of a square matrix. The observation is that, if $X=Y$ in distribution as random vectors, then $f(X)=f(Y)$ must also equal in distribution where $f$ can be a more or less arbitrary function, such as $\det$ or $\text{rank}$ or $\sum_i\lambda_i$ etc. The most difficult part is clearly finding a suitable $f$ which is somehow tractable under operations $\mathcal C,\mathcal R$ (unfortunately $\det$ doesn't fall under this category...), then spot a discrepancy between the distribution $f(\mathcal R(\mathcal C(\mathcal R(M))))$ and $f(M)$.
