We know,from generating function of Legendre polynomials,
$(1-2xz+z^2)^\frac{-1}{2}=\sum z^nP_n(x)$
Squaring both sides and writing the sum by separating the cross terms and "self-quadrature" terms,we get
$(1-2xz+z^2)^{-1}=\sum z^{2n}P_n^2(x)+2\sum z^{m+n}P_m(x)P_n(x)$
Integrating both sides between -1 and +1,we have,
$\int_{-1}^1\sum z^{2n}P_n^2(x)dx+\int_{-1}^1 2\sum z^{m+n}P_m(x)P_n(x)dx=\int_{-1}^{1}(1-2xz+z^2)^{-1}dx$
Now,
$\int_{-1}^{1}P_n(x)P_m(x)dx=0$ due to orthogonality of Legendre polynomials and the 2nd integral on LHS of the above equation is wrt x,the terms containing m and n can be taken out,and therefore the 2nd integral vanishes.
$\Rightarrow \int_{-1}^1\sum z^{2n}P_n^2(x)dx+0=\int_{-1}^{1}(1-2xz+z^2)^{-1}dx$
The integral on RHS,i.e.,
$
{\displaystyle\int}\dfrac{1}{-2zx+z^2+1}\,\mathrm{d}x $ can be solved as--
Substitute $u=-2zx+z^2+1$
$\dfrac{\mathrm{d}u}{\mathrm{d}x} = -2z$
$\Rightarrow \mathrm{d}x=-\dfrac{1}{2z}\,\mathrm{d}u$
The integral reduces to
$=-\class{steps-node}{\cssId{steps-node-1}{\dfrac{1}{2z}}}{\displaystyle\int}\dfrac{1}{u}\,\mathrm{d}u$
which is a standard integral.
Therefore the integral becomes
$=-\dfrac{\ln\left(u\right)}{2z}$
$=-\dfrac{\ln\left(-2zx+z^2+1\right)}{2z}$
$\Rightarrow \sum z^{2n}\int_{-1}^1P_n^2(x)dx=\frac{-1}{2z}[\ln(1-2xz+z^2)]_{-1}^1$
RHS=
$\frac{-1}{2z}\ln\frac{1-2z+z^2}{1+2z+z^2}=\frac{-1}{2z}\ln(\frac{1-z}{1+z})^2=\frac{1}{z}[\ln(1+z)-\ln(1-z)]$
Using Maclaurin Series of $\ln(1+z)$ and $\ln(1-z)$ and expanding,
RHS=
$\frac{1}{z}[(z-\frac{z^2}{2}+\frac{z^3}{3}-\frac{z^4}{4}+...)-(-z-\frac{z^2}{2}-\frac{z^3}{3}-\frac{z^4}{4}-...)]$
Observe that the terms containing even powers of x cancel and the odd terms gets doubled(due to presence of them in both the series indicated above).Taking the 1/z term inside the square brackets,
$\Rightarrow \int_{-1}^1\sum z^{2n}P_n^2(x)dx=2[1+\frac{z^2}{3}+\frac{z^4}{5}+..\frac{z^{2n}}{2n+1}.]$
Finally,equating the coefficients of $z^{2n}$ on both sides we have,
$\int_{-1}^{1} (P_n(x))^2 dx = \frac{2}{2n+1}$,which is the required relation.
Best Answer
After some thoughts, I couldn't find any reasonable way to do it directly, but at least now I understand how to fill the gap in the derivation of the genfunc of Legendre polynomials. First, let's expand the genfunc.
$$ G(t, x) = \frac{1}{\sqrt{1-2xt+t^2}} = \sum\limits_k \frac{t^k(2x-t)^k}{4^k}\binom{2k}{k}. $$ If we want to extract the coefficient near $t^n$, we get $$ [t^n] G(t, x) = \frac{1}{2^n}\sum\limits_k \binom{2k}{k} \binom{k}{n-k} (-1)^{n-k} x^{2k-n}, $$ where $k$ goes from $\lceil n/2 \rceil$ to $n$. Changing $k \to n-k$, we get $$ [t^n] G(t, x) = \frac{1}{2^n} \sum\limits_k (-1)^k \binom{2n-2k}{n-k} \binom{n-k}{k} x^{n-2k}. $$ On the other hand, $$ \binom{2n-2k}{n-k}\binom{n-k}{k}=\binom{2n-2k}{n-k,k,n-2k}=\binom{2n-2k}{n}\binom{n}{k}, $$ which gives one of the standard expressions $$ [t^n] G(t, x) =\boxed{ \frac{1}{2^n} \sum\limits_k (-1)^k \binom{n}{k}\binom{2n-2k}{n} x^{n-2k}} $$
Now, we need to connect the dots. To make it relevant to $\binom{n}{k}^2$, we define Legendre polynomials as $$ P_n(x) = \frac{1}{2^n} \sum\limits_{k=0}^n \binom{n}{k}^2(x-1)^{n-k} (x+1)^k. $$ It's easy to check that it compacts into $$ P_n(x) = [t^n] \frac{(t+x-1)^n(t+x+1)^n}{2^n} = [t^n] \frac{((t+x)^2-1)^n}{2^n}. $$
We can now expand it back to get the same expression: $$ [t^n]\frac{((x+t)^2-1)^n}{2^n} = [t^n]\frac{1}{2^n}\sum\limits_k (-1)^k \binom{n}{k}(x+t)^{2n-2k} =\boxed{ \frac{1}{2^n} \sum\limits_k (-1)^k \binom{n}{k}\binom{2n-2k}{n} x^{n-2k}} $$ I still don't like this approach, as I needed to know $G(t, x)$ in advance to justify it, and it's quite far from being direct in terms of the genfunc for $\binom{n}{k}^2$, but I suppose it's still nice to have proofs that are directly based on binomial identities, rather than physics...