Derive the recurrence relations

Solution 1:

Your polynomials are crying out to be put into a generating function, $$P = \sum_{m, n \geq 0} P_{m, n} \frac{u^m}{m!} \frac{v^n}{n!}.$$ Working backwards from your explicit formula, set $N = m+n-i-j$, $m' = m-2i, n' = n-2j$ and do routine simplifications to get $$ \begin{align*} P &= \sum_{N \geq 0} \sum_{(i, j, m', n') \vDash N} \binom{N}{i,j,m',n'} (2x)^{m'} (2y)^{n'} (-u^2)^i (-v^2)^j u^{m'} v^{n'} \\ &= \sum_{N \geq 0} (2xu + 2yv - u^2 - v^2)^N \\ &= \frac{1}{1 + u^2 + v^2 - 2xu - 2yv}. \end{align*} $$ That is, $$P_{m, n} = m! n! [u^m v^n] \frac{1}{1 + u^2 + v^2 - 2xu - 2yv}.$$

Sometimes in practice you'd be able to prove your interpretation satisfies the above generating function formula independently, which would prove the explicit formula by running the preceding argument backwards.

In any case, your initial recurrences can be packaged up into PDE's involving the $P$ generating function in a standard way. (See e.g. Wilf's generatingfunctionology for a nice, though unusual, treatment.) For example, $$\frac{\partial}{\partial u} P = \sum_{m, n \geq 0} P_{m+1,n} \frac{u^m}{m!} \frac{v^n}{n!}$$ and $$\frac{\partial}{\partial u} u P = \sum_{m, n \geq 0} (m+1) P_{m, n} \frac{u^m}{m!} \frac{v^n}{n!}.$$

Using these observations, your first recurrence relation, $$n \partial_x P_{m, n-1} = m \partial_y P_{m-1, n},$$ can be repackaged into $$\partial_{uvx} vP = \partial_{uvy} uP.$$ It's of course easy to have a computer verify this holds with the formula above.

Your five differential recurrence relations translate to $$ \begin{align*} \partial_{uvx} vP &= \partial_{uvy} uP \\ x \partial_{xu} P &= \partial_u u (\partial_u + \partial_x) P \\ y \partial_{xuv} P &= \partial_{uv^2} uP + \partial_{xuv} vP \\ y \partial_{vy} P &= \partial_v v (\partial_v + \partial_y) P \\ x \partial_{yuv} P &= \partial_v v \partial_{u^2} P + \partial_u u \partial_y. \end{align*} $$ Your "pure" differential equation translates to $$(1-x^2) \partial_{x^2} P - xy \partial_{xy} P - x \partial_{xv} vP - 2x \partial_x P + uy \partial_{uy} P + u \partial_u (\partial_u u + \partial_v v) P = 0.$$ Your "pure" recurrence relations translate to $$ \begin{align*} 2x (\partial_u u + \partial_v v + 2) \partial_{uv^2} P &= \partial_{u^2 v^2} P - (1+\partial_v v) \partial_{u^2} P + \partial_u u (\partial_u u + 2 \partial_v v + 3) \partial_{v^2} P \\ 2y (\partial_v v + \partial_u u + 2) \partial_{vu^2} P &= \partial_{v^2 u^2} P - (1+\partial_u u) \partial_{v^2} P + \partial_v v (\partial_v v + 2 \partial_u u + 3) \partial_{u^2} P. \end{align*} $$ All 8 of these relations hold for $P$ above.

These identities correspond to elements of the annihilator of $P$ in the Weyl algebra on 4 variables. These elements are, respectively, $$ \begin{align*} &\partial_{uvx} v - \partial_{uvy} u \\ &x \partial_{xu} - \partial_u u (\partial_u + \partial_x) \\ &y \partial_{xuv} - \partial_{uv^2} u - \partial_{xuv} v \\ &y \partial_{vy} - \partial_v v (\partial_v + \partial_y) \\ &x \partial_{yuv} - \partial_v v \partial_{u^2} - \partial_u u \partial_{yv} \\ &(1-x^2) \partial_{x^2} - xy \partial_{xy} - x \partial_{xv} v - 2x \partial_x + uy \partial_{uy} + u \partial_u (\partial_u u + \partial_v v) \\ &2x (\partial_u u + \partial_v v + 2) \partial_{uv^2} - \partial_{u^2 v^2} + (1+\partial_v v) \partial_v v \partial_{u^2} - \partial_u u (\partial_u u + 2 \partial_v v + 3) \partial_{v^2} \\ &2y (\partial_v v + \partial_u u + 2) \partial_{vu^2} - \partial_{v^2 u^2} + (1+\partial_u u) \partial_u u \partial_{v^2} - \partial_v v (\partial_v v + 2 \partial_u u + 3) \partial_{u^2}. \end{align*} $$

Your question--how to derive the last two "pure" recurrence relations from the previous six--is equivalent to asking if the last two elements are in the (left) ideal generated by the first six elements. We can check this in Macaulay2:

loadPackage "Dmodules"
W = QQ[x,y,u,v,Dx,Dy,Du,Dv, WeylAlgebra => {x=>Dx, y=>Dy, u=>Du, v=>Dv}]
rP = 1+u^2+v^2-2*x*u-2*y*v
I = RatAnn rP -- the ideal of relations from the explicit formula
J = ideal(Du*Dv*(Dx*v-Dy*u), x*Dx*Du - Du*u*(Du+Dx), y*Dx*Du*Dv-Du*Dv*Dv*u-Dx*Du*Dv*v, y*Dv*Dy-Dv*v*(Dv+Dy), x*Dy*Du*Dv-Dv*v*Du^2-Du*u*Dy*Dv, (1-x^2)*Dx^2-x*y*Dx*Dy-x*Dx*Dv*v-2*x*Dx+u*y*Du*Dy+u*Du*(Du*u+Dv*v))
isSubset(J,I) -- True: first 6 proposed relations hold for explicit formula
isSubset(I,J) -- False: there are more than these 6
f = 2*x*(Du*u+Dv*v+2)*Du*Dv^2-Du^2*Dv^2+(1+Dv*v)*Dv*v*Du^2-Du*u*(Du*u+2*Dv*v+3)*Dv^2
f % J -- 0: 7th relation is implied by first 6!
f // gens J -- write 7th relation in terms of first 6; pages of horrendous mess

This counts as a proof, though we have to trust the computer for the hard part. If you can't prove the geometric series formula directly, you could instead prove the other generators of "I" work in your context and presumably simplify the mess the computer had to do to get the 7th relation from the first 6. This does not seem worthwhile, however.