Evaluate $\int_{1}^{\infty}e^{-x}\ln^{2}\left(x\right)dx$

Evaluate :$$\int_{1}^{\infty}e^{-x}\ln^{2}\left(x\right)\mathrm{d}x$$

I've tried to solve this with some elegant substitutions like $t=e^x$ or $t=\ln\left(x\right)$ . I've also tried to integrate by parts without any success. any help would be good.

Starting the integral from $1$ is in some sense unnatural here, and as a result there is no closed form solution outside of those involving hypergeometric functions or the incomplete Gamma function. The integral should start from $0$, as starting from $0$ usually goes hand in hand with integrating $e^{-x}$ against another function, as the additive character $e^{-x}$ lives on the additive group of the reals $[0,\infty)$, not $[1,\infty)$. This is related to how the incomplete Gamma Function often has no nice closed form, and is more difficult to work with than the Gamma Function itself, and that for the Gamma function, we should start from $0$. It is not by chance that the Laplace Transform starts from $0$ instead of $1$, and we can view $\int_0^\infty e^{-x} \log^2x dx$ as the value of the Laplace transform of $\log^2x$ evaluated at a point. In particular the integral when started from $0$ behaves nicely, and we can prove that $$\int_0^\infty e^{-x}\log^2 x dx=\gamma^2+\frac{\pi^2}{6},$$ which I will now show.

Proof: Consider the exponential generating series $$F(z)=\sum_{k=0}^{\infty}\frac{z^{k}}{k!}\int_{0}^{\infty}e^{-x}\left(\log x\right)^{k}dx.$$ Rearranging, this equals $\int_{0}^{\infty}e^{-x}x^{z}dx=\Gamma(z+1),$ and so $$\int_{0}^{\infty}e^{-x}\left(\log x\right)^{2}dx=\frac{d^{2}}{dz^{2}}\Gamma(z+1)\biggr|_{z=0}=\Gamma''(1).$$ Taking the logarithm of the Weierstrass product for $\Gamma(z+1),$ and differentiating, we have that

$$\frac{\Gamma'(z+1)}{\Gamma(z+1)}=-\gamma+\sum_{n=1}^{\infty}\frac{z}{n(n+z)}.$$ Upon differentiating both sides of the above equation, we obtain

$$\frac{\Gamma''(z+1)}{\Gamma(z+1)}-\frac{\left(\Gamma'(z+1)\right)^{2}}{\Gamma^{2}(z+1)}=\sum_{n=1}^{\infty}\frac{1}{(n+z)^{2}},$$ and so $$\Gamma''(1)=\Gamma'(1)^{2}+\sum_{n=1}^{\infty}\frac{1}{n^{2}}=\gamma^{2}+\frac{\pi^{2}}{6},$$ which proves the result.

Let's define (for $\Re(a)>-1$) the function $$ f(a):=\int_{1}^{\infty}x^a\,e^{-x}\left(x\right)\,\mathrm{d}x$$ then this is the incomplete gamma function $\,f(a)=\Gamma(a+1,1)$.

But since $\,\displaystyle x^a=e^{a\ln(x)}\,$ we have : $$f''(a)=\int_{1}^{\infty}x^a\;\ln^{2}(x)\;e^{-x}\,\mathrm{d}x$$

So that your function is the second derivative of $\,\Gamma(a+1,1)\,$ at $\,a=0$ (so that $x^a$ disappears) or the second derivative relatively to $a$ of $\,\Gamma(a,1)\,$ at $\,a=1$ : $$f''(0)=\lim_{a\to 1}\Gamma(a,1)''$$ (this answer was also given by Mr.G in the comments)
Let's add that the derivatives are not much simpler so that it is probably better to keep it this way...

Of course as noted by Eric Naslund (+1) the answer is simpler when the lower bound is $0$ because in this case we get simply $\;\lim_{a\to 0}\Gamma(a+1,0)''=\Gamma''(1)$.