Class of integrals: $I(a)=\int_0^\infty \frac{dx}{e^x+ax}$

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*}

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.

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.

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}.$$

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.