New identity for sums of Bessel functions? I've come across the following proposed identity:
$$
2\sum_{q=1}^{\infty}\frac{q-1}{2q-1}\left(-\right)^{q}I_{q}^{2}\left(\alpha^{2}\right)+\alpha^{2}\sum_{q=-\infty}^{\infty}\frac{\left(-\right)^{q}}{2q+1}I_{q}\left(\alpha^{2}\right)I_{\left|q\right|+1}\left(\alpha^{2}\right)=\frac{1}{4}\sqrt{\frac{\left(1+e^{-2\alpha^{2}}\right)^{3}\left(e^{4\alpha^{2}}-1\right)}{1-e^{-2\alpha^{2}}}}-I_{0}^{2}\left(\alpha^{2}\right)
$$
In alternative form:
$$
\sum_{q=1}^{\infty}\frac{\left(-\right)^{q}I_{q}\left(\alpha^{2}\right)}{2q-1}\left[\left(q-1\right)I_{q}\left(\alpha^{2}\right)-\frac{\alpha^{2}I_{q+1}\left(\alpha^{2}\right)}{2q+1}\right]=\frac{1}{8}\sqrt{\frac{\left(1+e^{-2\alpha^{2}}\right)^{3}\left(e^{4\alpha^{2}}-1\right)}{1-e^{-2\alpha^{2}}}}-\frac{1}{2}I_{0}\left(\alpha^{2}\right)\left[I_{0}\left(\alpha^{2}\right)+\alpha^{2}I_{1}\left(\alpha^{2}\right)\right]
$$
Above, $I_q$ is the modified Bessel function of the first kind. I have verified this numerically for real $0\leq\alpha\leq7$. Could anyone provide some pointers as to how to derive this analytically?
 A: To prove this identity one needs the following formulas http://dlmf.nist.gov/10.29.E1 http://dlmf.nist.gov/10.32.E15
$$I_{q-1}(x)+I_{q+1}(x)=2I_q'(x),$$ $$I_{q-1}(x)-I_{q+1}(x)=2qI_q(x)/x,$$
$$\mathop{I_{\mu}}\nolimits\!\left(z\right)\mathop{I_{\nu}}\nolimits\!\left(%
z\right)=\frac{2}{\pi}\int_{0}^{\frac{1}{2}\pi}\mathop{I_{\mu+\nu}}\nolimits%
\!\left(2z\mathop{\cos}\nolimits\theta\right)\mathop{\cos}\nolimits((\mu-%
\nu)\theta)d\theta,$$
$$I_{1/2}(x)=\sqrt{\frac{2}{\pi x}}\sinh x,$$ $$I_{-1/2}(x)=\sqrt{\frac{2}{\pi x}}\cosh x,$$
and the addition theorem (equation 11.41.(6) on page 365, Treatise on the theory of Bessel funcions, by G.N. Watson)
$$
I_0(2a\sin{\theta}/{2})=\sum_{k=-\infty}^\infty (-1)^kI_k^2(a)e^{ik\theta}.
$$
First, by elementary simplifications one can write the initial identity in the form

$${a \sum _{n=-\infty}^{\infty} \frac{(-1)^n I_{n+1}(a) I_n(a)}{2 n+1}}={\sinh ^2(a)},$$

or equivalently
$${\sum _{n=-\infty}^{\infty} \frac{(-1)^n (n+1) I_{n+1}^2(a)}{2 n+1}}+a{\sum _{n=-\infty}^{\infty} \frac{(-1)^n I_{n+1}(a)I_{n+1}'(a)}{2 n+1}}={\sinh ^2(a)}.$$
By integrating the addition theorem with the help of the above identities and the fact $\int_0^{2\pi}e^{i\theta(k+1/2)}d\theta=4i/(2k+1),~k\in \mathbb{Z}$ one can show that
$$
{\sum _{n=-\infty}^{\infty} \frac{(-1)^n I_{n+1}^2(a)}{2 n+1}}=\int_0^{\pi/2}I_0(2a\sin\theta)\sin\theta d\theta,
$$
$$
{\sum _{n=-\infty}^{\infty} \frac{(-1)^n I_{n+1}(a)I_{n+1}'(a)}{2 n+1}}=\int_0^{\pi/2}I_1(2a\sin\theta)\sin^2\theta d\theta,
$$
$$
2{\sum _{n=-\infty}^{\infty} \frac{(-1)^n (n+1)I_{n+1}^2(a)}{2 n+1}}=\int_0^{\pi/2}I_0(2a\sin\theta)\sin\theta d\theta-1.
$$
The identity under consideration becomes 
$$
\frac{1}{2}\int_0^{\pi/2}I_0(2a\sin\theta)\sin\theta d\theta+a \int_0^{\pi/2}I_1(2a\sin\theta)\sin^2\theta d\theta=\sinh^2a+\frac{1}{2}=\frac{1}{2}\cosh 2a,
$$
or
$$
\frac{d}{da}\left( a\int_0^{\pi/2}I_0(a\sin\theta)\sin\theta d\theta\right)=\cosh a.
$$
This is a consequence of 
$$
\int_0^{\pi/2}I_0(a\sin\theta)\sin\theta d\theta=\frac{\sinh a}{a}
$$
which can be derived from the formulas given at the beginning of this post.
