# How to derive the $\frac{4\pi}{3}\vec{p}\delta^3(\vec{r})$ element for the dipole field, from its potential?

This might be a bit more general question about how to figure out what is the appropriate (delta) expression in singular points, but e.g. for the dipole, we can derive its potential by a taylor approiximation:
$$\Phi(\vec{r})=\frac{\vec{p}\cdot\vec{r}}{r^3}.$$
However, this was derived from $$\int d^3r' \rho(\vec{r'})-\vec{r'}\nabla\frac{1}{r},$$ for which the above potential expression is valid only for $$r\neq 0$$.
However, the field is said to be derived from the above potential (which seems to be true only for $$r\neq 0$$), and it is
$$\vec{E}(\vec{r})=\frac{3\hat{r}(\hat{r}\cdot\vec{p})-\vec{p}}{r^3}-\frac{4\pi}{3}\vec{p}\delta^3(\vec{r}).$$
The first term can be derived by calculating $$E=-\nabla\Phi$$ when assuming $$r\neq 0$$, but for the second element it's unclear how to derive it.
As this basically a result of $$\Delta \frac{1}{r}$$ I can see why $$4\pi\delta(\vec{r})$$ is involved, but how does one derive this exact term? and why how come the field is a sum of a delta function at 0, and a function the isn't defined/converges at 0?

Hint: Formally one should introduce testfunctions to deal with distributions. Another more physical approach is to regularize the dipole potential $$\Phi_{\varepsilon}~=~ \frac{\vec{p}\cdot\vec{r}}{(r^2+\varepsilon)^{3/2}}, \tag{1}$$ similar to my Phys.SE answer here. The regularized dipole potential $$\Phi_{\varepsilon}\in C^{\infty}(\mathbb{R}^3)$$ is infinitely many times differentiable. The regularized electric field then becomes: $$\vec{E}_{\varepsilon}~=~-\vec{\nabla}\Phi_{\varepsilon} ~=~\frac{3(\vec{p}\cdot\vec{r})\vec{r}-r^2\vec{p} }{(r^2+\varepsilon)^{5/2}} - \vec{p}\frac{\varepsilon}{(r^2+\varepsilon)^{5/2}}. \tag{2}$$ Next it is straightforward to check that the last term is a regularized 3D Dirac delta distribution $$\frac{\varepsilon}{(r^2+\varepsilon)^{5/2}}~\to~ \frac{4\pi}{3}\delta^3(\vec{r}) \quad\text{for}\quad\varepsilon\to 0^+. \tag{3}$$