Let $\mathbf{\mathcal{G}}$ denote the Oseen-Burgers tensor and $\mathbf{\Pi}$ denote the fundamental pressure vector in $\mathbb{R}^3$, i.e. on components we have $$\mathbf{\mathcal{G}}_{jk}(\mathbf{x}-\mathbf{x_0})=\frac{1}{8\pi}\left\{\frac{\delta_{jk}}{|\mathbf{x}-\mathbf{x_0}|}+\frac{\hat{x}_j\hat{x}_k}{|\mathbf{x}-\mathbf{x_0}|^3}\right\}, \mathbf{\Pi}_k(\mathbf{x}-\mathbf{x_0})=\frac{1}{8\pi}\frac{2\hat{x}_k}{|\mathbf{x}-\mathbf{x_0}|^3},$$ where $$\hat{\mathbf{x}}:=\mathbf{x}-\mathbf{x_0}=(\hat{x}_1,..., \hat{x}_n).$$
Then, the pair $(\mathbf{\mathcal{G}}, \mathbf{\Pi})$ satisfies the following equations: $$\Delta_{\mathbf{x}}\mathbf{\mathcal{G}}(\mathbf{x}-\mathbf{x_0})-\nabla_{\mathbf{x}}\mathbf{\Pi}(\mathbf{x}-\mathbf{x_0})=-\delta_{\mathbf{x_0}}(\mathbf{x})\mathbb{I}, \operatorname{div}_{\mathbf{x}}\mathbf{\mathcal{G}}(\mathbf{x}-\mathbf{x_0})=0, \forall \mathbf{x}\in \mathbb{R}^3.$$
This part I understand and it is all clear, but the following result I do not understand. Let $\mathbf{S}(S_{ijk})$ denote the fundamental stress tensor of the Stokes system in $\mathbb{R}^3$ and let $\mathbf{\Lambda}(\Lambda_{ik})$ denote the fundamental pressure tensor of the Stokes system in $\mathbb{R}^3$, i.e. $$S_{ijk}(\mathbf{x}-\mathbf{x_0}):=-\mathbf{\Pi}_j(\mathbf{x}-\mathbf{x_0})\delta_{ik}+\frac{\partial \mathbf{\mathcal{G}}_{ij}}{\partial x_k}(\mathbf{x}-\mathbf{x_0})+\frac{\partial \mathbf{\mathcal{G}}_{kj}}{\partial x_i}(\mathbf{x}-\mathbf{x_0})=-\frac{3}{4\pi}\frac{\hat{x}_i\hat{x}_j\hat{x}_k}{|\mathbf{x}-\mathbf{x_0}|^5},$$ $$\Lambda_{ik}(\mathbf{x}-\mathbf{x_0}):=-2\frac{\partial \Pi_i(\mathbf{x}-\mathbf{x_0})}{\partial x_k}=-2\frac{\Pi_k(\mathbf{x}-\mathbf{x_0})}{\partial x_i}=\frac{1}{2\pi}\left(-\frac{\delta_{ik}}{|\mathbf{x}-\mathbf{x_0}|^3}+3\frac{\hat{x}_i\hat{x}_k}{|\mathbf{x}-\mathbf{x_0}|^5}\right).$$
Then, I am told that the pair $(S_{ijk}, \Lambda_{ik})$ satisfies the Stokes system in $\mathbb{R}^3$ whenever $\mathbf{x}\ne\mathbf{x_0}$, i.e. $$\begin{cases}
-\Delta_{\mathbf{x_0}}S_{ijk}(\mathbf{x}-\mathbf{x_0})+\frac{\partial \Lambda_{ik}(\mathbf{x}-\mathbf{x_0})}{\partial x_{0;j}}=0, \\
\frac{\partial S_{ijk}(\mathbf{x}-\mathbf{x_0})}{\partial x_{0;j}}=0,
\end{cases}$$
where $\mathbf{x_0}=(x_{0;1}, ..., x_{0:n})$. Of course, the Einstein summation convention is used.
This claim is supposedly obvious using the relations I know about $\mathbf{\mathcal{G}}$ and $\mathbf{\Pi}$, but I cannot see how to obtain it. I tried writing on components the fact that $(\mathbf{\mathcal{G}}, \mathbf{\Pi})$ is a fundamental solution of the Stokes system and using the symmetry of these tensors, but I got nowhere.
WLOG let $\boldsymbol x_0=0$, because we can always just shift our coordinate system later.
The equations read
$$\mathbf G(\boldsymbol x) =\frac{1}{8\pi}\left[\frac{\mathbf I}{|\boldsymbol x|}+\frac{\boldsymbol x\otimes \boldsymbol x}{|\boldsymbol x|^3}\right] \tag{1}$$ $$ \boldsymbol \Pi (\boldsymbol x)=\frac{1}{8\pi}\frac{2\boldsymbol x}{|\boldsymbol x|^3}\tag{2}$$ $$\Delta \mathbf G-\nabla \boldsymbol \Pi=-\mathbf I~\delta^{\text{D}}\tag{3}$$ $$\nabla\cdot\mathbf G=0\tag{4}$$ Where I have used the superscript $\mathrm D$ to distinguish the Dirac delta from the Kronecker delta. From this point I will suppress the function argument $(\boldsymbol x)$.
We now write the fundamental stress and pressure tensors: $$S_{ijk}=-\delta_{ik}\Pi_j+\nabla_kG_{ij}+\nabla_iG_{jk}~~~ \\\left\{=-\frac{3}{4\pi}\frac{x_ix_jx_k}{|\boldsymbol x|^5} \right\}\tag{5A}$$ Note that we could also write this as $$S_{ijk}=-\delta_{ij}\Pi_k+\nabla_jG_{ik}+\nabla_i G_{jk}\tag{5B}$$ We also define $\Lambda$, the fundamental pressure tensor: $$ \Lambda_{ij}=-2\nabla_j\Pi_i \\~~\left\{=\frac{1}{2\pi}\left(-\frac{\delta_{ij}}{|\boldsymbol x|^3}+3\frac{x_ix_j}{|\boldsymbol x|^5}\right)\right\}\tag{6}$$
Let's show that
$$-\Delta S_{ijk}+\nabla_k \Lambda_{ij}=0\tag{7}$$ $$ \text{and}$$ $$ \nabla^kS_{ijk}=0\tag{8}$$
These are the same as the equations you wrote, because of the symmetry of the indices.
Let's start with $(8)$. This one is very easy simply using the above equations, but it is also not too difficult as a direct computation, so I will show it for demonstration purposes. $$\nabla_l S_{ijk}=\frac{-3}{4\pi}\frac{\partial}{\partial x^l}\left(\frac{x_ix_jx_k}{|x|^5}\right) \\ =-\frac{3}{4\pi}\left(\frac{\partial_l(x_ix_jx_k)~|x^5|-\partial_l\big(|x|^5\big)~x_ix_jx_k}{|x|^{10}}\right)$$ Recall that $$\nabla |\boldsymbol x|=\frac{\boldsymbol x}{|\boldsymbol x|}=\implies \partial_l\big(|x|\big)=\frac{x_l}{|x|} \\ \implies \partial_l\big(|x|^p\big)=p|x|^{p-1}\frac{x_l}{|x|}=p|x|^{p-2}x_l$$ We may also calculate $$\partial_l(x_ix_jx_k)=x_ix_j\partial_lx_k+x_k\partial_l(x_ix_j) \\ =x_ix_j\partial_lx_k+x_kx_i\partial_lx_j+x_kx_j\partial_lx_i \\ =x_ix_j\delta_{lk}+x_ix_k\delta_{lj}+x_jx_k\delta_{li}$$ Therefore, $$\nabla_l S_{ijk}=\frac{-3}{4\pi}\left(\frac{(x_ix_j\delta_{lk}+x_ix_k\delta_{lj}+x_jx_k\delta_{li})~|x|^5-5|x|^3x_ix_jx_kx_l}{|x|^{10}}\right)$$ Performing some Einstein summation, $$\nabla^kS_{ijk}=\frac{-3}{4\pi}\left(\frac{\big( x_ix_j\delta^k_k+x_ix_k\delta ^k_j+x_jx_k\delta^k_i \big)|x|^5-5|x|^3 x_ix_jx_kx^k}{|x|^{10}}\right)$$ Of course, $\delta^k_k=3$, and the other deltas have the effect of replacing the $k$ with an $i$ or $j$ instead. Also, since $x^kx_k=|x|^2$, this reads $$\nabla^kS_{ijk}=\frac{-3}{4\pi}\left(\frac{\big( 3x_ix_j+x_ix_j+x_jx_i \big)|x|^5-5|x|^3x_ix_j |x|^2}{|x|^{10}}\right) \\ =\frac{-3}{4\pi}\left(\frac{5x_ix_j|x|^5-5x_ix_j|x|^5}{|x|^{10}}\right)=0$$ You can also show this by observing, from $(5\mathrm A )$, that $$\nabla^kS_{ijk}=-\nabla^k\delta_{ik}\Pi_j+\nabla^k\nabla_kG_{ij}+\nabla^k\nabla_iG_{jk} \\ =\underbrace{-\nabla_i\Pi_j+\Delta G_{ij}}_{0,~\text{by}~(3A)}+\nabla_i(\underbrace{\nabla^kG_{jk}}_{=(\nabla\cdot G)_j=0}) \\ =0.$$ Or, another way, from $(5\mathrm B)$, $$\nabla^kS_{ijk}=-\delta_{ij}\nabla^k\Pi_k +\nabla_i((\nabla\cdot G)_j)+\nabla_j((\nabla\cdot G)_i)=-\delta_{ij}\nabla_k\Pi^k$$ I'll leave it as an exercise for you to show that $\nabla_k\Pi^k=0$, i.e $\Pi$ is divergence-free.
I will now assume $\boldsymbol x\neq 0$, so that we can take the Dirac delta as simply being zero.
We move on now to showing $(7)$. Using our definitions $(5B)$ and $(6)$, we get $$-\Delta S_{ijk}+\nabla_k\Lambda_{ij}=-\delta_{ij}\Delta \Pi_k+\Delta(\nabla_iG_{jk}+\nabla_jG_{ik})-\nabla_k (2\nabla_j\Pi_i)$$ In regards to the first term, see that, using $(3)$, and the fact that the Laplacian and the divergence commute, $$\Delta \boldsymbol \Pi=\nabla\cdot\nabla\boldsymbol \Pi=\nabla \cdot(\Delta\mathbf G)=\Delta(\nabla\cdot\mathbf G)=0$$ This leaves us with $$-\Delta S_{ijk}+\nabla_k\Lambda_{ij}=\nabla_i(\Delta G_{jk})+\nabla_j(\Delta G_{ik})-2\nabla_k\nabla_j \Pi_i $$
Now remember, due to the index order convention of the covariant derivative, i.e that $(\nabla v)_{ij}=\nabla_jv_i$ and NOT $\nabla_iv_j$, we have that $\nabla_i (\Delta G_{jk})=\nabla_i((\nabla \Pi)_{jk})=\nabla_i\nabla_k\Pi_j$ and hence $$-\Delta S_{ijk}+\nabla_k\Lambda_{ij}=\nabla_i\nabla_k \Pi_j+\nabla_j\nabla_k\Pi_i-2\nabla_k\nabla_j\Pi_i \\ =\nabla_i\nabla_k\Pi_j-\nabla_j\nabla_k\Pi_i \\=\nabla_i(\Delta G_{jk})-\nabla_j(\Delta G_{ik})$$ Almost there. Now all there is left to do is to use the symmetry of $\mathbf G$: $$-\Delta S_{ijk}+\nabla_k\Lambda_{ij}=\nabla_i(\Delta G_{jk})-\nabla_j(\Delta G_{ik})=\nabla_i(\Delta G_{kj})-\nabla_j(\Delta G_{ki}) \\ =\nabla_i(\nabla_j\Pi_k)-\nabla_j(\nabla _i\Pi_k) \\ =(\nabla_i\nabla_j-\nabla_j\nabla_i)\Pi_k \\ =0.$$
Hope this helps. Please send me a snippet from your course notes if you can as they look interesting.