Integral of Bessel function, power, and hyperbolic cosecant I have the integral
$$
I(a)=\int\limits_0^\infty dk \ \frac{k^3J_1(ak)}{\sinh(k )}
$$
Where $J_1$ is a Bessel function of the first kind. By plotting the integrand, it seems $I(a)$ exists as the integrand is finite and goes to zero for large $k$. I wonder if $I(a)$ may be expressed in `common' special functions? Mathematica does not evaluate it, and I have been unable to find it in Gradshtein.
I have been able to find an approximation to $I(a)$ as $a \to 0$ by replacing $J_1$ with its asymptotic form for small argument. Unfortunately, in my application the interesting case is $a\to \infty$. If $I(a)$ cannot be simplified, how can we develop an approximation to it for large $a$?
Background: The integral arises in a Laplacian boundary value problem for the induced charge on a grounded plate.
 A: By $(10.22.49)$ and $(15.4.19)$ in the DLMF, we find
\begin{align*}
& \int_0^{ + \infty } {\frac{{t^3 }}{{\sinh t}}J_1 (at)dt}  = 2\int_0^{ + \infty } {t^3 e^{ - t} \frac{1}{{1 - e^{ - 2t} }}J_1 (at)dt} 
\\ & = 2\int_0^{ + \infty } {t^3 \sum\limits_{n = 0}^\infty  {e^{ - (2n + 1)t} } J_1 (at)dt} \\ & = 2\sum\limits_{n = 0}^\infty  {\int_0^{ + \infty } {t^3 } e^{ - (2n + 1)t} J_1 (at)dt} 
\\ & = 24a\sum\limits_{n = 0}^\infty  {\frac{1}{{(2n + 1)^5 }}F\!\left( {\frac{5}{2},3;2; - \frac{{a^2 }}{{(2n + 1)^2 }}} \right)} 
\\ & = 24a\sum\limits_{n = 0}^\infty  {\left( {(2n + 1)^2  - \frac{{a^2 }}{4}} \right)\frac{1}{{((2n + 1)^2  + a^2 )^{7/2} }}} \\ &=
24a\sum\limits_{n = 0}^\infty  {\frac{1}{{((2n + 1)^2  + a^2 )^{5/2} }}}  - 30a^3 \sum\limits_{n = 0}^\infty  {\frac{1}{{((2n + 1)^2  + a^2 )^{7/2} }}} 
\end{align*}
provided $\Re a>0$. Here $F$ stands for the hypergeometric function.
Addendum. I shall give the asymptotics for $a\to +\infty$. Let us introduce the generalised Mathieu series via
$$
S_{\mu ,\gamma } (a;\lambda ) = \sum\limits_{n = 1}^\infty  {\frac{{n^\gamma  }}{{(n^\lambda   + a^\lambda  )^\mu  }}} \quad \quad (\mu  > 0,\quad \lambda  > 0,\quad \lambda \mu  - \gamma  > 1).
$$
With this notation
$$
\sum\limits_{n = 0}^\infty  {\frac{1}{{((2n + 1)^2  + a^2 )^\mu  }}}  = S_{\mu ,0} (a;2) - 4^{ - \mu } S_{\mu ,0} (a/2;2) 
$$
provided $\mu>\frac{1}{2}$. The precise asymptotics of the generalised Mathieu series was derived in this paper. In particular,
$$
S_{\mu ,0} (a;2) = \frac{{\sqrt \pi  \Gamma \left( {\mu  - \frac{1}{2}} \right)}}{{2\Gamma (\mu )a^{2\mu  - 1} }} - \frac{1}{{2a^{2\mu } }} + \frac{{\pi ^\mu  }}{{\Gamma (\mu )a^\mu  }}e^{ - 2\pi a} (1 + o(1))
$$
as $a\to +\infty$. Consequently,
$$
\sum\limits_{n = 0}^\infty  {\frac{1}{{((2n + 1)^2  + a^2 )^\mu  }}}  = \frac{{\sqrt \pi  \Gamma \left( {\mu  - \frac{1}{2}} \right)}}{{4\Gamma (\mu )a^{2\mu  - 1} }} - \frac{{\pi ^\mu  }}{{\Gamma (\mu )(2a)^\mu  }}e^{ - \pi a} (1 + o(1))
$$
as $a\to +\infty$. Combining this with the exact series representation above, we find
$$
\int_0^{ + \infty } {\frac{{t^3 }}{{\sinh t}}J_1 (at)dt}  =
\frac{{\sqrt 2 \,\pi ^3 }}{{a^{1/2} }}e^{ - \pi a} (1 + o(1))
$$
as $a\to +\infty$. More precise asymptotics can be derived by using more terms in the asymptotic expansion of the generalised Mathieu series which can be found in the paper cited above.
A: Here is a more self contained approach.
We rewrite the integral in terms of Hankel functions $H_1(x)$ of the first kind since they have the nice property of decaying exponentially in the upper half of the complex plane. By Lemma 1 proven below, we get:
$$
I(a)=\frac12\Re\int_{\mathbb{R}}\frac{k^3 H_{1}(ak)}{\sinh(k)}dk =\frac12\Re \int_{\mathbb{R}} g(x)
$$
Now due to the above mentioned decay, we close the contour of integration in the upper half plane and get by residue theorem
$$
I(a)=\Re(\pi i \sum_{k>0} \text{res}(g(z),z=i \pi k) )
$$
it is worth noting that since $x^3H_1(x)$ has no poles in the complex plane we only sum over the zeros of $\sinh(x)$.
In the asymptotic limit of large $a$ only the residue which is closest to the imaginary axis will contribute to lowest order (additional terms are suppressed by $O(e^{- k \pi a })$)
$$
I(a)\sim_{a\rightarrow \infty}\Re(\pi i  \text{res}(g(z),z=i \pi ) )=\Re(-\pi^4 H_{1}(i a \pi))
$$
by the standard asymptotic expansions for the Hankel function this is equivalent to

$$
I(a)\sim_{a\rightarrow \infty} \frac{\pi^3\sqrt{2}e^{-\pi a}}{\sqrt a}\left(1+\frac{3}{ 8 a \pi}+o(a^{-1})\right)
$$

which coincides with the other answer to leading order. I expect this technique to work for every $k^l H_{\nu}(a k)$ such that $l,\nu  \in 2\mathbb{N}+1$ or $l,\nu  \in 2\mathbb{N}$.
Lemma 1:
$\Re(H_{\nu}(x))$ is even for $\nu \in 2\mathbb{N}+1$ and $x \in \mathbb{R}$
Proof:
$$
\Re(H_{\nu}(x)) = \begin{cases}
  -J_{\nu}(x)  & x <0 \\
  J_{\nu}(x) & x>0
\end{cases}
$$
the first case follows from the connection formula 10.4.7 and $J_{\nu}(-x)=e^{i\pi \nu}J_{\nu}(x)$.
the second case is due to the very definiton of the Hankel function for $x>0$. Since $J_{\nu}(x)$ is odd, we are done. QED
A: Well, we are trying to find the following integral:
$$\mathcal{I}_\text{n}\left(\alpha\right):=\int_0^\infty\frac{x^\text{n}\mathscr{J}_1\left(\alpha x\right)}{\sinh\left(x\right)}\space\text{d}x=\int_0^\infty x^\text{n}\mathscr{J}_1\left(\alpha x\right)\text{csch}\left(x\right)\space\text{d}x\tag1$$
Using the 'evaluating integrals over the positive real axis' property of the Laplace transform, we can write:
$$\mathcal{I}_\text{n}\left(\alpha\right)=\int_0^\infty\mathscr{L}_x\left[x^\text{n}\mathscr{J}_1\left(\alpha x\right)\right]_{\left(\sigma\right)}\cdot\mathscr{L}_x^{-1}\left[\text{csch}\left(x\right)\right]_{\left(\sigma\right)}\space\text{d}\sigma\tag2$$
Using known results, we can write:
$$\mathcal{I}_\text{n}\left(\alpha\right)=\left(-1\right)^\text{n}\cdot\int_0^\infty\frac{\partial^\text{n}}{\partial\sigma^\text{n}}\left(\frac{\alpha}{\alpha^2+\sigma\left(\sigma+\sqrt{\alpha^2+\sigma^2}\right)}\right)\cdot\left(2\sum_{\text{k}\space\ge\space0}\delta\left(\sigma-2\text{k}-1\right)\right)\space\text{d}\sigma\tag3$$
We can rewrite this a bit:
$$\mathcal{I}_\text{n}\left(\alpha\right)=2\left(-1\right)^\text{n}\sum_{\text{k}\space\ge\space0}\int_0^\infty\frac{\partial^\text{n}}{\partial\sigma^\text{n}}\left(\frac{\alpha}{\alpha^2+\sigma\left(\sigma+\sqrt{\alpha^2+\sigma^2}\right)}\right)\cdot\delta\left(\sigma-2\text{k}-1\right)\space\text{d}\sigma\tag4$$
Using the property:
$$\int_0^\infty\text{y}\left(x\right)\delta\left(x-\text{p}\right)\space\text{d}x=\text{y}\left(\text{n}\right)\theta\left(\text{n}\right)\tag5$$
We can write:
$$\mathcal{I}_\text{n}\left(\alpha\right)=2\left(-1\right)^\text{n}\sum_{\text{k}\space\ge\space0}\left\{\theta\left(1+2\text{k}\right)\cdot\left.\frac{\partial^\text{n}}{\partial\sigma^\text{n}}\left(\frac{\alpha}{\alpha^2+\sigma\left(\sigma+\sqrt{\alpha^2+\sigma^2}\right)}\right)\right|_{\space\sigma\space=\space1+2\text{k}}\right\}\tag6$$
Now, using the fact that when $\text{k}\in\mathbb{N}$ we get $\theta\left(1+2\text{k}\right)=1$. So we can conclude:
$$\mathcal{I}_\text{n}\left(\alpha\right)=2\left(-1\right)^\text{n}\sum_{\text{k}\space\ge\space0}\left\{\left.\frac{\partial^\text{n}}{\partial\sigma^\text{n}}\left(\frac{\alpha}{\alpha^2+\sigma\left(\sigma+\sqrt{\alpha^2+\sigma^2}\right)}\right)\right|_{\space\sigma\space=\space1+2\text{k}}\right\}\tag7$$
