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-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!