According to WolframAlpha, the integral $$\int_0^{\infty} \log(x)^2 e^{-x} \, \mathrm{d}x$$ has closed form $\gamma^2 + \frac{\pi^2}{6}$, where $\gamma$ is the Euler-Mascheroni constant.
The term $\frac{\pi^2}{6}$ is $\zeta(2)$, and other integrals of the form $\int_0^{\infty} \log(x)^n e^{-x} \, \mathrm{d}x$ clearly have some relation to values of $\zeta(s)$ at integers. I do not know where to go with this. The usual trick I know for integrating something of the form $\log(x)^n f(x)$ is to compare integrals along the top and bottom of the slit in a keyhole contour, but this doesn't converge with the term $e^{-x}$.
Using $$\Gamma(z) = \int_0^\infty e^{-x}x^{z-1} dx$$ we obtain $$\int_0^\infty e^{-x}\ln^2 x dx = \Gamma''(1) = \gamma^2+\zeta(2)$$
To derive $\Gamma''(1)$, use the expansion and take expoential on both sides: $$\ln\Gamma(1+z) = -\gamma z + \sum_{k=2}^\infty \frac{(-1)^k\zeta(k)}{k}z^k$$