Show that in $\mathbb R^3$ the function $F(x)=\frac{-1}{4\pi|x|}e^{-|x|}$ is a fundamental solution of the operator $\Delta-I$ 
Show that in $\mathbb R^3$ the function $F(x)=\frac{-1}{4\pi|x|}e^{-|x|}$ is a fundamental solution of the operator $\Delta-I$

By the way,the function $F$ is the Yukawa potential.What I try is that if $F$ is any fundamental solution, taking Fourier transform about $(\Delta-I)F=\delta$ plus noting that the characteristic polynomial of the operator is $-(1+4\pi^2|\xi|^2)$ I can get $$F=\int _{\mathbb R^d}\frac{1}{-(1+4\pi^2|\xi|^2)}e^{2\pi ix\xi}d\xi=-\frac{1}{2}e^{-|x|}$$ in the sense of distribution.So in $\mathbb R^3$ I think $-\frac{1}{2}e^{-|x|}$ may be a fundamental solution which is different from the given function $\frac{-1}{4\pi|x|}e^{-|x|}$ (Is what I calculate right?this incidentally show that in general the fundamental solution is not unique?).But how to show $\frac{-1}{4\pi|x|}e^{-|x|}$ is also a fundamental solution?The hint says to use two identity $$\int_{|\xi|=1}e^{2\pi i\xi x}d\sigma(\xi)=\frac{2\sin(2\pi|x|)}{|x|},\hat{Q_y}(\xi)=e^{-2\pi y|\xi|}\frac{sign(\xi)}{i}$$ where $Q$ is the conjugate Possion kernel.But I can't apply them to solve the problem.Can anyone give me a helping hand,thank you
 A: Let $\displaystyle G(r)=-\frac{e^{-r}}{4\pi r}$. Then, for $r\ne 0$, the gradient and Laplacian of $G$ are given by
$$\begin{align}
\nabla G&=\frac{e^{-r}(1+r)}{4\pi r^2}\tag1\\\\
\nabla^2 G&=G\tag2
\end{align}$$
Now,  for any $\phi\in \mathbb{S}$ we have in distribution
$$\begin{align}
\langle \nabla^2G,\phi\rangle&=\langle G,\nabla^2\phi \rangle \tag3\\\\
&=\int_{\mathbb{R^3}} G(r) \nabla^2 \phi(\vec r)\,d^3\vec r \tag4\\\\
&=- \int_{\mathbb{R^3}} \nabla G(r)\cdot \nabla \phi(\vec r)\,d^3\vec r \tag5 \\\\
&= -\lim_{\varepsilon\to 0}\int_{\mathbb{R^3}\setminus B(0,\varepsilon)} \left(\nabla \cdot \left(\phi(\vec r)\nabla G(r) \right)-G(r)\phi(\vec r )\right)\,d^3\vec r \tag6\\\\
&=\lim_{\varepsilon\to0}\frac1{4\pi}\int_0^{2\pi}\int_0^{\pi} \phi(\varepsilon)e^{-\varepsilon}(1+\varepsilon)\,\sin(\theta)\,d\theta\,d\phi+\int_{\mathbb{R^3}}G(r)\phi(\vec r)\,d^3\vec r\tag7\\\\
&=\phi(0)+\int_{\mathbb{R^3}}G(r)\phi(\vec r)\,d^3\vec r \tag8
\end{align}$$
Using $(8)$ and $(2)$, it is easy to see that
$$\langle (\nabla^2-I)G, \phi \rangle =\phi(0)$$
from which we deduce that in distribution
$$(\nabla^2-I)G=\delta$$
where $\delta$ is the Dirac Delta distributon.


NOTES:
We applied the definition of the distributional derivative to arrive at the right-hand side of $(3)$.
In going from $(4)$ to $(5)$, we used the product rule, $\nabla \cdot (\psi\vec A)=\nabla \psi \cdot \vec A+\psi \nabla\cdot \vec A$ with $\psi =G$ and $\vec A=\nabla \phi$.  Then, we applied the Divergence Theorem.
In going from $(5)$ to $(6)$, we applied the product rule again, appealed to $(2)$, and wrote the integral as the limit of the integral that excludes the spherical regionn, $|\vec r|<\varepsilon$ ($B(0,\varepsilon)$).
In going from $(6)$ to $(7)$ we applied the Divergence Theorem and noted that the outer unit normal to $B(0,  \varepsilon)$ is $-\frac{\vec r}{|\vec r|}$.
