Class of integrals: $I(a)=\int_0^\infty \frac{dx}{e^x+ax}$ I'm investigating integrals in the form
$$I(a):=\int_0^\infty \frac{dx}{e^x+ax}$$
So far, I haven't been able to find any special values other than $I(0)=1$, and I've only managed to evaluate these similar indefinite integrals:
$$\int \frac{x-1}{e^x+ax}dx=-\frac{\ln(1+axe^{-x})}{a}+C$$
$$\int \frac{xdx}{e^x+x+1}=-\ln(1+e^{-x}(x+1))+C$$
I've also found the following series representation for $I(a)$:
$$I(a)=\sum_{n=0}^\infty \frac{(-a)^n n!}{(n+1)^{n+1}}$$
...which looks remarkably similar to the Maclaurin series for the Lambert-W function.
QUESTION: Can anyone find any non-trivial special values of this integral? I find this unlikely because of the weird series representation of $I(a)$, so if this isn't feasible, can anyone find any interesting properties or functional/differential equations for $I(a)$?
UPDATE: I've managed to show that
$$\lim_{a\to\infty }\frac{aI(a)}{\ln(a)}=1$$
 A: Here are my notes on this so far with $a=1$. I hope these are useful.
It seems that
$$
\int \frac{dx}{e^x + x} = \sum_{n=0}^\infty \left[\sum_{k=0}^n\frac{(n-k+1)^k}{k!}\right]\frac{(-1)^nx^{n+1}}{(n+1)}
$$
this doens't look like the Cauchy product of series. It seems that
$$
\sum_{k=0}^n\frac{(n-k+1)^k}{k!} \sim \kappa e^{W(1)n}
$$
where $W(x)$ is the Lambert-W function, and $W(1)=\Omega$. Taking a limit
$$
\lim_{n \to \infty} \left(e^{- W(1) n}\sum_{k=0}^n\frac{(n-k+1)^k}{k!}\right)=\kappa\approx 1.251190909867859\cdots
$$
if this was a valid series expansion, then the function is zero at $x=0$, so the infinite asymptotics may yield a result. (I may be wrong)
Edit: the more general series with $a$ seems to be
$$
\int \frac{dx}{e^x+ax} = \sum_{n=0}^\infty \left[\sum_{k=0}^n \frac{a^k(k+1)^{n-k}}{(n-k)!} \right] \frac{(-1)^n x^{n+1}}{n+1}
$$
a similar treatment seems to give 
$$
\sum_{k=0}^n \frac{a^k(k+1)^{n-k}}{(n-k)!} \sim e^{\left(\log(a)+W\left(\frac{1}{a}\right)\right)n}
$$
this came from guessing and using the inverse symbolic calculator.
Note:
If it assists with Sangchul Lee's integral representation, it appears that
$$
\Re\left(\int_0^{1/e} u^n W(-u) \; du \right)= \frac{n!-q(n)}{(n+1)^{n+2} e^{n+1}}, \;\; q(n) \in \mathbb{N}
$$
where the $q's$ go like $3,9,53,462,5319,76008,1296273,25679664,579336363,\cdots$, for $n=0,1,\cdots$ but it is not clear what these are. Further it seems that
$$
\Re\left(\int_0^{1/e} u^n W(-u) \; du \right)= (n+1)^{-n-2} \Gamma (n+1)+\left(\frac{1}{n+1}\right)^{n+3} ((n+1) \Gamma (n+2,n+1)-\Gamma (n+3,n+1))
$$
so then 
$$
q(n) = \left(\frac{1}{n+1}\right)^n (n+1)^{2 n+1} \left(e^{n+1} E_{-n-1}(n+1)+1\right)
$$
for exponential integral function.
A: Starting from the series that you already got
$$
\eqalign{
  & I(a) = \int_0^\infty  {{{dx} \over {e^{\,x}  + ax}}}  = \int_0^\infty  {{{e^{\, - x} dx} \over {\left( {1 + axe^{\, - x} } \right)}}}  =   \cr 
  &  = \int_0^\infty  {\sum\limits_{0\, \le \,k} {\left( { - 1} \right)^{\,k} \left( {a^{\,k} x^{\,k} e^{\, - \,\left( {k + 1} \right)\,x} } \right)} \;dx} 
  = \sum\limits_{0\, \le \,k} {\left( { - 1} \right)^{\,k} {{k!} \over {\left( {k + 1} \right)^{k + 1} }}a^{\,k} }  \cr} 
