How to evaluate integral $\int_{0}^{\infty} \left(\frac{1-e^{-x}}{x}\right)^n dx$. First, according to
\begin{align*}
\int_{0}^{\infty} x^{-m}(1-e^{-x})^{n} \, dx
=\frac{n}{1-m}\int_{0}^{\infty} x^{1-m}(1-e^{-x})^{n} \, dx
-\frac{n}{1-m}\int_{0}^{\infty} x^{1-m}(1-e^{-x})^{n-1} \, dx,
\end{align*}
which  can be denoted by 
\begin{align*}
I_{m,n} = \frac{n}{1-m}I_{m-1,n}-\frac{n}{1-m}I_{m-1,n-1}. \tag{1}
\end{align*}
I also find
\begin{align*}
\int_{0}^{\infty} e^{-nx}\ln x \, dx = -\frac{\gamma}{n}-\frac{\ln n}{n}.
\end{align*}
So 
\begin{align*}
\int_{0}^{\infty} x^{-2} \left(1 - e^{-x}\right)^{n} \, dx
&= -\int_{0}^{\infty} \left(1 - e^{-x} \right)^{n} \, d\left(\frac{1}{x}\right)
= n \int_{0}^{\infty} \frac{\left( 1 - e^{-x} \right)^{n - 1} e^{-x}}{x} \, dx \\
&= n \int_{0}^{\infty}  \frac{\sum_{k = 0}^{n - 1} \binom{n-1}{k} ( -1 )^{k} e^{-(k + 1)x} }{x} \, dx \\
&= n\sum_{k = 0}^{n - 1} \binom{n-1}{k} (-1)^{k} (k+1) \int_{0}^{\infty} e^{-(k+1)x} \ln x \, dx \\
&= n\sum_{k = 0}^{n - 1} \binom{n-1}{k} (-1)^{k} (k+1) \left( -\frac{\gamma}{k+1} - \frac{\ln (k + 1)}{k + 1} \right) \\
&= -n\sum_{k = 0}^{n - 1} \binom{n-1}{k} (-1)^{k} \left( \gamma + \ln (k + 1) \right) \\
&= -n\sum_{k = 0}^{n - 1} \binom{n-1}{k} (-1)^{k} \ln (k + 1).
\end{align*}
Hence, it seems that we can solve this integral by recursive relation $\text{(1)}$, but how to get a accurate result about the integral.

A basic idea is to use integration by parts, as you did for $I_{2,n}$.
Let $2 \leq m \leq n$ be integers and denote $f(x) = (1 - e^{-x})^{n}$. We have two observations:


*

*$f^{(k)}(x) = O(x^{n-k})$ near $x = 0$ for $0 \leq k \leq n$.

*$f^{(k)}(x) = O(e^{-x})$ near $x = \infty$ for $1 \leq k \leq n$.


Using this, applying integration by parts $m$ times, we obtain
\begin{align*}
I_{m,n}
= \int_{0}^{\infty} \frac{(1 - e^{-x})^{n}}{x^{m}} \, dx
&= -\frac{1}{(m-1)!} \int_{0}^{\infty} f^{(m)}(x) \log x \, dx \\
&= \frac{(-1)^{m-1}}{(m-1)!} \sum_{k=1}^{n} \binom{n}{k} (-1)^{k} k^{m} \int_{0}^{\infty} e^{-kx} \log x \, dx \\
&= \frac{(-1)^{m}}{(m-1)!} \sum_{k=1}^{n} \binom{n}{k} (-1)^{k} k^{m-1} (\gamma + \log k) \\
&= \frac{(-1)^{m}}{(m-1)!} \sum_{k=1}^{n} \binom{n}{k} (-1)^{k} k^{m-1} \log k.
\end{align*}
