The result for any positive integer n is the sum of two terms: (1) an Exp[x y] term multiplied by a polynomial in x and y and (2) a polynomial in x and `y'. The first term is

Exp[x  y] LaguerreL[n, (x - y)^2]/n!

but I haven't figured out a general form for the second polynomial (assuming there is a closed form for that polynomial). So far for that second polynomial the sum of the x and y exponents is always even and the coefficient for x^i y^j is the same as the coefficient for x^j y^i. That's all I've determined. Maybe someone else can see the pattern.

From here we have:

$L_p^l\left(a x^2\right)=\sum _{k=0}^p \frac{(-1)^k \binom{p+l}{p-k} \left(a x^2\right)^k}{k!}$

Then we have:

.

$\sum _{l=0}^{\infty } \frac{(c x y)^l L_p^l\left(a x^2\right) L_p^l\left(b y^2\right)}{(l+p)!}=\\\sum _{l=0}^{\infty } \frac{(c x y)^l \left(\sum _{k=0}^p \frac{(-1)^k \binom{p+l}{p-k} \left(a x^2\right)^k}{k!}\right) \sum _{j=0}^p \frac{(-1)^j \binom{p+l}{p-j} \left(b y^2\right)^j}{j!}}{(l+p)!}=\\\sum _{k=0}^p \sum _{j=0}^p \left(\sum _{l=0}^{\infty } \frac{(-1)^{j+k} \left(a x^2\right)^k (c x y)^l \left(b y^2\right)^j \binom{l+p}{-j+p} \binom{l+p}{-k+p}}{j! k! (l+p)!}\right)=\\\sum _{k=0}^p \sum _{j=0}^p \frac{(-1)^{j+k} \left(a x^2\right)^k \left(b y^2\right)^j \binom{p}{-j+p} \binom{p}{-k+p} \, _2F_2(1,1+p;1+j,1+k;c x y)}{j! k! p!}$

Going from a single sum to a triple sum will probably not help me. I would need a closed form solution because this will then go inside an integral.

I believe there probably will be a closed form. I saw the Hardy-Hille formula on Wikipedia:https://en.wikipedia.org/wiki/Laguerre_polynomials which is of the similar form, just with an extra factorial in the numerator. Further, when specifying the p values from 0 to 10, Mathematcia does give closed form expressions for those particular p values.

I need closed form expression because this is a small part of my project.