I am curious whether the integral
$$\int_1^\infty \frac{1}{x(e^{a x}-1)} dx, \qquad a>0$$
has an exact solution. The best I can do so far is to see that is has the lower bound
$$\int_1^\infty \frac{1}{x(e^{a x}-1)} dx \geqslant \int_1^\infty \frac{1}{xe^{a x}} dx = E_1(a),$$
where $E_1(a)$ is the exponential $E_n$ integral.
Mathematica can't solve it, although I've seen Mathematica fail to solve definite integrals that have exact solutions in terms of obscure special functions, even if Mathematica knows those special functions.
Any tips or ideas would be appreciated.
EDIT:
I've also found an upper bound,
$$\int_1^\infty \frac{1}{x(e^{a x}-1)} dx \leqslant \int_1^\infty \frac{1}{e^{a x}-1} dx = 1 - \frac{\ln(e^a-1)}{a},$$
so I at least have upper and lower bounds. But maybe there are stricter bounds, or even better, an exact solution?
You can write $$ \frac{1}{e^{ax}-1} = \sum_{j=1}^\infty e^{-ajx}$$ converging uniformly for $x \ge 1$. This leads to
$$ \sum_{j=1}^\infty \Gamma(0,ja) = \sum_{j=1}^\infty \text{Ei}_1(ja)$$
but still not a closed form.
EDIT: For $a=1$, for example, Maple gives the approximate numerical result $$0.2867924309492973407247190257941567908651386511054584788045817762268118195038714403648237791703892668$$ which produces no match in the Inverse Symbolic Calculator. Of course there is no way to prove that there is no closed form.