A class of generalized Integrals involving polygamma functions I recently came across some nice integrals and I have a few questions about them: You may have heard of the Euler-Mascheroni constant $\gamma$ and if you did so, you may know the following integral: 
$I_1=\int_0^\infty e^{-x}\ln(x) dx=-\gamma$
But what happens if you somehow change this integral a little bit? If you replace $\ln(x)$ with $\ln^2(x)$, you get an suprising answer. 
$I_2=\int_0^\infty e^{-x}\ln^2(x) dx=\gamma^2+\frac{\pi^2}{6}$
This already looks quite weird. Why would the square of $\gamma$ appear right here. And what has this to do with $\zeta(2)=\frac{\pi^2}{6}$. When I've tried to generalize this even more things got even weirder.
$I_m=\int_0^\infty e^{-x}\ln^m(x) dx$
If you now evaluate this with $m\in\mathbb{Z}^+$ you get answers involving integer powers of $\gamma$ and the Riemann Zeta function: $\gamma^m, \zeta(m)$. Intuitively, why would powers of $\gamma$ pop out right here? They seem to appear nowhere else. And what has this to do with the zeta function? Is there a general formula for $I_m$, where $m$ is an integer or even a real number? I really have no idea.
I'm thanking every one who replies.
 A: Starting from the integral representation of Gamma
$$
\int_0^\infty  {e^{\, - \,t} t^{\,x} \,dt}  = \Gamma (x + 1)
$$
and taking the derivative wrt $x$
$$
{\partial  \over {\partial x}}\int_0^\infty  {e^{\, - \,t} t^{\,x} \,dt}  = {\partial  \over {\partial x}}\int_0^\infty  {e^{\, - \,t} e^{\,x\ln \,t} \,dt}  = \int_0^\infty  {e^{\, - \,t} t^{\,x} \ln t\,dt}  = \Gamma '(x + 1)
$$
and the second derivative
$$
{{\partial ^{\,2} } \over {\partial x^{\,2} }}\int_0^\infty  {e^{\, - \,t} t^{\,x} \ln t\,dt}  = \int_0^\infty  {e^{\, - \,t} t^{\,x} \left( {\ln t} \right)^{\,2} \,dt}  = \Gamma ''(x + 1)
$$
and so forth.
So what you are looking for corresponds to $m!$ times the coefficient of $x^m$ in the power expansion of $\Gamma (x+1)$.
$$
\int_0^\infty  {e^{\, - \,t} t^{\,x} \left( {\ln t} \right)^{\,m} \,dt}  = \left. {\Gamma ^{\,\left( m \right)} (x + 1)\;} \right|_{\,x\, = \,0}  = m!\left[ {x^{\,m} } \right]\Gamma (x + 1)
$$
These can be expressed as linear combination of the coefficients of the expansion of $\ln \Gamma (1+x)$, i.e. in terms of $\psi ^{\,k} (1)$, and thus in terms 
of $\gamma ,\; \zeta(k)$.
