I am interested in reasonable to obtain upper and lower bounds, as well as techniques used in doing so, for $\int_0^\infty e^{-u(u^\epsilon -1)}\,du$ for epsilon in some right-neighborhood of zero, i.e. valid for all $0<\epsilon<\delta$ with $\delta>0$. I believe I was able to prove the integral is $\mathcal{O}(1/\epsilon)$ as $\epsilon\to0^+$ by using the fact that $u^\epsilon -1\ge \epsilon\log u$ and then $\int_e^\infty e^{-\epsilon u\log u}\,du\le\int_e^\infty e^{-\epsilon u}\,du$ but didn't get much further.
Edit: So it seems that thanks River Li's comment, we now have the bounds $$\frac{c}{\sqrt{\epsilon}}<\int_0^\infty e^{-u(u^\epsilon -1)}\,du < \frac{C}{\epsilon}$$ for some constant $C, c$ in some right neighborhood of zero. Can we do better?
The integral is asymptotically $-1/(\epsilon\log\epsilon)$ as $\epsilon\to0^+$; that is, $$\int_0^\infty e^{-u(u^\epsilon-1)}\,du=-\frac{f(\epsilon)}{\epsilon\log\epsilon},\qquad\lim_{\epsilon\to0^+}f(\epsilon)=1.$$ Indeed, put $u=-x/(\epsilon\log\epsilon)$, then $f(\epsilon)=\int_0^\infty e^{g(x,\epsilon)}\,dx$ with $\lim\limits_{\epsilon\to0^+}g(x,\epsilon)=-x$ and $$g(x,\epsilon)\leqslant\frac{x}{\log\epsilon}\log\frac{-x}{\epsilon\log\epsilon}\leqslant-x+\frac1e$$ because of your $u^\epsilon-1\geqslant\epsilon\log u$ and $y\log y\geqslant-1/e$. Thus (by DCT) $$\lim_{\epsilon\to0^+}f(\epsilon)=\int_0^\infty e^{-x}\,dx=1.$$ A by-product is an upper bound $\color{blue}{-e^{1/e}/(\epsilon\log\epsilon)}$ for the integral, and a lower bound of the same scale should exist (for a bounded range of $\epsilon$), although an explicit one may be harder to obtain.