The proof relies on a representation formula for the Green function
on compact manifolds, given in [1]}.
For \((x, \, y)\in M\times M\)  with \(x\ne y\) , we define
\( G_0(x, \, y) := -\frac{1}{2\pi } \log \operatorname{dist}_g(x, \, y) \, f(\operatorname{dist}_g(x, \, y)),\) 
where \(f\in C^\infty _{\mathrm {c}}[0, \, +\infty )\)  is a cut-off
function, such that \(f(t)=1\)  for \(t\)  sufficiently small
(say, \(t\)  is smaller than half the injectivity radius of \(M\) )
and \(f(t)=0\)  for \(t\)  larger than the injectivity radius of \(M\) .
We define
\( \Gamma _1(x, \, y) := -\Delta _x G_0(x, \, y) + \delta _y\) 
By working in geodesic coordinates centered at \(y\) ,
we obtain (see [1]})
\( \Gamma _1(x, \, y) = \frac{1}{2\pi \operatorname{dist}_g(x, \, y)}\partial _\nu \left(\log \sqrt{\det g(x)}\right)\) 
where \(\partial _\nu \)  denotes the derivative in the radial direction.
The metric tensor in geodesic coordinates satisfies
\(g(x) = \operatorname{Id}+ \mathrm {O}(\operatorname{dist}_g(x, \, y)^2)\) ,
so \(\Gamma _1\)  is bounded. We define
inductively the sequence of functions \(\Gamma _k\)  as
\(\Gamma _{k+1}(x, \, y) := \int _M \Gamma _k(x, \, z) \, \Gamma _1(z, \, y) \, \mathrm {vol}_g(z)\qquad \textrm {for } k\ge 1, \ (x, \, y)\in M\times M, \ x\ne y.\) 
An argument by induction, based on [1]},
shows that \(\Gamma _k\)  extends to a function of class \(C^{k-1}(M\times M)\) ,
for \(k\ge 2\) . Now, \(H\)  can be written in the form
\( H(x, \, y) = G(x, \, y) - G_0(x, \, y)= \sum _{k=1}^2 \int _M \Gamma _k(x, \, z) \, G_0(z, \, y) \, \mathrm {vol}_g(z)+ F(x, \, y),\) 
for some function \(F\)  that satisfies
\( \Delta _x F(x, \, y) = \Gamma _3(x, \, y) + \frac{1}{\mathrm {Vol}(M)}\) 
(see Equations (17) and (18) in [1]}).
By differentiating (REF ) with respect to \(y\) , we obtain
\(\Delta _x \nabla _y F(x, \, y) = \nabla _y \Gamma _3(x, \, y)\) 
Since \(\Gamma _3\)  is of class \(C^2\) , Schauder's estimates
implies that \(\nabla _y F(\cdot , \, y)\in C^{2, \beta }(M)\) 
for any \(\beta \in (0, \, 1)\)  and any \(y\in M\) .
Moreover, by differentiating (REF ) and (REF ), we see that
\( \left| \nabla _x \Gamma _1(x, \, y) \right|_g + \left| \nabla _y G_0(x, \, y) \right|_g\lesssim \operatorname{dist}_g(x, \, y)^{-1}\) 
while on the other hand \(\left| \nabla _x \Gamma _2(x, \, y) \right|_g \le C\) 
because \(\Gamma _2\)  is of class \(C^1\) . Equation (REF ) gives
\( \begin{split}\nabla _x\nabla _y H(x, \, y)= \sum _{k=1}^2 \int _M \nabla _x\Gamma _k(x, \, z) \, \nabla _y G_0(z, \, y)\, \mathrm {vol}_g(z) + \mathrm {O}(1)\end{split}\) 
The estimate (REF ) follows from (REF ) and (REF ) —
see e.g. [1]} for the details.
