Elliptic regularity for the Neumann problem I'm trying to understand how to establish regularity for elliptic equations on bounded domains with Neumann data.
For simplicity, let's presume we are focusing on $-\Delta u = f$ in $\Omega$ and $\frac{\partial u}{\partial \nu} = 0$ on $\partial \Omega$. 
Interior regularity works the same as always.
When proving boundary regularity, for the dirichlet boundary case we first consider some ball $B(0,1) \cap \mathbb{R}_+^n$ and let $\xi = 1$ on $B(0,1/2)$, $\xi = 0$ on $\mathbb{R}^n - B(0,1)$ and then estimate all derivatives $\frac{\partial^2 u}{\partial x_i \partial x_j}$ except $\partial^2 u/\partial x_n^2$. Two main points are needed
1) $\xi$ vanishes on the curved part of $B(0,1) \cap \mathbb{R}_+^n$\
2) $u=0$ on $\{x_n=0\}$.
This allows us to let $-\partial_{x_i} (\xi \partial_{x_j} u)$ (with derivatvies replaced by difference quotients) be an admissible test function for our weak definitoin of a solution.
I presume the main difficulty in neumann boundary data is making your test function be admissible. In other words, we would need $\int v = 0$ since our existence was established on $H^1(\Omega)$ restricted to mean value zero functions. 
So in order to proceed, can we just subtract off a constant from our original $-\partial_{x_i}(\xi \partial_{x_j}u)$?
Is there some more natural way to establish regularity in this case? I do not want to take advantage of the fact that we have a green's function in this case however as I only chose the Laplace equation for simplicity.
 A: It's true that to prove existence for the Neumann problem you work in the space $X=\{v\in H^1(\Omega):\, \int_\Omega v\,dx=0\}$ and so a weak solution $u$ satisfies 
$$\int_\Omega\nabla u\cdot\nabla v\,dx=\int_\Omega fv\,dx$$
for all $v\in X$. However, if you take $w\in H^1(\Omega)$ and consider the function $v=w-\frac1{|\Omega|}\int_\Omega v\,dy$ you get 
$$\int_\Omega\nabla u\cdot\nabla w\,dx=\int_\Omega f\left(w-\frac1{|\Omega|}\int_\Omega v\,dy\right)\,dx=\int_\Omega fw\,dx-0$$
since $\int_\Omega f\,dx=0$. So You can actually work in the entire space $H^1(\Omega)$ rather than in the smaller space $X$.
Difference quotients work fine since the constraint is gone. If $\Omega=B^+(0,1)$ you can do difference quotients in the tangential directions $e_i$, $i=1,\cdots, n-1$ (so no need to assume $u=0$ on $\{x_n=0\}$) to prove that all second derivatives but $\partial^2 u/\partial x_n^2$ are in $L^2 (B^+(0,1/2))$ and then you use the equation (which is satisfied pointwise a.e. by interior regularity) to conclude that also $\partial^2 u/\partial x_n^2$ is  in $L^2 (B^+(0,1/2))$.
A: In the case that you mentioned, we want to avoid this cut-off/difference quotients approach, since it could be hard to prove that $\partial_{x_i} (\xi \partial_{x_j} u)$ is a valid test function.
In general, when working with regularity theory, another standard approach is to use an 'approximated problem'. However, the kind of the approximated problem, of course, depends on  the PDE.
For the Neumann like problem I suggest the following approximation:
First observe that since $\int_\Omega f = 0$, we can easily construct a sequence $\{f_n\} \subset C^\infty(\Omega)$ such that $f_n \to f$ in $L^2(\Omega)$ and $\int_\Omega f_n=0, \ \forall n \in \mathbb{N}$.
Then we consider  $u_n \in C^\infty(\Omega)$ such that
$(*) -\Delta u_n + \frac{1}{n} u_n = f_n, \mbox{ in } \Omega \ \ $  and $\dfrac{\partial u_n}{\partial \nu}=0 \ , \mbox{on }\partial \Omega.$
The sequence $\{ u_n \}$ can be obtained by the use of Theorem 2.2.2.5, p.91 and  Theorem 2.5.1.1, p. 121 of Grisvard's Book p. 91.
In fact you just need to use a boostrap argument to $-\Delta u_n = -\frac{1}{n} u_n + f_n$.
Notice that $\int_\Omega u_n=n\int_\Omega f_n = 0$.
Now, you use $u_n$  as your test functions and obtain the following estimate:
$(**) \|\nabla u_n\|_{L^2}^2 \leq \|f_n\|_{L^2}^2$, $\forall \ n \in \mathbb{N}$
Now you use $-\Delta u_n $, as a test function to your PDE ( observe that $-\Delta u_n$ is a valid test function, anyway we don't need to worry about it since the approximated equation holds everywhere).
After integrating by parts, by using $(**)$  with some standard manipulations with your boundary terms you end up with
$\|D^2 u_n \|_{L^2}^2 \leq C(\partial \Omega)\|f_n\|_{L^2}^2$, $\forall \ n \in \mathbb{N}$.
(For instance, see Grisvard's book p.132-138, in particular eq. 3.1.1.5)
The key point for the above estimation is to control the boundary elements in terms of the mean curvature of $\partial \Omega$.
Now, since $\int_\Omega u_n =0$ we conclude that 
$\|u_n \|_{H^2}^2 \leq C(\partial \Omega)\|f_n \|_{L^2}^2 $
so that 
$\|u_n \|_{H^2}^2 \leq C(\partial \Omega) \|f\|^2$, the $L^2$ norm of $f$.
In this way we obtain $u\in H^2$ such that $u_n \to u$ weakly in $H^2$ and strongly in $H^1$.
Observe that the latter convergence is sufficient to handle the term $\dfrac{1}{n}u_n$.
Then, we can pass to the limit in equation (*) so that $u$ is a strong solution of
$-\Delta u = f$ in $\Omega$
$\dfrac{\partial u}{\partial \nu}=0$ on $\partial \Omega$ 
with 
$\|u\|_{H^2} \leq C(\partial \Omega) \|f\|$, the $L^2$ norm of $f$.
