Relationship between Green's function and geodesic distance?

Question

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!

Answer 1

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.