How to approximate the integral $\int_{\mu}^{\infty} \frac{\sqrt{x^2-\mu^2}}{e^x-1} dx$ for a constant $\mu \gg 1$ I think this integral has no exact solution, but I've found that you can approximate it for large $\mu$ as:
$$\int_{\mu}^{\infty} \frac{\sqrt{x^2-\mu^2}}{e^x-1} dx \approx \sqrt{\frac{\pi \mu}{2}} e^{-\mu}$$
This is used a couple of times in a book I'm reading. I've checked it numerically and it works, but I wanted to get the symbolic expression.
I've tried expanding the square root:
$$\sqrt{x^2-\mu^2}=\mu\sqrt{-(1-x^2/\mu^2)}\approx i\mu \left(1-\frac{1}{2}\frac{x^2}{\mu^2}-\frac{1}{8}\frac{x^4}{\mu^4}-...\right)$$
That means that I get an imaginary solution as far as I know, because the first term is a simple logarithmic integral:
$$\int_{\mu}^{\infty} \frac{i\mu}{e^x-1}dx=-i\mu\ln{|e^{-\mu}-1|}$$
And the following terms are integrals that I don't know how to solve. So I guess that this is not the way to proceed.
Any ideas on how to do this approximation, or a more exact solution? Appreciate your help.
 A: With the substitution $z=x-\mu$, the integral becomes
$$\int_{\mu}^{\infty} \frac{\sqrt{x^2-\mu^2}}{e^x-1}\,dx=\int_0^\infty \frac{\sqrt{z^2+2\mu z}}{e^{\mu+z}-1}\,dz=\frac{\sqrt{2\mu}}{e^\mu }\int_0^\infty\frac{\sqrt{z+z^2/(2\mu)}}{e^z-e^{-\mu}}\,dz$$
Taking the limit $\mu\to\infty$ on the integrand and substituting $w=\sqrt{z}$, the integral becomes  $$\int_0^\infty \sqrt{z}\,e^{-z}\,dz=\int_0^\infty 2w^2 e^{-w^2}\,dx=\int_{-\infty}^\infty w^2 e^{-w^2}\,dw.$$
This is a standard variation on the Gaussian integral and evaluates to $\sqrt{\pi}/2$. Hence $$\int_{\mu}^{\infty} \frac{\sqrt{x^2-\mu^2}}{e^x-1}\,dx\sim \frac{\sqrt{2\mu}}{e^\mu}\frac{\sqrt{\pi}}{2}=\sqrt{\frac{\pi u}{2}}e^{-\mu}$$ in the $\mu\to\infty$ limit as claimed.
A: Having been stuck with this integral, now using @Semiclassical's elegant approach, I think that we can go further.
$$I=\frac{\sqrt{2\mu}}{e^\mu }\int_0^\infty\frac{\sqrt{z+\frac {z^2}{2\mu}}}{e^z-e^{-\mu}}\,dz\tag1$$
Expanding $\sqrt{z+\frac {z^2}{2\mu}}$ as a series, we have
$$\color{blue}{I=2\sqrt{\pi}\,\,\sum_{n=0}^\infty \frac{\Gamma \left(n+\frac{3}{2}\right)\,\,\text{Li}_{n+\frac{3}{2}}\left(e^{-\mu }\right)}{(2\mu)^{n-\frac{1}{2}}\,\,\Gamma \left(\frac{3}{2}-n\right)\,\, \Gamma
   (n+1)} }\tag2$$
and the first term is then
$$a_0=\sqrt{\frac{\pi \mu}{2}}  \text{Li}_{\frac{3}{2}}\left(e^{-\mu }\right)=\sqrt{\frac{\pi \mu}{2}} \sum_{k=1}^\infty \frac {e^{-k\mu}}{k^{\frac 32}}$$ and then the required approximation using the first term of the last summation. Doing the same (that is to say using $\text{Li}_{n+\frac{3}{2}}\sim e^{-\mu}~~\forall n$) a better approximation would be
$$I=\sqrt{\frac{\pi\mu}{2}} e^{-\mu } \left(1+t-\frac{5 t^2}{6}+O\left(t^3\right) \right) \qquad \text{with} \qquad t=\frac{3}{8 \mu }$$
Edit
To better see the approximation, we could rewrite
$$I=\sqrt{\frac{\pi \mu}{2}} e^{-\mu } \big[e^{\mu } \,\text{Li}_{\frac{3}{2}}\left(e^{-\mu }\right) \big]\Bigg[1+\sum_{n=1} ^\infty B_n\Bigg]$$ with
$$B_n=\frac{ \Gamma \left(n+\frac{3}{2}\right)}{(2\mu)^n\,\,\Gamma   \left(\frac{3}{2}-n\right)\,\,\Gamma (n+1)}\,\,\frac{\text{Li}_{n+\frac{3}{2}}\left(e^{-\mu }\right)  } {\text{Li}_{\frac{3}{2}}\left(e^{-\mu }\right)  }$$
$$e^{\mu } \,\text{Li}_{\frac{3}{2}}\left(e^{-\mu }\right)=1+\frac{e^{-\mu }}{2 \sqrt{2}}+\frac{e^{-2 \mu }}{3 \sqrt{3}}+\cdots$$
$$\frac{\text{Li}_{n+\frac{3}{2}}\left(e^{-\mu }\right)  } {\text{Li}_{\frac{3}{2}}\left(e^{-\mu }\right)  }=1- 2^{-n-\frac{3}{2}} \left(2^n-1\right)e^{-\mu }+\cdots$$
