Systematic procedure for calculating derivatives involving Dirac delta terms As we know, $\nabla \cdot (\hat{r}/r^2) = 4\pi \delta^3(\vec{r})$. This particular instance is easily justified by appealing to the divergence theorem, but there are also more complicated examples such as this one. My question is actually inspired by a post over at physics SE (link) where it seems that a delta function term has been left unaccounted for somehow, but I'm not sure where.
In general, derivatives are easy to calculate systematically using the chain rule, product rule, etc. but it seems that when the result contains Dirac delta terms, some ad hoc method is used to derive them. Is there a systematic way to do it?
 A: I don't know what you mean specifically by a "systematic way to do it."  But to begin, the Dirac Delta and its derivatives are Generalized Functions, also known as Distributions.
Distributions are linear Functionals that map test functions from $C_C^\infty$ onto the reals.  For the Dirac Delta, the functional definition is given as
$$\langle \delta_a, \phi \rangle =\phi(a) \tag 1$$
where $\phi$ is a suitable test function.

DERIVATIVE OF A DISTRIBUTION:
Suppose that $T$ is a distribution.  We define the $n$'th derivative $T^{(n)}$ of $T$ such that for any $\phi\in C_C^\infty$, we have
$$\langle T^{(n)},\phi \rangle =(-1)^n\langle T,\phi^{(n)}\rangle \tag2$$
It is important to understand the difference between a function and its distributional counterpart.  And it is, therefore, important to understand that differentiation of a function is not the same as differentiation of a distribution.

We can apply "systematic" rules such as the product rule for differention to distributions.  Suppose that $T_1$ and $T_2$ are distributions and that the product $T_1T_2$ is a well-defined distribution (Note that the product of distributions might not be well-defined.  For example, the product of the Dirac Delta and the Heaviside function is not defined.).  Then, we have
$$\begin{align}
\langle (T_1T_2)',\phi\rangle&=-\langle T_1T_2,\phi'\rangle\\\\
&=-\langle T_1, T_2\phi'\rangle\\\\
&=-\langle T_1, (T_2\phi)'-T_2'\phi\rangle\\\\
&=\langle T_1',T_2\phi\rangle+\langle T_1,T_2'\phi\rangle\\\\
&=\langle T_1'T_2,\phi\rangle+\langle T_1T_2',\phi\rangle\\\\
&=\langle T_1T_2'+T_1T_2',\phi\rangle
\end{align}$$
Hence, in distribution, we assert that $(T_1T_2)'=T_1T_2'+T_1T_2'$.

In the following examples, we will systematically apply $(2)$ to determine the distributional interpretations of the (i) derivative of the Heaviside function $H$ and (ii) second mixed partials of the scalar potential function $T=\frac1{|\vec x|}$.


EXAMPLE $(1)$:
Let $H(t)$ denote the Heaviside (or unit step) function.  The derivative of $H(t)$ is $0$ for $t\ne 0$ and fails to exist at $0$.
As a distribution, the derivative, $H'$ of $H$ can be evaluated as follows.  Let $\phi \in C_C^\infty$.  Then, applying $(2)$ we have
$$\begin{align}
\langle H',\phi\rangle&=-\langle H, \phi'\rangle\\\\
&=-\int_{-\infty}^\infty H(t)\phi'(t)\,dt\\\\
&=-\int_0^\infty \phi'(t)\,dt\\\\
&=\phi(0)\tag3
\end{align}$$
whence we see from $(3)$ that in distribution $H'=\delta_0$


EXAMPLE $2$:
As an example, define the function $T$ as $T(|\vec x|)=\frac1{|\vec x|}$ for $|\vec x|\ne 0$.  The second mixed partial derivatives for $|\vec x|\ne 0$ are given by
$$\frac{\partial^2}{\partial x_i \partial x_j}\left(\frac1{|\vec x|}\right)=\frac{3x_ix_j-|\vec x|^2\delta_{ij}}{|\vec x|^5}$$
where $\delta_{ij}$ is the Kronecker Delta.  Irrespective of how $T(0)$ is defined (or not defined), the second mixed partials of $T(|\vec x|)$ fail to exist at $0$.
Now, let $T(|\vec x|)$ be a distribution.  Using $(2)$, we see that the second partial derivatives of $T$ are given by
$$\begin{align}
\langle \partial_{ij}T, \phi\rangle &=\langle T, \partial_{ij}\phi\rangle\\\\
&=\int_{\mathbb{R}^3}\frac1{|\vec x|} \frac{\partial^2 \phi}{\partial x_i \partial x_j}\,d^3\vec x\\\\
&=-\int_{\mathbb{R}^3}\frac{\partial}{\partial x_i}\left( \frac1{|\vec x|}\right)\frac{\partial \phi}{ \partial x_j}\,d^3\vec x\\\\
&=-\lim_{\varepsilon\to 0^+}\int_{\mathbb{R}^3\setminus B_\varepsilon(0)}\frac{\partial}{\partial x_i}\left( \frac1{|\vec x|} \right)\frac{\partial \phi}{ \partial x_j}\,d^3\vec x\tag4
\end{align}$$
where $B_\varepsilon(0)$ is a sphere of radius $\varepsilon$ centered at $0$.

Then integrating by parts the integral on the right-hand side of $(4)$ reveals
$$\begin{align}
\langle \partial_{ij}T, \phi\rangle &=-\phi(0)\int_0^{2\pi}\int_0^\pi n_i n_j \,\sin(\theta)d\theta\,d\phi+\lim_{\varepsilon\to 0^+}\int_{\mathbb{R}^3\setminus B_\varepsilon(0)}\frac{\partial^2}{\partial x_i\partial x_j}\left(\frac1{|\vec x|}\right)\phi\,d^3\vec x\\\\
&=-\frac{4\pi}{3}\phi(0)\delta_{ij}+\lim_{\varepsilon\to 0^+}\int_{\mathbb{R}^3\setminus B_\varepsilon(0)}\frac{3x_ix_j-|\vec x|^2\delta_{ij}}{|\vec x|^5}\phi\,d^3\vec x\tag5
\end{align}$$
where $n_i$ is the $i$'th component of the outer unit normal to $B_\varepsilon(0)$.

Therefore, we see from $(5)$ and $(1)$ that in distribution we have
$$\frac{\partial^2}{\partial x_i \partial x_j}\left(\frac1r\right)=-\frac{4\pi}{3}\delta_{ij}\delta_0(\vec x)+\text{PV}\left(\frac{3x_ix_j-|x|^2\delta_{ij}}{|\vec x|^5}\right)$$

SPECIAL CASES:
For $i=j$, we see that the second derivative of $\frac1{|\vec x|}$ with respect to $x$ is
$$\frac{\partial^2}{\partial x^2}\left(\frac1{|\vec x|}\right)=-\frac{4\pi}{3}\delta_0(\vec x)+\text{PV}\left(\frac{2x^2-y^2-z^2}{|\vec x|^5}\right)$$
And the Laplacian is $-4\pi \delta_0(\vec x)$ as expected.