$$
and which converges for
$$
\left| a \right|x/e^{\,x}  < \left| a \right|1/e < 1\quad  \Rightarrow \quad \left| a \right| < e
$$
From this related post we get
$$
\sum\limits_{1\, \le \,\,n} {{1 \over {n^{\,n} }}x^{\,n} }
  = x\sum\limits_{0\, \le \,\,n} {{1 \over {\left( {n + 1} \right)^{\,n + 1} }}x^{\,n} }
  = x\int_{\,0}^{\,1} {t^{\, - \,x\,t} dt} 
$$
and since
$$
A(z) = \sum\limits_{0\, \le \,n} {a_n \,z^n } \quad  \Leftrightarrow \quad \int_{\;t\, = \,0}^\infty  {e^{\, - \,t} A(z\,t)\,d\,t}
  = \sum\limits_{0\, \le \,n} {n!a_n z^{\,n} } 
$$
we get another integral representation
$$ \bbox[lightyellow] {  
\eqalign{
  & I( - x) = \sum\limits_{0\, \le \,\,n} {{{n!} \over {\left( {n + 1} \right)^{\,n + 1} }}x^{\,n} }
  = \int_0^\infty  {{{e^{\, - u} du} \over {1 - x\,u\,e^{\, - u} }}}  =   \cr 
  &  = \int_{\,u\, = \,0}^{\,\infty } {e^{\, - \,u} \int_{\,t\, = \,0}^{\,\,1} {t^{\, - \,x\,u\,t} dt\,} du}
  = \int_{\,t\, = \,0}^{\,1} {\int_{\,u\, = \,0}^{\,\infty } {e^{ - \,u\left( {1 + x\,t\ln t} \right)}\, dt\,} du}  =   \cr 
  &  = \int_{\,t\, = \,0}^{\,1} {{{dt} \over {\left( {1 + x\,t\ln t} \right)}}}  \cr} 
}$$
Now the second line tells us that
$$
I( - 1/s) = \int_{\,u\, = \,0}^{\,\infty } {e^{\, - \,u} \int_{\,t\, = \,0}^{\,\,1} {t^{\, - \,\,\left( {u/s} \right)\,t} dt\,} du} 
$$
i.e.
$$ \bbox[lightyellow] {  
\eqalign{
  & {1 \over s}I( - 1/s)
 = \int_{\,\alpha \, = \,0}^{\,\infty } {e^{\, - \,s\,\alpha } \left( {\int_{\,t\, = \,0}^{\,\,1} {t^{\, - \,\,\alpha \,t} dt\,} } \right)d\alpha }  =   \cr 
  &  = \int_0^\infty  {{{e^{\, - u} } \over {s - \,u\,e^{\, - u} }}du}  = \int_{\,t\, = \,0}^{\,1} {{{dt} \over {\left( {s + \,t\ln t} \right)}}}  \cr} 
}$$
so that our integral is tied to the Laplace transform of the interesting function
$\int_{0}^1 {t^{-xt}}dt=\sum_{n=1}^\infty \frac{x^{n-1}}{n^n} =  $ Sphd$(-x;1)$
cited by JJacquelin in his answer to the already cited post.
A: Notice first that
$$ I(a) = \int_{0}^{\infty} \frac{xe^{-x}}{1+axe^{-x}} \, dx. $$
Indeed, this follows from $\int_{0}^{\infty} \frac{x-1}{e^x + ax} \, dx = 0$ using OP's computation. Since the graph of $x \mapsto xe^{-x}$ is unimodal, for each $u$ in the range we may define the 'width' $l(u)$ of the graph of $xe^{-x}$ at height $u$.
$\hspace{8em}$ 
To be precise, we define $l(u)$ as the Lebesgue measure of the set $\{ x > 0 : xe^{-x} > u \}$. Then
$$ I(a)
= \int_{0}^{\infty} \left( \int_{xe^{-x}}^{\infty} \frac{du}{(1+au)^2} \right) \, dx
\stackrel{\text{(Fubini)}}{=} \int_{0}^{\infty} \frac{l(u)}{(1+au)^2} \, du. $$
Now $l$ can be written explicitly in terms of the Lambert W-function:
$$ l(u) = \begin{cases}
W(-u) - W_{-1}(-u), & \text{if } u \leq \frac{1}{e} \\
0, & \text{if } u > \frac{1}{e}
\end{cases} $$
So it follows that
$$ I(a)
= \int_{0}^{\frac{1}{e}} \frac{l(u)}{(1+au)^2} \, du
= \int_{0}^{\frac{1}{e}} \frac{W(-u) - W_{-1}(-u)}{(1+au)^2} \, du. \tag{1} $$
This suggests that the asymptotic behavior of $I(a)$ as $a\to\infty$ is intimately related to the asymptotic behavior of $W_{-1}(u)$ as $u\to 0$. For instance, using the fact that
$$ l(u) = -W_{-1}(-u) + \mathcal{O}(1) = -\log u + \log\log(1/u) + \mathcal{O}(1) $$
on $(0, 1/e]$ as $u\to0$, we obtain
$$ I(a) = \frac{\log a}{a} + \frac{\log\log a}{a} + \mathcal{O}\left(\frac{1}{a}\right) \quad \text{as } a \to \infty. \tag{2} $$

