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!