Asymptotic Expansion of Bessel Function Integral I have an integral:
$$I(y) = \int_0^\infty \frac{xJ_1(yx)^2}{\sinh(x)^2}\ dx $$
and would like to asymptotically expand it as a series in $1/y$. Does anyone know how to do this? By numerically computing the integral it appears that
$$I(y) = \frac 12 - \frac 1 {\pi y}+ \frac {3\zeta(3)}{4y^3\pi^3} + O(y^{-5}) $$
but this is just (high precision) guesswork and I would like to understand the series analytically.
 A: Inserting the Mellin-Barnes representation for the square of the Bessel function (DLMF),
\begin{equation}
 J_{1}^2\left(xy\right)=\frac{1}{2\pi i}\int_{c-i\infty}^{c+i
\infty}\frac{\Gamma\left(-t\right)\Gamma\left(2t+3\right)}{\Gamma^2\left(t+2\right)\Gamma%
\left(t+3\right)}\left(\frac{xy}{2}\right)^{2t+2}\,dt
\end{equation} 
where $-3/2<\Re (c)<0$, and changing the order of integration, one obtains
\begin{equation}
 I(y)=\frac{1}{2\pi i}\int_{c-i\infty}^{c+i
\infty}\frac{\Gamma\left(-t\right)\Gamma\left(2t+3\right)}{\Gamma^2\left(t+2\right)\Gamma%
\left(t+3\right)}\left(\frac{y}{2}%
\right)^{2t+2}\,dt\int_0^\infty \frac{x^{2t+3}}{\sinh^2x}\,dx
\end{equation} 
From G. & R. (3.527.1)
\begin{equation}
\int_0^\infty \frac{x^{2t+3}}{\sinh^2x}\,dx=\frac{1}{2^{2t+2}}\Gamma\left( 2t+4 \right)\zeta\left( 2t+3 \right)
\end{equation} 
valid for $t>-1$. Thus we choose $-1<\Re(c)<0$
and thus
\begin{equation}
 I(y)=\frac{1}{2\pi i}\int_{c-i\infty}^{c+i
\infty}\frac{\Gamma\left(-t\right)\Gamma\left(2t+3\right)\Gamma\left( 2t+4 \right)\zeta\left( 2t+3 \right)}{\Gamma^2\left(t+2\right)\Gamma\left(t+3\right)}\left(\frac{y}{4}\right)^{2t+2}\,dt
\end{equation} 
To evaluate asymptotically this integral, one can close the contour by the large left half-circle. Poles are situated at $t=-1$ and $t=-\frac{2n+1}{2}$, with $n=1,2,3\ldots$. With the help of a CAS, the first corresponding residues are:
\begin{equation}
 R_{-1}=\frac{1}{2}\quad ;\quad R_{-3/2}=-\frac{1}{4\pi}\quad ;\quad R_{-5/2}=-\frac{3}{64\pi}\zeta'(-2)\quad ;\quad R_{-7/2}=\frac{15}{8192\pi}\zeta'(-4)
\end{equation} 
(General expressions can probably be found, if necessary). The derivative of the Riemann Zeta function at even integer values are involved and can be simply expressed. We obtain finally
\begin{equation}
 I(y)=\frac{1}{2}-\frac{1}{\pi}y^{-1}+\frac{3\zeta(3)}{4\pi^3}y^{-3}+\frac{45\zeta(5)}{32\pi^5}y^{-5}+O\left( y^{-7} \right)
\end{equation} 
