Let me just observe that while rigorous notation is a must for any
serious mathematician it can sometimes block the path of the beginning
reader. What is happening here is very simple. We need the cycle index
of the permutations from $S_2\times S_3$ acting on the slots of the
matrix. While the present case can be computed by inspection I will
explain the general method.
We require the cycle index $Z(M_{2,3})$ of the pairs $(\sigma,\tau)$
of row and column permutations acting simultaneously on the slots of
the matrix ($6$ of them, and with $12$ permutations total). Start
from the two cycle indices for the row and column permutations
$$Z(S_2) = \frac{1}{2} a_1^2 + \frac{1}{2} a_2$$
and
$$Z(S_3) =\frac{1}{6} b_1^3 + \frac{1}{2} b_1 b_2 + \frac{1}{3} b_3.$$
The method here is as follows. We draw a diagram of the cycle types of
$\sigma$ and $\tau,$ perhaps one beneath the other. Now for the cycle
index of the cartesian product of $S_2$ and $S_3$ we must factor the
combined action of the two permutations on the slots into cycles. We
represent row-column pairs i.e. slots by marking say the row and the
column and connecting them with a edge, not to be confused with the
directed edges of the two cycles. This edge travels in parallel along
the two cycles it is on and returns to its initial position after
$\mathrm{lcm}(l_1, l_2)$ steps, where $l_{1,2}$ are the lengths of the
cycles. As the pair of cycles contributes $l_1\times l_2$ slot
identifiers we get for the contribition to the combined cycle index
the term
$$a_{\mathrm{lcm}(l_1, l_2)}^{l_1 l_2/\mathrm{lcm}(l_1, l_2)}.$$
We now do the computation. There are thus six possible combinations of
cycle types that combine to form $Z(M_{2,3})$. We process these in
turn, leaving the most difficult part for last.
First, combining $a_1^2$ and $b_1^3.$ This fixes all six slots for a
contribution of
$$\frac{1}{12} a_1^6.$$
Second, combining $a_1^2$ and $b_3.$ This partitions the pairs into
three-cycles for a contribution of
$$\frac{1}{6} a_3^2.$$
Third, combining $a_2$ and $b_1^3.$ Here we have everything on
two-cycles to get
$$\frac{1}{12} a_2^3.$$
Fourth, combining $a_2$ and $b_3$. This will produce a six-cycle to get
$$\frac{1}{6} a_6.$$
Now for the tricky part. Fifth, combining $a_1^2$ and $b_1 b_2.$ There
are two pairs that are fixed by these, and two pairs on two-cycles
and we obtain
$$\frac{1}{4} a_1^2 a_2^2.$$
Finally, sixth, combining $a_2$ and $b_1 b_2.$ We have everything on
two-cycles and obtain
$$\frac{1}{4} a_2^3.$$
Adding everything we now have the cycle index
$$Z(M_{2,3})
= \frac{1}{12} a_1^6
+ \frac{1}{3} a_2^3
+ \frac{1}{6} a_3^2
+ \frac{1}{4} a_1^2 a_2^2
+ \frac{1}{6} a_6.$$
In order to apply Burnside and ask about colorings with at most $N$
colors we have that the assignment of the colors must be constant on
each cycle and we obtain
$$\frac{1}{12} N^6
+ \frac{1}{3} N^3
+ \frac{1}{6} N^2
+ \frac{1}{4} N^4
+ \frac{1}{6} N.$$
This yields the sequence
$$M_n = 1, 13, 92, 430, 1505, 4291, 10528, 23052, \ldots$$
which is OEIS A027670. In particular we
find there are $92$ colorings using at most three distinct colors. We
could apply PET at this point since we have the cycle index. E.g. for
three colors we obtain
$$1/12\, \left( R+G+B \right) ^{6}
\\ +1/4\, \left( R+G+B \right) ^{2} \left( {B}^{2}+{G}^{2}+{R
}^{2} \right) ^{2}+1/6\, \left( {B}^{3}+{G}^{3}+{R}^{3} \right) ^{2}
\\ +1/3\, \left( {B}^{2}+
{G}^{2}+{R}^{2} \right) ^{3}+1/6\,{B}^{6}
+1/6\,{G}^{6}+1/6\,{R}^{6}$$
which expands to
$${B}^{6}+{B}^{5}G+{B}^{5}R+3\,{B}^{4}{G}^{2}+3\,{B}^{4}GR
\\ +3\,{B}^{4}{R}^{2}+3\,{B}^{3}{G}^{3}+6\,{B}^{3}{G}^{2}R
\\ +6\,{B}^{3}G{R}^{2}+3\,{B}^{3}{R}^{3}+3\,{B}^{2}{G}^{4}
\\ +6\,{B}^{2}{G}^{3}R+11\,{B}^{2}{G}^{2}{R}^{2}+6\,{B}^{2}G{R}^{3}
\\ +3\,{B}^{2}{R}^{4}+B{G}^{5}+3\,B{G}^{4}R+6\,B{G}^{3}{R}^{2}
\\ +6\,B{G}^{2}{R}^{3}+3\,BG{R}^{4}+B{R}^{5}+{G}^{6}+{G}^{5}R
\\ +3\,{G}^{4}{R}^{2}+3\,{G}^{3}{R}^{3}+3\,{G}^{2}{R}^{4}
\\+G{R}^{5}+{R}^{6}.$$
The reader might want to verify some of these with pen and paper.
We may also apply inclusion-exclusion to obtain the count of colorings
of the matrix with exactly $N$ colors. The nodes of the poset
correspond to sets $P\subseteq [N]$ of colors which include colorings
that use some subset of these colors, with the top node representing
at most $N$ colors, which is the only $P$ that includes colorings with
exactly $N$ colors. Colorings with exactly $p \lt N$ colors where
$p\ge 1$ are included at all nodes that are supersets of the $p$
colors. We thus obtain for the total weight
$$\sum_{q=p}^N {N-p\choose q-p} (-1)^{N-q}
= (-1)^{N-p} \sum_{q=0}^{N-p} {N-p\choose q} (-1)^q = 0$$
since $N-p\ge 1.$ Colorings with less than $N$ colors have total
weight of zero in the poset. We thus obtain
$$\sum_{q=1}^N {N\choose q} (-1)^{N-q} M_q.$$
We get a finite sequence since with six slots it is not possible to
have a coloring with more than six distinct colors:
$$1, 11, 56, 136, 150, 60$$
The last term represents colorings with exactly six colors. This means
all slots in the matrix are distinct. Therefore all orbits have the
same size, the number of permutations in the group, which is twelve,
and indeed $6!/12 = 60.$ The term for two colors indicates that the
two monochrome colorings have been excluded.
This MSE link
has the Maple code for the general case.
Addendum. Here is what we mean when we say in the introduction
that the cycle index can be computed by inspection. This refers to the
isomorphism between $M_{2,3} = S_2\times S_3$ and $D_6$, the dihedral
group (reflections and rotations of regular polygons) acting on six
slots in this case. The cycle indices $Z(D_p)$ are tabulated and have
simple closed forms, consult
e.g. Wikipedia. Label the
vertices of a hexagon in clockwise order with the labels $(0,0),
(1,1), (0,2), (1,0), (0,1)$ and $(1,2).$ Then it is not difficult to
see that the rotations of the hexagon are in a bijection with the
pairs of cycles from $C_2 \times C_3$ embedded in $M_{2,3}.$ E.g. the
rotation that takes the top vertex to its clockwise neighbor
corresponds to the two cycles $(0,1)$ and $(0,1,2)$. The reflections
in an axis passing through opposite vertices preserve parity
(permutation $(0)(1)$ from $S_2$) and fix one element from $0,1,2$
while permuting the other two in two-cycles. The reflections in an
axis passing through opposite edges flip parity (permutation $(0,1)$
from $S_2$) and fix one of three elements from $0,1,2$ while
exchanging the other two. In this way we have bjectively accounted for
all permutations and the proof of the isomorphism is complete.
Best Answer
TL;DR (or for those unfamiliar with group theory) A brief summary of the results, without proofs.
Any pair of such $n\times n$-matrices with the same unordered set of rows, columns, diagonals and antidiagonals differ by a permutation of their entries, and hence by a permutation of the index set $$\{(i,j):\ 1\leq i,j\leq n\}.$$ It turns out that any permutation that preserves the set of rows, columns, diagonals and antidiagonals is in fact an affine map, i.e. it is some combination of a rotation, reflection, stretching and translation${}^{\ast}$, but in general not all affine maps preserve the set of rows, columns, diagonals and antidiagonals.
You had already found that the translations, reflections and rotations by a quarter turn preserve the set of rows, columns, diagonals and antidiagonals; that's almost everything! What remains are sometimes more rotations, and a lot of stretchings. The stretching of a matrix $M=(m_{ij})_{1\leq i,j\leq n}$ by a factor $1\leq c\leq n$ is the matrix $$N=(n_{i,j})_{1\leq i,j\leq n}\qquad\text{ with }\qquad n_{i,j}:=m_{ci,cj},$$ where we take $ci$ and $cj$ modulo $n$. For example, for $n=5$ and $c=2$ the matrix $$\begin{pmatrix} 1&2&3&4&5\\ 6&7&8&9&10\\ 11&12&13&14&15\\ 16&17&18&19&20\\ 21&22&23&24&25 \end{pmatrix} \qquad\text{ becomes }\qquad \begin{pmatrix} 13&11&14&12&15\\ 3&1&4&2&5\\ 18&16&19&17&20\\ 8&6&9&7&10\\ 23&21&24&22&25 \end{pmatrix}. $$ Such a stretching is well-defined only if $c$ is coprime to $n$, i.e. if $\gcd(c,n)=1$. It turns out that all such stretching with $\gcd(c,n)$ preserve the set of rows, columns, diagonals and antidiagonals.
The question then remains; which rotations preserve the set of rows, columns, diagonals and antidiagonals? The results are fundamentally different depending on whether $n$ is even or odd.
While it is clear what a quarter turn of a matrix is, I should clarify a one-eighth turn. The one-eighth turn of a matrix $M=(m_{ij})_{1\leq i,j\leq n}$ by a factor $1\leq c\leq n$ is the matrix $$N=(n_{i,j})_{1\leq i,j\leq n}\qquad\text{ with }\qquad n_{i,j}:=m_{i+j,-i+j},$$ where again we take $i+j$ and $-i+j$ modulo $n$. For example, the one-eighth turn of the matrix $$\begin{pmatrix} 1&2&3&4&5\\ 6&7&8&9&10\\ 11&12&13&14&15\\ 16&17&18&19&20\\ 21&22&23&24&25 \end{pmatrix} \qquad\text{ is }\qquad \begin{pmatrix} 21&9&17&5&13\\ 14&22&10&18&1\\ 2&15&23&6&19\\ 20&3&11&24&7\\ 8&16&4&12&25 \end{pmatrix}. $$ Combining all these symmetries quickly becomes confusing, some group theory is essential to get a clear picture. But it turns out the symmetries go together quite well; only rotations and stretching get mixed up, because a half-turn is the same as stretching by a factor $-1$.
To generate all matrices with the same set of rows, columns, diagonals and antidiagonals as a given matrix $M$, simply apply all the symmetries as follows:
$\ast$ There are also the shear maps, but these preserve the set of rows, columns, diagonals and antidiagonals only if these are all lines in the $n\times n$-plane, which is the case if and only $n\leq3$.
Here follows the original proof for odd $n$.
Observation: Any two rows, columns, diagonals or antidiagonals that are not in the same category (e.g. that are not both rows) meet in precisely one point if and only if $n$ is odd.
I will use this observation throughout the proof without mention, as the proof concerns only odd $n$. Because this fails for even $n$ the geometry is significantly less intuitive, so I will not consider this case for now.
Observation: The set of $n\times n$-matrices whose entries are precisely the integers from $1$ to $n^2$ inclusive, is naturally in bijection group of symmetries of the affine $n$-plane $\Bbb{A}_2(n):=(\Bbb{Z}/n\Bbb{Z})^2$.
The question then becomes
Proposition 1. For odd $n>3$ the set of horizontal, vertical, diagonal and antidiagonal lines is preserved by all symmetries of the form $$\Bbb{A}_2(n)\ \longrightarrow\ \Bbb{A}_2(n):\ \binom{i}{j}\ \longmapsto\ uA\binom{i}{j}+B,$$ with $u\in(\Bbb{Z}/n\Bbb{Z})^{\times}$, $B\in(\Bbb{Z}/n\Bbb{Z})^2$ and $A\in\langle\rho,\mu\rangle$ where $$\rho:=\begin{pmatrix} \hphantom{-}1&1\\ -1&1 \end{pmatrix} \qquad\textit{ and }\qquad \mu:=\begin{pmatrix} 1&\hphantom{u}0\\ 0&-1 \end{pmatrix}. $$
Proof. It is clear that translation by any element of $(\Bbb{Z}/n\Bbb{Z})^2$ sends rows to rows, columns to columns, diagonals to diagonals and antidiagonals to antidiagonals. To see that the same is true for multiplication by any $u\in(\Bbb{Z}/n\Bbb{Z})^{\times}$, note that any line in $\Bbb{A}_2(n)$ is of the form $$L_{a,b,c}:=\{ax+by=c:\ x,y\in\Bbb{Z}/n\Bbb{Z}\},$$ and so multiplication by $u\in(\Bbb{Z}/n\Bbb{Z})^{\times}$ maps this line to the line $L_{a,b,uc}$; in particular the slope of the line is preserved, so this scaling also sends rows to rows, columns to columns, diagonals to diagonals and antidiagonals to antidiagonals. It remains to show that $\rho$ and $\mu$ preserve the set of rows, columns, diagonals and antidiagonals.
Geometrically, multiplication by $\mu$ is a reflection in the first row. This sends rows to rows and columns to columns, and swaps diagonals and antidiagonals. Similarly, multiplication by $\rho$ is a rotation around the point $(0,0)$ over an angle of $-\tfrac{\pi}{4}$. This sends rows to diagonals, diagonals to columns, columns to antidiagonals, and antidiagonals to rows. $\hspace{20pt}\square$
Theorem 2. For odd $n>3$, every symmetry of $\Bbb{A}_2(n)$ that preserves the set of rows, columns, diagonals and antidiagonals is of the form described in Proposition 1.
Proof. See below.
Corollary 3. For any $n\times n$-matrix whose entries are precisely the integers from $1$ to $n^2$ inclusive, there are precisely $8\varphi(n)n^2$ such matrices with the same set of rows, columns, diagonals and antidiagonals, where $\varphi(n)$ denotes Eulers totient function.
Lemma 1: Let $n$ be an odd positive integer. Let $\sigma$ be a symmetry of $\Bbb{A}_2(n)$ that preserves the set of rows, columns, diagonals and antidiagonals such that $\sigma((0,0))=(0,0)$ and $\sigma((0,1))$ is on the first row. Then there exists some $d\in(\Bbb{Z}/n\Bbb{Z})^{\times}$ such that $$\forall k\in\Bbb{Z}/n\Bbb{Z}:\qquad\sigma((0,k))=(0,dk).$$
Proof. Let $d\in\Bbb{Z}/n\Bbb{Z}$ be such that $\sigma((0,1))=(0,d)$, and note that $d\neq0$. Because $\sigma$ preserves the set of rows, columns, diagonals and antidiagonals, there are at most six candidates for $\sigma((1,1))$. After all, the point $(1,1)$ is in the same column as $(0,1)$ and on the same diagonal as $(0,0)$, and $\sigma$ maps these to a column, diagonal or antidiagonal through $(0,d)$ and $(0,0)$. The following picture clarifies the situation:
This picture shows the points $(0,0)$ and $(0,d)$ in black, their columns, diagonals and antidiagonals in blue, and the possible positions of $\sigma((1,1))$ also in blue.
The point $(0,2)$ is the unique point on the line through $(0,0)$ and $(0,1)$ that is on a line with $(1,1)$, and that is distinct from $(0,0)$ and $(0,1)$. For each candidate point for $\sigma((1,1))$ this unique point and corresponding line are shown in red in the picture above. Some elementary algebra shows that these points are $(0,-d)$, $(0,\tfrac{d}{2})$ and $(0,2d)$, just as the picture suggests, where $(0,\tfrac{d}{2})=(0,\frac{n+1}{2}d)$ and $(0,-d)=(0,(n-1)d)$.
We will prove that $\sigma((0,2))=(0,2d)$. Note that for $n=3$ we have $2d=-d=\frac{n+1}{2}d$ and so there is nothing to prove. So suppose $n>3$.
Consider the possible positions for $\sigma((1,2))$. Note that $(1,2)$ is the unique point that is on a line with $(0,1)$, $(0,2)$ and $(1,1)$ but not on a line with any pair of them. Also note that $(1,2)$ is not on a line with $(0,0)$ because $n>3$, and it is on the unique line through $(1,1)$ that does not contain $(0,0)$, $(0,1)$ or $(0,2)$.
If $\sigma((1,1))=(d,0)$, which is the bottom left blue point in the picture above, then $\sigma((0,2))=(0,-d)$. Then $\sigma((1,2))$ is on the unique line through $\sigma((1,1))=(d,0)$ that does not contain $(0,0)$, $(0,d)$ and $(0,-d)$, and it is on one of the two lines through $(0,-d)$ that does not contain $(0,0)$ , $(0,d)$ or $(d,0)$. These lines are shown in green in the following picture:
Because $\sigma((1,2))$ is not on a line with $\sigma((0,0))=(0,0)$, we must have $\sigma((1,2))=(-2d,d)$. This point is shown in green in the picture above. But for $n>3$ this point is not on a line with $\sigma((0,1))=(0,d)$, a contradiction. This proves that $\sigma((1,1))\neq(d,0)$, and by symmetry also that $\sigma((1,1))\neq((-d,0))$. This implies that $\sigma((0,2))\neq(0,-d)$, as the first picture shows.
A similar geometrical argument shows that if $\sigma((1,1))=(\pm\tfrac{d}{2},\tfrac{d}{2})$ then the point $(-1,1)$ has nowhere to go, from which it follows that $\sigma((0,2))\neq(0,\tfrac{d}{2})$. I omit the details.
This proves that $\sigma((0,2))=(0,2d)$, and a simple induction argument then shows that $\sigma((0,k))=(0,dk)$ for all $k\in\Bbb{Z}/n\Bbb{Z}$. Because $\sigma$ is injective we have $d\in(\Bbb{Z}/n\Bbb{Z})^{\times}$.$\hspace{20pt}\square$
Corollary: In the situation above we have $\sigma((1,1))=(\pm d,d)$.$\hspace{20pt}\square$
Lemma 2: Let $n$ be an odd positive integer. Let $\sigma$ be a symmetry of $\Bbb{A}_2(n)$ that preserves the set of rows, columns, diagonals and antidiagonals that is the identity on $(0,0)$, $(0,1)$ and $(1,1)$. Then $\sigma$ is the identity on $\Bbb{A}_2(n)$.
Proof. To be supplied; similar geometrical arguments as before.$\hspace{20pt}\square$
Proof of Theorem 1. Let $\sigma$ be a symmetry of $\Bbb{A}_2(n)$ that preserves the set of rows, columns, diagonals and antidiagonals. Let $(i,j)=\sigma((0,0))$ and let $\tau$ denote the translation by $(-i,-j)$. Then the permutation $\tau\sigma$ fixes the point $(0,0)$, and $\tau$ is the unique translation for which this holds.
Because $\tau\sigma$ preserves lines we know that $\tau\sigma((0,1))$ is on a horizontal, vertical or diagonal line through $\tau\sigma((0,0))=(0,0)$. Then there is a unique $k\in\{0,1,2,3\}$ such that $\rho^k\tau\sigma((0,1))$ is on the horizontal line through $(0,0)$, and hence by Lemma 1 the restriction of $\rho^k\tau\sigma$ to the first row is given by multiplication by some $d\in(\Bbb{Z}/n\Bbb{Z})^{\times}$. Let $\gamma_d(i,j):=(di,dj)$ denote the scaling by $d$, which preserves the sets of rows, columns, diagonals and antidiagonals. Then $\gamma_d\rho^k\tau\sigma$ is the identity on the first row, and by the Corollary we have $$\gamma_d\rho^k\tau\sigma((1,1))=(\pm1,1).$$ Then there is a unique $l\in\{0,1\}$ such that $$\mu^l\gamma_d\rho^k\tau\sigma((1,1))=(1,1).$$ Now it follows from Lemma 2 that this map is the identity, and hence $$\sigma((x,y))=(\tau^{-1}\rho^{-k}\gamma_d^{-1}\mu^{-l})(x,y)=d^{-1}\rho^{-k}\mu^l\binom{x}{y}-\binom{i}{j},$$ for all $(x,y)\in\Bbb{A}_2(n)$. This shows that $\sigma$ is a symmetry of the form described in Theorem 1.$\hspace{20pt}\square$