I want to find a closed form to the bivariate generating function $$ G(x, y) = \sum\limits_{i, j} \binom{i+j}{i}^2 x^i y^j. $$ Ideally, I would prefer a direct approach that is based on the definition above.
I know that there is a closed form here, as one can reduce the summation to Legendre polynomials: $$ G(x, y) = \sum\limits_n y^n \sum\limits_k \binom{n}{k}^2 \left(\frac{x}{y}\right)^k = \sum\limits_n (y-x)^n P_n \left(\frac{y+x}{y-x}\right), $$ where $P_n(x)$ is the $n$-th Legendre polynomial. From this, using the generating function formula $$ \sum\limits_n P_n(x) t^n = \frac{1}{\sqrt{1-2xt+t^2}}, $$ we get the closed-form expression for $G(x, y)$ as $$ \boxed{G(x, y) = \frac{1}{\sqrt{1-2(y+x)+(y-x)^2}}} $$ But I totally fail to see any meaningful way to derive it in a more direct and self-contained way. Any hints? And while we're at it, are there similar closed-form expressions for higher powers of $\binom{n}{k}$?