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.
The integral can be studied with Laplace's method. One has
$$
I = \int_1^{\infty}\left(\frac{\ln x}{x} \right)^{\lambda} dx = \int_1^{\infty} \exp\left(-\lambda \left(\ln x -\ln \ln x\right)\right)dx, \quad \lambda \rightarrow \infty
$$
Let $p(x) = \ln x -\ln \ln x$, a minimum occurs around $x=e$ and Taylor expanding around this point, one has
$$
p(x) = 1 + \frac{1}{2e^2} (x-e)^2 - \frac{5}{6e^3}(x-e)^3 +...
$$
Split the integral into two parts, the first one from $1$ to $e$ ($I_1$) and the second from $e$ to infinity ($I_2$). For the second integral, one may identify from the link the parameters $p(a) = 1,P=(2e^2)^{-1},\mu=2,Q=1,\lambda=1$, at the point $a=e$, which gives the dominating term
$$
I_2 = \int_1^{\infty}\left(\frac{\ln x}{x} \right)^{\lambda} dx \sim \frac{Q}{\mu}\Gamma\left(\frac{\lambda}{\mu} \right) \frac{e^{-\lambda p(a)}}{(Px)^{\lambda/\mu}} = e\sqrt{\frac{{\pi}}{2}} \frac{e^{-\lambda}}{\lambda^{1/2}}
$$
One may similarly show that $I_1 = I_2$ which then gives
$$
I = \int_1^{\infty}\left(\frac{\ln x}{x} \right)^{\lambda} dx \sim \sqrt{2\pi}e\lambda^{-1/2}e^{-\lambda}, \quad \lambda \rightarrow \infty
$$
to leading order. From the same link, it is shown how to include additional terms if necessary.
Best Answer
We have
\begin{align*} I(\lambda) &= 2 \int_{0}^{\infty} e^{-x^2}(\cosh x)^{\lambda} \, \mathrm{d}x \\ &= 2^{1-\lambda} \int_{0}^{\infty} e^{-x^2+\lambda x}(1 + e^{-2x})^{\lambda} \, \mathrm{d}x \\ &= e^{\lambda^2/4} 2^{1-\lambda} \int_{-\lambda/2}^{\infty} e^{-y^2}(1 + e^{-2y-\lambda})^{\lambda} \, \mathrm{d}y. \end{align*}
Using this and invoking dominated convergence theorem, it is easy to check that
$$ I(\lambda) \sim e^{\lambda^2/4}2^{1-\lambda}\sqrt{\pi} $$
as $\lambda \to \infty$. Furthermore, it is possible to show that the relative error of this asymptotic formula,
$$ \left| \frac{I(\lambda)}{e^{\lambda/4}2^{1-\lambda}\sqrt{\pi}} - 1 \right|, $$
decays exponentially fast as $\lambda \to \infty$: