Need good reference or a proof on regularity of solution to Neumann problem

This question has already a bounty winning answer: however, since I found it interesting, upvoted it the very first day it was posted, and found unusual difficulties in finding a nice answer, I nevertheless want to share my own researches and thoughts.

The only reference I am aware of which gives a proof of this result, is the textbook by Mikhailov [1]. Precisely, he deals with the regularity problem for the first boundary problem (Dirichlet problem) and the second boundary problem (Neumann problem) for the Poisson equation in §2.3, pp. 216-226. Mikhailov proves the result directly: he first solves the problem in the case of homogeneous boundary conditions, proving the following theorem:

Theorem 4 ([1], p. 217). If $f\in H^k(\Omega)$ and $\partial\Omega\in C^{k+2}$ for certain $k\ge 0$, then the generalized solutions $u(x)$ of the first and second boundary-value problems with homogeneous boundary condition for the Poisson equation belong to $H^{k+2}(\Omega)$ and satisfy (in the case of the second boundary problem it is assumed that $\int_\Omega u \mathrm{d}x=0$) the inequality $$ \Vert u\Vert_{ H^{k+2}(\Omega)}\le C\Vert f\Vert_{ H^{k}(\Omega)} $$ where the constant $C>0$ does not depend on f.

By using this result, he extends the regularity theorem 4 above to the case of non homogeneous boundary conditions ([2], p. 226) by reducing non homogeneous boundary conditions to homogeneous ones. He explicitly does it only for the Dirichlet problem but the same method nevertheless works for the Neumann problem: let's see this. Consider a solution $u(x)$ of the problem \eqref{eqlocalvar-Neumann} above and a function $\Phi\in H^{k+2}(\Omega)$ whose normal derivative on $\partial\Omega$ is $g\in H^{k+3/2}(\partial\Omega)$, i.e. $$ \frac{\partial\Phi}{\partial\nu}=g\;\text{ on }\;\partial\Omega\label{1}\tag{1} $$ Define $w=u-\Phi$: then, for every $v\in H^{1}(\Omega)$, $$\label{GeneralizedNeumann}\tag{GN} \begin{split} \int_{\Omega} \nabla w \cdot \nabla v \, \mathrm{d}x&=\int_{\Omega} \nabla u \cdot \nabla v \, \mathrm{d}x - \int_{\Omega} \nabla \Phi \cdot \nabla v \, \mathrm{d}x \\ &= \int_\Omega f v \, \mathrm{d}x + \int_{\partial \Omega}gv \, \mathrm{d}\sigma(x) - \int_{\Omega} \nabla \Phi \cdot \nabla v \, \mathrm{d}x \\ &=\int_\Omega f v \, \mathrm{d}x + \int_{\partial \Omega}gv \, \mathrm{d}\sigma(x) - \int_{\Omega} \nabla\cdot(v\nabla\Phi) \, \mathrm{d}x + \int_{\Omega} v\Delta\Phi \, \mathrm{d}x \\ &=\int_\Omega \big[\,f + \Delta\Phi\big]v\, \mathrm{d}x + \int_{\partial \Omega}gv \, \mathrm{d}\sigma(x) - \int_{\partial\Omega} \nu \frac{\partial\Phi}{\partial\nu} \, \mathrm{d}\sigma(x), \end{split}\label{2}\tag{2} $$ and substituting \eqref{1} in \eqref{2} we get $$ \int_{\Omega} \nabla w \cdot \nabla v \, \mathrm{d}x=\int_\Omega f_1 v \, \mathrm{d}x $$ where $f_1=f+\Delta\Phi\in H^k(\Omega)$. This means that $w$ solves the Neumann problem \eqref{eqlocalvar-Neumann} with homogeneous boundary conditions (i.e. $g\equiv 0$), and by applying theorem 4 we have $$ w=u-\Phi\in H^{k+2}(\Omega) \iff u=w+\Phi\in H^{k+2}(\Omega) $$

Notes

• In Mikhailov's notation, $H^0=L^2$: also he does not use the standard notation for traces, for example expressing that $\varphi\in H^{1/2}(\partial G)$ by saying that $\varphi$ is a function in $L^2(\partial G)$ which is the trace of a function $\Phi\in H^1(G)$.
• Regarding the third boundary problem (the so called Robin problem), Mikhailov ([1], footnote *, p. 217) remarks that analysis of the regularity of solutions can be dealt as done in theorem 4 above for the solutions to the first and second boundary value problems, provided certain conditions are assumed.

[1] V. P. Mikhailov (1978), Partial differential equations, Translated from the Russian by P.C. Sinha. Revised from the 1976 Russian ed., Moscow: Mir Publishers, p. 396 MR0601389, Zbl 0388.3500.

A standard reference is the book of Grisvard Elliptic Problems in Nonsmooth Domains Grisvard. Chapter 2 is about regularity of elliptic equations in smooth domains. The theorem you want is Theorem 2.2.2.5. As for proofs, Grisvard mainly deals with the Dirichlet case and then in a short remark explains how to modify it for Neumann. It might not be what you want. If you want all the proofs, take a look at Theorems 130- 136 in these lectures notes. lecture notes