En primer lugar, de acuerdo a
\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*}
que puede ser denotado por
\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*}
También me parece
\begin{align*} \int_{0}^{\infty} e^{-nx}\ln x \, dx = -\frac{\gamma}{n}-\frac{\ln n}{n}. \end{align*}
Por lo tanto
\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*}
Por lo tanto, parece que podemos resolver esta integral por la recurrente relación $\text{(1)}$, pero cómo obtener un resultado exacto de la integral.