Idea of the proof: Fix a standard radial [mollifier](/page/Mollifier) $\eta \in C_c^\infty(\mathbb{R}^n)$ with $\int_{\mathbb{R}^n}\eta\, d\mathcal{L}^n = 1$, put $\eta_\varepsilon(x) := \varepsilon^{-n}\eta(x/\varepsilon)$, and define $u^\varepsilon := \eta_\varepsilon * u$ on $U_\varepsilon$. We compute $u^\varepsilon(x)$ as a chain of equalities, explicitly marking each use of coarea, constancy on spheres, and the normalization.

Let $x \in U_\varepsilon$; then \begin{align*} u^\varepsilon(x) &= \int_U \eta_\varepsilon(x-y)\, u(y)\, d\mathcal{L}^n(y) && \text{(definition of convolution on $U_\varepsilon$)} \\[3pt] &= \frac{1}{\varepsilon^n} \int_U \eta\!\left(\frac{x-y}{\varepsilon}\right) u(y)\, d\mathcal{L}^n(y) && \text{(scaling of the mollifier)} \\[3pt] &= \frac{1}{\varepsilon^n} \int_{B(x,\varepsilon)} \eta\!\left(\frac{|x-y|}{\varepsilon}\right) u(y)\, d\mathcal{L}^n(y) && \text{($\mathrm{supp}\,\eta \subset B(0,1)$ $\Rightarrow$ $|x-y|<\varepsilon$)} \\[3pt] &= \frac{1}{\varepsilon^n} \int_{0}^{\varepsilon} \left( \int_{\partial B(x,r)} \eta\!\left(\frac{|x-y|}{\varepsilon}\right) u(y)\, d\mathcal{H}^{n-1}(y) \right) dr && \text{(coarea: $y \mapsto r=|x-y|$, $|\nabla r|=1$ a.e.)} \\[3pt] &= \frac{1}{\varepsilon^n} \int_{0}^{\varepsilon} \eta\!\left(\frac{r}{\varepsilon}\right) \left( \int_{\partial B(x,r)} u(y)\, d\mathcal{H}^{n-1}(y) \right) dr && \text{(on $\partial B(x,r)$, $|x-y|=r$ so $\eta(\tfrac{|x-y|}{\varepsilon})$ is constant in $y$)} \\[3pt] &= \frac{1}{\varepsilon^n} \int_{0}^{\varepsilon} \eta\!\left(\frac{r}{\varepsilon}\right)\, \mathcal{H}^{n-1}(\partial B(x,r))\, u(x)\, dr && \text{(apply (MV) to $\partial B(x,r) \subset U$)} \\[3pt] &= \int_{0}^{\varepsilon} \left( \int_{\partial B(x,r)} \eta_\varepsilon(x-y)\, d\mathcal{H}^{n-1}(y) \right) u(x)\, dr && \text{(rewrite: $\tfrac{1}{\varepsilon^n}\eta(\tfrac{r}{\varepsilon})\, \mathcal{H}^{n-1}(\partial B(x,r))$}\\[-2pt] &&& \text{equals $\int_{\partial B(x,r)} \eta_\varepsilon(x-y)\, d\mathcal{H}^{n-1}(y)$ since the integrand is constant)} \\[3pt] &= u(x) \int_{0}^{\varepsilon} \left( \int_{\partial B(0,r)} \eta_\varepsilon(z)\, d\mathcal{H}^{n-1}(z) \right) dr && \text{(change variables $z := x-y$; translation invariance of $\mathcal{H}^{n-1}$)} \\[3pt] &= u(x) \int_{B(0,\varepsilon)} \eta_\varepsilon(z)\, d\mathcal{L}^n(z) && \text{(coarea again: surface-to-volume for $\eta_\varepsilon$)} \\[3pt] &= u(x) \cdot 1 && \text{(normalization $\int_{\mathbb{R}^n} \eta_\varepsilon\, d\mathcal{L}^n = 1$)} \\[3pt] &= u(x). \end{align*}

Thus $u^\varepsilon \equiv u$ on $U_\varepsilon$. Since $u^\varepsilon \in C^\infty(U_\varepsilon)$ (convolution with $C_c^\infty$ kernel), $u$ is $C^\infty$ at every interior point (choose $\varepsilon < \mathrm{dist}(x,\partial U)$). Therefore $u \in C^\infty(U)$.