Divergence of electric field of point dipole

As @Ján Lalinský mentions, the addition of the $delta$ contribution comes out when discussing a volume integral of the field. A finer point that I think is worth mentioning is that the addition of the $delta$ term comes as relevant once we set the volume integral of the singular non-$delta$ term to $0$ (you can see this in the discussion at the end of section $4.1$ in Jackson's Classical Electrodynamics book). As mentioned there, this makes sense when looking at the angular integration, but it is less obvious when one looks at the radial integration.

When asking about the divergence of the field, it is somehow more straightforward to talk about the Laplacian of the point dipole potential, which you can find in this answer. When you compute the electric field via $-nabla Phi$ from that potential, you find only the non-$delta$ term, as long as you ignore all subtleties with the singular point at the location of the dipole. Interestingly enough, when computing $-nabla Phi$, it is possible to somehow motivate the $delta$ term, by splitting up double partial derivatives $partial_i partial_j$ into a trace and a traceless terms, in which the trace brings out a Laplacian of a $frac{1}{r}$. In the end, this is an attempt to fully capture the singular behavior of the expression.

You can also address the issue of the divergence with a tangential computation: by definition of the dipole moment, $ mathbf{p} = int mathbf{x'} rho(mathbf{x'}) d^3x'$. Demanding $ mathbf{p} = epsilon_0 int mathbf{x'} ,nabla'cdot mathbf{E}(mathbf{x'}) d^3x'$ does give you another angle to tackle the problem and test your expressions for the field and charge density. This is another integral that is ultimately related to the volume integral of the field. It also touches on the importance of preserving the contributions from the boundary terms. By imposing this constraint and using the prescription set for the volume integral, you can find some consistency in the end.