This can indeed be done with the Polya Enumeration Theorem. There is a
very simple yet powerful observation that can be made here: if we have
the cycle index of the edge permutation group induced by the action of
the symmetric group on the vertices, the cycle index of the edge
permutation group including self-loops can be obtained by multiplying
the factorizations into disjoint cycles of the edge permutation and
the original vertex permutation and sum these contributions for all
vertex permutations. (Of course we work directly with the cycle index
of the symmetric group which is much more efficient than iterating
over all permutations.) But the cycle index of the edge permutation
group is not difficult to compute, the code can be found at this MSE
Link.
This gives the following cycle indices for the enumeration of graphs
that may have self-loops: for $n=2$ we have
$$1/2\,{a_{{1}}}^{3}+1/2\,a_{{1}}a_{{2}},$$
for $n=3$ we get
$$1/6\,{a_{{1}}}^{6}+1/2\,{a_{{1}}}^{2}{a_{{2}}}^{2}
+1/3\,{a_{{3}}}^{2},$$
and finally for $n=4$ we obtain
$$1/24\,{a_{{1}}}^{10}+1/4\,{a_{{1}}}^{4}{a_{{2}}}^{3}+
1/3\,{a_{{3}}}^{3}a_{{1}}+1/8\,{a_{{1}}}^{2}{a_{{2}}}^{4}
+1/4\,a_{{2}}{a_{{4}}}^{2}.$$
As an example we substitute $1+z$ into the cycle index for $n=3$ to get
$${z}^{6}+2\,{z}^{5}+4\,{z}^{4}+6\,{z}^{3}+4\,{z}^{2}+2\,z+1.$$
There are indeed two graphs on three vertices with $1$ edge including
loops, namely a single edge or a single loop. Similarly there are four
graphs with $2$ edges, namely two loops, two edges, an edge with a
loop on one end and an edge that is opposed to a vertex with a loop on
it. Furthermore the six graphs with three edges are: a triangle with
no loops, a fork with a loop on the handle, a fork with a loop on one
of its ends, a line segment with two loops at either end, a line
segment with a loop on one end and a loop not incident on the line
segment, and finally, three loops. It goes on like this.
These cycle indices produce the following total counts of
non-isomorphic undirected graphs with possible self-loops present:
1
6
20
90
544
5096
79264
2208612
113743760
10926227136
1956363435360
652335084592096
405402273420996800
470568642161119963904
1023063423471189431054720
4178849203082023236058229792
32168008290073542372004082199424
468053896898117580623237189908861696
12908865186281154904787051023018388368640
676599895346467962508983189158040730734206464
67557268911383303274368343026941404659790140175872
12878644470123443279570180725329554086442149175832124416
4696891990987444795146988875290693218960182452250871238964224
3283275444870880220030739747435965751837707129958501790119044590592
4406671511641658245922265014648856899986657018190959680004287879813376000
11374011362523188765310166602160674880959112891073505609667787021125321629659136
This sequence of values points us to OEIS
A000666 where a large amount of additional
material can be found.
For those who are curious to see the details, here is the Maple code
to compute these values:
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_edg_loops :=
proc(n)
option remember;
local all, term, termvars, res, l1, l2, inst1, u, v,
uidx, vidx;
if n=0 then return 1; fi;
if n=1 then return a[1]; fi;
all := 0:
for term in pet_cycleind_symm(n) do
termvars := indets(term); res := 1;
# edges on different cycles of different sizes
for uidx to nops(termvars) do
u := op(uidx, termvars);
l1 := op(1, u);
for vidx from uidx+1 to nops(termvars) do
v := op(vidx, termvars);
l2 := op(1, v);
res := res *
a[lcm(l1, l2)]
^((l1*l2/lcm(l1, l2))*
degree(term, u)*degree(term, v));
od;
od;
# edges on different cycles of the same size
for u in termvars do
l1 := op(1, u); inst1 := degree(term, u);
# a[l1]^(1/2*inst1*(inst1-1)*l1*l1/l1)
res := res *
a[l1]^(1/2*inst1*(inst1-1)*l1);
od;
# edges on identical cycles of some size
for u in termvars do
l1 := op(1, u); inst1 := degree(term, u);
if type(l1, odd) then
# a[l1]^(1/2*l1*(l1-1)/l1);
res := res *
(a[l1]^(1/2*(l1-1)))^inst1;
else
# a[l1/2]^(l1/2/(l1/2))*a[l1]^(1/2*l1*(l1-2)/l1)
res := res *
(a[l1/2]*a[l1]^(1/2*(l1-2)))^inst1;
fi;
od;
all := all + term*res;
od;
all;
end;
pet_varinto_cind :=
proc(poly, ind)
local subs1, subs2, polyvars, indvars, v, pot, res, k;
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;
VGF :=
proc(n)
option remember;
local ind;
ind := pet_cycleind_edg_loops(n);
expand(pet_varinto_cind(1+z, ind));
end;
COUNT :=
proc(n)
option remember;
local q;
subs([seq(a[q]=2, q=1..n*(n+1)/2)],
pet_cycleind_edg_loops(n));
end;
