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?
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.
