Are Neumann Laplacian eigenfunctions in $C(\overline{\Omega})$? Consider that $u\in H^1(\Omega)$ with $\Delta u\in L^2(\Omega)$ (in the distributional sense) such that for some $\lambda>0$ we have that:
$$\begin{cases} \Delta u(x)=\lambda u(x), & x\in\Omega\\ \dfrac{\partial u}{\partial\nu} (x)=0, & x\in\partial\Omega\end{cases}$$
We assume $\Omega\subseteq\mathbb{R}^2$ to be an open, connected, bounded and has a uniform Lipschitz boundary.
1) Is it true that $u\in H^2(\Omega)$? If this is not true:
2) Is it true that $u\in C(\overline{\Omega})$? If this is not true:
3) How can we prove that $u\in L^{\infty}(\Omega)$?
I know that (3) is valid from the inequality posted here: Contractivity of Neumann Laplacean
But I do not know how to prove that inequality. Maybe it can be done in an easier way...
I wonder if there is an estimate of the form:
$\Vert u\Vert_{\infty}\leq c\lambda^{\alpha}$, where $c$ is a constant depending on $\Omega$?
I know that such estimates hold for for the Dirichlet laplacian.
I found some references about the problem here:
https://math.stackexchange.com/questions/2309436/regularity-of-laplacian-eigenfunctions-in-convex-polygon
but they do not represent counter-examples for any of my questions.
P.S. I found that (2) might also be true from that post: Is the linear span of the Neumann eigenfunctions dense in $C(\overline{D})$ but I did not understand the argument. Is there any clear way of getting more information about the regularity of the eigenfunction $u$ of $\lambda$?
 A: For (3), this follows from De-Giorgi-Nash-Moser so long as you are OK with dependence of the constant on the $C^{0, 1}$ character of the domain, not just e.g. the size.
This can be checked by running through the De Giorgi proof (all the energy inequalities are valid for the Neumann problem), or, with some extra care, flattening the boundary. There is a book by Gary Lieberman, "Oblique Derivative Problems for Elliptic Equations," which is a good reference for the regularity theory in general.
There are some tricks that will give a power dependence on $\lambda$ like stated in the question for some large, explicit $\alpha$. First treat $\lambda u$ as a zero-order term and run a finite number of Moser iterations (can be skipped in 2D to get $\alpha = 1$) to show that $\|u\|_{L^p} \leq C \lambda^\alpha \|u\|_{L^2}$ for a $p > \frac{n}{2}$, then treat $\lambda u = f$ in $- \Delta u = f$ and show that $\|u\|_{L^\infty} \leq C \|f\|_{L^p}$. This will generally not give the optimal value $\alpha$.
This also answers (2); De Giorgi-Nash-Moser gives that $u \in C^{0, \beta}$ for some $\beta$ which will depend on the Lipschitz constant of the domain.
(1) is false, as can be checked on non-convex polygons in the plane by blowing up at the corners. The blow-up limit is essentially explicit (like $z^\alpha$ for some $\alpha > \frac{1}{2}$, but tending to $\frac{1}{2}$ as the corner becomes more oblique), so locally this has $|D^2 u|^2 \approx |z|^{2\alpha - 4}$, which fails to be integrable. While obviously not a proof, this is the correct heuristic here.
It's unclear to me whether or not (3) remains true with $c=c(|\Omega|)$ only, at least specifically in 2D.
