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\"ormander'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.