How can we obtain the $-\frac{4\pi}3\mu(x)$ term? Given the expression 
$$K_{ik} := \frac{\partial}{\partial x_k} \int_{\mathscr X} \frac{y_i-x_i}{|y-x|^3} \mu(y) dy,$$
where $\mathscr X=\mathbb R^3$, how does one derive the expression
\begin{align}
K_{kk} &= -\frac{4\pi}3\mu(x) + \int_{y : |y-x| \le d} (\mu(y)-\mu(x))\left(\frac{3(y_k-x_k)^2}{r^5} -\frac 1{r^3} \right) \\
&\qquad+\left(\int_{y : d \le |y-x| \le R} + \int_{y:|y-x| \ge R} \right) \mu(y) \left(\frac{3(y_k-x_k)^2}{r^5} - \frac 1{r^3} \right) dy
\end{align}
My question in particular is how can we obtain the $-\frac{4\pi}3\mu(x)$ term?
My original approach is passing the derivative operator $\frac{\partial}{\partial x_k}$ inside the integral for the expression of $K_{ik}$, carry out the differentiation, substitute $i=k$, and add/subtract $\mu(x)$ in the $|y-x|\le d$ integral.
This comes from pages 344-345 of the journal that contains the paper "Global symmetric solutions of the initial value problem of stellar dynamics" by Jurgen Batt.
 A: $$K_{ik}=-\frac{\partial}{\partial x_i}\int\frac{\partial}{\partial y_k}\frac{1}{|y-x|}\mu(y)dy$$
$$\quad=\int\frac{\partial}{\partial y_i}\frac{\partial}{\partial y_k}\frac{1}{|y-x|}\mu(y)dy$$
now use that
$$\frac{\partial}{\partial y_i}\frac{\partial}{\partial y_k}\frac{1}{|y-x|}=-\frac{4\pi}{3}\delta(y-x)\delta_{ik}-\delta_{ik}|y-x|^{-3}+3(y_i-x_i)(y_k-x_k)|y-x|^{-5}$$
the Dirac delta function appears because of the Laplacian identity $\Delta |y|^{-1}=-4\pi\delta(y)$ (which incidentally holds only in three dimensions)

some more steps, at the request of the OP; note the identity
$$\int_{|y-x|<d}\frac{\partial}{\partial y_i}\frac{\partial}{\partial y_k}\frac{1}{|y-x|}dy=\frac{1}{3}\delta_{ik}\int_{|y|<d}\Delta |y|^{-1}\,dy=-\frac{4\pi}{3}\delta_{ik}$$
then decompose the integrals in the definition of $K$ into
$$K_{ik}=\int_{|y-x|<d}\frac{\partial}{\partial y_i}\frac{\partial}{\partial y_k}\frac{1}{|y-x|}\mu(y)\,dy+\int_{|y-x|>d}\frac{\partial}{\partial y_i}\frac{\partial}{\partial y_k}\frac{1}{|y-x|}\mu(y)\,dy$$
$$\quad=\mu(x)\int_{|y-x|<d}\frac{\partial}{\partial y_i}\frac{\partial}{\partial y_k}\frac{1}{|y-x|}\,dy+\int_{|y-x|<d}\frac{\partial}{\partial y_i}\frac{\partial}{\partial y_k}\frac{1}{|y-x|}[\mu(y)-\mu(x)]\,dy+\int_{|y-x|>d}\frac{\partial}{\partial y_i}\frac{\partial}{\partial y_k}\frac{1}{|y-x|}\mu(y)dy$$
$$\quad=-\frac{4\pi}{3}\delta_{ik}\mu(x)+\int_{|y-x|<d}\frac{\partial}{\partial y_i}\frac{\partial}{\partial y_k}\frac{1}{|y-x|}[\mu(y)-\mu(x)]\,dy+\int_{|y-x|>d}\frac{\partial}{\partial y_i}\frac{\partial}{\partial y_k}\frac{1}{|y-x|}\mu(y)dy$$

higher dimensional generalizations:
the OP asks in a comment for the more general integral
$$K_{ik}^{(n)}=\int_{\mathbb{R}^n}\frac{\partial}{\partial y_i}\frac{\partial}{\partial y_k}\frac{1}{|y-x|^{n-2}}\mu(y)d y$$
where now $x,y\in\mathbb{R}^n$ and $n\geq 3$. This can be evaluated in the same way, using the $n$-dimensional generalization of the Laplacian identity,
$$\Delta_n|r|^{2-n}=-(n-2)S_n\delta_n(r),\;\;n\geq 3$$
with $S_{n}$ the surface area of the $n$-dimensional unit sphere.
A: This is just a further detalization of the Carlo Beenakker's answer. Let us write $$K_{ik}=\int\left (\frac{\partial^2}{\partial r_i \partial r_k}\frac{1}{r}\right )\mu(\vec{y})d\vec{y},$$ where $r_i=y_i-x_i$. But, as Carlo indicated, $$\frac{\partial^2}{\partial r_i \partial r_k}\frac{1}{r}=\frac{3r_ir_k-r^2\delta_{ik}}{r^5}-\frac{4\pi}{3}\,\delta(\vec{r})\delta_{ik}.$$ The derivation of this formula can be found, for example, in http://arxiv.org/abs/1009.2480 (Generalized second-order partial derivatives of 1/r, by V. Hnizdo) and, as Carlo indicated, it is valid in three dimensions. Then $i=k$, we get $$K_{kk}=-\frac{4\pi}{3}\,\mu({\vec{x}})+\int\frac{3(y_k-x_k)^2-r^2}{r^5}\mu(\vec{y}) d\vec{y}.$$
This differs from the expression given in the question by only $$-\mu(\vec{x})\int\limits_{r\le d}\frac{3(y_k-x_k)^2-r^2}{r^5}\,d\vec{y}$$
term. But this last integral equals to zero, because by symmetry $$\int\limits_{r\le d}\frac{r_1^2}{r^5}\,d\vec{y}=\int\limits_{r\le d}\frac{r_2^2}{r^5}\,d\vec{y}=\int\limits_{r\le d}\frac{r_3^2}{r^5}\,d\vec{y}=\frac{1}{3}\int\limits_{r\le d}\frac{r_1^2+r_2^2+r_3^2}{r^5}=\frac{1}{3}\int\limits_{r\le d}\frac{1}{r^3}d\vec{y}.$$
