On the domain of the Neumann Laplacian Let $U$ be a bounded domain of $\mathbb{R}^d$,  and write $m$ for the Lebesgue measure on $U$. For $k=1,2$, we denote by $H^k(U)$ the set of all locally $m$-integrable functions $u\colon U \to \mathbb{R}$
such that for any multi-index $\alpha$ with $|\alpha|\le k$, the weak derivative $D^\alpha u$ exists and belongs to $L^2(U,m)$.
Define the Neumann Laplacian $(L,\text{Dom}(L))$ on $U$ by
\begin{eqnarray*}
\text{Dom}(L) & = & \{u \in H^1(U) :  H^1(U) \ni v \mapsto \int_{U}\nabla u\cdot \nabla v\,dm  \\ 
& & \text{ is continuous on $L^2(U,m)$}\} \\
-\int_{U}v Lu\,dm &= &\int_{U} \nabla u\cdot \nabla v\,dm,\quad u \in \text{Dom}(L),\,v \in H^1(U).
\end{eqnarray*}
I know that if $U$ is a bounded $C^2$ domain, $\text{Dom}(L) \subset H^2(U)$ (see Section~10.6.2 in [1] ). Even if $U$ is a bounded Lipschitz domain, does this inclusion hold ? I don't think this is correct in general.
For example, if $U$ is a $C^{1,1}$ domain or convex domain, does this holds?
I would like to know various conditions for $U$ such that the inclusion  $\text{Dom}(L) \subset H^2(U)$ holds.
 A: This is a partial (positive) answer for the convex case only but not every detail has been worked out.
Let first $U$ be convex and smooth and all functions be in $C^3$ up to the boundary. Integrating by parts one obtains
$$
\int_U |\Delta u|^2=\int_U |D^2 u|^2+\int_{\partial U}\sum_{i,j}(D_{jj}uD_i u-D_{ij}uD_j u)\nu_i d\sigma.
$$
Here $|D^2 u|^2 =\sum_{i,j}|D_{ij}u|^2$. The first term in the boundary integral vanishes because $\langle \nabla u, \nu \rangle =0$. Differentiating this equality along any tangent vector $h$ to $\partial U$ one gets $\langle D^2 u\, h, \nu \rangle +\langle \nabla u , D_h \nu \rangle=0$.
If $h=\nabla u$, the convexity of $U$ gives $\langle h, D_h \nu \rangle \geq 0$ and then $\langle D^2 u\, \nabla u, \nu\rangle \leq 0$ and finally by the equality displayed $\|D^2 u\|_2 \leq \|\Delta u\|_2$. This is an a priori estimates for smooth functions and smooth convex sets, with a constant independent on the smoothness of $U$ and an approximation argument should give the result.
