Let $\Omega$ be a bounded domain of $\mathbb{R}^d$. We define the Sobolev space $H^1(\Omega)$ by
\begin{align*}
H^1(\Omega)=\{f \in L^2(\Omega,m) \mid |\nabla f| \in L^2(\Omega)\}.
\end{align*}
Here, $m$ denotes the Lebesgue measure on $\Omega$, and $\nabla f$ is the distributional gradient of $f$. The Neumann Laplacian $(L,\text{Dom}(L))$ on $\Omega$ is defined as
\begin{align*}
\text{Dom}(L)&=\left\{f \in L^2(\Omega,m) : H^1(\Omega) \ni g \mapsto \int_{\Omega}\nabla f\cdot \nabla g\,dm \text{ is continuous on $L^2(\Omega,m)$}\right\} \\
-\int_{\Omega}g Lf\,dm&=\int_{\Omega} \nabla f\cdot \nabla g\,dm,\quad f \in \text{Dom}(L),\,g \in H^1(\Omega).
\end{align*}
Question. When $\partial \Omega$ is Lipschitz, we define
\begin{align*}
C=\left\{g \in C^2(\overline{\Omega}) :\frac{\partial g}{\partial n}(x)=0,\quad \sigma\text{-a.e.} \right\}
\end{align*}
Here, we denote by $n$ the inward unit normal vector on $\partial \Omega$, and $\sigma$ is the surface measure on $\partial D$. Can we show that
$\{(g,Lg) \mid g \in C\}$ is a dense subspace of $\{(f,Lf) \mid f \in \text{Dom}(L)\}$ with respect to the $L^2$-norm? This means that $C$ is a core of $(L,\text{Dom}(L))$.
Although this fact should be valid, I do not know the proof.
Please let me know if you have any references.
Best Answer
What helped me understand how things work with cores is the following consideration:
Now take your question for example. Let $$E=\{e\ : \ e \text{ is a finite linear combination of eigenfunctions of L}\}\subseteq \mathrm{D}(L).$$ Now, $$\{(e,L e) \ : \ e \in E\}$$ is a dense subset of $$\{(f,Lf) \ : \ f \in \mathrm{D}(L)\},$$ because of the density of the eigenfunctions in the domain. In particular, your question is now whether $$E \subseteq C =C^2_\mathrm{neu}(\Omega).$$ Now we have reduced this fairly abstract problem to a very concrete one. Pick any $\lambda\in(-\infty,0]$ (the Laplacian is a negative operator) and suppose that $e_\lambda \in L^2(\Omega), \|e_\lambda\|=1$ solves $$\Delta e_\lambda = \lambda e_\lambda.$$ What is the regularity of $e_\lambda$?
This kind of question is what Laplacians are build for, and maybe you already know that eigenfunctions are smooth. If not, to see this we can restate the problem as $$e_\lambda=(-\Delta +1)^{-1}\big( (1-\lambda) e_\lambda\big).$$ Now one can estimate $$\|e_\lambda\|_{H^2_{\mathrm{neu}}} = \bigg\|(-\Delta +1)^{-1} \big( (1-\lambda) e_\lambda\big)\bigg\|_{H^2_{\mathrm{neu}}} \leq C(1-\lambda) \|e_\lambda\|_{L^2} =C(1-\lambda) ,$$ where the constant $C>0$ depends on how you define the $H^2_{\mathrm{neu}}$ norm (choose your favourite version of Sobolev space). Iterating this inequality you find that $e_\lambda \in H^{2m}_{\mathrm{neu}}$ for any $m\in\mathbb{N}$. Then by Sobolev embedding you find that $e_\lambda \in C^{\infty}_{\mathrm{neu}}(\Omega)\subseteq C$.
A final note: this proof uses a lot of nontrivial tools, such as Sobolev embeddings, spectral decompositions etc.. With it I am not aiming for a simple direct proof (no doubt possible with some approximation via $C^2$ functions), but I hope it gives the idea of what's going on.