If $T_f$ is a distribution, i.e. a linear functional, continuous according to the convergence defined here, defined on the space $K$ of the functions of class $C^\infty$ that are null outside a bounded interval (which is not the same for all functions), its derivative is defined as $$\frac{dT_f}{dx}(\varphi):=-T_f(\varphi')$$where $\varphi'$ is the derivative of $\varphi$. The symbolic writing $T_f(\varphi)=\int_{-\infty}^\infty f(x)\varphi(x)dx$ is often used to write such a functional, since, if $g$ is (Riemann or Lebesgue) integrable on every bounded interval, then $\int_{-\infty}^\infty g(x)\varphi(x)dx$ indeed is such a continuous functional. In this context, we can symbolically define the "derivative" $f'$, for any $T_f$, even if the symbolic writing $f$ does not refer to an integrable function, according to the expression$$\int_{-\infty}^\infty f'(x)\varphi(x)dx:=-\int_{-\infty}^\infty f(x)\varphi'(x)dx=:\frac{dT_f}{dx}(\varphi).$$
Let us come to my question. While studying physics, in particular the theory of electromagnetism and the derivation of the Biot-Savart law from Ampère's law, I always find the equality$$\nabla^2\left(\frac{1}{\|\boldsymbol{x}-\boldsymbol{x}_0\|}\right)=-4\pi\delta(\boldsymbol{x}-\boldsymbol{x}_0)$$where $\nabla^2$ is the Laplacian$^1$. I suppose that, in the tridimensional case, with $\phi:\mathbb{R}^3\to\mathbb{R}$, $\int_{\mathbb{R}^3}\frac{\partial f(\boldsymbol{x})}{\partial x_i}\phi(\boldsymbol{x}) dx_1dx_2dx_3$ is analogously defined as $\frac{\partial T_f}{\partial x_i}(\varphi)$, which I suppose to be analogously defined, in turn, as $-\int_{\mathbb{R}^3}f(\boldsymbol{x})\frac{\partial \phi(\boldsymbol{x})}{\partial x_i} dx_1dx_2dx_3$, although I say I suppose because I have not found a rigourous definition of such derivatives on line nor in cartaceous texts; as to mathematical resources, I have studied Kolmogorov-Fomin's Элементы теории функций и функционального анализа, which only focuses on the monodimensional $\varphi:\mathbb{R}\to\mathbb{R}$ case. Once fixed a proper definition of such derivatives, how can it be proved that $\nabla^2(\|\boldsymbol{x}-\boldsymbol{x}_0\|^{-1})=-4\pi\delta(\boldsymbol{x}-\boldsymbol{x}_0)$?
$^1$ The link contains a derivation of $\nabla^2(\|\boldsymbol{x}-\boldsymbol{x}_0\|^{-1})=-4\pi\delta(\boldsymbol{x}-\boldsymbol{x}_0)$ (equations (18)-(24)), but I do not understand it: I would understand it if we could apply Gauss's divergence theorem at (20), but I know it for functions of class $C^1(\mathring{A})$, $\overline{V}\subset\mathring{A}$, only, while $\nabla\left(\frac{1}{\|\boldsymbol{x}-\boldsymbol{x}_0\|}\right)$ is not even defined for $\boldsymbol{x}=\boldsymbol{x}_0$; the other derivation of the identity $\nabla^2(\|\boldsymbol{x}-\boldsymbol{x}_0\|^{-1})$ $=-4\pi\delta(\boldsymbol{x}-\boldsymbol{x}_0)$ that I have found uses a "weak limit", but it does not use the formal definiton of derivative of a distribution that I have written above. These are the two only references addressing what I am asking that I have managed to find.