Are mixed derivatives of $1/r$ a distribution? 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$?
 A: OP's differentiation formulas 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.
We regularize $1/r$ as a smooth function
$$ u_{\varepsilon}(r)~:=~\frac{1}{(r^2+\varepsilon)^{1/2}} 
~\rightarrow~ {\rm P.V.}\frac{1}{r}
\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 u_{\varepsilon}(r)}{\partial x_i}~=~-\frac{x_i}{(r^2+\varepsilon)^{3/2}},\tag{B} $$
$$ \frac{\partial^2 u_{\varepsilon}(r)}{\partial x_i\partial x_j}~=~3\frac{x_ix_j}{(r^2+\varepsilon)^{5/2}}-\frac{\delta_{ij}}{(r^2+\varepsilon)^{3/2}}~\rightarrow~ {\rm P.V.}\left(\frac{3x_ix_j}{r^5} -\frac{\delta_{ij} }{r^3}\right)
\quad\text{for}\quad\varepsilon\to 0^+, \tag{C} $$
$$\nabla^2u_{\varepsilon}(r)
~=~-\frac{3\varepsilon}{(r^2+\varepsilon)^{5/2}}~\rightarrow~ -4\pi\delta^3({\bf r})
\quad\text{for}\quad\varepsilon\to 0^+. \tag{D}$$
In order to make sense of eq. (C) [which OP is inquiring about] we apparently need the principal value distributions 
$${\rm P.V.} \frac{1}{r^p}[f]~:=~\lim_{\varepsilon\to 0^+}\int_{\mathbb{R}^3} \mathrm{d}^3{\bf r}\frac{f({\bf r})}{(r^2+\varepsilon)^{p/2}}, \qquad p~\leq~3,\tag{E}$$
$${\rm P.V.} \frac{x_ix_j}{r^p}[f]~:=~\lim_{\varepsilon\to 0^+} \int_{\mathbb{R}^3} \mathrm{d}^3{\bf r}\frac{x_ix_jf({\bf r})}{(r^2+\varepsilon)^{p/2}}, \qquad p~\leq~5.\tag{F}$$
On one hand, eqs. (E) & (F) do not make sense for smooth test functions $f\in C^{\infty}_c(\mathbb{R}^3)$ with compact support but they do make sense if the test functions $f$ are restricted to vanish $f({\bf 0})=0$ at the origin ${\bf r}={\bf 0}$, because then the singularity is removable. On the other hand, applying this restriction $f({\bf 0})=0$, we are not able to detect Dirac delta contributions in eq. (C), which seems to be OP's main motivation to start with. 
This issue does not affect eq. (D), which is a well-known representation for the 3D Dirac delta distribution.
A: 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*}