We also notice that for $n \geq 1$,
\begin{align*}
\left( \frac{d}{da} \right)^n (aI(a))
&= (-1)^{n-1} n! \int_{0}^{\infty} \frac{x^{n-1}e^{-nx}}{(1 + axe^{-x})^{n+1}} \, dx \\
&= \frac{(-1)^{n-1} n!}{a^n} \int_{0}^{\infty} \frac{u^{n-1}e^{-nu/a}}{(1 + ue^{-u/a})^{n+1}} \, du \\
&\sim \frac{(-1)^{n-1} (n-1)!}{a^n} \quad \text{as } a \to \infty.
\end{align*}
A: Not an answer, but an observation (to express my interest in your question myself). Another integral representation of $I(a)$ is
$$I(a)=\int_0^{+\infty} \frac{x\,dx}{e^x+ax}$$
(follows from the first of indefinite integrals in your question). Also
$$I'(a)=-\int_0^{+\infty} \frac{x\,dx}{(e^x+ax)^2}=-\int_0^{+\infty} \frac{x^2\,dx}{(e^x+ax)^2}.$$
A: Too long for comment.
Using integration by parts, easy to obtain for $m\ge0,\ n\ge1:$
$$J(m,n)= \int\limits_0^1 t^m \log^n t\,\mathrm dt = \dfrac {t^{m+1}}{m+1}\log^n t\Bigg|_0^1 - \dfrac n{m+1}\int\limits_0^1 t^m\log^{n-1}t\,\mathrm dt= -\dfrac{n}{m+1}J(m, n-1),$$
$$J(m,n)=(-1)^n\dfrac{n!}{(m+1)^{n+1}}.\tag1$$
This allow to calculate Taylor series for the integral
$$I(a,m,n) = \int\limits_0^\infty \dfrac{e^{-mx}x^n}{e^x+ax}\,\mathrm dx = \int\limits_0^1\dfrac{t^m\log^n t}{1-at\log t}\,\mathrm dt = \sum_{k=0}^\infty J(m+k,n+k)a^k,$$
$$I(a,m,n) = \sum_{k=0}^\infty(-1)^{n+k}\dfrac{(n+k)!}{(m+k+1)^{n+k+1}}a^k.\tag2$$
Formula $(2)$ can be useful for the further investigations.
