How to evaluate integral $\int_{0}^{\infty} \left(\frac{1-e^{-x}}{x}\right)^n dx$.

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