Asymptotics of $\int_{0}^{+\infty}\!\!\frac{dx}{\sinh^2(\epsilon \sqrt{x^2+1}) } $ for $\epsilon$ near $0$ 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.
 A: 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*}
A: 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)}.$$
