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.

• Are you sure there is not a typo in your approximation? I can fairly quickly get $I(y)\sim1/2-1/y.$ – skbmoore Nov 14 '18 at 2:43
• It’s possible, a colleague did the numerics and he may have made a typo in the notes. Did you get the 1/y term analytically? – djbinder Nov 14 '18 at 2:51
• Even my comment had a mistake. the $\pi$ is in the denominator, as in Paul Enta's answer. – skbmoore Nov 14 '18 at 16:53
• Thanks for pointing this out, I've fixed the location of the $\pi$s – djbinder Nov 14 '18 at 23:11

Inserting the Mellin-Barnes representation for the square of the Bessel function (DLMF), $$$$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$$$$ where $$-3/2<\Re (c)<0$$, and changing the order of integration, one obtains $$$$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$$$$ From G. & R. (3.527.1) $$$$\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)$$$$ valid for $$t>-1$$. Thus we choose $$-1<\Re(c)<0$$ and thus $$$$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$$$$ 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: $$$$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)$$$$ (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 $$$$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)$$$$