How to prove a closed expression for $\int_0^{\infty } \left(\frac{1-e^{-q}}{q}\right)^n \, dq$? Recently, in the context of https://math.stackexchange.com/a/3392284/198592 , I stumbled over the integral
$$f(n) = \int_0^{\infty } \left(\frac{1-e^{-q}}{q}\right)^n \, dq\tag{1}$$
and wanted to find a closed expession for it.
The integral is convergent $n\gt 1$. Indeed, close to $q=0$ the integrand behaves as $1-\frac{q}{2}+\frac{q^2}{6}+O(q^3)$ so that there is no songularity there. For $q\to \infty$ the integrand becomes $\frac{1}{q^n}$. Hence the integral in convergent under the conditions stated.
Here we confine ourselves to the case of integer powers, i.e. $n=2, 3, ...$.
The first few values can be easily calculated
$$f(n=2..5) = \left\{\log (4),\log \left(\frac{81 \sqrt{3}}{64}\right),\frac{88 \log (2)}{3}-18 \log (3),\frac{5}{24} (-544 \log (2)+162 \log (3)+125 \log (5))\right\}\tag{2}$$
A general formula seemed hard to guess from these cases. 
Alternatively, a direct attack on the integral seems to suffer from the singular negative power of $q$ after a binomial expansion of the integrand.
I suggest you try for yourself to find it.
Here is my result
$$f(n) = \frac{1}{(n-1)!}\sum _{k=1}^n (-1)^{k+n} k^{n-1} \log (k) \binom{n}{k}\tag{3}$$
I did not know it in advance. For instance I haven't found it in Gradshteiyn/Ryshik. (EDIT After finishing the OP I found that 3.411.19 is related).
But now we have it it should not be difficult for the reader to prove it.
 A: You may use integration by parts and the Laplace transform/Frullani's theorem. For instance, in the $n=2$ case we have
$$ \int_{0}^{+\infty}\left(\frac{1-e^{-x}}{x}\right)^2\,dx \stackrel{\text{Binomial}}{=} \int_{0}^{+\infty}\frac{1-2e^{-x}+e^{-2x}}{x^2}\,dx \stackrel{\text{IBP}}{=}2\int_{0}^{+\infty}\frac{e^{-x}-e^{-2x}}{x}\,dx = 2\log(2) $$
while in general (I hope not to mess signs here)
$$\begin{eqnarray*}\int_{0}^{+\infty}\left(\frac{1-e^{-x}}{x}\right)^n&=&\int_{0}^{+\infty}\sum_{k=0}^{n}\binom{n}{k}(-1)^k e^{-kx}\frac{dx}{x^n}\\&=&\frac{(-1)^{n-1}}{(n-1)!}\int_{0}^{+\infty}\sum_{k=1}^{n}\binom{n}{k}(-1)^k k^{n-1}e^{-kx}\frac{dx}{x}\\&=&\frac{(-1)^{n-1}}{(n-1)!}\int_{0}^{+\infty}\sum_{k=1}^{n}\binom{n}{k}(-1)^k k^{n-1}\frac{ds}{k+s}\end{eqnarray*} $$
by using $\int_{0}^{+\infty}f(x)\frac{dx}{x}=\int_{0}^{+\infty}(\mathcal{L}f)(s)\,ds$. Of course $\sum_{k=1}^{n}\binom{n}{k}(-1)^k k^{n-1}=0$ (which follows by applying $n$ times the forward difference operator to the polynomial $x^{n-1}$) ensures the convergence of the last integral and provides the following closed form:
$$\int_{0}^{+\infty}\left(\frac{1-e^{-x}}{x}\right)^n\,dx= \frac{(-1)^n}{(n-1)!}\sum_{k=1}^{n}\binom{n}{k}(-1)^k k^{n-1}\log(k).$$
A: For integers $m,n$ ($2\leqslant m\leqslant n$), consider $$I(m,n)=\int_0^\infty x^{-m}(1-e^{-x})^n\,dx.$$
Integration by parts gives $$I(m,n)=\frac{n}{m-1}\int_0^\infty x^{-(m-1)}(1-e^{-x})^{n-1}e^{-x}\,dx,\tag{1}\label{reduction}$$
hence for $m>2$ we have a recurrence $$I(m,n)=\frac{n}{m-1}\big(I(m-1,n-1)-I(m-1,n)\big).\tag{2}\label{recurrence}$$
To handle $m=2$, recall the following extension of Frullani integrals: $$\int_0^\infty\left(\sum_{k=1}^{n}b_k e^{-a_k x}\right)\frac{dx}{x}=-\sum_{k=1}^{n}b_k\ln a_k$$ for $a_1,\ldots,a_n>0$ and $b_1+\ldots+b_n=0$. (Proven by induction.) Thus, from $\eqref{reduction}$, $$I(2,n)=-n\sum_{k=0}^{n-1}\binom{n-1}{k}(-1)^k\ln(k+1)=\sum_{k=1}^{n}\binom{n}{k}(-1)^k k\ln k.$$
This is the base case of $$\boxed{I(m,n)=\displaystyle\frac{1}{(m-1)!}\sum_{k=1}^{n}\binom{n}{k}(-1)^{m+k}k^{m-1}\ln k}$$ which is proven using induction on $m$ and the recurrence $\eqref{recurrence}$.
A: Notice that
$$\frac{1-e^{-q}}{q} = \int_0^1 e^{-qx}dx$$
So for $n>1$ an integer, we can use Fubini's theorem to rewrite the integral in the following way:
$$\int_0^\infty \left(\int_0^1 e^{-qx}dx\right)^n dq = \int_{[0,1]^n} \int_0^\infty e^{-(x_1+\cdots+x_n)q}dq dx_1\cdots dx_n $$
$$= \int_{[0,1]^n}\frac{1}{x_1+\cdots+x_n}dx_1\cdots dx_n$$
which seems like the kind of integral we should ask a statistician friend about. 

We can use $x_1+\cdots+x_n = n\left(1-\left(1-\frac{x_1+\cdots+x_n}{n}\right)\right)$ to get a geometric series in terms of $\left(1-\frac{x_1+\cdots+x_n}{n}\right)$. Using multinomial expansion to expand the geometric series powers and integrating the subsequent product of monomial terms leaves us with the following double sum:
$$\frac{1}{n}\sum_{k=0}^\infty \sum_{m_0+\cdots +m_n=k} {k \choose m_0,\cdots,m_n} \left(-\frac{1}{n}\right)^{k-m_0} \prod_{j=1}^n \frac{1}{m_j+1}$$
