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
Pick a base vertex $v$. For any other vertex $w$, pick a path $P(v,w)$ from $v$ to $w$. $P(v,w)$ exists because $G$ is connected.
For a spanning subgraph $H\subseteq_{sp} G$, we speak of "flipping an edge $e$" to mean removing $e$ from $H$ if $e\in H$, and adding $e$ to $H$ if $e$ is not in $H$. This transformation yields another spanning subgraph. Note that flipping all edges along a path changes the degree of the end vertices of the path by one (provided the path is not a loop), but not any other vertices.
Define an equivalence relation $\sim$ on the set of all spanning subgraphs of $G$ as follows: $H_1\sim H_2$ if $H_2$ is obtained from $H_1$ by flipping all the edges along the path $P(v,w)$ for some $w$. Thus each equivalence class contains $2^{n-1}$ spanning subgraphs, as we have $n-1$ paths $P(v,w)$ to flip.
For any spanning subgraph $H$, there is a unique equivalent graph with even degree for all vertices $w\neq v$ (just flip all paths $P(v,w)$ for which $w$ has odd degree). Call this equivalent subgraph $H'$. In $H'$, $v$ must also have even degree, because the sum of degrees for all vertices is even (being equal to $2\times\text{number of edges}$), and the sum of degrees of all vertices $w\neq v$ is also even by hypothesis. Hence, each equivalence class has a unique even spanning subgraph.
There are $2^e$ spanning subgraphs, and thus $2^{e-n+1}$ equivalence classes, and thus $2^{e-n+1}$ even spanning subgraphs.