What's the divergence of Oseen tensor for the Zimm model?

Question

While studying the Zimm model on the Doi & Edwards book "The Theory of Polymer Dynamics", I faced equation 4.41 which states that: $$ \frac{\partial}{\partial\boldsymbol{R}_j}\cdot\boldsymbol{H}_{ij}=0 \ , \ \ \text{with} \ \ \ i,j=1,...,N \ \ \ \text{and} \ \ i≠j, $$ where $\boldsymbol{R}_j$ is the position vector of the $i$-th particle and $\boldsymbol{H}_{ij}$ is the Oseen tensor defined as $$ \boldsymbol{H}_{ij}=\frac{1}{8\pi\eta_s\lvert\boldsymbol{r}_{ij}\rvert}\left[\mathbb{I}+\widehat{\boldsymbol{r}_{ij}}\widehat{\boldsymbol{r}_{ij}}\right] \ \ \ \text{where} \ \ \boldsymbol{r}_{ij}= \boldsymbol{R}_i-\boldsymbol{R}_j \ \ , \ \widehat{\boldsymbol{r}_{ij}}=\frac{\boldsymbol{r}_{ij}}{\lvert\boldsymbol{r}_{ij}\rvert} \ \ \ \text{and} \ \ \mathbb{I} \ \ \text{the} \ \ 3\times3 \ \ \text{identity matrix.} $$ I tried to prove equation 4.41 of the book using mixed coordinates and defining the Oseen tensor in the following way using index notation $$ H_{ij}=\frac{1}{8\pi\eta_sr}\left[\delta_{ij}+\frac{x_ix_j}{r^2}\right] $$ and taking the derivatives differentiating sequentially and applying index rules $$ \frac{\partial{H_{ij}}}{\partial{x_k}}=\frac{\delta_{ik}x_j+\delta_{jk}x_i-\delta_{ij}x_k}{r^3}-3\frac{x_ix_jx_k}{r^5} $$ (neglecting the constant factors in which appears the viscosity of the fluid $\eta_s$) but this never becomes zero even if I put $k=j$. Could you tell me if I'm reasoning in the wrong way and if there is a more formal way to take the divergence of a dyadic tensor such as the Oseen tensor?