The function has a curious property: for any linear function , and any point , the integral evaluates to . This is easy to check using the fact that odd powers of integrate to zero:

More generally, for any integer there exists a unique symmetric polynomial that has degree in and separately and satisfies for all polynomials of degree at most . For example, (obviously) and

The formula is not really intuitive, and a 3d plot would not help the matter much. To visualize , I plotted , , and below (green, red, blue respectively).

For and large , the function approaches the Dirac delta at , although the convergence is slow, especially when is close to . I don’t think there is anything good to be said about the case .

The existence and uniqueness of are a consequence of the Riesz representation of linear functionals on an inner product space. Indeed, polynomials of degree at most form such a space with inner product , and the functional is linear for any fixed . Hence, this functional can be written as for some . The function is a *reproducing kernel* for this space. Its symmetry is not immediately obvious.

The Legendre polynomials are an orthogonal basis of ; more precisely, form an orthonormal basis. It’s a general fact about reproducing kernels that

(which, incidentally, proves the symmetry ). Indeed, taking this sum as the definition of and writing , we find

This is the Sage code used for the above plots.

```
n = 20
k = sum([(j+1/2)*legendre_P(j,x)*legendre_P(j,y) for j in range(0,n+1)])
plot(k(x,y=-3/4),(x,-1,1),color='green') + plot(k(x,y=0),(x,-1,1),color='red') + plot(k(x,y=1/2),(x,-1,1),color='blue')
```

Higher degrees cause some numerical issues…

Post motivated by Math.SE question

Shouldn’t the multiplying factor for \(\tilde{P}_j\) be \((j+1/2)^{1/2}\)?

You are right; I corrected this. It was a typo; the code with (j+1/2)*legendre_P(j,x)*legendre_P(j,y) was based on sqrt(j+1/2) being the factor, which becomes (j+1/2) when applied to both polynomials.