$\Delta u$ is bounded. Can we say $u\in C^1$?

Let $\Omega\subset\mathbb{R}^n$ be a bounded open set. Let us say it has a Lipschitz boundary.

Consider the Laplacian $\Delta$ in the classical sense. Suppose $\Delta u=\frac{\partial^2}{\partial x_1^2}u+\dotsb+\frac{\partial^2}{\partial x_n^2}u$ is bounded.

Q: Can we say $u\in C^1(\Omega)$? Does it depend on the dimension $n$? Can we claim the smoothness recursively, i.e., if $\Delta^m u$: bounded, then,...etc?

I was pondering about the relations of partial differentiability and continuity, and got confused.

Bounded partial derivatives imply continuity says "If all partial derivatives of f are bounded, then f is continuous on E.", but we cannot apply this argument recursively as we do not have the "cross term" $\frac{\partial^2}{\partial x_i\partial x_j}$.

We have $\Delta u\in L^2(\Omega)$ as $\Delta u$ is bounded on a bounded region $\Omega$. However, we cannot use this result Sobolov Space $W^{2,2}\cap W^{1,2}_0$ norm equivalence and say $u\in H^2(\Omega)$, because 1. $u$ does not necessarily vanish on the boundary, and 2. we are not sure if $\frac{\partial^2}{\partial x_1^2}u+\dotsb+\frac{\partial^2}{\partial x_n^2}u+(\text{partial derivatives of cross terms})u$ are bounded.

Aha, from
Equivalent Norms on Sobolev Spaces, $\Delta u\in L^2(\Omega)$ is enough to say $u\in H^2(\Omega)$.

But one thing is I do not know if we have the same kind of equivalence for $m>3$, and another thing is resorting to Sobolev embedding does not sound like a good idea as it depends on the dimension heavily.

I wonder I could show this directly.

Do you know the regularity of Elliptic equation? The answer of your question is not just a Sobolev space problem, it is an Elliptic PDE problem. The answer is by the regularity of Laplace equation, you can boost up the regularity of solution based on the smoothness of your boundary. Check this book, Section 6.3.1 Theorem 1 and Theorem 2. In those theorems, $b$ and $c$ are only $0$ and $a_{ii}=1$ the identity matrix. $f$ is $\Delta u$ itself. Then use the regularity again and again.

One approach is via the Newtonian potential. You'll need the following ingredients.

1. The solution of Laplace's equation are $C^\infty$ smooth, even when the equation holds in the weak (distributional) sense.
2. Let $v = \Gamma*(\Delta u)$; this is a function such that $\Delta v=\Delta u$
3. Since $\nabla \Gamma$ is locally integrable and $\Delta u$ is bounded, we can differentiate under the integral sign: $\nabla (\Gamma*(\Delta u)) = (\nabla \Gamma)*(\Delta u)$
4. The convolution of an $L^1$ function and an $L^\infty$ function is continuous.

The above argument works the same in all dimensions.

If $\Delta^2 u$ is bounded, then $\Delta u\in C^1$ by the above. Then follow the same steps again, #4 being "The convolution of an $L^1$ function and a $C^1$ function is $C^1$". And so on.

Note the above argument isn't getting you the best regularity possible; you only get one derivative of $u$ for each iteration of Laplacian. The approach via Sobolev embeddings would give a more precise result.