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).





Global regularity for Neumann problem
Proposition 7.7 of chapter 5 in the first volume of "Partial Differential Equations" by M. E. Taylor (Springer 1996) proves the regularity result for $C^\infty$ boundaries. (Notice that any two solutions differ at most by an additive constant.) The argument of the proof should go through for boundaries with $C^{2,1}$ regularity, because then the Neumann trace $\gamma\_1:H^2(\Omega)\to H^{1/2}(\partial \Omega)$ is known to be surjective, implying that it suffices to consider only the case $g=0$. I am not sure whether the assumption $\partial\Omega\in C^2$ is sufficient.


Among others the following sources contain detailed results about elliptic boundary problems which in particular apply to the Neumann boundary problem for the Laplace operator: J. Wloka, Partial Differential Equations, CUP, 1987; Chazarain and Piriou, Introduction to the theory of linear partial differential equations, North-Holland 1982, Chapter 5. The latter source uses the method of the Calderon projector to reduce an elliptic boundary value problem to a pseudo-differential equation on the boundary.


