This can be done using Polya's theorem, but we need the cycle index of the face permutation group $F$ (once we have this index we could also use Burnside).
We proceed to enumerate the permutations of this group. There is the identity, which contributes $$a_1^6.$$ There are three rotations for each pair of opposite faces that fix those faces (rotate about the axis passing through the center of the two faces). They contribute
$$3\times (2 a_1^2 a_4 + a_1^2 a_2^2).$$
There are rotations about an axis passing through opposite vertices, of which there are four pairs, giving
$$4\times 2 a_3^2.$$
Finally we may rotate about an axis passing through the centers of opposite edges and there are six of these, giving
$$6\times a_2^3.$$
It follows that the cycle index of $F$ is given by
$$Z(F) = \frac{1}{24}
\left(a_1^6 + 6a_1^2a_4 + 3a_1^2a_2^2 + 8a_3^2 + 6a_2^3\right).$$
With four colors we have
$$Z(F)(W+X+Y+Z) \\=
1/24\, \left( W+X+Y+Z \right) ^{6}\\+1/4\, \left( W+X+Y+Z \right) ^{2} \left( {W}^
{4}+{X}^{4}+{Y}^{4}+{Z}^{4} \right) \\+1/8\, \left( W+X+Y+Z \right) ^{2} \left( {W
}^{2}+{X}^{2}+{Y}^{2}+{Z}^{2} \right) ^{2}\\+1/3\, \left( {W}^{3}+{X}^{3}+{Y}^{3}+
{Z}^{3} \right) ^{2}\\+1/4\, \left( {W}^{2}+{X}^{2}+{Y}^{2}+{Z}^{2} \right) ^{3}.$$
This ordinary generating function provides more information than Burnside alone, namely through its coeffcients. We have, for example, that
$$[WXYZ^3] Z(F)(W+X+Y+Z) = 5,$$
which says that there are five colorings where colors $X, Y$ and $Z$ ocurr exactly once and color $Z$ three times. Obviously we also have
$$[WXY^3Z] Z(F)(W+X+Y+Z) = 5,$$
by symmetry. Furthermore, we have
$$[W^2X^2Y^2] Z(F)(W+X+Y+Z) = 6,$$
hence there are six cubes with colors $W,X,Y$ each ocurring twice. It is a very useful exercise to pick out one of these coefficients and verify its value using pen and paper.
Evaluating the substituted cycle index at $W=1, X=1, Y=1, Z=1$ we get
$$Z(F)(W+X+Y+Z)_{W=1, X=1, Y=1, Z=1} = 240.$$
The sequence of colorings when there are $N$ colors is as follows:
$$1, 10, 57, 240, 800, 2226, 5390, 11712, 23355, 43450, \ldots.$$
This is sequence A047780 from the OEIS.
Another related cycle index is the cycle index of the permutation group $G$ obtained when we add reflections to the admissible transformations of the cube. This yields the full symmetry group of the regular octrahedron.
$$ Z(G) =
1/48\,{a_{{1}}}^{6}+1/16\,{a_{{1}}}^{4}a_{{2}}+3/16\,{a_{{1}}}^{2}{a_{{2}}}^{2}+
1/8\,{a_{{1}}}^{2}a_{{4}}\\+{\frac {7}{48}}\,{a_{{2}}}^{3}+1/6\,{a_{{3}}}^{2}+1/6
\,a_{{6}}+1/8\,a_{{4}}a_{{2}}.$$
The corresponding sequence is
$$1, 10, 56, 220, 680, 1771, 4060, 8436, 16215, 29260,\ldots$$
which is A198833 from the OEIS.
This cycle index can be computed algorithmically without the need to classify the different types of symmetries. It suffices to encode the adjacencies of the faces and select those permutations of the symmetric group that turn out to be automorphisms.
This is the code. The reader may want to attempt the classification without a CAS.
As these groups contain more and more elements an algorithm can be useful.
with(numtheory);
with(group):
with(combinat):
pet_autom2cycles :=
proc(src, aut)
local numa, numsubs;
local marks, pos, cycs, cpos, clen;
numsubs := [seq(src[k]=k, k=1..nops(src))];
numa := subs(numsubs, aut);
marks := [seq(true, pos=1..nops(aut))];
cycs := []; pos := 1;
while pos <= nops(aut) do
if marks[pos] then
clen := 0; cpos := pos;
while marks[cpos] do
marks[cpos] := false;
cpos := numa[cpos];
clen := clen+1;
od;
cycs := [op(cycs), clen];
fi;
pos := pos+1;
od;
return mul(a[cycs[k]], k=1..nops(cycs));
end;
f := {{1,2},{1,3},{1,4},{1,5},
{2,3},{2,5},{2,6},
{3,4},{3,6},
{4,5},{4,6},
{5,6}};
cube_cind :=
proc()
option remember;
local ff, p, res, count, term;
count := 0; res := 0;
for p in permute(6) do
ff := subs([seq(k=p[k], k=1..6)], f);
if f = ff then
count := count+1;
term := pet_autom2cycles([seq(k,k=1..6)], p);
res := res+ term;
fi;
od;
print(count);
res/count;
end;
pet_varinto_cind :=
proc(poly, ind)
local subs1, subs2, polyvars, indvars, v, pot, res;
res := ind;
polyvars := indets(poly);
indvars := indets(ind);
for v in indvars do
pot := op(1, v);
subs1 :=
[seq(polyvars[k]=polyvars[k]^pot,
k=1..nops(polyvars))];
subs2 := [v=subs(subs1, poly)];
res := subs(subs2, res);
od;
res;
end;
v :=
proc(n)
option remember;
local p, k, gf;
p := add(cat(q, k), k=1..n);
gf := expand(pet_varinto_cind(p, cube_cind()));
subs({seq(cat(q, k)=1, k=1..n)}, gf);
end;
Note that we can also compute the cycle index of $F$ algorithmically. In order to accomplish this we work with the sets of faces incident on each other at a vertex (there are eight of these sets) and orient them by putting a clockwise spin on every vertex, thereby obtaining exactly one cycle of three faces per vertex. This spin is inverted by reflections, effectively removing them from the set of symmetries. (This is the moment to get pen and paper and verify the set of oriented vertex/face cycles.) The code is mostly like what we saw above, except for the following differences.
f := {[2,1,5], [3,1,2], [4,1,3],
[5,1,4], [6,2,5], [6,3,2],
[6,4,3], [6,5,4]
};
minfirst :=
proc(l)
local minpos, minval, pos;
minval := l[1]; minpos := 1;
for pos from 2 to nops(l) do
if l[pos]<minval then
minval := l[pos];
minpos := pos;
fi;
od;
[seq(l[k], k=minpos..nops(l)),
seq(l[k], k=1..minpos-1)];
end;
cube_cind :=
proc()
option remember;
local fA, fB, p, res, count, term;
fA := map(minfirst, f);
count := 0; res := 0;
for p in permute(6) do
fB := map(minfirst, subs([seq(k=p[k], k=1..6)], fA));
if fA = fB then
count := count+1;
term := pet_autom2cycles([seq(k,k=1..6)], p);
res := res+ term;
fi;
od;
print(count);
res/count;
end;
