I am trying to derive $$\mathbf E(\mathbf r) = \dfrac1{4\pi\varepsilon_0} \left(3(\bf p \cdot \mathbf {\hat {r}} ) \mathbf{\hat{r}} - \bf p \right)$$
From $$V(\mathbf r) =\dfrac{ \bf \hat r\cdot p}{4\pi\varepsilon_0 r^2} $$
Where $\bf p$ is the dipole vector defined by $$\mathbf p = \int_\text{Volume} \mathbf r^\prime {\rho}{(}\mathbf r^\prime{)}{d\tau^\prime}$$
Where $d\tau^\prime$ is the volume element,
Using $\mathbf E = {-\nabla V}$,
$$\mathbf E = -\dfrac1{4\pi\varepsilon_0}\left( \dfrac{\mathbf{ \hat{r}}}{r^2} \times (\nabla\times \mathbf p) + \mathbf p \times \left(\nabla\times \dfrac{\mathbf{ \hat{r}}}{r^2}\right) + (\mathbf p\cdot\nabla)\dfrac{\mathbf{ \hat{r}}}{r^2} + \left(\dfrac{\mathbf{ \hat{r}}}{r^2} \cdot\nabla\right) \mathbf p\right)$$.
Now I don't know how to simplify this expression without using a specific coordinate system. Any hints ?
You don't want to start pulling $\mathbf{p}$ apart: it's just a constant vector. Instead, find $$ -\nabla \left( \frac{\mathbf{r} \cdot \mathbf{p}}{r^3} \right) = -\left(\nabla \frac{1}{r^3} \right) \mathbf{r} \cdot \mathbf{p} - \frac{1}{r^3}\nabla (\mathbf{r} \cdot \mathbf{p}) = +3\frac{\mathbf{r}}{r^5} (\mathbf{r} \cdot \mathbf{p}) - \frac{\mathbf{p}}{r^3} = \frac{1}{r^3}(3 (\hat{\mathbf{r}} \cdot \mathbf{p})\hat{\mathbf{r}}-\mathbf{p}) $$ using $\nabla f(r) = f'(r)\mathbf{r}/r$ and $\partial_i r_j p_j = \delta_{ij} p_j = p_i $.