Why are functions with vanishing normal derivative dense in smooth functions?

Question

Let $M$ be a compact Riemannian manifold with piecewise smooth boundary. Why are smooth functions with vanishing normal derivative dense in $C^\infty(M)$ in the $H^1$ norm?

Here I define $C^\infty(M)$ to be those functions which have all orders of derivative continuous on $M$ and smooth in its interior. For example, $(x\mapsto \sin(\pi x))\in C^\infty([0,1])$ but $(x\mapsto \sqrt{x})\notin C^\infty([0,1])$.

I have crossposted this to MathOverflow.

Background

This is inspired by my belief that the form domain of the Friedrichs extension of the Neumann Laplacian on $M$ is equal to $H^1(M)$. If my belief is wrong, I would certainly accept as answer a counterexample, preferably with some discussion/references.

Here are the approaches I'm exploring.

• Given $u\in C^\infty(M)$ construct a function that agrees away from an $\epsilon$ neighborhood of $\partial M$, but has been modified to have zero normal derivative. This is described below, but runs into some trouble at corners and with smoothing.
• Given $u\in C^\infty(M)$, find some $\eta$ supported on an $\epsilon$ neighborhood of $\partial M$ such that $\nabla\eta|_{\partial M}$ is equal to the projection of $\nabla u|_{\partial M}$ onto the normal direction and $\eta$ and $|\nabla\eta|$ uniformly bounded in $\epsilon$. Then $u - \eta$ will be the desired approximation and uniform boundedness will imply $\|u - (u-\eta)\|_{H^1}\to 0$. I'm not sure if this is a search for an integrable harmonic vector field or if it's a constrained optimization problem.
• Simply show that any $u\in C^\infty(M)$ that is perpendicular to all smooth functions with vanishing normal derivative must be zero. In order to do this, I think it still runs into the same fundamental difficulty as the others, which is constructing functions with vanishing normal derivative supported on an $\epsilon$-neighborhood of the boundary.

