I'm looking to compute the integral $$I = \int_0^\infty dy \, \sinh^{d-2 \epsilon -1} y $$ as a Laurent series in $\epsilon$. (For those who are familiar with it, this is related to dimensional regularization.)
My attempt so far has been to use the familiar identity for the $\Gamma$ function, $$\alpha^{-s} \Gamma(s) = \int_0^\infty dx x^{s-1} e^{-\alpha x}.$$
From here, I can do the typical thing where I write $$I = \frac{1}{\Gamma(1+2\epsilon -d)} \int_0^\infty dy \, \int_0^\infty dx\, x^{2\epsilon-1} e^{-x\sinh(y)}.$$
So far, these expressions should agree in the limit where $\epsilon \rightarrow 0$ for all positive integer $d$. Now, the trouble comes when I do the typical trick of reordering the integral to instead integrate with respect to $y$ first, yielding
$$I = \frac{1}{\Gamma(1+2\epsilon -d)} \int_0^\infty dx \, x^{2\epsilon-1} \int_0^\infty dy\, e^{-x\sinh(y)}.$$
I've used Maple to compute the $y$ integral, obtaining $$\int_0^\infty dy\, e^{-x\sinh(y)} = \frac{\pi}{2} \left[H_0(x) - Y_0(x) \right]$$ where $H_0(x)$ is the Struve function of order $0$ and $Y_0(x)$ is the Bessel function of the second kind (of order $0$).
From here, Maple can (again) integrate directly the next integral. I won't bore you with the details, but the bottom line is that for even $d$, $$\lim_{\epsilon \rightarrow 0} I = \text{finite!}$$.
Now one can easily see that the original integral is clearly divergent for any integer $d>1$.
So what went wrong here? I know that I am playing with fire by working with divergent quantities and rearranging integrals all over the place (don't hate me, I'm a physicist!), but this seems bizarre to say the least.
Does anyone have any intuition as to why this may be going wrong, and how I can still do this type of manipulation legally?
Alternatively, does anyone have any idea how to explicitly compute the integral $I$ without resorting to the $\Gamma$ function trick?