What is the origin of the Dirac delta term in the dipole electric field?
I'll give you the derivation from my book which includes a nice way to see how the delta functions arise: ..............................................................................................................................................................
We can derive the potential field $\vec{A}$ and the electromagnetic fields $\vec{E}$ and $\vec{B}$ of a vector point dipole and an axial point dipole in the same way as we derive these in the case of a the point monopole. We start with a static point charge $\delta(\vec{r})$ and derive the dipole charge/current densities with the help of differential operators.
We apply the same differential operators for $\vec{A}$, $\vec{E}$ and $\vec{B}$ to obtain the dipole fields from the monopole fields. First we recall the fields of the monopole. The point charge obtains (over time) a potential field given by.
$ \mbox{field}\Big\{\,\delta(\vec{r})\,\Big\} ~~=~~ \frac{1}{4\pi r} $
The reversed operator which derives the source from the field is just the Laplacian operator. $ -\nabla^2\Big\{\,\frac{1}{4\pi r}\,\Big\} ~~=~~ \delta(\vec{r}) $
Integrating the delta function over space shows equal contributions from the three spatial components.
$ \int \delta(\vec{r})~d\vec{r} ~=~ -\int\left[\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}\right] \frac{d\vec{r}}{4\pi r} ~=~ \frac13+\frac13+\frac13 ~=~ 1 $
We can define vector dipole and axial dipole sources by using differential operators on the monopole $\delta(\vec{r})$ and derive their potential and electromagnetic fields.
The vector dipole is obtained by differentiating the monopole along the direction that we want the dipole to have. The result is a combination of a positive and a negative delta function. The axial dipole is defined by the curl of the monopole so that it gives a circular point current in the direction of the dipole.
charge/current densities, potentials and fields of dipoles:
$ \begin{array}{|lcll|} \hline &&& \\ j^o &=& ~~~\mbox{div}\,\left(~\vec{\mu} \,\delta(\vec{r})~\right) & \mbox{vector dipole charge density} \\ \vec{j} &=& ~~~\mbox{curl}\left(~\vec{\mu}\, \delta(\vec{r})~\right) & \mbox{axial dipole current density} \\ &&& \\ A^o &=& ~~~\mbox{div}\,\left(~\vec{\mu} \,\frac{1}{4\pi r}\,\right) & \mbox{vector dipole electric potential} \\ \vec{A} &=& ~~~\mbox{curl}\left(~\vec{\mu} \,\frac{1}{4\pi r}\,\right) & \mbox{axial dipole magnetic potential} \\ &&& \\ \mathsf{E}&=& -~\mbox{grad}\left(\mbox{div}\,\left(~\vec{\mu} \,\frac{1}{4\pi r}\,\right)\right)& \mbox{vector dipole electric field} \\ \mathsf{B}&=& +~\,\mbox{curl}\left(\mbox{curl} \left(~\vec{\mu}\,\frac{1}{4\pi r}\,\right)\right)& \mbox{axial dipole magnetic field} \\ &&& \\ \hline \end{array} $
The potential fields are obtained by applying the same differential operators on the field $1/4\pi r$ of the monopole rather than on the charge distribution. For a dipole moment in the $z$-direction we obtain in SI-units.
Potential fields of the electric and magnetic dipoles:
$\Phi ~=~ \frac{z}{4\pi\epsilon_o r^3}\,\mu~~~~~~~~ \vec{A} ~=~ \Big( -y ~,~ x~ ~,~ 0 \Big)~ \frac{\mu_o}{4\pi r^3}\,\mu $
The expressions for the electromagnetic fields do implicitly contain delta functions at the center with the right magnitude. These delta functions are easily lost if the derivation isn't careful enough. The E and B fields are related to each other by the standard vector identity.
$\mbox{curl}(\mbox{curl}\vec{X}) ~=~ \mbox{grad}(\mbox{div}\vec{X})-\nabla^2\vec{X} $
We have seen that the last term (the Laplacian) yields $\delta(\vec{r})$. The two therefor differ only at the center by a delta function in the $\vec{\mu}$ direction. We can see this explicitly if we align the dipoles with the z-axis, so that $\vec{\mu}$ has only a z-component, and write out the fields.
The only difference is in the z-components. The total difference between the two is the Laplacian and thus $\delta(\vec{r})$. The vector dipole gets -1/3 while the axial dipole gets +2/3.
$\begin{aligned} \mathsf{E} &= &\frac{1}{\epsilon_o}\bigg[~ \mathbf{\hat{x}}\, \frac{\partial}{\partial x}\frac{\partial}{\partial z} ~+~ \mathbf{\hat{y}}\, \frac{\partial}{\partial y}\frac{\partial}{\partial z} ~+~~~ \mathbf{\hat{z}}\, \frac{\partial^2}{\partial z^2} ~~~~~~~~~~~~~ &\bigg] ~ \frac{\mu}{4\pi r}\\ \\ \mathsf{B} &= &\mu_o\bigg[~ \mathbf{\hat{x}}\,\frac{\partial}{\partial x}\frac{\partial}{\partial z} ~+~ \mathbf{\hat{y}}\,\frac{\partial}{\partial y}\frac{\partial}{\partial z} ~-~ \mathbf{\hat{z}}\left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}\right)~ &\bigg] ~ \frac{\mu}{4\pi r}\\ \end{aligned} $
Which we can see from the third equation with the three 1/3 parts. In the general case, with arbitrary dipole direction $\vec{\mu}$ we get for the electromagnetic fields in vector form.
Electromagnetic dipole fields:
$ \begin{aligned} \mathsf{E} &= &\frac{1}{\epsilon_o}\left(~\frac{3\left(\vec{\mu}\cdot\hat{r}\right)\hat{r}-\vec{\mu}}{4\pi r^3} ~-~\frac13\vec{\mu}\,\delta(\vec{r})~\right) \\ \\ \mathsf{B} &= &\mu_o\left(~\frac{3\left(\vec{\mu}\cdot\hat{r}\right)\hat{r}-\vec{\mu}}{4\pi r^3} ~+~\frac23\vec{\mu}\,\delta(\vec{r})~\right) \end{aligned} $
Alternatively we can look at the fields explicitly expressed in the individual $x$, $y$ and $z$ components which gives. (with the dipole moment in the $z$-direction)
$\begin{aligned} \mathsf{E} &= &\frac{\mu}{4\pi\epsilon_o}\bigg[~ \mathbf{\hat{x}}\, \frac{3xz}{r^5} ~+~ \mathbf{\hat{y}}\, \frac{3yz}{r^5} ~+~~~ \mathbf{\hat{z}}\, \left(\frac{3zz}{r^5}-\frac{1}{r^3} - \frac{4\pi}{3}\delta(r)\right)~ &\bigg] \\ \\ \mathsf{B} &= &\frac{\mu_o\mu}{~4\pi~}\bigg[~ \mathbf{\hat{x}}\, \frac{3xz}{r^5} ~+~ \mathbf{\hat{y}}\, \frac{3yz}{r^5} ~+~~~ \mathbf{\hat{z}}\, \left(\frac{3zz}{r^5}-\frac{1}{r^3} + \frac{8\pi}{3}\delta(r)\right)~ &\bigg] \\ \end{aligned} $
In general the fields decrease with the third order of $r$ compared to $1/r^2$ for the charge which makes it possible not many effects at larger scale. At the other hand, when we go to smaller scales, the magnetic dipole fields become more powerful relative to the charge.
There are two common types of dipole vector fields in physics:
The "divergence-free" dipole field:
$\vec{V}(\vec{r}) \propto \frac{3\left(\vec{\mu}\cdot\hat{r}\right)\hat{r}-\vec{\mu}}{4\pi r^3} + \frac23\vec{\mu}\,\delta(\vec{r})$
The "curl-free" dipole field:
$\vec{V}(\vec{r}) \propto \frac{3\left(\vec{\mu}\cdot\hat{r}\right)\hat{r}-\vec{\mu}}{4\pi r^3} - \frac13\vec{\mu}\,\delta(\vec{r})$
The divergence of the first one is zero everywhere, even at the origin; the curl of the second one is zero everywhere, even at the origin. (Proof outline in the appendix.)
The electric field of an ideal electric dipole has zero curl, so it has to be the second one. The magnetic B field of an ideal magnetic dipole has zero divergence, so it has to be the first one. The magnetic H field of an ideal magnetic dipole has zero curl, so it has to be the second one. And so on!
It's certainly important to include the delta-functions: If you leave it out, Maxwell's equations are violated at the origin! For example, if you leave out the delta-function, then the electric field has nonzero curl, which violates Faraday's law.
A nice example is the field inside a sphere which is uniformly filled with dipoles. It's just proportional to the delta-function term. So the E-field inside a uniformly polarized sphere points antiparallel to the dipole moments, the B-field inside a uniformly magnetized sphere points parallel, the H-field points antiparallel.
Appendix: Proof outline that the "divergence-free dipole field" is really divergence-free everywhere (including the origin), and the "curl-free dipole field" is really curl-free everywhere (including the origin).
This is only tricky because of the origin. Away from the origin, the fields are the same, and they have both zero divergence and zero curl. You can prove this by straightforward calculus and algebra. To include the origin requires a different approach.
One approach -- a kind of physical approach -- is to write the divergence-free dipole field as a limit of finite, smooth, divergence-free fields, and to write the curl-free dipole field as a limit of finite, smooth, curl-free fields. For the first one, you can use the magnetic field of a circular wire loop carrying a current that gets larger and larger while the loop simultaneously gets smaller and smaller, keeping the product of current and area constant. For the second one, you can use the electric field of a pair of equal and opposite charges that have larger and larger magnitude as their separation gets smaller and smaller, keeping the product of charge and separation constant. This approach requires considerable work and care to make sure you're keeping track of exactly what's happening at and around the origin, to find out what multiple of the delta-function it's approaching there.
An alternative approach -- more mathematical -- is to write the divergence-free one as the curl of something, and the curl-free one as the gradient of something. (You can look up the "something" easily: It's the electric dipole potential in the first case, the magnetic dipole vector potential in the second case.) To prove the equalities away from the origin, it's just straightforward calculus. To prove the equalities at the origin, you do an integral transform (e.g., Kelvin-Stokes theorem) to relate what's happening at the origin to the well-behaved stuff that's happening away from the origin.