(I'm also happy to simply have a reference to follow up. A reference for domains in $\mathbb{R}^n$ should be fine, too, and just a partition of unity away from a Riemannian manifold.)

Current Approach

The approach I'm considering right now is an elaboration on possibility (2) from my list. The sketch is: set the boundary condition that $\nabla\eta$ on the boundary be equal to the normal component of $\nabla u$ and find an integrable harmonic vector field $E$ satisfying that boundary condition. Then cut off $E$ so it is only supported in a neighborhood of $\partial M$, and let $\eta$ be so that $E = \nabla\eta$.

Intuition is that the maximum principle will control $|E|$ and $|\eta|$, so that $\eta$ will be bounded above in the $H^1$ norm by a constant times the volume of the $\epsilon$-neighborhood of $\partial M$.

Another approach is inspired by zhw's comment below: approximate $\nabla u$ among $L^2$ vector fields, multiply it by a cutoff function so it is supported in the interior of $M$, and then integrate it to approximate $u$. This should work with some tweaking in $[0,1]$ but I'm not sure how well it will work in general.

Older work

My approach has been as follows. The intuition is to take an arbitrary smooth function, restrict it to the complement of a collar neighborhood, then extend the restriction to the collar neighborhood so that the value is constant on inward-normal geodesics. However I'm running into issues at corners.

Let $\epsilon > 0$ be such that $\{p\in M\ |\ d(p,\partial M) < \epsilon\}$ is a collar neighborhood of $\partial M$. Let $e_\epsilon(p)$ for $p\in\partial M$ be the smaller of $\epsilon$ or the greatest time parameter such that the inward normal geodesic collides with no other inward normal geodesic. Let $N$ be the set of inward normal vectors on $\partial M$ whose lengths are no greater than $e_\epsilon$. The interior of the set $\operatorname{exp}(N)$ is foliated by geodesics. Edit- Note this is not necessarily true if the manifold has inward corners.

Given a smooth, continuous function $u$ on $M$, define $\bar{u}$ to be the restriction of $u$ in the complement of $\operatorname{exp}(N)$ and on each geodesic leaf of the interior of $\operatorname{exp}(N)$ define $\bar{u}$ to take the value that $u$ takes on the inward limit of that leaf.

Here's trouble. As $\bar{u}$ need not be smooth, I'd like to mollify it. Edit- I had a detail incorrect. A standard mollifier produces a function defined on compact subsets of the interior of $M$. So the mollifier has to be modified. One idea I'm following up on is varying the support of the mollifying function based on distance to $\partial M$. I'm skeptical, as varying the mollifying function will add another component to the gradient, but if it works I'll post as an answer.

I believe the following construction works. Analogous to the $n=1$ case from the comments, the idea is to approximate $\nabla u$ by compactly supported $\mathbf{g}_n$ and construct a suitable sequence $u_n$ such that each $\nabla u_n = \mathbf{g}_n.$ I will use some results from the following reference.

Maz’ya, Vladimir G., Sobolev spaces. With applications to elliptic partial differential equations. Transl. from the Russian by T. O. Shaposhnikova, Grundlehren der Mathematischen Wissenschaften 342. Berlin: Springer (ISBN 978-3-642-15563-5/hbk; 978-3-642-15564-2/ebook). xxviii, 866 p. (2011). ZBL1217.46002.

Definition (1.1.8): We say $\Omega \subset \Bbb R^n$ is starshapred with respect to a ball $B \subset \Omega$ if for all $x \in B$ and $y \in \Omega,$ we have the segment $[x,y] \subset \Omega.$

It is shown in 1.1.9 (Lemma 1) that if $\Omega$ satisfies the cone property, then it can be written as a union of finitely many domains which are startshared with respect to a ball. Note this also holds for $C^{0,1}$ domains, and it should include any reasonable definition of a piecewise smooth boundary.

Theorem (1.1.10): If $\Omega$ is bounded and starshapred with respect to a ball $B_{\delta}$ and $u \in W^{1,p}(\Omega),$ then we have $$ u(x) = \delta^{-n} \int_{B_{\delta}} \varphi\left(\frac{y}{\delta}\right) u(y) \,\mathrm{d}y + \int_{\Omega} \frac{\mathbf{f}(x;r,\theta)}{r^{n-1}} \cdot \nabla u(y)\,\mathrm{d}y $$ for almost every $x \in \Omega,$ where $r = |y-x|,$ $\theta = (y-x)r^{-1},$ $\varphi \in C^{\infty}_c(B_1),$ and $\mathbf{f}$ is smooth and $\Bbb R^n$-valued such that $$ |\mathbf{f}(x;r,\theta)| \leq c \left(\frac{\operatorname{diam}(\Omega)}{\delta}\right)^{n-1}.$$

I will omit a detailed proof (which is in the text), but the idea is to fix a mollifier $\varphi \in C^{\infty}_c(B_1)$ and set $\varphi_{\delta}(y) = \delta^{-n} \varphi(y/\delta).$ Then for $u$ smooth we have the Taylor representation $$ u(x) = u(y) + \int_0^1 (1-t) \sum_{i=1}^n D^iu(y+t(x-y))(x_i-y_i) \,\mathrm{d}t, $$ for all $x \in \Omega$ and $y \in B_{\delta}.$ Multiplying by $\varphi_{\delta}(y)$ and integrating over $y \in B_{\delta}$ gives $$ u(x) = \delta^{-n} \int_{B_{\delta}} \varphi\left(\frac{y}{\delta}\right) u(y) \,\mathrm{d}y + \int_{B_{\delta}} \varphi_{\delta}(y)\int_0^1 (1-t) \sum_{i=1}^n D^iu(y+t(x-y))(x_i-y_i) \,\mathrm{d}t \,\mathrm{d}y, $$ and the second term can be made to be of the right form by switching to $(r,\theta)$ coordinates. If $\delta=1$ one obtains $$ f_i(x;r,\theta) = - \theta_i \int_r^{\infty} \varphi(x+ \rho\theta) \rho^{n-1} \,\mathrm{d}\rho, $$ which one can check satisfies the desired properties.

Using the above we can prove the reslt in the case when $M = \Omega$ is starsharped with respect to some ball $B_{\delta}.$ Let $\mathbf{g}_n \in C^{\infty}_c(\Omega;\Bbb R^n)$ such that $\mathbf{g}_n \to \nabla u$ in $L^2.$ Then define $$ u_n(x) = \delta^{-n} \int_{B_{\delta}} \varphi\left(\frac{y}{\delta}\right) u(y) \,\mathrm{d}y + \int_{\Omega} \frac{\mathbf{f}(x;r,\theta)}{r^{n-1}} \cdot \mathbf{g}_n(y)\,\mathrm{d}y. $$ A calculation gives $\nabla u_n = \mathbf{g}_n,$ so we have $L^2$ convergence of the derivatives. For $L^2$ convergence this is essentially proved in 1.1.11, but since the text is light on the details I will include a proof. Applying Hölder and Fubini we have $$\begin{split} \int_{\Omega} |u(x)-u_n(x)|^2 \,\mathrm{d} x &=\int_{\Omega} \left| \int_{\Omega} \frac{\mathbf{f}(x;r,\theta)}{r^{n-1}} \cdot (\nabla u(y) - \mathbf{g}_n(y))\,\mathrm{d}y \right|^2 \,\mathrm{d} x \\ &\leq\int_{\Omega}\int_{\Omega} \frac1{r^{1-n}} \,\mathrm{d} y \int_{\Omega} \left|\mathbf{f}(x;r,\theta) \cdot (\nabla u(y) - \mathbf{g}_n(y))\right|^2 \frac1{r^{n-1}}\,\mathrm{d}y \,\mathrm{d} x \\ &\leq C \int_{\Omega} \left(\int_{\Omega} |x-y|^{1-n} \,\mathrm{d} x\right)^2 |\nabla u(y) - \mathbf{g}_n(y)|^2 \,\mathrm{d}y, \end{split}$$ which converges to zero as $n \to \infty$ noting that $\int_{\Omega} |x-y|^{1-n} \,\mathrm{d} x \leq C(\Omega)$ uniformly in $y.$ Hence it follows that $u_n \to u$ in $H^1(\Omega).$

For the general case we can using a patching argument; let $\{\Omega_i\}$ be a finite covering of $M$ by subsets which can be identified by means of a diffeomorphism with a starshaped domain in $\Bbb R^n.$ Let $\{\psi_i\}$ be a subordinate partition of unity, approximate $u$ on each $\Omega_i$ by $u_n^{(i)}$ as above and set $u_n = \sum_i \psi_i u_n^{(i)}.$ Note this is only a sketch as it depends on the precise formulation of a piecewise smooth boundary, but I will leave you to fill in the details.