Confusion regarding the solution to $\nabla ^2 \phi(\textbf{r}) = \rho(\textbf{r})$ using Green's function We know we can solve one of the Maxwell's equation using Green's function. More specifically, we can solve
$$\nabla ^2 \phi(\textbf{r}) = \rho(\textbf{r})$$
using
$$\phi(\textbf{r}) = \int d\textbf{r}' G(\textbf{r},\textbf{r}') \rho(\textbf{r})
\qquad\text{where}\qquad
\nabla ^2 G(\textbf{r},\textbf{r}') = \delta (\textbf{r}-\textbf{r}'),$$
and
$$G(\textbf{r},\textbf{r}')=\frac{1}{|\textbf{r}-\textbf{r}'|}$$ follows.
However, I cannot understand how we obtain the delta function $\delta (\textbf{r}-\textbf{r}')$ by having the Laplacian act on $G(\textbf{r},\textbf{r}')$. Any thoughts on how can I obtain that?
 A: Define a vector field
$$
\mathbf{F} 
= \nabla \frac{1}{|\mathbf{r}|}
= -\frac{\mathbf{r}}{|\mathbf{r}|^3}
.
$$
If we take the divergence of this, we see that it vanishes:
$$
\nabla \cdot \mathbf{F}
= - \frac{(\nabla\cdot\mathbf{r})|\mathbf{r}|^3 - \mathbf{r}\cdot3|\mathbf{r}|^2\mathbf{r}/|\mathbf{r}|}{|\mathbf{r}|^6}
= -\frac{3|\mathbf{r}|-3|\mathbf{r}|^3}{|\mathbf{r}|^6}
= 0.
$$
But this calculation is only defined for $\mathbf{r} \neq \mathbf{0}.$ To cover origin we will use the divergence theorem:
$$
\iiint_\Omega \nabla\cdot\mathbf{F} \, dV = \iint_{\partial\Omega} \mathbf{F}\cdot\mathbf{n}\,dS,
$$
where $\Omega$ is some region with a smooth enough boundary $\partial\Omega.$ If $\Omega$ does not contain origin, both sides of the equality vanish. Now take $\Omega=B_r(\mathbf{0}),$ i.e. a ball with radius $r$ and center in origin. Then $\partial\Omega$ is the sphere $S_r(\mathbf{0})$ with radius $r$ and center in origin, and the right hand side becomes
$$
\iint_{\partial\Omega} \mathbf{F}\cdot\mathbf{n}\,dS
= \iint_{S_r(\mathbf{0})} \left(-\frac{\mathbf{r}}{|\mathbf{r}|^3}\right)\cdot\frac{\mathbf{r}}{|\mathbf{r}|} \,|\mathbf{r}|^2 d\omega
= -\iint_{S_r(\mathbf{0})} d\omega
= -4\pi
.
$$
(Here $\omega$ is the solid angle measure.)
Thus,
$$
\iiint_\Omega \nabla\cdot\mathbf{F} \, dV
= \begin{cases}
0, & \text{ if } \mathbf{0} \not\in \Omega \\
-4\pi, & \text{ if } \mathbf{0} \in \Omega \\
\end{cases}
$$
Therefore, $\nabla\cdot\mathbf{F}(\mathbf{r}) = -4\pi\,\delta(\mathbf{r}).$
A: REGULARIZING THE DIRAC DELTA:
As I showed in This Answer, we can show that $\nabla \cdot \left(\frac{\hat r}{r^2}\right)=4\pi \delta (\vec r)$ by using a regularization of the Dirac Delta.  To begin, let $\vec \psi$ be the function given by
$$\vec \psi(\vec r;a)=\frac{\vec r}{(r^2+a^2)^{3/2}} \tag 1$$
where we note that $\psi(\vec r;0)=\frac{\hat r}{r^2}$ for $\vec r\ne0$.
The divergence of $(1)$ is
$$\nabla \cdot \vec \psi(\vec r; a)=\frac{3a^2}{(r^2+a^2)^{5/2}}$$
And as shown in the referenced answer, for any test function $\phi$ we have
$$\begin{align}
\lim_{a \to 0}\int_V \nabla \cdot \vec \psi(\vec r; a)\phi(\vec r)\,dV=\begin{cases}4\pi \phi(0)&, \{0\}\in V\\\\
0&,\{0\}\notin V
\end{cases}
\end{align}$$
and it is in this sense that
$$\lim_{a\to 0} \nabla \cdot \vec \psi(\vec r;a)\sim 4\pi \delta(\vec r)$$
Enforcing the translation $\vec r\mapsto \vec r-\vec r'$ yields the coveted result
$$\bbox[5px,border:2px solid #C0A000]{\lim_{a\to 0} \nabla \cdot \vec \psi(\vec r-\vec r';a)\sim 4\pi \delta(\vec r-\vec r')}$$


CLASSICAL ANALYSIS:

We need not use the Dirac Delta to prove that $\nabla^2\int_{V}\rho(\vec r')G(\vec r,\vec r')\,dV'=\rho(\vec r)$.

For $\rho(\vec r)\in C^\infty_C$ the gradient of $\phi(\vec r)$ can be written
$$\begin{align}
\nabla \int_{V}\rho(\vec r')G(\vec r,\vec r')\,dV'&=\int_{V}\rho(\vec r')\nabla G(\vec r,\vec r')\,dV'\\\\
&=-\int_{V}\rho(\vec r')\nabla' G(\vec r,\vec r')\,dV'\\\\
&=-\oint_{\partial V}\rho(\vec r') G(\vec r,\vec r')\hat n'\,dS'+\int_{V}\nabla' \rho(\vec r')G(\vec r,\vec r')\,dV'\tag2
\end{align}$$
Taking the divergence of $(2)$ reveals
$$\begin{align}
\nabla^2 \int_{V}\rho(\vec r')G(\vec r,\vec r')\,dV'&=\oint_{\partial V}\rho(\vec r') \frac{\partial G(\vec r,\vec r')}{\partial n'}\,dS'-\int_{V}\nabla' \rho(\vec r')\cdot \nabla 'G(\vec r,\vec r')\,dV'\tag3
\end{align}$$
We may write the integrand of the integral on the right-hand side of $(3)$ as
$$\nabla' \rho(\vec r')\cdot \nabla 'G(\vec r,\vec r')=\nabla' \cdot (\rho(\vec r')\nabla' G(\vec r,\vec r'))-\rho(\vec r')\nabla'^2 G(\vec r,\vec r')$$
but cannot apply the Divergence Theorem since $\nabla'G(\vec r,\vec r')$is not continuously differentiable in $V$.  We can work around this issue by proceeding as follows.
We exclude the singularity at $\vec r'=\vec r$ from $V$ with a spherical volume $V_\varepsilon$ with center at $\vec r$ and with radius $\varepsilon$.  Then, using $\nabla'^2 G(\vec r,\vec r')=0$ for $\vec r'\in V-V\varepsilon$, we can write
$$\begin{align}
\int_{V}\nabla' \rho(\vec r')\cdot \nabla 'G(\vec r,\vec r')\,dV'&=\lim_{\varepsilon\to 0^+}\int_{V-V_\varepsilon}\nabla' \rho(\vec r')\cdot \nabla 'G(\vec r,\vec r')\,dV'\\\\
&=\lim_{\varepsilon\to 0}\int_{\partial V+\partial V_\epsilon}\rho(\vec r') \frac{\partial G(\vec r,\vec r')}{\partial n'}\,dS'\\\\
&=\int_{\partial V}\rho(\vec r') \frac{\partial G(\vec r,\vec r')}{\partial n'}\,dS'\\\\
&+\lim_{\varepsilon\to 0^+}\int_0^{2\pi}\int_0^\pi \rho(\vec r')\frac{\vec r-\vec r'}{\varepsilon^3}\cdot \frac{\vec r'-\vec r}{\varepsilon}\,\varepsilon^2\,\sin(\theta)\,d\theta\,d\phi\\\\
&=\int_{\partial V}\rho(\vec r') \frac{\partial G(\vec r,\vec r')}{\partial n'}\,dS'-4\pi \rho(\vec r)\tag4
\end{align}$$
Substituting $(4)$ into $(3)$, we find that
$$\bbox[5px,border:2px solid #C0A000]{\nabla^2 \int_{V}\rho(\vec r')G(\vec r,\vec r')\,dV'=4\pi \rho(\vec r)}\tag5$$


Inasmuch as $(5)$ is true for any test function $\rho(\vec r)$, we see that in the sense of distributions
$$\bbox[5px,border:2px solid #C0A000]{\nabla^2 G(\vec r,\vec r')=4\pi \delta(\vec r-\vec r')}$$
