OP's differentiation formulas (including OP's well-known first identity) 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.
If we regularize OP's tensor (5.28) as a smooth function$^1$
$$ T_{ij}~=~\frac{\delta_{ij}}{(r^2+\varepsilon)^{1/2}} + \frac{x_ix_j}{(r^2+\varepsilon)^{3/2}}
~\rightarrow~ {\rm P.V.}\left(\frac{\delta_{ij}}{r} + \frac{x_ix_j}{r^3}\right)
\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 T_{ij}}{\partial x^k}
~=~\frac{\delta_{ik}x_j +\delta_{jk}x_i- \delta_{ij}x_k}{(r^2+\varepsilon)^{3/2}} - 3\frac{x_ix_jx_k}{(r^2+\varepsilon)^{5/2}}, \tag{B}$$
$$ \begin{align}\nabla^2T_{ij}&~=~-\frac{\delta_{ij}}{(r^2+\varepsilon)^{3/2}} + 3\frac{\delta_{ij}r^2-7x_ix_j}{(r^2+\varepsilon)^{5/2}} +15\frac{x_ix_jr^2}{(r^2+\varepsilon)^{7/2}} \cr
&~=~\delta_{ij}\left(\frac{2}{(r^2+\varepsilon)^{3/2}}
-\color{red}{\frac{3\varepsilon}{(r^2+\varepsilon)^{5/2}}}\right)
- 3x_ix_j\left(\frac{2}{(r^2+\varepsilon)^{5/2}}
+\color{red}{\frac{5\varepsilon}{(r^2+\varepsilon)^{7/2}}}\right)\cr
&~\rightarrow~ \delta_{ij}\left({\rm P.V.}\frac{2}{r^3} - \color{red}{8\pi\delta^3({\bf r})}\right) - {\rm P.V.} \frac{6x_ix_j}{r^5} \quad\text{for}\quad\varepsilon\to 0^+ .\end{align} \tag{C}$$
In eq. (C) the 1st & 3rd terms reproduce the terms in eq. (5.118). The 2nd & 4th term are indeed 3D Dirac delta distributions (of same size). Interestingly, there are no Dirac delta contributions in the off-diagonal sectors $i\!\neq\!j$. $\Box$
--
$^1$ Let us for simplicity drop an overall trivial multiplicative factor $\frac{1}{8\pi\eta_0}$. Here the principal value distributions
$${\rm P.V.} \frac{1}{r^p}~:=~\lim_{\varepsilon\to 0^+}\frac{1}{(r^2+\varepsilon)^{p/2}}, \qquad p~\leq~3,\tag{D}$$
$${\rm P.V.} \frac{x_ix_j}{r^p}~:=~\lim_{\varepsilon\to 0^+} \frac{x_ix_j}{(r^2+\varepsilon)^{p/2}}, \qquad p~\leq~5,\tag{E}$$
are well-defined for test functions $\in C^{\infty}_c(\mathbb{R}^3)$ that vanish at $r=0$ for the observed powers $p$. See also this related math.SE post.