# Series applied to an infinite sum does not work

When I have

Sum[(-1)^n x^(n^2) y^n, {n, 0, ∞}]


and I try evaluating

Series[Sum[(-1)^n x^(n^2) y^{n}, {n, 0, ∞}], {x, 0, 3}, {y, 0, 2}]


it does not return the terms in the series that I want.

I am not interested in changing my Sum expression to be finite as in, say, {n, 0, 3} instead of {n, 0, ∞}, because I will be making many such functions defined as infinite series and will want to constantly manipulate them using algebraic operations, so I will be generating multiple series. This also means I need to keep the operation Series, since I will be feeding, e.g., products of such sums with other functions.

It does not work , because MMA dosen't know a closed-form solution for the Sum. Workaround:

Only for: -1 < x < 1, -1 < y < 1

func = Sum[MellinTransform[(-1)^n*y^n*x^(a n^2), a, s] // PowerExpand , {n,
0, \[Infinity]}, Assumptions -> {-1 < x < 1, -1 < y < 1, s > 0}]
sol = Simplify[Series[func, {x, 0, 8}, {y, 0, 8}] // Normal] // Expand;
sol2 = 1 + InverseMellinTransform[Total@Table[sol[[n]], {n, 2, Length[sol]}], s, a] /. a -> 1

(*1 + x y (-1 + x^3 y - x^8 y^2 + x^15 y^3 - x^24 y^4 + x^35 y^5 -
x^48 y^6 + x^63 y^7)*)


Check:

 N[sol2 /. x -> -1/2 /. y -> 1/2, 20]
(*1.2658700952308663684*)
NSum[(-1)^n*y^n*x^(n^2) /. x -> -1/2 /. y -> 1/2 , {n, 0, \[Infinity]}, WorkingPrecision -> 20]
(*1.2658700952308663684*)