This is a case of Power Group Enumeration with the group permuting
the slots being the eight symmetries $G_N$ of the $N\times N$ square
and the group acting on the $Q$ colors being the symmetric group
$S_Q$.
The cycle indices for $G_N$ were carefully documented and computed at
the following MSE link I.
The cycle index of the symmetric group can be computed from the
classical recurrence by Lovasz.
It then remains to apply the Power Group Enumeration formula /
algorithm as documented at the following
MSE link II.
We get for the case of coloring a square with at most two
interchangeable colors
$$1, 4, 51, 4324, 2105872, 4295327872, 35184441295872,
\\ 1152921514807410688,\ldots$$
which is [OEIS A182044](https://oeis.org/A182044)
where a closed formula can be found.
For at most three colors we get the sequence
$$1, 6, 490, 901012, 17653123147, 3126972103187548,
\\ 4985402694248850150928, 71535079589434063488162675274,
\\ 9238051838396620455005158025879826301,\ldots $$
Finally for at most four colors we get the sequence
$$1, 7, 1555, 22396971, 5864091026091, 24595658827938966187,
\\ 1650586719048786316922366635, 1772303994379887884341412962742479531,
\\ 30447950777727144129269702965544605972525918891,\ldots$$
The sequence of colorings using any number of available interchangeable colors is also quite interesting and starts
$$1, 7, 2966, 1310397193,
\\ 579823814813639193, 477464341236883456112705749206,
\\ 1340767144321669800049265230788088027597024910,
\\ 21516767919669856366796245458194341840929552797762722429430679631,
\ldots$$
An implementation of this algorithm is included below.
with(combinat);
pet_cycleind_symm :=
proc(n)
option remember;
local l;
if n=0 then return 1; fi;
expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;
pet_cycleind_grid :=
proc(N)
option remember;
local cind;
if type(N, even) then
cind :=
1/8*(a[1]^(N^2)+3*a[2]^(N^2/2)+
2*a[1]^N*a[2]^((N^2-N)/2) + 2*a[4]^(N^2/4));
else
cind :=
1/8*(a[1]^(N^2)+4*a[1]^N*a[2]^((N^2-N)/2)+
a[1]*a[2]^((N^2-1)/2) + 2*a[1]*a[4]^((N^2-1)/4));
fi;
cind;
end;
gridcols :=
proc(N, Q)
option remember;
local idx_slots, idx_cols, res, term_a, term_b,
v_a, v_b, inst_a, inst_b, len_a, len_b, p, q;
if N > 1 then
idx_slots := pet_cycleind_grid(N);
else
idx_slots := [a[1]];
fi;
if Q > 1 then
idx_cols := pet_cycleind_symm(Q);
else
idx_cols := [a[1]];
fi;
res := 0;
for term_a in idx_slots do
for term_b in idx_cols do
p := 1;
for v_a in indets(term_a) do
len_a := op(1, v_a);
inst_a := degree(term_a, v_a);
q := 0;
for v_b in indets(term_b) do
len_b := op(1, v_b);
inst_b := degree(term_b, v_b);
if len_a mod len_b = 0 then
q := q + len_b*inst_b;
fi;
od;
p := p*q^inst_a;
od;
res := res +
lcoeff(term_a)*lcoeff(term_b)*p;
od;
od;
res;
end;