Relationship between Green's function and geodesic distance? I am interested in showing that a certain Green's function can be used to approximate the distance function on a Riemannian manifold in the following sense.  Let $(M,g)$ be a Riemannian manifold and consider a ball $B \subset M$ centered at a distinguished point $p \in M$ whose radius is no larger than the injectivity radius at $p$.  Let $d: B \rightarrow \mathbb{R}$ be the geodesic distance to $p$.  (The reason for considering the ball $B$ instead of the entire domain is simply to avoid issues concerning the cut locus, where the distance function fails to be smooth.)  If $\Delta$ is the negative-(semi)definite Laplace-Beltrami operator on $M$, then there is a 1-parameter family of Green's functions $u_t$ defined as solutions to
$$(\mathrm{id}-t\Delta) u_t = \delta_p$$
where the parameter $t$ is positive and $\delta_p$ is a Dirac delta at $p$.
Question: does $\nabla u_t / |\nabla u_t|$ approach $-\nabla d$ as $t \rightarrow 0$?  Alternatively, do the level sets of $u_t$ approach geodesic circles as $t \rightarrow 0$?
Several closely related results suggest that the answer is likely positive, mostly related to analysis of the heat kernel.  In particular, Varadhan in his classic paper ("On the behavior of the fundamental solution of the heat equation with variable coefficients") considers a similar boundary value problem $(\mathrm{id} - t\Delta)v_t = 0$ on a domain $\Omega$ with $v_t |\partial\Omega = 1$ and shows that $\lim_{t \rightarrow 0} -\sqrt{t}/2 \log v_t = d$, i.e., the function itself converges to the distance function.  However, this result does not (as far as I know) explicitly guarantee convergence of the gradients.  In a similar vein, Malliavin and Stroock ("Short time behavior of the heat kernel and its logarithmic derivatives") essentially show that the gradient of the heat kernel converges to the gradient of the distance function.  However, the heat kernel is a solution to the parabolic problem $\dot{u} = \Delta u$ for some duration $t>0$ with initial conditions $u_0 = \delta_p$ -- i.e., it is not the same as the elliptic problem described above.  I am also aware of some results by Bardi ("An Asymptotic Formula for the Green's Function of an Elliptic Operator"), but they do not seem relevant in this case because they do not consider operators with a constant component ($\mathrm{id}$) as in the case above.  Basically what I'm saying here is that the result is almost certainly true (and not a major departure from what's already mentioned in this paragraph), but I'm having trouble nailing down a concrete reference to cite.
Thanks!
 A: First observe that on a compact  Riemann manifold $(M, g)$ the operator   $1+t \Delta$, $t>0$. $\Delta = d^*d: C^\infty(M)\to C^\infty(M)$ has a unique fundamental solution. Jacques Hadamard  has constructed   very explicit  asymptotic expansions for this fundamental solution which  lead to convergent series in the case of real analytic manifolds and metrics.
A modern description of  Hadamard's  construction can be found in volume  3, Sec. 17.4  of  L. Hörmander's  four volume on linear partial differential operators.   (Hadamard's original memoir is also very useful, but harder to penetrate.) I will give a brief description  of the      fundamental solution $S_r$ of $(r+\Delta)$, $r>0$.
For $\nu=0,1,2,\dotsc $  and $r>0$ denote by  $F_{\nu,r}(x)$ the generalized function (a.k.a. distribution) on $\mathbb{R}^n$   described  as $\newcommand{\ii}{\boldsymbol{i}}$ a Fourier transform of a temperate distribution.
$$F_{\nu,r}(x)= \nu! (2\pi)^{-n} \int_{\mathbb{R}^n} e^{ \ii\langle x,\xi\rangle} (|\xi|^2+r)^{-\nu-1} d\xi. $$
The function $F_\nu$ can be expressed explicitly in terms of Bessel functions.  Note that
$$F_{\nu,r}(x)= r^{\frac{n-\nu-1}{2}}  F_\nu(\sqrt{r} x),\;\;F_\nu(x):=F_{\nu,r=1}(x). $$
If $$\Delta=-\sum_k\partial^2_{x_k}$$
denotes  the (geometers') Laplacian in $\mathbb{R}^n$, then
$$(r+\Delta)F_{0,r}=\delta_0,\;\;(r+\Delta)F_{\nu,r}=\nu F_{\nu-1, r},\;\;\forall \nu>0. $$
One can show that the  generalized function  $F_\nu$ depends only on  the distance $|x|$.
Going back  to the Riemann manifold  $(M,g)$ we denote by $d: M\times M\to \mathbb{R}$ the geodesic  distance function.
The Green function $G(x,y)$ then  has an asymptotic expansion
$$ G(x,y)\sim \sum_{\nu=0}^\infty  U_\nu(x,y) F_{\nu,r}( \; d(x,y)\;)$$
valid  for $d(x,y)$ sufficiently small, where the  functions $U_\nu(x,y)$ are  explicitly described in the above reference.   If $(M,g)$  is real analytic, then the above  series converges in an appropriate sense.
This asymptotic expansion ought to be enough to   investigate your question.
Update I want to add a   "philosophical" comment. The question you asked  is a special case of the following more general question.
Suppose that $f:\mathbb{R}\to\mathbb{R}$ is a continuous function. For simplicity, let us assume it is also bounded. We can define the bounded symmetric operator $f(\sqrt{\Delta})$ where $\Delta$ is the Laplacian on an $m$-dimensional manifold $M$.  Investigate the  behavior of $f(\varepsilon\Delta)$ as $\varepsilon \to 0$.Your case  corresponds  to $f(x)=(1+x^2)^{-1}$. The  heat equation problems  correspond to $f(x)=e^{-x^2}$.     Suppose that $f$   is a  symbol  of order $k$, where $k$ could  be $-\infty$.   For example $(1+x^2)^{-k}$ is a symbol of order  $-2k$, while $e^{-x^2}$ is a symbol of order $-\infty$.
In any case, when $f$ is a symbol,  then $f(\Delta)$  is a pseudodifferential operator, and as such it has a Schwartz kernel which is a distribution on $M\times M$. $\newcommand{\ve}{\varepsilon}$ Your question is about the  behavior as $\ve\to 0$ of the  Schwartz kernel of $f(\ve \Delta)$ along normal directions to the diagonal of $M\times M$.
If  $f$ is rapidly decaying  at $\infty$, say $f(x) < (1+x^2)^{-m}$, $m=\dim M$,  then the Schwartz kernel of $f(\Delta)$  is given  by a continuous function and  one can  be quite precise about the behavior of the  kernel of $f(\ve \Delta)$. In fact,  the faster the  decay of $f$ at $\infty$, the more accurate one can be about the behavior of  the Schwartz kernel of $f(\ve \Delta)$. The radial symmetry you are talking about is then a simple consequence if $f$ decays faster than  $|x|^{-N}$, $N$ sufficiently large. (I believe that $N>2m$ ought to do it but I don't want to be too firm.) If $f$ has exponential decay at $\infty$ one can be remarkable accurate and recover the radial symmetry you are mentioning.  Your  question  involves  the symbol $(1+x^2)^{-1}$  that isn't decaying fast enough at $\infty$. Translation:  your problem requires a bit of care.
