So, in three-dimensions we famously have the result that the Laplacian acting on $1/r$ is a distribution: $$\vec{\nabla}^2\frac{1}{4\pi r}=-\delta^3(\vec{r})$$ where $\delta^3(\vec{r})$ is the Dirac-delta function.
My question: how should one think of the mixed derivative $\partial_{i}\partial_{j}\frac{1}{4\pi r}=?$.
Naively, taking derivatives, one gets $$\partial_{i}\partial_{j}\frac{1}{4\pi r}=\frac{1}{4\pi}\left(\frac{3 r_i r_j}{r^5}-\frac{\delta_{ij}}{r^3}\right)$$ but, tracing over indices does not reproduce the $\delta$-function piece, of course. So instead, it naively seems that we should have something like $$\partial_{i}\partial_{j}\frac{1}{4\pi r}\stackrel{?}{=}\frac{1}{4\pi}\left(\frac{3 r_i r_j}{r^5}-\frac{\delta_{ij}}{r^3}\right)-\left(\frac{\delta_{ij}}{3}+c(\delta_{ij}/3-r_ir_j/r^2)\right)\delta^3(\vec{r})$$ which reproduces the original relation for any value of $c$ upon contracting indices. So, is something like the above correct? If so, is there a unique way of fixing $c$?
Probably you should make an approach something like this. Let \begin{equation*} \Phi(x)=\frac{1}{4\pi r}\quad\text{with}\quad r(x)=\left(x_{1}^{2}+x_{2}^{2}+x_{3}^{2}\right)^{\frac{1}{2}}. \end{equation*} We hope to describe $\partial_{i} \partial_{j} \Phi(x)$ as some sort of distribution. A distribution must be integrated against a smooth function, so let's do that: \begin{equation*} \int_{\mathbb{R}^{3}}\partial_{i}\partial_{j}\Phi(x)f(x)dx =\int_{\mathbb{R}^{3}\setminus B_{\epsilon}} \partial_{i}\partial_{j}\Phi(x)f(x)dx + \int_{B_{\epsilon}} \partial_{i}\partial_{j}\Phi(x)f(x)dx. \end{equation*} We are isolating the singularity inside a small ball. I guess you are happy with the first term there (it can be evaluated using the expression you derived) so let's focus on the second. Integrating by parts: \begin{equation*} \int_{B_{\epsilon}} \partial_{i}\partial_{j}\Phi(x)f(x)dx =-\int_{B_{\epsilon}}\partial_{i}\Phi(x) \partial_{j}f(x)dx +\frac{1}{\epsilon}\int_{\partial B_{\epsilon}}\partial_{i}\Phi(x) f(x)x_{j}dS(x) \end{equation*}
The first term integrate by parts again \begin{equation*} -\int_{B_{\epsilon}}\partial_{i}\Phi(x) \partial_{j}f(x)dx= \int_{B_{\epsilon}}\Phi(x) \partial_{i}\partial_{j}f(x)dx -\frac{1}{\epsilon}\int_{\partial B_{\epsilon}}\Phi(x) \partial_{j}f(x)x_{i}dS(x). \end{equation*} But we can ignore all of that as $\epsilon\to 0$ because \begin{equation*} \left\lvert \int_{B_{\epsilon}}\Phi(x) \partial_{i}\partial_{j}f(x)dx\right\rvert \leq \left\lVert \partial_{i}\partial_{j}f \right\rVert_{L^{\infty}} \int_{B_{\epsilon}}\lvert\Phi(x)\rvert dx\leq C \epsilon^{2} \end{equation*} and \begin{equation*} \left\lvert \frac{1}{\epsilon}\int_{\partial B_{\epsilon}}\Phi(x) \partial_{j}f(x) x_{i} dS(x)\right\rvert \leq \left\lVert \partial_{j}f \right\rVert_{L^{\infty}} \int_{\partial B_{\epsilon}}\lvert\Phi(x)\rvert dS(x)\leq C \epsilon \end{equation*} We are left with \begin{equation*} \frac{1}{\epsilon}\int_{\partial B_{\epsilon}}\partial_{i}\Phi(x) f(x)x_{j} dS(x)= -\frac{1}{4\pi\epsilon^{4}}\int_{\partial B_{\epsilon}}x_{i}x_{j}f(x)dS(x). \end{equation*} I had to look it up but it seems that \begin{equation*} \int_{\partial B_{\epsilon}}x_{i}x_{j}dS(x)=\frac{4\pi}{3}\epsilon^{4}\delta_{ij}. \end{equation*} Therefore as $\epsilon\to 0$ \begin{equation*} \int_{B_{\epsilon}} \partial_{i}\partial_{j}\Phi(x)f(x)dx\to -\frac{1}{4\pi\epsilon^{4}}\int_{\partial B_{\epsilon}}x_{i}x_{j}f(x)dS(x)\to -\frac{1}{3}\delta_{ij}f(0). \end{equation*} So in conclusion one might write \begin{equation*} \partial_{i}\partial_{j}\Phi(x)= \begin{cases} \frac{1}{4\pi}\left(\frac{3 x_{i}x_{j}}{r^{5}}-\frac{\delta_{ij}}{r^{3}}\right)&\text{for }x\neq 0\\ -\frac{1}{3}\delta_{ij}\delta(x)&\text{for } x = 0. \end{cases} \end{equation*}