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.
You are right, the partial derivatives of distributions on higher-dimensional spaces are defined as you surmised,
$$\frac{\partial T}{\partial x_i} \colon \varphi \mapsto -T\biggl[\frac{\partial \varphi}{\partial x_i}\biggr],$$
more generally
$$D^{\alpha} T \colon \varphi \mapsto (-1)^{\lvert\alpha\rvert} T[D^{\alpha}\varphi]$$
for higher derivatives.
From that you obtain $(\nabla^2 T)[\varphi] = T[\nabla^2\varphi]$, and for the locally integrable function $f(x) = \lVert x-x_0\rVert^{-1}$ we therefore have
$$(\nabla^2 T_f)[\varphi] = T[\nabla^2\varphi] = \int_{\mathbb{R}^3} f(x)\cdot \nabla^2\varphi(x)\,dx_1\,dx_2\,dx_3.\tag{1}$$
Now choose $R > 0$ so large that $\operatorname{supp} \varphi \subset B_R(x_0)$, and choose $0 < \varepsilon < R$. Then
$$\int_{\mathbb{R}^3} f(x)\cdot \nabla^2\varphi(x)\,dx_1\,dx_2\,dx_3 = \lim_{\varepsilon \searrow 0} \int_{\varepsilon < \lVert x-x_0\rVert < R} f(x)\cdot \nabla^2\varphi(x)\,dx_1\,dx_2\,dx_3$$
since $f$ is locally integrable and $\nabla^2\varphi$ is continuous.
Away from $x_0$, $f$ is smooth, and a slightly tedious computation shows that $\nabla^2 f \equiv 0$ on $\mathbb{R}^3\setminus \{x_0\}$, so
$$\int_{\mathbb{R}^3} f(x)\cdot \nabla^2\varphi(x)\,dx_1\,dx_2\,dx_3 = \lim_{\varepsilon \searrow 0} \int_{\varepsilon < \lVert x-x_0\rVert < R} f(x)\cdot \nabla^2\varphi(x) - \varphi(x)\cdot \nabla^2 f(x)\,dx_1\,dx_2\,dx_3,$$
and to that integral you can apply Green's formula to obtain
\begin{align} \int_{\mathbb{R}^3} f(x)\cdot \nabla^2\varphi(x)\,dx_1\,dx_2\,dx_3 &= \lim_{\varepsilon \searrow 0}\; \Biggl(\int_{\lVert x-x_0\rVert = R} f(x)\frac{\partial \varphi}{\partial \nu}(x) - \varphi(x)\frac{\partial f}{\partial\nu}(x)\,dS\\ &\qquad\qquad + \int_{\lVert x-x_0\rVert = \varepsilon} f(x)\frac{\partial \varphi}{\partial \nu}(x) - \varphi(x)\frac{\partial f}{\partial\nu}(x)\,dS\Biggr)\\ &= \lim_{\varepsilon \searrow 0} \int_{\lVert x-x_0\rVert = \varepsilon} f(x)\frac{\partial \varphi}{\partial \nu}(x) - \varphi(x)\frac{\partial f}{\partial\nu}(x)\,dS, \end{align}
where $dS$ denotes the surface measure of the sphere, $\frac{\partial}{\partial\nu}$ the directional derivative in direction of the outer normal, and the first integral on the right hand side vanishes since $\varphi$ vanishes in a neighbourhood of the outer sphere. The outer normal on the sphere $\lVert x-x_0\rVert = \varepsilon$ is $-\frac{x-x_0}{\lVert x-x_0\rVert}$,
so we are left with
$$\frac{1}{\varepsilon}\int_{\lVert x-x_0\rVert = \varepsilon} f(x)\langle \nabla\varphi(x),x-x_0\rangle - \varphi(x)\langle \nabla f, x-x_0\rangle\,dS.$$
An easy estimate shows
$$\frac{1}{\varepsilon}\int_{\lVert x-x_0\rVert = \varepsilon} f(x)\langle\nabla\varphi(x),x-x_0\rangle\,dS \leqslant \frac{1}{\varepsilon^2}\int_{\lVert x-x_0\rVert = \varepsilon} \lVert\nabla\varphi\rVert\cdot\varepsilon\,dS \leqslant C\frac{4\pi\varepsilon^2}{\varepsilon} \xrightarrow{\varepsilon \searrow 0} 0.$$
For the other part of the integral, computing $\frac{\partial f}{\partial\nu}$ gives
$$\int_{\lVert x-x_0\rVert = \varepsilon} \varphi(x)\cdot \frac{1}{\varepsilon^2}\,dS,$$
and since $\varphi$ is continuous at $x_0$ we have
$$\lim_{\varepsilon\searrow 0} \int_{\lVert x-x_0\rVert = \varepsilon} \varphi(x)\cdot \frac{1}{\varepsilon^2}\,dS = 4\pi\varphi(x_0).$$
Collecting everything without mucking up the signs, we get
$$(\nabla^2 f)[\varphi] = \int_{\mathbb{R}^3} f(x)\cdot \nabla^2\varphi(x)\,dx_1\,dx_2\,dx_3 = -4\pi \varphi(x_0) = -4\pi\delta_{x_0}[\varphi].$$
That holds for all test functions $\varphi$, hence $\nabla^2 f = -4\pi\delta_{x_0}$.