Finding The Laplacian of $\frac{1}{|\vec x|^2}$ in the sense of distributions I am currently trying to find the Laplacian of $u(x)=\frac{1}{|\vec x|^2}$ in $3$ Dimensions in the sense of distributions. Now, let us have a test function $\phi(x)$. We consider two cases, when $0\notin \text{supp}\>\phi$, and when $0\in \text{supp }\phi$. Now in first case, it is obvious that the Laplacian in the sense of distributions is just the pointwise Laplacian. Now I considered when $0\in\text{supp }\phi$. Then, let $B(0,\epsilon)$ be some ball centered at $0$ with radius $\epsilon$, and consider $\Omega$ the support of $\phi$. Then, consider $\Omega_\epsilon=\Omega\setminus B(0,\epsilon)$. Integrating over this using Green's Second Identity yielded:$$\iiint_{\Omega_{\epsilon}}\frac{1}{|x|^2}\Delta\phi dx = \iiint_{\Omega_{\epsilon}}\Delta\left(\frac{1}{|x|^2}\right)\phi dx+\iint_{\partial \Omega_{\epsilon}}\frac{1}{|x|^2} \frac{\partial\phi}{\partial\vec n} dS_{x}-\iint_{\partial \Omega_{\epsilon}}\frac{\partial}{\partial\vec n}\left(\frac{1}{|x|^2}\right)\phi dS_x\>\>\>\>(1)$$
Now, I focused on the last integral on the right hand side above:
$$\iint_{\partial \Omega_{\epsilon}}\frac{\partial}{\partial\vec n}\left(\frac{1}{|x|^2}\right)\phi dS_x=8\pi\epsilon\left(\frac{1}{4\pi\epsilon^2}\iint_{\partial B(0,\epsilon)}\phi dS_x\right)$$
By the Mean Value Property, the value in the brackets above gave $\phi(0)$ and if we take $\epsilon\rightarrow 0$, we get the above tends to $0$. Now I looked at the second integral on the right hand side of $(1)$. 
$$\iint_{\partial B(0,\epsilon)}\frac{1}{|x|^2}\frac{\partial\phi}{\partial \vec n} dS_x = \frac{1}{\epsilon^2}\iint_{\partial B(0,\epsilon)}\frac{\partial\phi}{\partial\vec n}dS_x$$
Now since $\phi$ is bounded as it is compact, its derivatives hold the same characteristics, and thus:
$$\frac{1}{\epsilon^2}\iint_{\partial B(0,\epsilon)}\frac{\partial\phi}{\partial\vec n}dS_x\leq C\cdot \frac{1}{\epsilon^2}\cdot 4\pi\epsilon^2$$
The above gives a constant. Thus, am I supposed to conclude that the distributional Laplacian when the support of $\phi$ includes $0$ is simply the pointwise Laplacian plus some arbitrary constant?
 A: No need for "two cases": the second case is the general case, since we are never going to use the assumption $0\in \operatorname{supp}\phi$.  
Incorrect

By the Mean Value Property, the value in the brackets above gave $ϕ(0)$

The function $\phi$ is not harmonic, so it does not satisfy the Mean Value Property. But it is true that 
$$\frac{1}{4\pi\epsilon^2}\iint_{\partial B(0,\epsilon)}\phi dS_x \to \phi(0)$$
because $\phi$ is continuous. However, there is more to say about this term...
Incorrect
$$\iint_{\partial \Omega_{\epsilon}}\frac{\partial}{\partial\vec n}\left(\frac{1}{|x|^2}\right)\phi dS_x=8\pi\epsilon\left(\frac{1}{4\pi\epsilon^2}\iint_{\partial B(0,\epsilon)}\phi dS_x\right)$$
The normal derivative of $1/|x|^2$ on the sphere of radius $\epsilon$ is $2\epsilon^{-3}$, not $2\epsilon^{-1}$. As a result, this term tends to infinity, which is troublesome (but see below). 
Correct but incomplete
$$\frac{1}{\epsilon^2}\iint_{\partial B(0,\epsilon)}\frac{\partial\phi}{\partial\vec n}dS_x\leq C\cdot \frac{1}{\epsilon^2}\cdot 4\pi\epsilon^2$$
This term actually tends to $0$, because the divergence theorem implies
$$
\iint_{\partial B(0,\epsilon)}\frac{\partial\phi}{\partial\vec n}dS_x
= \iiint_{B(0, \epsilon)} \Delta \phi(x)\,dx \le C \epsilon^3
$$
Conclusion
We need the limit as $\epsilon\to 0$ of 
$$
\iiint_{\Omega_{\epsilon}}\Delta\left(\frac{1}{|x|^2}\right)\phi dx
-2\epsilon^{-3}\iint_{\partial \Omega_{\epsilon}} \phi dS_x
\tag{A}$$
The pointwise Laplacian should not be casually dismissed: it is not locally  integrable, hence does not define a distribution. The formula (A) represents its regularization; the subtracted term tempers the singularity at $0$.
Using the divergence formula, one can compute
$$
\iiint_{\Omega_{\epsilon}}\Delta\left(\frac{1}{|x|^2}\right) dx = 2\epsilon^{-3} (4\pi \epsilon^2)$$
which allows us to rewrite (A) as 
$$
\iiint_{\Omega_{\epsilon}}\Delta\left(\frac{1}{|x|^2}\right)(\phi-\phi(0))\,dx
-2\epsilon^{-3}\iint_{\partial \Omega_{\epsilon}}( \phi-\phi(0)) dS_x
\tag{B}$$
Here the second term tends to zero because $\phi(x) - \phi(0 ) = \nabla \phi(0) x + O(|x|^2)$, and the integral of $x$ over a sphere is zero. Final answer: the Laplacian is
$$
\phi\mapsto \lim_{\epsilon\to 0} \iiint_{\Omega_{\epsilon}}\Delta\left(\frac{1}{|x|^2}\right)(\phi-\phi(0))\,dx
$$ 
where you can put the formula for pointwise $\Delta\left(\frac{1}{|x|^2}\right)$ (some constant times $|x|^{-4}$) to make it more explicit. 
This is not "pointwise Laplacian plus a constant", this is pointwise Laplacian regularized.
