Let $P_n(x)$ be the $n$-th degree Legendre polynomial. Let $k$ be a nonnegative integer less than or equal to both $n,m$. How to prove that $$ \int_{-1}^1 (1-x^2)^k D^kP_n(x) D^kP_m(x)\,dx = \frac{2(n+k)!}{(n-k)!(2n+1)}\delta_{mn}, $$ where $D^k$ stands for derivative of order $k$, and $\delta_{mn}$ stands for Kronecker delta.
I think somehow we will use the orthogonality of $P_n, P_m$ and the norm of $P_n$, but I don't see how.
Let $I(k,n,m)$ be the considered integral.
(1) Since $\deg P_n=n$, we see immediately that $I(k,n,m)=0$ if $k>n$ or $k>m$.
(2) Since $I(k,n,m)=I(k,m,n)$, we see that the remaining case to be considered is $0\leq m\leq n$.
(3) Assume that $0\leq m\leq n$. Noting that $-1$ and $1$ are multiple zeros of $(x^2-1)^kD^kP_m$ of multiplicity $k$, we obtain by using k-successive integrations by parts $$ I(k,n,m)=\int_{-1}^1P_n(x)D^k\left((x^2-1)^kD^{k}P_m(x)\right)dx $$
As a conclusion, the only cases where $I(k,n,m)$ is different from $0$ is when $n=m$ and $0\leq k\leq n$, and in this case $$ I(k,n,n)=\frac{ (n+k)!}{(n-k)!}\frac{2}{2n+1}. $$