Let $\Omega\subset \mathbb{R}^d $ be a bounded open subset ($d\in \mathbb{N}$) and denote $\partial\Omega$ its boundary which we assume to be Lipschitz. The classical inhomogeneous Neumann problem for Laplace operator associate to data $f:\Omega\to\mathbb{R}$ and $g: \partial\Omega \to \mathbb{R}$ (measurable functions) consists on finding a function $u:\Omega\to \mathbb{R}$ satisfying

\begin{equation}\label{eqlocal-Neumann}\tag{$N_1$} -\Delta u = f \quad\text{in}~~~ \Omega \quad\quad\text{ and } \quad\quad \frac{\partial u}{\partial \nu}= g ~~~ \text{on}~~~ \partial \Omega. \end{equation}

In the standard setting one usually choose $ f$ in $L^2(\Omega)$ or in the dual space of $H^1(\Omega)$ and $g$ can be choose in the trace spaces of $H^1(\Omega)$ denote by $H^{\frac{1}{2}}(\partial\Omega)$ or in its dual $H^{-\frac{1}{2}}(\partial\Omega)$.

Assume $f\in L^{2}(\Omega)$ and $g \in H^{1/2}(\partial\Omega)$. We have the following Green-Gauss formula

$$\label{eqgreen-Gauss} \int_{\Omega} (-\Delta) u v \, \mathrm{d}x = \int_{\Omega} \nabla u \cdot \nabla v \, \mathrm{d}x- \int_{\partial \Omega} \gamma_{1} u \gamma_{0}v \, \mathrm{d}\sigma(x), \quad u\in H^{2}(\Omega) ~\hbox{and}~v\in H^{1}(\Omega). $$

Henceforth, on $\partial \Omega$ we merely write $\gamma_0 v= v$ and $ \displaystyle\gamma_1 v=\displaystyle \frac{\partial v}{\partial \nu} $.

Clearly from this Green-Gauss formula, if $u\in H^{2}(\Omega)$ and solves \eqref{eqlocal-Neumann} then $u $ satisfies the variational problem

$$\label{eqlocalvar-Neumann}\tag{$V_1$} \int_{\Omega} \nabla u \cdot \nabla v \, \mathrm{d}x= \int_{\partial \Omega} f v \, \mathrm{d}x + \int_{\partial \Omega}gv \, \mathrm{d}\sigma(x), \qquad \hbox{for all } ~~v\in H^{1}(\Omega). $$ In particular if we put $v=1$ the above formulation becomes $$\label{eqlocalcompatible-Neumann}\tag{$C_1$} \int_{\Omega}f\mathrm{d}x+ \int_{\partial \Omega}g\mathrm{d}\sigma(x)=0 $$ which is the compatibility condition.

Vice versa we have the following **global regularity result**.

**Theorem**

Assume $\Omega\subset \mathbb{R}^{d}$ is bounded open with $C^2$-boundary. If a function $ u \in H^{1}(\Omega)$ is solution to \eqref{eqlocal-Neumann} with $f\in L^{2}(\Omega)$ and $g \in H^{1/2}(\partial\Omega)$ then it belongs to $ H^{2}(\Omega).$

Question:In which book or recommendable reference can I find the proof of the above theorem? I know that the proof of this Theorem for the corresponding Dirichlet problem has been done in the book by Brezis or by Evans. Patently, both references avoid the inhomogeneous Neumann problem.

**Remark**
Moreover, observe that if $u$ solves \eqref{eqlocal-Neumann} or \eqref{eqlocalvar-Neumann} so does $\tilde{u} = u+c$ for every $c\in\mathbb{R} $ (that is invariant under additive constant).