Variation of the Green function with respect to the metric Consider a (closed) Riemann surface and let $G(x,y)$ be the Green function of the Laplace-Beltrami operator. We can informally identify $G$ with the two-point correlation function for the Gaussian random field:
$$G(x,y)=\left<\phi(x)\phi(y)\right>=\frac{1}{Z} \int \mathcal{D}\phi\;\phi(x)\phi(y) \exp\left(-\frac{1}{2}\int \left| \nabla\phi(z)\right|^2 dV_g(z)\right).$$
This is easy to vary with respect to the metric $g_{\mu\nu}$:
$$\delta G(x,y)=-\frac{1}{2}\int dV_g(z)\;\delta g^{ij}(z)\left<\phi(x)\phi(y) \nabla_i\phi(z)\nabla_j\phi(z)\right>_c,$$
where $\left<\right>_c$ is the connected correlation function (the disconnected diagram $x\leftrightarrow y,z\leftrightarrow z$ is cancelled by the variation of the partition function). We can use Wick's theorem to compute this, getting 
$$\boxed{\frac{\delta G(x,y)}{\delta g^{ij}(z)}=-\nabla_{[i} G(x,z)\nabla_{j]} G(z,y)}$$
(derivatives w.r.t. $z$). This formula looks wickedly similar to the Hadamard variation formula for the variation of the boundary of the domain in flat space. Yet I haven't been able to find any mentions of this in the mathematical literature.
Furthermore, if we define the regularized Green's function (a.k.a. the Robin function) by
$$G^R(x)=\lim_{y\to x} \left(G(x,y)-\frac{1}{2\pi}\ln d(x,y)\right),$$
where $d(x,y)$ is the local geodesic distance, then I'm conjecturing the following variational formula:
$$\frac{\delta G^R(x)}{\delta g^{ij}(z)}=-\nabla_{i} G(x,z)\nabla_{j} G(x,z)-\frac{1}{4\pi}\nabla_i \nabla_j G(x,z).$$
The second term is motivated by the well-known formula for conformal variations (where it is $\frac{1}{4\pi}\delta_x(z)$) and seems to be necessary to cancel the second order pole in this variation. Edit: this guess turned out to be wrong, see my answer below.
 A: It seems that the naive derivation from the path integral only picks up the term coming from quasiconformal variations. Combining the known result for the quasiconformal variation with the much more well known conformal variation, I arrived at the expression
$$\frac{\delta G(x,y)}{\delta g^{\mu\nu}(z)}=(-\nabla_{(\mu}G(z,x)\nabla_{\nu)}G(z,y)+\frac{1}{2}g_{\mu\nu}(z)\nabla_\rho G(z,x)\nabla^\rho G(z,y))+\frac{1}{2V}(G(x,z)+G(z,y))g_{\mu\nu}(z).$$
The traceless part in the brackets represents the quasiconformal variation (which in complex conformal coordinates reduces to just $8\partial_z G(z,x)\partial_z G(z,y)$, as found in the literature). The last "trace" term is the conformal variation, which is not probed by the path integral because it only shows up at finite volume.
I believe this formula to be correct based on the independence of conformal and quasiconformal variations. The resulting tensor is also divergence-free away from $x$ and $y$, as required by general covariance. The more difficult part is proving the following conjecture for the Robin function:
$$\frac{\delta G^{R}\left(x\right)}{\delta g^{\mu\nu}\left(z\right)}=-\nabla_{(\mu}G\left(z,x\right)\nabla_{\nu)}G^R\left(z\right)+\frac{1}{2}g_{\mu\nu}(z)\nabla_\rho G(z,x)\nabla^\rho G^R(z)+\frac{1}{V}G\left(z,x\right)g_{\mu\nu}\left(z\right)-\frac{1}{4\pi}\nabla_{\mu}\nabla_{\nu}G\left(x,z\right)-\frac{\delta(x,z)}{4\pi}g_{\mu\nu}(z)+\frac{1}{8\pi V}g_{\mu\nu}(z)+\frac{1}{2V}f_{\mu\nu},$$
where $f_{\mu\nu}$ is an unknown symmetric traceless tensor defined by
$$\nabla^\mu f_{\mu\nu}=\nabla_\nu G^R,$$
as required by the divergence-free condition on the metric variation of $G^R$.
This result is based on the exact formula
$$G^R(x)=\frac{1}{4\pi}\int G\left(x,y\right)R\left(y\right)\mathrm{d}V(y)+\frac{1}{V}\zeta^{R}\left(1\right)-c,$$
where $\zeta^{R}(s)=\zeta(s)-\frac{V}{4\pi}\frac{1}{s-1}$ and $\zeta$ is the spectral zeta function of the positive Laplace-Beltrami operator, $c=\frac{1}{2\pi}\left(\gamma-\ln2\right)$ (reference below).
Steiner, Jean, A geometrical mass and its extremal properties for metrics on $S^2$, Duke Math. J. 129, No. 1, 63-86 (2005). ZBL1144.53055.
