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).
$$
