Try to prove these beautiful formulas without using the Euler or the Weirstrass definitions of the Gamma function : $$ \ \ \displaystyle\int_{0}^{+\infty}{\mathrm{e}^{-x}\ln{x}\,\mathrm{d}x}=-\gamma \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left(1\right) $$ $$ \ \displaystyle\int_{0}^{+\infty}{\mathrm{e}^{-x}\ln^{2}{x}\,\mathrm{d}x}=\gamma^{2}+\displaystyle\frac{\pi^{2}}{6} \ \ \ \ \ \ \ \ \ \left(2\right) $$
Is there a generalisation giving the value of $ \int\limits_{0}^{+\infty}{\mathrm{e}^{-x}\ln^{n}{x}\,\mathrm{d}x} $, for any integer $ n $ ?
One can take derivatives of the gamma function to find what you want: $$ \int_0^\infty e^{-x} \log^n(x) dx = \frac{d^n}{ds^n} \int_0^\infty e^{-x} x^s dx\Big|_{s=0} = \frac{d^n}{ds^n}\Gamma(s+1)\Big|_{s=0}.$$ It's easy to program this procedure in a CAS. However, a closed-form formula is possible if one uses the known expansion $$ \log{\Gamma(s+1)} = -\gamma\,s + \sum_{k=2}^\infty \frac{\zeta(k)}{k} (-s)^k$$ You'd have to exponentiate, and that leads to the partial Bell polynomials.