How to prove that $\int_0^{\infty} \log^2(x) e^{-kx}dx = \dfrac{\pi^2}{6k} + \dfrac{(\gamma+ \ln(k))^2}{k}$?

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)!} $$

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.