Sum of squared Negative Binomial probability masses Let $(p_k)_{k=0, \dots, \infty}$ denote the probability masses of a Negative Binomial distribution with parameters $r>0$ and $p\in]0,1[$. I'm looking for the sum of their squares,
$$\sum_{k=0}^\infty p_k^2$$
as a function of $r$ and $p$ (other parameterizations are also fine). In other words I am interested in (the exponential of) the second-order Renyi entropy of a Negative Binomial distribution. 
Background:


*

*I'd like to use this to evaluate Brier and spherical scores.

*Czado et al. (2009) write that this expression can be evaluated analytically, but don't give further information, and I'm kind of stuck.

*No, this is not homework, although I imagine it could be ;-) Any hints or pointers to literature would be almost as much appreciated as a full solution.

*Here is the analogous question for the Poisson distribution.
 A: Take two independent Poisson random variables $X$ and $X'$,  with means $\lambda$ and $\lambda'$. 
The formula answering your question in the Poisson case is a particular case of the identity $$\Pr(X =  X' \mid \lambda, \lambda') = \exp(-\lambda)\exp(-\lambda')I_0(2\sqrt{\lambda}{\sqrt{\lambda'}}).$$
> lambda <- 1; lambdaa <- 2
> sum(dpois(0:100,lambda)*dpois(0:100,lambdaa))
[1] 0.2117121
> gsl::bessel_I0(2*sqrt(lambda*lambdaa)) * exp(-lambda-lambdaa)
[1] 0.2117121

Your problem is the same as calculating $\int\int\Pr(X= X' \mid \lambda, \lambda') f_{a,b}(\lambda)f_{a,b}(\lambda')d\lambda d\lambda'$ where $f_{a,b}$ is the $\Gamma(a,b)$ pdf, setting $r=a$ and $p=\frac{b}{1+b}$ with your notations, because of the link between the negative binomial distribution and the Poisson-Gamma distribution.
Let's start by $$\int \exp(-\lambda)I_0(2\sqrt{\lambda}{\sqrt{\lambda'}})f_{a,b}(\lambda)d\lambda = \frac{b^a}{\Gamma(a)}\int \lambda^{a-1}\exp\bigl(-(b+1)\lambda\bigr)I_0(2\sqrt{\lambda}{\sqrt{\lambda'}})d\lambda.$$
According to Mathematica this is equal to 
$$
{\left(\frac{b}{1+b}\right)}^a {}_1\!F_1\left(a, 1, \frac{\lambda'}{b+1}\right)
$$
where ${}_1\!F_1$ is the Kummer hypergeometric function. 
Now we can even get something for 
$$
\begin{multline}
\int  \exp(-\lambda'){}_1\!F_1\left(a, 1, \frac{\lambda'}{b+1}\right)f_{a',b'}(\lambda')d\lambda' \\
=  \frac{{b'}^{a'}}{\Gamma(a')}\int {\lambda'}^{a'-1}\exp\bigl(-(b'+1)\lambda'\bigr) {{}_1\!F_1}\left(a, 1, \frac{\lambda'}{b+1}\right)d\lambda'.
\end{multline}
$$
Indeed, Mathematica gives 
$$
{\left(\frac{b'}{1+b'}\right)}^{a'}
{}_2\!F_1\left(a, a', 1, \frac{1}{(b+1)(b'+1)}\right)
$$
where ${}_2\!F_1$ is the Gauss hypergeometric function.
The final result is beautiful:
$$
{\left(\frac{b}{1+b}\right)}^{a}{\left(\frac{b'}{1+b'}\right)}^{a'}{}_2\!F_1\left(a, a', 1,  \frac{1}{(b+1)(b'+1)}\right),
$$
and even a bit more beautiful with your notations:
$$
p^{a}{p'}^{a'}{}_2\!F_1\left(a, a', 1,  (1-p)(1-p')\right)
$$
Check:
> a <- 2; A <- 3; b <- 5; B <- 8
> (b/(1+b))^a*(B/(1+B))^A*gsl::hyperg_2F1(a,A,1,1/(b+1)/(B+1))
[1] 0.5450618
> sum(dnbinom(0:100, a, b/(1+b))*dnbinom(0:100, A, B/(1+B)))
[1] 0.5450618

In the special case you are interested in, the sum of squares is
$$
p^{2a}{}_2\!F_1\left(a, a, 1, {(1-p)}^2\right),
$$
and the second-order Renyi entropy is 
$$
-2a \log p - \log {}_2\!F_1\left(a, a, 1, {(1-p)}^2\right).
$$
