Solution of the integral $\int_1^\infty \frac{1}{x(e^{a x}-1)} dx$? 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?
 A: This is the best I could come up with: by the formula for the sum of the geometric series, valid for the common ratio $e^{-ax}$ which has magnitude smaller than $1$ for all $x \in [1,\infty)$, and by Weierstrass' theorem on total convergence,
$$\begin{split}
I(a) &= \int_1^\infty \frac{1}{x} \cdot \frac{1}{e^{ax} -1} dx = \int_1^\infty \frac 1 x \frac{e^{-ax}}{1 - e^{-ax}} dx = \int_1^\infty \frac{e^{-ax}}{x} \sum_{n=0}^\infty (e^{-ax})^n dx \\ &= \sum_{n=1}^\infty \int_1^\infty \frac{e^{-anx}}{x}dx
= \sum_{n=1}^\infty \int_{an}^\infty e^{-y} y^{-1} dy = \sum_{n=1}^\infty \Gamma(0,an)
\end{split}$$
where $\Gamma(s,t)$ is the incomplete Gamma function. Since for all $\xi > 0$ one has that $\Gamma(0,\xi) = -\operatorname{Ei}(-\xi) = E_1(\xi)$, we get
$$I(a) = -\sum_{n=1}^\infty \operatorname{Ei}(-an) = \sum_{n=1}^\infty E_1(an).$$

Edit. We may also employ Feynman's trick:
$$\begin{split}
I'(a) &= \int_1^\infty \frac{\partial}{\partial a}\left[\frac{1}{x(e^{ax} - 1)} \right] dx = \int_1^\infty \frac{-e^{ax}}{(e^{ax} - 1)^2} dx \\ &= \int_1^\infty \frac{dx}{2-2\cosh(ax)} = \frac{1}{a} \int_a^\infty \frac{dy}{1-\cosh y} \\
&= \frac 1 {a(1-e^a)} = \frac 1{2a} \left[1 - \coth\left(\frac a 2\right) \right]
\end{split}$$
Therefore, since $I(\infty) = 0$,
$$I(a) = I(\alpha)\Big|_\infty^a = \int_{\infty}^a I'(\alpha) d\alpha = \int_\infty^a \frac {d\alpha} {\alpha(1-e^\alpha)} = \int_\infty^a \frac{1}{2\alpha}\left[1 - \coth\left(\frac \alpha 2 \right) \right]d\alpha. $$
But I don't know whether these integrals are solvable in terms of known functions.
A: 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.
A: I have been working on polynomial expansions for this integral form. In general we have
$$\begin{align}
\int_{x}^{\infty}\frac{e^{-at}}{t\left(e^{-t}-1\right)}dt=\space
&\frac{\ln2\pi}{2}-\left(a-\frac{1}{2}\right)\gamma-\ln\Gamma{\left(a\right)}
\\
&-\frac1x-\left(a-\frac12\right)\ln x+\sum_{n=2}^{\infty}\left(-1\right)^{n}B_n\left(a\right)\frac{x^{n-1}}{n!\left(n-1\right)}
\end{align}\tag{1}$$
and
$$\begin{align}
-\int_{x}^{\infty}\frac{2e^{-at}}{t\left(e^{-t}+1\right)}dt=\space
&\ln2+\gamma+2\ln\Gamma\left(\frac{a+1}{2}\right)-2\ln\Gamma\left(\frac{a}{2}\right)
\\
&+\ln x+\sum_{n=1}^{\infty}\left(-1\right)^{n}E_n\left(a\right)\frac{x^{n}}{n!n}
\end{align}\tag{2}$$
For this problem we can make use of the following transformations.
$$\int_{x}^{\infty}\frac{e^{-at}}{t\left(e^{-bt}-1\right)}dt=\int_{bx}^{\infty}\frac{e^{-\frac{a}{b}t}}{t\left(e^{-t}-1\right)}dt$$
for $b\in \mathbb{R}$ and
$$\int_{x}^{\infty}\frac{e^{-at}}{t\left(e^{bt}-1\right)}dt=-\int_{bx}^{\infty}\frac{e^{-\left(\frac{a}{b}+1\right)t}}{t\left(e^{-t}-1\right)}dt$$
for $b\ge 0$.
Using this second transformation with equation $(1)$ and substituting the values from the original post, we get
$$\int_{1}^{\infty}\frac{1}{t\left(e^{bt}-1\right)}dt=-\frac{\ln2\pi}{2}+\frac{\gamma}{2}+\frac{1}{b}+\frac{\ln b}{2}-\sum_{n=2}^{\infty}\frac{B_n}{n!\left(n-1\right)}b^{n-1}$$
for $0\le b \le 2\pi$.
