I've only found two definitions of the charge density of a pure electric dipole:
\begin{equation}
\rho(\vec{r})=\vec{p}\cdot\nabla\delta(\vec{r}-\vec{r}_0)
\end{equation}
and
\begin{equation}
\rho(\vec{r})=-\vec{p}\cdot\nabla\delta(\vec{r}-\vec{r}_0).
\end{equation}
Neither of these came with reasons or derivations, so I'm not really sure which to use?
Best Answer
In electrostatics, when you write the multipole expansion of the potential you find \begin{equation} \Phi(\vec{x})=\frac{1}{4 \pi \epsilon_0} \left[ \frac{q}{r} + \frac{\vec{p} \cdot \vec{x}}{r^3} + ... \right] \, , \end{equation} where $r=|\vec{x}|$ and $...$ indicates the higher order multipole terms. We can find the effective charge distribution with \begin{equation} \nabla^2 \Phi = - \frac{\rho}{\epsilon_0} \, . \end{equation}
Now we just need appropriate expressions for the laplacian of the terms we have. Remember that \begin{equation} \frac{\vec{x}}{r^3} = - \nabla \frac{1}{r} \quad \mathrm{and} \quad \nabla^2 \frac{1}{r} = - 4 \pi \, \delta(\vec{x}) \, . \end{equation} Now we will be able to reproduce the result if we notice that
\begin{eqnarray} \nabla^2 \left[ \frac{\vec{p} \cdot \vec{x}}{r^3}\right] &=& \nabla \cdot \left[ \nabla \left( \frac{\vec{p} \cdot \vec{x}}{r^3}\right) \right] = \nabla \cdot \left[ (\vec{p} \cdot \nabla) \left( \frac{\vec{x}}{r^3}\right) \right] = \vec{p} \cdot \nabla \left[ \nabla \cdot \left(\frac{\vec{x}}{r^3}\right) \right] \\ &=& - \vec{p} \cdot \nabla \left[ \nabla \cdot \left( \nabla \frac{1}{r}\right) \right] = - \vec{p} \cdot \nabla \left[ \nabla^2 \frac{1}{r} \right] = 4 \pi \, \vec{p} \cdot \nabla \delta (\vec{x}) \, . \end{eqnarray}
Check each equality carefully. Hence we find the effective charge density \begin{equation} \rho(\vec{x}) = q \, \delta(\vec{x}) - \vec{p} \cdot \nabla \delta(\vec{x}) \, + \, ... \, . \end{equation}
EDIT:
Alternatively you can say \begin{eqnarray} -4\pi\nabla \delta(\vec{x}) &=&\nabla \left[ \nabla \cdot \left(\nabla \frac{1}{r}\right)\right]=\nabla \left[ \nabla \cdot \left(- \frac{\vec{x}}{r^3}\right)\right] = \nabla \times \left[ \nabla \times \left(- \frac{\vec{x}}{r^3}\right) \right] + \nabla^2 \left(- \frac{\vec{x}}{r^3}\right) \\ &=& - \nabla^2 \left( \frac{\vec{x}}{r^3}\right) \, . \end{eqnarray} Now I think everything is consistent, sorry for the confusion.