How to prove that $\int_0^{\infty} \log^2(x) e^{-kx}dx = \dfrac{\pi^2}{6k} + \dfrac{(\gamma+ \ln(k))^2}{k}$? I was answering this question: $\int_0^\infty(\log x)^2(\mathrm{sech}\,x)^2\mathrm dx$ and in my answer, I encountered the integral
$$\int_0^{\infty} \log^2(x) e^{-kx}dx$$
which according to WolframAlpha for $k=1,2,3$ and thereby generalizing gives us

$$\int_0^{\infty} \log^2(x) e^{-kx}dx = \dfrac{\pi^2}{6k} + \dfrac{(\gamma+ \ln(k))^2}{k} \tag{$\star$}$$

Also, in general, is there a nice form for

$$\int_0^{\infty} \log^m(x) e^{-kx} dx \tag{$\perp$}$$

EDIT
For $(\star)$, if we let $kx = t$, we then get
$$I = \int_0^{\infty} \log^2(t/k) e^{-t} \dfrac{dt}k = \dfrac1k\int_0^{\infty}\log^2(t) e^{-t}dt-\dfrac{2\log(k)}{k}\int_0^{\infty}\log(t) e^{-t} dt + \dfrac{\log^2(k)}k$$
Considering the above, $(\star)$ and $(\perp)$ boil down to evaluating the integral

$$\color{red}{I(m) = \int_0^{\infty} \log^m(x) e^{-x} dx} \tag{$\spadesuit$}$$

Some values of $I(m)$.
\begin{align}
I(1) & = -\gamma\\
I(2) & = \zeta(2) + \gamma^2\\
I(3) & = -2\zeta(3) - 3 \gamma \zeta(2) - \gamma^3\\
I(4) & = \dfrac{27}2 \zeta(4) +8 \gamma \zeta(3) + 6 \gamma^2 \zeta(2) + \gamma^4\\
I(5) & = -24 \zeta(5) - \dfrac{135}2 \gamma \zeta(4) - 20 \gamma^2 \zeta(3) - 20\zeta(2) \zeta(3) - 10 \gamma^3 \zeta(2) -\gamma^5
\end{align}
 A: As mentioned in this answer
$$
\int_0^\infty\log(t)\,e^{-t}\,\mathrm{d}t=\Gamma'(1)=-\gamma
$$
Using the definition
$$
\Gamma(x)=\int_0^\infty t^{x-1}\,e^{-t}\,\mathrm{d}t
$$
and taking the derivative $n$ times, we get
$$
\Gamma^{(n)}(1)=\int_0^\infty\log(t)^n\,e^{-t}\,\mathrm{d}t
$$
In the aforementioned answer, we also have
$$
\Gamma'(x+1)=\Gamma(x+1)\left(-\gamma+\sum_{k=1}^\infty\left(\frac1k-\frac1{k+x}\right)\right)\tag{$\ast$}
$$
Taking the derivative of $(\ast)$ at $x=0$ yields
$$
\begin{align}
\Gamma''(1)
&=\Gamma'(1)\left(-\gamma\right)+\Gamma(1)\zeta(2)\\
&=\gamma^2+\zeta(2)
\end{align}
$$
Taking the second derivative of $(\ast)$ at $x=0$ yields
$$
\begin{align}
\Gamma'''(1)
&=\Gamma''(1)(-\gamma)+2\Gamma'(1)\zeta(2)+\Gamma(1)(-2\zeta(3))\\
&=-\gamma^3-3\gamma\zeta(2)-2\zeta(3)
\end{align}
$$
We can use Leibniz rule and
$$
\frac{\mathrm{d}^n}{\mathrm{d}x^n}\left(-\gamma+\sum_{k=1}^\infty\left(\frac1k-\frac1{k+x}\right)\right)=(-1)^{n+1}n!\,\zeta(n+1)\quad\text{ for }n\ge1
$$
applied to $(\ast)$, to get the recursion
$$
\Gamma^{(n)}(1)=-\gamma\,\Gamma^{(n-1)}(1)+(n-1)!\sum_{k=2}^n(-1)^k\zeta(k)\frac{\Gamma^{(n-k)}(1)}{(n-k)!}
$$
A: Another method involving Mellin Transforms.
[This is my first time using this tactic so please correct me if I make a mistake, I am answering this question in part as practice]
Consider $F(s)$ as the Mellin Transform of $f(x)$
$$F(s)=\int_0^\infty x^{s-1} f(x)\text{ d}x$$
And thus
$$F''(s)=\int_0^\infty x^{s-1}\ln^2 x \space f(x) \text{ d}x$$
We assume that $f(x)=e^{-kx}$ and we are left with a nice known identity:
$$F(s)=\int_0^\infty x^{s-1}e^{-kx}=\frac{\Gamma(s)}{k^{s+1}}$$
Thus the answer to our integral is truly:
$$\lim_{s\to 3} \left(\left[\frac{\Gamma(s)}{k^{s-1}}\right]''\right)$$
Unfortunately, I am unable to calculate this limit, but to me, if you wanted to take a route similar to this, it seems much cleaner in the long run.
