I have been interested in the gamma function and its derivatives recently and was surprised to see these constants appear in the values of the derivatives of the gamma function: $$\Gamma'(1) = -\gamma$$ $$\Gamma''(1) = \gamma^2 + \zeta(2)$$ $$\Gamma'''(1) = -2\zeta(3)-\gamma^3-\frac{\gamma\pi^2}{2}$$
I have read that $$-\int_0^\infty e^{-t}\ln(t)dt$$ is a definition of $\gamma$, which also happens to be $-\Gamma'(1)$, so this isn't very surprising to see showing up in higher derivatives of $\Gamma(z)$. I am more curious about where the zeta values are coming from. I should clarify that I am currently taking calculus 2 and the basic integration techniques I know don't seem to get me anywhere in solving these by hand, so I got these closed form values from Wolfram Alpha. Maybe there is some more advanced technique which would easily explain where the zetas are coming from?
This is a consequence of applying repeatedly the Leibniz integral rule for differentiation under the integral sign (the case with constant limits of integration) which you might do in Calculus II: $$\begin{aligned} \Gamma'(1)=\left[\dfrac{\mathrm d}{\mathrm dz}\displaystyle\int_0^\infty x^{z-1}e^{-x}\mathrm dx\right|_{z=1}&=\left[\displaystyle\int_0^\infty \dfrac{\partial e^{(z-1)\ln x}}{\partial z}e^{-x}\mathrm dx\right|_{z=1}\\ &=\left[\displaystyle\int_0^\infty \ln xe^{(z-1)\ln x}e^{-x}\mathrm dx\right|_{z=1}\\ &=\displaystyle\int_0^\infty \ln x\ e^{-x}\mathrm dx=-\gamma \end{aligned} $$ $$\begin{aligned} \Gamma''(1)=\left[\dfrac{\mathrm d^2}{\mathrm dz^2}\displaystyle\int_0^\infty x^{z-1}e^{-x}\mathrm dx\right|_{z=1}&=\left[\dfrac{\mathrm d}{\mathrm dz}\displaystyle\int_0^\infty \dfrac{\partial e^{(z-1)\ln x}}{\partial z}e^{-x}\mathrm dx\right|_{z=1}\\ &=\left[\dfrac{\mathrm d}{\mathrm dz}\displaystyle\int_0^\infty \ln xe^{(z-1)\ln x}e^{-x}\mathrm dx\right|_{z=1}\\ &=...\\ &=\left[\displaystyle\int_0^\infty \ln^2 xe^{(z-1)\ln x}e^{-x}\mathrm dx\right|_{z=1}\\ &=\displaystyle\int_0^\infty \ln^2 x\ e^{-x}\mathrm dx=\gamma^2+\zeta(2) \end{aligned} $$ Idem for the third derivative.
In general, $$\Gamma^{(n)}(1)=\left[\dfrac{\mathrm d^n}{\mathrm dz^n}\displaystyle\int_0^\infty x^{z-1}e^{-x}\mathrm dx\right|_{z=1}=\displaystyle\int_0^\infty \ln^n x\ e^{-x}\mathrm dx$$
If you wanted the explicit calculations for each integral there might be some MSE post covering this or some online paper talking about these types of integrals.