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