We have $$ \Gamma(z)=\int_0^\infty x^{z-1}e^{-x}\;dx \tag{1} $$ We also have $$ \frac{d^n}{dz^n}\, x^{z-1}e^{-x}=\log(x)^n e^{-x}x^{z-1}, \;\; z>1 \tag{2} $$ If we create an operator $$ \hat{O}=\sum_{n=0}^\infty \frac{1}{n!}\frac{d^n}{dz^n} = e^{\frac{d}{dz}} $$ and apply it to both sides of $(1)$, we get $$ \hat{O}\Gamma(z)=\hat{O}\int_0^\infty x^{z-1}e^{-x}\;dx $$ if we allow this under the integral sign, then $$ \sum_{n=0}^\infty \frac{1}{n!}\frac{d^n}{dz^n}\Gamma(z)=\int_0^\infty \left(\sum_{n=0}^\infty \frac{1}{n!}\frac{d^n}{dz^n}x^{z-1}e^{-x}\right)\;dx $$ $$ \sum_{n=0}^\infty \frac{1}{n!}\frac{d^n}{dz^n}\Gamma(z)=\int_0^\infty \left(\sum_{n=0}^\infty \frac{\log(x)^n}{n!}\right)x^{z-1}e^{-x}\;dx $$ $$ \sum_{n=0}^\infty \frac{1}{n!}\frac{d^n}{dz^n}\Gamma(z)=\int_0^\infty x^{z}e^{-x}\;dx $$ $$ \sum_{n=0}^\infty \frac{1}{n!}\frac{d^n}{dz^n}\Gamma(z)=\Gamma(z+1) $$ this seems to hold numerically for $z>1$, for $z=1$, the value jumps between $1/2$ and $3/2$, which is curious as the average is right, but then $(2)$ doesn't hold anymore. The $n^{th}$ derivative of $\Gamma(z)$ seems to get increasingly complicated with large $n$.There is always a factor of $\Gamma(z)$ as well so it seems more to be a statement that $$ \frac{1}{\Gamma(z)}\sum_{n=0}^\infty \frac{1}{n!}\frac{d^n}{dz^n}\Gamma(z)=z \tag{3} $$
Can this be proven in another way?
Is there any information about the derivative of $\Gamma(z)$ that makes this identity more obvious?
I notice that if we replace $\frac{d^n}{dz^n}\Gamma(z)$ with $\Gamma(z) \log(z)^n$, then the equation would still hold. Is that related at all?
I've noticed that if we replace $\psi^{(t)}(z)$ as $z^t$ then $$ \frac{1}{\Gamma(z)}\frac{d^n}{dz^n} \Gamma(z) \big|_{\psi^{(t)}(z) \to z^t} = \sum_{k=0}^n \left\{\begin{matrix} n \\ k \end{matrix}\right\}z^{n-k} $$ with $\left\{\begin{matrix} n \\ k \end{matrix}\right\}$ as Stirling numbers of the second kind.
We end up with a kind of ladder operator $$ e^{\frac{d}{dz}}\Gamma(z) = \Gamma(z+1) $$ $$ e^{-\frac{d}{dz}}\Gamma(z) = \Gamma(z-1) $$
Indeed, you have discovered the shift operator. Note that $$ e^{\partial_x}x^a=\sum_{k=0}^\infty\frac{1}{k!}\partial_x x^a=x^a\sum_{k=0}^\infty\frac{(-a)_k}{k!}(-x)^{-k}=x^a{_1F}_0(-a;-;-x^{-1})=(x+1)^a, $$ where $(s)_k$ is the Pochhammer symbol and ${_pF}_q(\mathbf a;\mathbf b;z)$ is the generalized hypergeometric function. Given this identity we can write $t^x$ as a power series to find $$ e^{\partial_x}t^x=\sum_{k=0}^\infty\frac{\log^kt}{k!}e^{\partial_x}x^k=\sum_{k=0}^\infty\frac{\log^kt}{k!}(x+1)^k=t^{x+1} $$ and so, $$ e^{\partial_x}\Gamma(x)=\int_0^\infty e^{\partial_x}t^{x-1}e^{-t}\,\mathrm dt=\int_0^\infty t^xe^{-t}\,\mathrm dt=\Gamma(x+1). $$ With a little work you can easily extend this result to the more general operator shift operator $e^{t\partial_x}$.