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 ?
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.$$