Let $(M^n,g)$ be a closed Riemannian manifold with Laplace–Beltrami $\Delta_g$. It is known that, for $s>0$, the resolvent $\mathcal G_s$ of $(-\Delta_g)^s$ has integral kernel, say $G_s$, satisfying $$\left| G_{n/2}(x,y)+ c_n \log d(x,y) \right| \leq C_g$$ for some (unique) universal constant $c_n>0$ and some constant $C_g=C_g(M)$. Is there some reference to the fact that the desingularized kernel $$(x,y)\longmapsto G_{n/2}(x,y)+c_n \log d(x,y)$$ is smooth around $x=y$?
Note: it seems to me that this should follow from some facts in [1], but their arguments are somewhat sketchy, and I would rather like a textbook reference with a solid proof.
[1] S. Minakshisundaram, Å. Pleijel, Some Properties of the Eigenfunctions of The Laplace-Operator on Riemannian Manifolds, Canadian J. Math. 1949