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}.\tag{1}$$
However, this was derived from $$\int d^3r' \rho(\vec{r'})-\vec{r'}\nabla\frac{1}{r},\tag{2}$$ 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}).\tag{3}$$
The first term can be derived by calculating $$\vec{E}=-\nabla\Phi\tag{4}$$ 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?