Series of products of Bessel functions of the first kind

Question

I am a physicist and encountered the following series in my research: \begin{equation} \sum_{n=-\infty}^\infty n^p J_n^2(z), \end{equation} where $J_n(z)$ is a Bessel function of the first kind and $p$ a nonnegative integer. I would like to know if there exists a general expression of this series in terms of a Taylor series in $z$? From the symmetry properties of the Bessel function I already know that the sum vanishes for odd $p$ and that it has to be an even function of $z$. Hence, we can restrict ourselves to the sum \begin{equation} \sum_{n=-\infty}^\infty n^{2m} J_n^2(z), \end{equation} where $m$ is a nonnegative integer. By using the generating function of Bessel functions, I have written a program in Mathematica that essentially finds the series expansion for a given $m$. By looking at a few cases, I have established the following hypothesis: \begin{equation} \sum_{n=-\infty}^\infty n^{2m} J_n^2(z) = \sum_{l=0}^m c_{m,l} \frac{z^{2l}}{(2l)!}. \end{equation} So essentially I would like to know if there exists a general expression for the coefficients $c_{m,l}$. For example, by putting $z=0$ one finds $c_{m,0}=\delta_{m,0}$ and by taking derivatives of the series and setting $z=0$, one finds $c_{m,1}=1$ for $m\geq1$ and $c_{m,2}=3(4^{m-1}-1)$ for $m\geq2$ and zero otherwise. One can keep on taking derivatives, i.e., the $2l$th derivative at the origin will be nonzero only for $m\geq l$, but this would require a general expression for the $2l$th derivative of the series at the origin which seems kind of cumbersome.

I would also be interested to know if there exist similar expressions for the more general series, \begin{equation} \sum_{n=-\infty}^\infty n^p J_n(z) J_{n+m}(z), \end{equation} for example \begin{equation} \sum_{n=-\infty}^\infty n J_n(z) J_{n+m}(z) = \frac{z}{2} \delta_{m,\pm1}. \end{equation}

Answer 1

A possible approach is to use the generating function $$ \sum_{n\in\mathbb{Z}}J_n(z)J_{n+m}(z)t^{2n+m}=I_m\big(z(t-1/t)\big)\qquad(z,t\in\mathbb{C}) $$ which can be deduced from known generating functions $$ \sum_{n\in\mathbb{Z}}J_n(z)t^n=e^{\frac z2\left(t-\frac 1t\right)}, \quad \sum_{n\in\mathbb{Z}}I_n(z)t^n=e^{\frac z2\left(t+\frac 1t\right)}. $$

In detail, Cauchy's integral formula yields $$J_n(z)=\frac1{2\pi i}\oint e^{\frac z2\left(w-\frac1w\right)}\frac{dw}{w^{n+1}}\implies J_{n}(z)t^n=\frac1{2\pi i}\oint e^{\frac z2\left(tw-\frac1{tw}\right)}\frac{dw}{w^{n+1}},$$ so that \begin{align}\sum_{n\in\mathbb{Z}}J_n(z)J_{n+m}(z)t^{2n+m}&=\frac1{2\pi i}\sum_{n\in\mathbb{Z}}J_n(z)t^n\oint e^{\frac z2\left(tw-\frac1{tw}\right)}\frac{dw}{w^{n+m+1}}\\&=\frac1{2\pi i}\oint e^{\frac z2\left(tw-\frac1{tw}\right)}e^{\frac z2\left(\frac tw-\frac wt\right)}\frac{dw}{w^{m+1}},\end{align} which equals $I_m\big(z(t-1/t)\big)$, again by Cauchy.

Assume $m\geqslant 0$ (the other case is treated similarly). Replace $t$ by $e^{t/2}$ to get $$ \sum_{n\in\mathbb{Z}}J_n(z)J_{n+m}(z)e^{nt}=e^{-mt/2}I_m\big(2z\sinh(t/2)\big). $$ Now take $d^p/dt^p$ and put $t=0$. Using the power series for $I_m$, we get $$ \sum_{n\in\mathbb{Z}}n^p J_n(z)J_{n+m}(z)=\sum_{\substack{n\ \geqslant\ 0\\2n+m\leqslant p}}\frac{(z/2)^{2n+m}}{n!(n+m)!}\sum_{k=0}^{2n+m}(-1)^k\binom{2n+m}{k}(n-k)^p. $$ I don't know whether the inner sum can be simplified (or expressed using known/named numbers).