Esto no es una respuesta completa, en lugar de dos observaciones:
Deje $\mathcal{I}_{k,n}$ denotar la integral. Integración por partes, asumiendo $k \geqslant 2$ $n \geqslant 1$ da la ecuación de recurrencia:
$$
\mathcal{I}_{k+1,n} = n \mathcal{I}_{k,n-1} + k \mathcal{I}_{k,n} \etiqueta{1}
$$
que se pueden utilizar para reducir el problema a $k=1$. Deje $\mathcal{J}_{k,n} = \frac{1}{n!}\mathcal{I}_{k,n}$. Es fácil ver que $\mathcal{J}_{k,n}$ satisface la relación de recurrencia de la unsigned los números de Stirling de primera especie:
$$
\mathcal{J}_{k+1,n} = \mathcal{J}_{k,n-1} + k \mathcal{J}_{k,n} \etiqueta{2}
$$
Por el otro lado:
$$
\frac{\mathrm{d}^n}{\mathrm{d}^n} x^{m-1} = x^{m-1} \log(x)^n
$$
Por lo tanto:
$$
\int_0^\infty \mathrm{e}^{-x} x^{k-1} \log(x)^n \mathrm{d} x = \left. \frac{\mathrm{d}^n}{\mathrm{d}^n} x^{m-1} \int_0^\infty \mathrm{e}^{-x} x^{m-1} \mathrm{d} x \right|_{s=k} = \left. \frac{\mathrm{d}^n}{\mathrm{d}^n} \Gamma(s) \right|_{s=k} = \frac{\mathrm{d}^n}{\mathrm{d} k^n} \exp(\log\Gamma(k) )
$$
Ahora uso Faà di Bruno fórmula:
$$
\frac{\mathrm{d}^n}{\mathrm{d} k^n} \exp(\log\Gamma(k) ) = \Gamma(k) \sum_{m=0}^n B_{n,m}\left( \psi(k), \psi^{(1)}(k), \ldots, \psi^{(n-1)}(k)\right)
$$
donde $\psi(x)$ es la función digamma, y $\psi^{(m)}(x)$ es el polygamma función, y donde $B_{n,k}$ stand para la Faà di Bruno polinomios, que son homogéneos:
$$
B_{n,k}\left(\lambda x_1, \lambda x_2, \ldots, \lambda x_n\right) = \lambda^k B_{n,k}\left(x_1, x_2, \ldots, x_n\right)
$$
$$
B_{n,k}\left(\lambda x_1, \lambda^2 x_2, \ldots, \lambda^n x_n\right) = \lambda^n B_{n,k}\left(x_1, x_2, \ldots, x_n\right)
$$
Suponga que la recurrencia de la relación $(1)$ fue utilizado para reducir a $k=1$ de los casos. Para esto $k$, los valores de polygammas están relacionados con la $\zeta$-valores:
$$
\psi(1) = -\gamma \qquad \psi^{(m)}(1) = \left. \frac{\mathrm{d}^m}{\mathrm{d}^m} \psi(s) \right|_{s=1} = (-1)^{m+1} m! \sum_{n=1}^\infty \frac{1}{n^{m+1}} = (-1)^{m+1} m! \zeta(m+1)
$$
Así
$$
\int_0^\infty \mathrm{e}^{s} \log(x)^n \mathrm{d} x = (-1)^n \sum_{m=0}^n B_{n,m} \left(\gamma \zeta(2), 2! \zeta(3), \ldots, n! \zeta(n+1)\right)
$$
y esto explica la clasificación.