OP's differentiation formulas can of course be understood pointwise on $\mathbb{R}^3\backslash\{0\}$ where the functions are smooth. The interesting non-trivial question is whether they can be promoted to distributions on the full space $\mathbb{R}^3$? Well, let's see.
We regularize $1/r$ as a smooth function
$$ u_{\varepsilon}(r)~:=~\frac{1}{(r^2+\varepsilon)^{1/2}}
~\rightarrow~ {\rm P.V.}\frac{1}{r}
\quad\text{for}\quad\varepsilon\to 0^+ \tag{A}$$
in $C^{\infty}(\mathbb{R}^3)$, in the sense of generalized functions. Then the derivatives are well-defined:
$$ \frac{\partial u_{\varepsilon}(r)}{\partial x_i}~=~-\frac{x_i}{(r^2+\varepsilon)^{3/2}},\tag{B} $$
$$ \frac{\partial^2 u_{\varepsilon}(r)}{\partial x_i\partial x_j}~=~3\frac{x_ix_j}{(r^2+\varepsilon)^{5/2}}-\frac{\delta_{ij}}{(r^2+\varepsilon)^{3/2}}~\rightarrow~ {\rm P.V.}\left(\frac{3x_ix_j}{r^5} -\frac{\delta_{ij} }{r^3}\right)
\quad\text{for}\quad\varepsilon\to 0^+, \tag{C} $$
$$\nabla^2u_{\varepsilon}(r)
~=~-\frac{3\varepsilon}{(r^2+\varepsilon)^{5/2}}~\rightarrow~ -4\pi\delta^3({\bf r})
\quad\text{for}\quad\varepsilon\to 0^+. \tag{D}$$
In order to make sense of eq. (C) [which OP is inquiring about] we apparently need the principal value distributions
$${\rm P.V.} \frac{1}{r^p}[f]~:=~\lim_{\varepsilon\to 0^+}\int_{\mathbb{R}^3} \mathrm{d}^3{\bf r}\frac{f({\bf r})}{(r^2+\varepsilon)^{p/2}}, \qquad p~\leq~3,\tag{E}$$
$${\rm P.V.} \frac{x_ix_j}{r^p}[f]~:=~\lim_{\varepsilon\to 0^+} \int_{\mathbb{R}^3} \mathrm{d}^3{\bf r}\frac{x_ix_jf({\bf r})}{(r^2+\varepsilon)^{p/2}}, \qquad p~\leq~5.\tag{F}$$
On one hand, eqs. (E) & (F) do not make sense for smooth test functions $f\in C^{\infty}_c(\mathbb{R}^3)$ with compact support but they do make sense if the test functions $f$ are restricted to vanish $f({\bf 0})=0$ at the origin ${\bf r}={\bf 0}$, because then the singularity is removable. On the other hand, applying this restriction $f({\bf 0})=0$, we are not able to detect Dirac delta contributions in eq. (C), which seems to be OP's main motivation to start with.
This issue does not affect eq. (D), which is a well-known representation for the 3D Dirac delta distribution.