$\frac{\mathrm d^2 \log(\Gamma (z))}{\mathrm dz^2} = \sum\limits_{n = 0}^{\infty} \frac{1}{(z+n)^2}$
Solution 1:
Use the hadamard product formula
$\Gamma(z) = \frac{e^{-\gamma z}}{z} \prod_{k=1}^\infty \left( 1 + \frac{z}{k} \right)^{-1} e^{z/k} $
Then, note that
$\frac{d \log(\Gamma(z))} {dz} = \frac{\Gamma'(z)}{\Gamma(z)} $
For an infinite product, there is an easy way to compute this expression. If
$ f(z) = \prod f_n(z)$
then it is not hard to prove that
$ \frac{f'(z)}{f(z)} = \sum \frac{f_n'(z)}{f_n(z)} $
applying this to the Gamma function gives
$ \frac{\Gamma'(z)}{\Gamma(z)} = -\gamma - \frac{1}{z} + \sum_{k=1}^\infty \frac{-1}{k(1 + z/k)} + \frac{1}{k} $
Then we have to take one more derivative to get
$ \frac{d^2 \log(\Gamma(z))} {dz^2} = \frac{1}{z^2} + \sum_{k=1}^\infty \frac{1}{(k + z)^2} = \sum_{n=0}^\infty \frac{1}{(z+n)^2} $
Solution 2:
Let's assume the Gauss formula $$\frac{\Gamma'(a)}{\Gamma(a)}+\gamma=\int_{0}^{1}\frac{1-t^{a-1}}{1-t}dt$$ holds (where $\gamma$ is the Euler–Mascheroni constant). Integrating the identity $$\frac{1-t^{a-1}}{1-t}=\sum_{k=0}^{\infty}(t^{k}-t^{a+k-1})$$ yields the series $$\frac{d\ln \Gamma(a)}{da}+\gamma=\sum_{k=0}^{\infty}\left(\frac{1}{k+1}-\frac{1}{a+k}\right)$$ which converges uniformly on finite intervals $a\in[0,A]$. Now we can differentiate the latter series in $a$ to obtain that $$\frac{d^2\ln\Gamma(a)}{da^2}=\sum_{k=0}^{\infty}\frac{1}{(a+k)^2}.$$ The differentiation is valid since the resulting series converges uniformly for $a\geq 0.$
Derivation of the Gauss formula.
Using the basic properties of the beta function, we get $$\Gamma(b)-B(a,b)=\Gamma(b)-\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}=\frac{b\Gamma(b)(\Gamma(a+b)-\Gamma(a))}{b\Gamma(a+b)}$$ $$=\frac{\Gamma(b+1)}{\Gamma(a+b)}\cdot\frac{\Gamma(a+b)-\Gamma(a)}{b}.$$ Passing to the limit $b\to0$ yields $$\frac{d\ln \Gamma(a)}{da}=\frac{\Gamma'(a)}{\Gamma(a)}=\lim\limits_{b\to 0}(\Gamma(b)-B(a,b)),$$ or $$\frac{\Gamma'(a)}{\Gamma(a)}=\lim\limits_{b\to 0}\int_{0}^{\infty}x^{b-1}\left(e^{-x}-\frac{1}{(1+x)^{a+b}}\right)dx=\int_{0}^{\infty}\left(e^{-x}-\frac{1}{(1+x)^{a}}\right)\frac{dx}{x}.\qquad(1)$$ Identity (1) can be used to define the Euler constant $\gamma$ $$\qquad\qquad\qquad\qquad\qquad\qquad-\gamma:=\frac{\Gamma'(1)}{\Gamma(1)}=\int_{0}^{\infty}\left(e^{-x}-\frac{1}{1+x}\right)\frac{dx}{x}.\qquad\qquad\qquad\qquad\qquad\qquad(2)$$ Subtracting (2) from (1) and using the substitution $t=\frac{1}{1+x}$ we obtain that $$\frac{\Gamma'(a)}{\Gamma(a)}-\frac{\Gamma'(1)}{\Gamma(1)}=\int_{0}^{1}\frac{1-t^{a-1}}{1-t}dt.$$