The following data augment those in the first answer, namely we treat
the problem of computing $Z(Q_{n})$ for $K_{n,n}$ where we include
reflections. The algorithm here is slightly different from Harary and
Palmer. Clearly all permutations from the first answer represented in
$Z(Q_{n,n})$ contribute as well. Now for the reflections, which are
simple, fortunately. These (the vertex permutations) can be obtained
from the permutations in $S_n$ by placing vertices from component $A$
alternating with component $B$ on a cycle with twice the length of a
cycle contained in a permutation $\gamma$ of $S_n.$ The induced action
on the set of edges is the contribution to $Z(Q_n).$ In the following
we have used an algorithmic approach rather than the formula from
Harary and Palmer, namely we construct a canonical representative of
each permutation shape from the cycle index $Z(S_n)$ that doubles the
length of every cycle and alternates elements from the two
components. We let that act on the set of edges and factor the result.
This gives the following cycle indices, the first of which matches the
result given by Harary and Palmer.
For $n=3$,
$$Z(Q_3) =
{\frac {{a_{{1}}}^{9}}{72}}+1/6\,{a_{{1}}}^{3}{a_{{2}}}^{3}
\\+1/8\,a_{{1}}{a_{{2}}}^{4}+1/4\,a_{{1}}{a_{{4}}}^{2}+1/9\,{
a_{{3}}}^{3}+1/3\,a_{{3}}a_{{6}}
,$$
for $n=4$,
$$Z(Q_4) =
{\frac {{a_{{1}}}^{16}}{1152}}+{\frac {{a_{{1}}}^{8}{a_{{2}
}}^{4}}{96}}+{\frac {5\,{a_{{1}}}^{4}{a_{{2}}}^{6}}{96}}\\+{
\frac {{a_{{1}}}^{4}{a_{{3}}}^{4}}{72}}+{\frac {17\,{a_{{2}
}}^{8}}{384}}+1/12\,{a_{{1}}}^{2}a_{{2}}{a_{{3}}}^{2}a_{{6}
}+1/8\,{a_{{1}}}^{2}a_{{2}}{a_{{4}}}^{3}\\+1/18\,a_{{1}}{a_{{
3}}}^{5}+1/6\,a_{{1}}a_{{3}}{a_{{6}}}^{2}+1/24\,{a_{{2}}}^{
2}{a_{{6}}}^{2}+{\frac {19\,{a_{{4}}}^{4}}{96}}+1/12\,a_{{4
}}a_{{12}}+1/8\,{a_{{8}}}^{2}
,$$
and for $n=5$,
$$Z(Q_5) =
{\frac {{a_{{1}}}^{25}}{28800}}+{\frac {{a_{{1}}}^{15}{a_{{
2}}}^{5}}{1440}}+{\frac {{a_{{1}}}^{9}{a_{{2}}}^{8}}{288}}+
{\frac {{a_{{1}}}^{10}{a_{{3}}}^{5}}{720}}+{\frac {{a_{{1}}
}^{5}{a_{{2}}}^{10}}{192}}\\+{\frac {{a_{{1}}}^{3}{a_{{2}}}^{
11}}{96}}+{\frac {a_{{1}}{a_{{2}}}^{12}}{128}}+{\frac {{a_{
{1}}}^{6}{a_{{2}}}^{2}{a_{{3}}}^{3}a_{{6}}}{72}}+{\frac {{a
_{{1}}}^{4}{a_{{3}}}^{7}}{72}}\\+{\frac {{a_{{1}}}^{5}{a_{{4}
}}^{5}}{480}}+1/24\,{a_{{1}}}^{3}{a_{{2}}}^{3}{a_{{4}}}^{4}
+{\frac {{a_{{2}}}^{5}{a_{{3}}}^{5}}{720}}+1/48\,{a_{{1}}}^
{3}a_{{2}}{a_{{4}}}^{5}\\+1/48\,{a_{{1}}}^{2}{a_{{2}}}^{4}a_{
{3}}{a_{{6}}}^{2}+{\frac {{a_{{2}}}^{5}{a_{{3}}}^{3}a_{{6}}
}{72}}+1/32\,a_{{1}}{a_{{2}}}^{2}{a_{{4}}}^{5}\\+1/48\,{a_{{2
}}}^{5}a_{{3}}{a_{{6}}}^{2}+1/36\,{a_{{2}}}^{2}{a_{{3}}}^{5
}a_{{6}}+1/12\,{a_{{1}}}^{2}a_{{2}}a_{{3}}{a_{{6}}}^{3}\\+{
\frac {3\,a_{{1}}{a_{{4}}}^{6}}{32}}+{\frac {{a_{{2}}}^{2}{
a_{{3}}}^{3}{a_{{6}}}^{2}}{72}}+1/24\,{a_{{1}}}^{2}a_{{3}}{
a_{{4}}}^{2}a_{{12}}+1/24\,a_{{2}}a_{{3}}{a_{{4}}}^{2}a_{{
12}}\\+{\frac {13\,{a_{{5}}}^{5}}{600}}+1/8\,a_{{1}}{a_{{8}}}
^{3}+1/12\,a_{{3}}a_{{4}}a_{{6}}a_{{12}}+{\frac {{a_{{5}}}^
{3}a_{{10}}}{60}}+1/30\,{a_{{5}}}^{2}a_{{15}}\\+1/8\,a_{{5}}{
a_{{10}}}^{2}+1/20\,a_{{5}}a_{{20}}+1/30\,a_{{10}}a_{{15}}
.$$
These cycle indices can be calculated for large $n$ but the pattern
should be clear. The count of these graphs gives the sequence
$$2, 6, 26, 192, 3014, 127757, 16853750, \\ 7343780765, 10733574184956,
52867617324773592,\\882178116079222400788, 50227997322550920824045262,
\ldots $$
which points us to OEIS A007139 where
we find confirmation of this calculation.
This was the Maple code used to obtain these values.
with(numtheory);
pet_cycleind_symm :=
proc(n)
local p, s;
option remember;
if n=0 then return 1; fi;
expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;
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 := Array([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;
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;
pet_flatten_term :=
proc(varp)
local terml, d, cf, v;
terml := [];
cf := varp;
for v in indets(varp) do
d := degree(varp, v);
terml := [op(terml), seq(v, k=1..d)];
cf := cf/v^d;
od;
[cf, terml];
end;
pet_flat2rep_cyc :=
proc(f)
local p, q, res, cyc, t, len;
q := 1; res := [];
for t in f do
len := op(1, t);
cyc := [seq(p, p=q+1..q+len-1), q];
res := [op(res), cyc];
q := q+len;
od;
res;
end;
pet_cycs2table :=
proc(cycs)
local pairs, cyc, p, ent;
pairs := [];
for cyc in cycs do
pairs :=
[op(pairs),
seq([cyc[p], cyc[1 + (p mod nops(cyc))]],
p = 1 .. nops(cyc))];
od;
map(ent->ent[2], sort(pairs, (a,b)->a[1] < b[1]));
end;
pet_cycleind_knn :=
proc(n)
option remember;
local cindA, cindB, sind, t1, t2, p, q, cyc1, cyc2,
flat, len, len1, len2, v1, v2,
edges, edgeperm, rep, cycsA, cycsB, cycsrc, cyc;
if n=1 then
sind := [a[1]];
else
sind := pet_cycleind_symm(n);
fi;
cindA := 0;
for t1 in sind do
for t2 in sind do
q := 1;
for v1 in indets(t1) do
len1 := op(1, v1);
for v2 in indets(t2) do
len2 := op(1, v2);
len := lcm(len1, len2);
q := q *
a[len]^((len1*len2/len) *
degree(t1, v1)*degree(t2, v2));
od;
od;
cindA := cindA +
lcoeff(t1)*lcoeff(t2)*q;
od;
od;
edges := [seq(seq({p, q+n}, q=1..n), p=1..n)];
cindB := 0;
for t1 in sind do
flat := pet_flatten_term(t1);
cycsA := pet_flat2rep_cyc(flat[2]);
cycsB := [];
for cycsrc in cycsA do
cyc := [];
for q in cycsrc do
cyc := [op(cyc), q, q+n];
od;
cycsB := [op(cycsB), cyc];
od;
rep := pet_cycs2table(cycsB);
edgeperm :=
subs([seq(q=rep[q], q=1..2*n)], edges);
cindB := cindB + flat[1]*
pet_autom2cycles(edges, edgeperm);
od;
(cindA+cindB)/2;
end;
Q :=
proc(n)
option remember;
local cind;
cind := pet_cycleind_knn(n);
subs([seq(a[p]=2, p=1..n*n)], cind);
end;
Remarks, per request, Dec 2020. Here are some observations
concerning reflections. We need the factorization into cycles of the
vertices in order to compute the action on the edges. Now a reflection
when factored into cycles must alternate between left and right
vertices. That means left and right vertices are interleaved on those
cycles. Hence we can obtain the factorizations of all reflections by
taking a permutation from $S_n$ (which gives the left vertices) and
doubling all its cycles in length, inserting a permutation of the
right vertices, for a total of $n! \times n!$ reflections. The first
factorial comes from the permutations of the left vertices which can
be interleaved with the right ones in $n!$ ways, yielding the second
factorial. Now consider the factorization into cycles of the action on
the edges. Left and right vertices are disjoint, so no matter the
order in which we interleave the right vertices, we get the same shape
of factorizations into cycles (just permute the right vertices, which
maintains the cycle structure of the edges). Therefore we can take one
representative for each term in $Z(S_n)$, interleave it with some
permutation of the right, and factor that (i.e. factor the action on
the edges) to get the contribution to $Z(Q_n)$, which must be
multiplied by $n!$ for the number of reflections covered by this
term. This is what the algorithm does. The $n!$ is canceled when we
divide by the total number of permutations in the group.
Best Answer
Hint. An induced subgraph is uniquely determined by its vertices, that is $G$ has as many induced subgraphs as $V$ has subsets. Likewise the spanning subgraphs correspond to subsets of $E$.