How to find an asymptotic expansion, for $\epsilon$ near $0$, of the following integral $$ I(\epsilon):=\int_{0}^{+\infty}\frac 1{\sinh^2 (\epsilon \sqrt{x^2+1}) } {\rm d}x. $$ As $\epsilon \rightarrow 0$, I have easily obtained $$ I(\epsilon) \sim \frac{\pi}{2\:\epsilon^2}. $$
Some clear steps leading to an extended expansion would be appreciated.
My calculation shows that
$$I(\epsilon) := \int_{0}^{\infty} \frac{dx}{\sinh^{2}(\epsilon\sqrt{x^{2}+1})} = \frac{\pi}{2\epsilon^{2}} - \frac{1}{\epsilon} + \pi \epsilon \sum_{n=1}^{\infty} \frac{1}{(\pi^{2}n^{2} + \epsilon^{2})^{3/2}} \tag{1}. $$
In particular, if we expand the infinite sum on the RHS, we get
$$ I(\epsilon) = \frac{\pi}{2\epsilon^{2}} - \frac{1}{\epsilon} + \sum_{n=0}^{\infty} \binom{-3/2}{n} \frac{\zeta(2n+3)}{\pi^{2n+2}} \epsilon^{2n+1}. $$
Indeed, from the following expansion
$$ \frac{1}{\sinh^{2}z} = \sum_{n=-\infty}^{\infty} \frac{1}{(z-i\pi n)^{2}} = \frac{1}{z^{2}} + 2 \sum_{n=1}^{\infty} \frac{z^{2} - \pi^{2}n^{2}}{(z^{2} + \pi^{2} n^{2})^{2}}, $$
it follows from term-wise integration that
\begin{align*} \int_{0}^{R} \frac{dx}{\sinh^{2}(\epsilon\sqrt{x^{2}+1})} = \frac{\arctan R}{\epsilon^{2}} + \sum_{n=1}^{\infty} &\Bigg( \frac{2\epsilon \arctan \left( R\epsilon \big/ \sqrt{\pi^{2}n^{2} + \epsilon^{2}} \right)}{(\pi^{2}n^{2} + \epsilon^{2})^{3/2}} \\ &\quad + \frac{2}{R}\frac{1}{\pi^{2}n^{2}+\epsilon^{2}} \\ &\quad - \frac{2(R^{2}+1)}{R(\pi^{2}n^{2} + (R^{2}+1)\epsilon^{2})} \Bigg). \tag{2} \end{align*}
Here, term-wise integration is possible from Fubini's theorem together with the following estimate:
$$ \int_{0}^{R} \left| \frac{\epsilon^{2}(x^{2}+1) - \pi^{2}n^{2}}{(\epsilon^{2}(x^{2}+1) + \pi^{2} n^{2})^{2}} \right| = \frac{\arctan \left( R\epsilon \big/ \sqrt{\pi^{2}n^{2} + \epsilon^{2}} \right)}{\epsilon\sqrt{\pi^{2}n^{2}+\epsilon^{2}}} \lesssim_{\epsilon, R} \frac{1}{n^{2}}. $$
(Notice that the arctan term also contributes to order $n^{-1}$. It means that this argument fails if we consider $R = \infty$. This is why we consider proper integral first.)
Finally, taking $R \to \infty$ to (2) yields (1). When doing this, the only non-trivial calculation is to check that
$$ \lim_{R\to\infty} \sum_{n=1}^{\infty} \frac{2(R^{2}+1)}{R(\pi^{2}n^{2} + (R^{2}+1)\epsilon^{2})} = \frac{1}{\epsilon}. $$
But this follows from the squeezing lemma combined with the following inequality
$$ C(1, R) \leq \sum_{n=1}^{\infty} \frac{2(R^{2}+1)}{R(\pi^{2}n^{2} + (R^{2}+1)\epsilon^{2})} \leq C(0, R), $$
where
\begin{align*} C(a, R) &:= \int_{a}^{\infty} \frac{2(R^{2}+1)}{R(\pi^{2}x^{2} + (R^{2}+1)\epsilon^{2})} \, dx \\ &= \frac{\sqrt{R^{2}+1}}{R\epsilon} \frac{\arctan\left( \epsilon\sqrt{R^{2}+1} \big/ a\pi \right)}{\pi/2}. \end{align*}
What about considering more terms of the Weierstrass product of $\sinh$?
The first term of the asymptotics is easily recovered from: $$ I = \int_{0}^{\pi/2}\frac{d\theta}{\cos^2\theta\sinh^2\frac{\epsilon}{\cos\theta}}$$ by approximating $\sinh t$ with $t$. If we approximate $\sinh t$ with $t\left(1+\frac{t^2}{\pi^2}\right)$ or, even better, with $ t\, e^{t^2/6}$, we get: $$ I \approx \frac{1}{\varepsilon^2}\int_{0}^{\pi/2}\exp\left(-\frac{\epsilon^2}{3\cos^2\theta}\right)\,d\theta=\frac{\pi}{2\epsilon^2}\operatorname{Erfc}\left(\frac{\epsilon}{\sqrt{3}}\right)\approx\color{red}{\frac{\pi}{2\epsilon^2}\left(1-\frac{2\epsilon}{\sqrt{3\pi}}\right)}.$$
