Laplacian of a distribution Here is a small result annoying me : 

Let $u$ be a distribution ($u \in \mathcal{D}'(\mathbb{R}^n)$) such that $\Delta u$ is continuous. Then $u$ is continuous. 

I am not able to prove this, and have no idea of where to start. Anoyone with a link, a proof or hints ?
 A: Step 1 :
Let us define
$$\Gamma(x,y):=\begin{cases}\frac{1}{(2-n)\omega_{n-1}}|x-y|^{2-n} & \text{if }n\geq3 \\
\frac{1}{2\pi}\ln|x-y| & \text{if }n=2
\end{cases}$$
for all $x,y\in\mathbb{R}^n$, $x\neq y$, where $\omega_{n-1}$ is the Lebesgue measure of the unit sphere.
We show that
$$\Delta_x\Gamma(x,.)=\Delta_y\Gamma(x,.)=\delta_x.$$
Indeed, for all $1\leq i,j\leq n$, we have
\begin{align}
\partial_{x_i}\Gamma(x,y) &= \frac{1}{\omega_{n-1}}\frac{x_i-y_i}{|x-y|^n}, \\
\partial_{y_i}\Gamma(x,y) &= \frac{1}{\omega_{n-1}}\frac{y_i-x_i}{|x-y|^n}, \\
\partial_{x_j}\partial_{x_i}\Gamma(x,y) &= \partial_{y_j}\partial_{y_i}\Gamma(x,y) = \frac{1}{\omega_{n-1}}\frac{|x-y|^2\delta_{ij}-n(x_i-y_i)(x_j-y_j)}{|x-y|^{n+2}}.
\end{align}
Remark the symmetry in $x$ and $y$ for the second derivatives.
Next fix an arbitrary $R>0$ such that $x\in B_0(R)$, and take $\psi\in\mathcal{C}^{\infty}_c\left(\mathbb{R}^n;\mathbb{R}\right)$ - in fact, we only need that $\psi\in\mathcal{C}^{0,\alpha}_c\left(\mathbb{R}^n;\mathbb{R}\right)$ with $\alpha>0$; we then write
\begin{align} \int_{B_0(R)}\partial_{x_j}\partial_{x_i}\Gamma(x,y)\psi(y)\mathrm{d}y & = \frac{1}{\omega_{n-1}}\int_{B_0(R)}\frac{|x-y|^2\delta_{ij}-n(x_i-y_i)(x_j-y_j)}{|x-y|^{n+2}}\psi(y)\mathrm{d}y \\
&= \frac{1}{\omega_{n-1}}\int_{B_0(R)}\frac{|x-y|^2\delta_{ij}-n(x_i-y_i)(x_j-y_j)}{|x-y|^{n+2}}\left(\psi(y)-\psi(x)\right)\mathrm{d}y \\
& -\frac{\psi(x)}{\omega_{n-1}}\int_{\partial B_0(R)}\frac{(x_i-y_i)}{|x-y|^{n}}\frac{y_j}{R}\mathrm{d}\sigma_y
\end{align}
where we use the Green's formula to integrate by part to get the second term in the last expression (observe that taking the partial derivative of $\Gamma$ with respect to $x_i$ and then integrating with respect to $y_i$ gives us a minus sign). We can check now that these integrals are finite, thanks to the Hölder continuity of $\psi$. Then
\begin{align}\int_{B_0(R)}\Delta_x\Gamma(x,y)\psi(y)\mathrm{d}y & = -\frac{\psi(x)}{\omega_{n-1}}\int_{\partial B_0(R)}\sum_{i=1}^{n}\frac{(x_i-y_i)}{|x-y|^{n}}\frac{y_i}{R}\mathrm{d}\sigma_y \\
& = -\frac{\psi(x)}{\omega_{n-1}}\int_{\partial B_0(R)}\frac{(x-y)}{|x-y|^{n}}\cdot\frac{y}{R}\mathrm{d}\sigma_y \\
& = -\frac{\psi(x)}{\omega_{n-1}}\int_{\partial B_0(1)}\frac{(x-R\xi)}{|x-R\xi|^{n}}\cdot\frac{R\xi}{R}R^{n-1}\mathrm{d}\sigma_\xi \\
& = -\frac{\psi(x)}{\omega_{n-1}}\int_{\partial B_0(1)}\frac{(\frac{x}{R}-\xi)}{|\frac{x}{R}-\xi|^{n}}\cdot\xi\mathrm{d}\sigma_\xi \\
\end{align}
and we get the desired result by letting $R\to+\infty$.
Step 2 : Choose $x,y\in\mathbb{R}^n$, $x\neq y$ and let $\delta:=|x-y|$ and $\xi:=\frac{x+y}{2}$. Next, for $m\geq 1$, let $\rho_m\in\mathcal{C}^{\infty}\left(B_\xi(2\delta);\mathbb{R}\right)$ a bump function in $B_\xi(2\delta)$ such that
$$\rho_m(t)=\begin{cases} 0 & \text{if }0\leq t<\frac{1}{m} \\
1 & \text{if }\frac{2}{m}\leq t\leq\delta-\frac{2}{m} \\
0 & \text{if }\delta-\frac{1}{m}<t\leq\delta
\end{cases}.$$
Define $\varphi_m^x(z):=\Gamma(x,z)\rho_m(|x-z|)$ and $\varphi_m^x(z):=\Gamma(y,z)\rho_m(|y-z|)$ for all $z\in B_\xi(2\delta)$. As the cut-off $\rho_m$ kills the singularity of $\Gamma$, these functions are smooth in $B_\xi(2\delta)$. Furthermore,
\begin{align} \|\Gamma(x,.)-\varphi_m^x\|_{L^1\left(B_\xi(2\delta);\mathbb{R}\right)} &= \frac{1}{(n-2)\omega_{n-1}}\int_{B_\xi(2\delta)}\frac{1-\rho_m(|x-z|)}{|x-z|^{n-2}}\mathrm{d}z \\
&\leq \frac{1}{(n-2)\omega_{n-1}}\left(\int_{B_\xi(2\delta)\cap B_x(\frac{2}{m})}\frac{1-\rho_m(|x-z|)}{|x-z|^{n-2}}\mathrm{d}z+\int_{B_\xi(2\delta)\setminus B_x(\delta-\frac{2}{m})}\frac{1-\rho_m(|x-z|)}{|x-z|^{n-2}}\mathrm{d}z\right) \\
&\leq \frac{1}{(n-2)\omega_{n-1}}\left(\int_{B_\xi(2\delta)\cap B_x(\frac{2}{m})}\frac{\mathrm{d}z}{|x-z|^{n-2}}+\int_{B_\xi(2\delta)\setminus B_x(\delta-\frac{2}{m})}\frac{\mathrm{d}z}{|x-z|^{n-2}}\right)
\end{align}
and for all $R>0$ we have
\begin{align}\int_{B_x(R)}\frac{\mathrm{d}z}{|x-z|^{n-2}} &= \int_{B_0(R)}\frac{\mathrm{d}\eta}{|\eta|^{n-2}} \\
&= \int_{s=0}^R\frac{1}{s^{n-2}}\int_{\partial B_0(s)}\mathrm{d}\sigma\mathrm{d}s \\
&= \omega_{n-1}\int_{s=0}^R s\mathrm{d}s \\
&= O\left(R^2\right)
\end{align}
whereas
\begin{align}\int_{B_\xi(2\delta)\setminus B_x(\delta-\frac{2}{m})}\frac{\mathrm{d}z}{|x-z|^{n-2}} & \leq \int_{B_\xi(2\delta)\setminus B_x(\delta-\frac{2}{m})}\frac{\mathrm{d}z}{\left(\delta-\frac{2}{m}\right)^{n-2}} \\
&\leq \frac{\frac{2}{m}}{\left(\delta-\frac{2}{m}\right)^{n-2}}
\end{align}
and thus we get
$$\|\Gamma(x,.)-\varphi_m^x\|_{L^1\left(B_\xi(2\delta);\mathbb{R}\right)} = O\left(\frac{1}{m}\right)$$
that is $\rho^x_m\to\Gamma(x,.)$ in $L^1$-norm. Idem for $\rho_m^y$ and $\Gamma(y,.)$.
Step 3 : Fix $\varepsilon>0$ ; because $\Gamma(.,z)$ is continuous out off a neighborhood of $z$, there exists $\delta>0$ such that, for all $x,y\in B_\xi(2\delta)$ with $|x-y|<\delta$, we have $|\varphi^x_m(z)-\varphi^y_m(z)|<\varepsilon$. Hence, as $\Delta u\in L^{\infty}\left(B_\xi(2\delta);\mathbb{R}\right)$ (because it is supposed to be continuous), it comes
\begin{align} 2\sup_{B_\xi(2\delta)}|\Delta u|\delta\varepsilon &\geq \int_{B_\xi(2\delta)}\Delta u(z)\left(\varphi^x_m(z)-\varphi^y_m(z)\right)\mathrm{d}z \\
&= \int_{B_\xi(2\delta)}\Delta u(z)\left(\Gamma(x,z)-\Gamma(y,z)\right)\mathrm{d}z+O\left(\frac{1}{m}\right)
\end{align}
provided that $|x-y|<\delta$. To conclude, recall the symmetry $\Delta_x\Gamma(x,.)=\Delta_y\Gamma(x,.)$ seen in the Step $1$, so that
\begin{align} 2\sup_{B_\xi(2\delta)}|\Delta u|\delta\varepsilon &\geq  \int_{B_\xi(2\delta)}u(z)\left(\Delta_x\Gamma(x,z)-\Delta_y\Gamma(y,z)\right)\mathrm{d}z+O\left(\frac{1}{m}\right) \\
&= \int_{B_\xi(2\delta)}u(z)\left(\delta_x(z)-\delta_y(z)\right)\mathrm{d}z+O\left(\frac{1}{m}\right) \\
&= u(x)-u(y)+O\left(\frac{1}{m}\right).
\end{align}
Finally, for $m\gg 1$, we have
$$u(x)-u(y)\leq C\varepsilon\quad\quad\quad\quad C>0.$$
A: If $E$ is a distribution atisfying $\Delta E = \delta _0$ then $T:=u-E* \Delta u$ satisfies $\Delta T=0$ so that $T$ is $\mathcal{C}^{\infty}$ and then $u$ is continuous. 
But the problem is with the fact that I can not take $E* \Delta u$ because $E$ is not compactly supported, but still if $\phi$ denotes a test function such that $\phi = 1$ in an open ball centered in $x$ then $u-(\phi E) *u \in \mathcal{C}^{\infty}$ and the result is proved. 
