An infinite series involving Legendre polynomials For $x \in [-1,1]$ and $0 \le g < 1$, consider the convergent series
$$
H(x,g) = \sum_{k = 0}^\infty (2k+1) g^k P_k(x)^2
$$
where $P_k$ is the $k$-th Legendre polynomial. Then $H(1,g) = \frac{1+g}{(1-g)^2}$, and numerical evidence suggests that $H(0,g) \sim \frac{2}{\pi(1-g)}$ as $g \to 1$. 
Is this correct? Is this known? Is anything known about the function $H$?
 A: This is correct and known. In fact one has (formula 5.10.2.1 from Vol. II of Prudnikov-Brychkov-Marychev)
$$\mathcal{F}(t,y,z):=\sum_{k=0}^{\infty}t^k P_k(\cos y)P_k(\cos z)=\frac{4}{\pi(u_++u_-)}\mathbf{K}\left(\frac{u_+-u_-}{u_++u_-}\right),$$
where $u_{\pm}=\sqrt{1-2t\cos(y\pm z)+t^2}$ and $\displaystyle\mathbf{K}(k)=\int_0^{\pi/2}\frac{dx}{\sqrt{1-k^2\sin^2 x}}$ denotes the complete elliptic integral of the first kind.
The sum to be calculated is then given by
\begin{align}H(x,g)&=\left[\left(2t\frac{d}{dt}+1\right)\mathcal{F}(t,y,y)\right]_{t=g,\cos y=x}=\\
&=\frac4{\pi}\left(2g\frac{d}{dg}+1\right)\frac{1}{v(x,g)+1-g}\mathbf{K}\left(\frac{v(x,g)-1+g}{v(x,g)+1-g}\right),
\end{align}
where $v(x,g)=\sqrt{(1+g)^2-4gx^2}$. After differentiation we obtain somewhat cumbersome expression
\begin{align}
H(x,g)=&\frac{2(1+g)(1-g+v(x,g))}{\pi (1-g)v^2(x,g)}\mathbf{E}\left(\frac{v(x,g)-1+g}{v(x,g)+1-g}\right)+\\
+&\frac{(1+g)(1-g-v(x,g))}{\pi g(1-x^2)v(x,g)}\mathbf{K}\left(\frac{v(x,g)-1+g}{v(x,g)+1-g}\right),
\end{align}
where $\displaystyle\mathbf{E}(k)=\int_0^{\pi/2}\sqrt{1-k^2\sin^2 x}\,dx$ denotes complete elliptic integral of the second kind. 
Now the exact expression
$$H(0,g)=\frac{4\mathbf{E}(g)}{\pi(1-g^2)}-\frac{2\mathbf{K}(g)}{\pi},$$
and $\mathbf{E}(1)=1$ imply the asymptotics conjectured in the question.
A: There is a useful identity for products of Legendre polynomials (see DLMF (34.3.19)) which in the case of interest states
$P_k^2(x) = \sum\limits _{l=0}^\infty C_{kl} P_{2l}$ where $C_{kl}=\begin{pmatrix}k & k & 2l \\ 0 & 0 & 0 \end{pmatrix}^2$ and $\begin{pmatrix}l_1 & l_2 & l_3 \\ 0 & 0 & 0 \end{pmatrix}$ is the Wigner-$3j$ symbol. (We have taken the third parameter to be even since this symbol will vanish if the sum of parameters is odd.) Thus $H(x,t)=\sum_{k,l\geq 0}C_{kl} (2k+1)g^{k}P_{2l}(x)$.
We can make this more explicit by observing that (DLMF 34.3.5)
$$C_{kl}=\begin{pmatrix}k & k & 2l \\ 0 & 0 & 0 \end{pmatrix}^2
=(-1)^{k+l}\frac{(2l)!^2(2k-2l)!}{(2k+2l+1)!}\frac{(k+l)!^2}{k!^4\,2!}\\
\implies \quad
H(x,t)=\sum_{k,l\geq 0} (-1)^{k+l}\frac{(2l)!^2(2k-2l)!}{(2k+2l+1)!}\frac{(k+l)!^2}{k!^4\,l!^2}(2k+1)g^k P_{2l}.$$ Thus we have traded $P_k^2$'s for even Legendre polynomials and combinatorial coefficients. (I'll see if I can develop this answer further, particularly the cases $x=0,1$.)

