Calculation of a vector by taking the gradient of the integral of its divergence I have encountered several times of a special way of calculating a vector from the divergence of the vector. It has at least appeared in the theories of elasticity and electrodynamics.
If I define the operator of gradient is $\mathbf{\nabla}$ and the operator of a divergence is $\mathbf{\nabla} \cdot$, then a vector, $\mathbf{u}$, can be calculated by the following way:
$$\mathbf{u}(x,y,z) = - \frac{1}{4\pi} \mathbf{\nabla} \int{\frac{\mathbf{\nabla} \cdot \mathbf{u}(x',y',z') }{r} d V'}$$
in which $x$, $y$ and $z$ represent coordinates, $r = \sqrt{(x-x')^2 + (y-y')^2 + (z-z')^2}$, $\mathbf{\nabla} \cdot \mathbf{u} $ is a given function defined in all space and vanishing at infinity and the curl of $\mathbf{u}$ is zero. (Note that the integral works on the primed coordinates.)
I would like to know is there a name for this method? and how to prove it?
 A: What is going on is you're solving Poisson's equation in 3 dimensions. So, the problem is you're given some vector field $\mathbf{u}$ on $\Bbb{R}^3$, and you have specified its divergence $\nabla\cdot \mathbf{u}=-\rho$, for some known function $\rho$. You also impose the condition that it has no curl (electrostatics). As a result, we can introduce a scalar potential $\phi$ such that $\mathbf{u}=-\nabla \phi$ (a very special case of Poincare's lemma). As a result, we seek a function $\phi$ on $\Bbb{R}^3$ such that
\begin{align}
\nabla^2\phi&=\rho
\end{align}
(because we want $\nabla^2\phi=\nabla\cdot (\nabla\phi)=\nabla\cdot (-\mathbf{u})=-\nabla\cdot\mathbf{u}=-(-\rho)=\rho$). This equation is known as Poisson's equation, and is pretty much the starting point of serious electrostatics.
Now, the question becomes how might one go about solving such a PDE? Well, talking about it with full mathematical detail would require us to talk about distribution theory, which I certainly won't do now. So, let's discuss heuristically. Notice that the Laplacian is a linear operator, so in particular suppose someone has specified "charge densities" $\rho_1,\rho_2$ and a number $c\in\Bbb{R}$, and suppose by some miracle you manage to find functions $\phi_1,\phi_2$ such that $\nabla^2\phi_1=\rho_1$ and $\nabla^2\phi_2=\rho_2$. Now, suppose I tell you to find a function $\phi$ such that
\begin{align}
\nabla^2\phi=c\rho_1+\rho_2.
\end{align}
What would you do? Well, the most obvious thing would be to take $\phi=c\phi_1+\phi_2$, because then we clearly have
\begin{align}
\nabla^2(\phi)&=\nabla^2(c\phi_1+\phi_2)=c\nabla^2\phi_1+\nabla^2\phi_2=c\rho_1+\rho_2,
\end{align}
exactly as desired. This is also known as the principle of superposition. By the way, the function $\phi=c\phi_1+\phi_2$ is not the only possible solution; there are infinitely many others (to ensure uniqueness, we need to specify boundary conditions, such as vanishing at infinity).
Now, the next step is to take this further. Let us consider the problem of having a charge density equal to that of a unit positive point charge located at the origin, i.e a Dirac delta $\delta_0$ (it is from here that we need distribution theory to properly discuss things). This means we're looking for some function $G_0$ such that $\nabla^2G_0=\delta_0$; such a function is known as the fundamental solution to the Laplacian operator, or as the Green's function for the Laplacian operator. Let us for now leave aside the question of how to actually find $G_0$. Supposing the existence of $G_0$, let us try to reconstruct solutions to the general problem.

*

*First, what happens if the point charge is located at some point $p\in\Bbb{R}^3$, not necessarily at the origin, i.e the charge density is $\delta_p$. How can we find a solution in this case? Well, it is fairly easy to convince yourself that the function $\phi(\xi)=G_0(\xi-p)$ solves the equation $\nabla^2\phi=\delta_p$. In other, words, if you translate the location of the point charge, then we had better translate the argument of the function $G_0$ as well.


*What if we had several point charges at several locations, each of varying charges? Say you had point charges $q_1,\dots, q_k$ located at positions $p_1,\dots, p_k$. So, the charge density is $\rho=q_1\delta_{p_1}+\cdots q_k\delta_{p_k}$ (recall that Dirac delta in 3 spatial dimensions has units of 1/length$^3$, so this $\rho$ indeed has units of charge density). How can we find a solution in this case? Well, by the linearity and translation arguments I just stated, we have that $\phi(\xi)=\sum_{i=1}^kq_iG_{p_i}(\xi)=\sum_{i=1}^kq_iG_0(\xi-p_i)$ will satisfy Poisson's equation for this specific charge distribution. So, we're taking a "weighted sum" of functions in order to match the correct charge density.


*Finally, suppose we don't just have a finite set of point charges, and that we have an arbitrary (say compactly supported) charge distribution $\rho$. Then, we should expect to replace sums by integrals to get that $\phi(\xi)=\int_{\Bbb{R}^3}\rho(p)G_0(\xi-p)\,dV(p)$, i.e we integrate the charge density multiplied by the fundamental solution over all space (integrals with respect to $p$). Or if you prefer the primed notation, $\phi(\xi)=\int_{\Bbb{R}^3}\rho(\xi')G_0(\xi-\xi')\,dV'$.
It turns out that $G_0(\xi):=\frac{1}{4\pi |\xi|}$ satisfies $\nabla^2G_0=\delta_0$, and that this is the unique solution which is a function vanishing at infinity, so the above integral expression becomes $\phi(\xi)=\frac{1}{4\pi}\int_{\Bbb{R}^3}\frac{\rho(\xi')}{|\xi-\xi'|}\,dV'$.
In more fancy mathematical language, this is a convolution of the Green's function/fundamental solution with the charge density, written $\phi=G_0*\rho$, and hence (in the distributional sense) $\nabla^2\phi=(\nabla^2(G_0))*\rho=\delta_0*\rho=\rho$.
As a fun fact, if you were to solve Poisson's equation in $n$-spatial dimensions, then the fundamental solution is $G_0(\xi)=\frac{1}{(n-2)A_{n-1}}\frac{1}{|\xi|^{n-2}}$, where $A_{n-1}$ is the surface area of the $(n-1)$-dimensional sphere (if $n=3$, this is the surface area of the 2-dimensional sphere, which is $4\pi$; and this is exactly the $4\pi$ appearing in your formula). If you want the full gory details I have previously written up an answer regarding this matter (or you can simply google it).

To emphasize this again: the crux is not really that we have represented a vector field using its divergence. In E&M/Newtonian gravity, one should think of this as a solution obtained by superposition. This approach is useful for thinking because it applies to many other PDEs: the heat equation, the wave equation (obviously important if you talk about classical theory of light, i.e electromagnetic waves), Helmholtz's equation etc, all of which are obviously very important in physics.
