Let me post a pointer. The following MSE
link from eleven
months ago features orbital chromatic polynomials, which count
proper colorings of a graph under the symmetries of its
automorphisms. There is extensive documentation at that link. The code
that was posted there is easy to apply here: the underlying graph with
edges for adjacency this the three-by-three grid graph. We encode it
as follows:
SQUARE3BY3 :=
proc()
option remember;
return
[9,
{{1, 2}, {2, 3},
{4, 5}, {5, 6},
{7, 8}, {8, 9},
{1, 4}, {2, 5}, {3, 6},
{4, 7}, {5, 8}, {6, 9}},
[[1,2,3,4,5,6,7,8,9], # identity
[3,6,9,2,5,8,1,4,7], # 90 degrees
[7,4,1,8,5,2,9,6,3], # -90 degrees
[9,8,7,6,5,4,3,2,1], # 180 degrees
[7,8,9,4,5,6,1,2,3], # horizontal flip
[3,2,1,6,5,4,9,8,7], # vertical flip
[1,4,7,2,5,8,3,6,9], # falling diagonal
[9,6,3,8,5,2,7,4,1]]]; # rising diagonal
end;
The Maple command OCP(SQUARE3BY3()); then immediately yields the
OCP:
$$P(k) = 1/8\,{k}^{9}+8\,k-{\frac {133\,{k}^{2}}{4}}-3/2\,{k}^{8}
+{\frac {33\,{k}^{7}}{4}}-{\frac {53\,{k}^{6}}{2}}
+{\frac {217\,{k}^{5}}{4}}-{\frac {291\,{k}^{4}}{4}}
+{\frac {507\,{k}^{3}}{8}}.$$
This yields for up to twelve colors the sequence
$$0, 2, 69, 1572, 19865, 153480, 830802, 3476144, 12003462,
\\35757630, 94780235, 228579252, \ldots$$
which confirms the value for three colors that was first to appear.
Remark, as per comments. The value $P(k)$ of this OCP counts
the number of colorings using at most $k$ colors. We can compute
$P'(k)$ which gives the proper colorings using exactly $k$ colors by
inclusion-exclusion. Here the nodes of the poset are subsets $Q$ of
$[k]$ representing proper colorings using some subset of the colors in
$Q$, which is counted by $P(|Q|).$ A coloring using exactly the colors
from some set $R$ is represented by all nodes corresponding to
supersets $Q$ of $R.$ With the weight being $(-1)^{k-|Q|}$ those
colorings using exactly $k$ colors only occur at $Q=[k]$ with weight
$(-1)^{k-|Q|} = (-1)^{0} = 1.$ A coloring using exactly $R\subset [k]$
colors is represented by all $Q$ such that $R \subseteq Q \subseteq
[k]$, with total weight
$$\sum_{R'\subseteq [k] \setminus R} (-1)^{k-|R\cup R'|}
= (-1)^{k-|R|}
\sum_{r=0}^{|[k]\setminus R|} {|[k]\setminus R| \choose r} (-1)^{-r}
= 0.$$
Hence only the colorings with exactly $k$ colors contribute and we
find
$$P'(k) = \sum_{Q\subseteq [k]} (-1)^{k-|Q|} P(|Q|) =
\sum_{q=0}^k {k\choose q} (-1)^{k-q} P(q).$$
This gives the finite sequence
$$0, 2, 63, 1308, 12675, 56520, 120960, 120960, 45360, 0, \ldots$$
because it is clearly impossible to color the grid using more than
nine distinct colors. Observe also the entry for three colors, which
is $P(3) - {3\choose 2} P(2) = 69 - 3\times 2$ i.e. we have subtracted
the colorings using two colors (there are no coloring using one color
and hence $P(2)$ counts colorings with exactly two colors). Also note
that with nine colors all orbits have the same size, namely eight, and
indeed we obtain $9!/8 = 45360.$ As to what happens when there are
more than nine colors we can recover $P(k)$ as follows:
$$\sum_{q=0}^9 {k\choose q} P'(q).$$
Addendum. The reader might be interested to know that we can
compute the OCP for larger grids, using the following code in
conjunction with the quoted link:
SQUARE :=
proc(n)
option remember;
local src, rot, automs, edges, v2n;
src := [seq(seq([p, q], q=0..n-1), p=0..n-1)];
edges :=
{seq(seq({[p, q], [p+1, q]},
p=0..n-2), q=0..n-1),
seq(seq({[p, q], [p, q+1]},
p=0..n-1), q=0..n-2)};
rot := v -> [v[2], n-1-v[1]];
automs :=
[src, # identity
map(rot, src), # 90 degrees
map(v -> rot(rot(v)), src), # 180 degrees
map(v -> rot(rot(rot(v))), src), # 270 degrees
map(v -> [n-1-v[1], v[2]], src), # horizontal flip
map(v -> [v[1], n-1-v[2]], src), # vertical flip
map(v -> rot([n-1-v[1], v[2]]),
src), # rising diagonal
map(v -> rot(rot(rot([n-1-v[1], v[2]]))),
src)]; # falling diagonal
v2n :=
[seq(seq([p, q] = 1 + p*n + q, q=0..n-1), p=0..n-1)];
[n*n, subs(v2n, edges), subs(v2n, automs)];
end;
We obtain for a four-by-four the OCP
$$1/8\,{k}^{16}-3\,{k}^{15}+{\frac {69\,{k}^{14}}{2}}
-{\frac {2015\,{k}^{13}}{8}}+{\frac {10437\,{k}^{12}}{8}}
\\-{\frac {20307\,{k}^{11}}{4}}+15333\,{k}^{10}-{\frac {292907\,{k}^{9}}{8}}
-{\frac {848501\,{k}^{7}}{8}}+{\frac {1023195\,{k}^{6}}{8}}
\\-{\frac {240539\,{k}^{5}}{2}}+{\frac {557915\,{k}^{8}}{8}}
-{\frac {8807\,k}{4}}+{\frac {112831\,{k}^{2}}{8}}
+{\frac {683997\,{k}^{4}}{8}}-{\frac {347485\,{k}^{3}}{8}}$$
with the sequence
$$0, 1, 1155, 759759, 103786510, 4767856260, 107118740001, \ldots$$
We get for a five-by-five the OCP
$$1/8\,{k}^{25}+{\frac {69997383\,{k}^{17}}{8}}-5\,{k}^{24}
+{\frac {195\,{k}^{23}}{2}}-1233\,{k}^{22}+{\frac {45399\,{k}^{21}}{4}}
\\-80919\,{k}^{20}+{\frac {928545\,{k}^{19}}{2}}
-{\frac {17590911\,{k}^{18}}{8}}-{\frac {118477969\,{k}^{16}}{4}}
+{\frac {172111059\,{k}^{15}}{2}}
\\-{\frac {1726958987\,{k}^{14}}{8}}
+{\frac {3754019329\,{k}^{13}}{8}}-{\frac {1770719251\,{k}^{12}}{2}}
\\+{\frac {5797425049\,{k}^{11}}{4}}-2053661272\,{k}^{10}
+{\frac {20055169857\,{k}^{9}}{8}}+{\frac {9236896437\,{k}^{7}}{4}}
\\-{\frac {6780818949\,{k}^{6}}{4}}+{\frac {8083053959\,{k}^{5}}{8}}
-{\frac {20932696169\,{k}^{8}}{8}}+4017958\,k
\\-{\frac {145271789\,{k}^{2}}{4}}-{\frac {3768579695\,{k}^{4}}{8}}
+{\frac {1292510453\,{k}^{3}}{8}}$$
with the sequence
$$0, 2, 76332, 2557101612, 6352711134515, 2747239197568620,
\\378972203462839707, \ldots$$