[proofplan] We prove boundary $H^2$ regularity by localizing near a boundary point, flattening the boundary via a $C^2$-diffeomorphism, and then separately establishing tangential and normal regularity of the transformed solution $v$ on the half-ball $B_1^+$. Tangential regularity is obtained by testing the weak formulation with difference quotients $v = -D_k^{-h}(\zeta^2 D_k^h v)$ in tangential directions $k < n$, which is valid because the flat boundary condition $v(y', 0) = 0$ is preserved by tangential translations. Normal regularity follows by solving the PDE for $\partial_{y_n y_n} v$, expressing it in terms of already-controlled tangential and mixed derivatives. A partition of unity and interior regularity complete the global result. [/proofplan] [step:Flatten the boundary near $x_0 \in \partial U$ via a $C^2$-diffeomorphism] Fix $x_0 \in \partial U$. Since $\partial U$ is $C^2$, there exist a radius $r > 0$ and a $C^2$-diffeomorphism: \begin{align*} \Phi: U \cap B(x_0, r) &\to B_1^+ := \{y \in B_1(0) : y_n > 0\} \end{align*} with $\Phi(x_0) = 0$ and $\Phi(\partial U \cap B(x_0, r)) \subset \{y_n = 0\}$. Let $\Psi := \Phi^{-1}$. Define the transformed function $v(y) := u(\Psi(y))$. Since $u \in H^1_0(U)$, we have $v \in H^1(B_1^+)$ with $v(y', 0) = 0$ on the flat boundary. [/step] [step:Transform the equation and verify ellipticity is preserved] By the chain rule, the weak formulation for $u$ transforms into a weak formulation for $v$ on $B_1^+$: \begin{align*} \int_{B_1^+} \sum_{k,l=1}^n A_{kl}(y)\, \partial_{y_k} v\, \partial_{y_l} \phi \, d\mathcal{L}^n(y) = \int_{B_1^+} \tilde{f}\, \phi \, d\mathcal{L}^n(y) + \text{lower-order terms}, \end{align*} for all $\phi \in H^1_0(B_1^+)$, where $\tilde{f} \in L^2(B_1^+)$ and the new coefficients are: \begin{align*} A_{kl}(y) = \sum_{i,j=1}^n a_{ij}(\Psi(y))\, \frac{\partial \Phi_k}{\partial x_i}(\Psi(y))\, \frac{\partial \Phi_l}{\partial x_j}(\Psi(y))\, J(y), \end{align*} with $J(y) = |\det D\Psi(y)|$. These coefficients belong to $C^1(\overline{B_1^+})$ since $a_{ij} \in C^1(\bar{U})$ and $\Phi \in C^2$. For any $\xi \in \mathbb{R}^n$, setting $\eta := (D\Phi)^\top \xi$: \begin{align*} \sum_{k,l=1}^n A_{kl}(y)\, \xi_k\, \xi_l = J(y) \sum_{i,j=1}^n a_{ij}(\Psi(y))\, \eta_i\, \eta_j \ge J(y)\, \theta\, |\eta|^2. \end{align*} Since $\Phi$ is a diffeomorphism, $D\Phi$ is non-singular, so $|\eta|^2 \ge c\, |\xi|^2$ for some $c > 0$. Hence the transformed operator is uniformly elliptic with constant $\theta' := \theta\, c\, \inf_{B_1^+} J > 0$. [/step] [step:Establish tangential regularity via difference quotients] For $k \in \{1, \dots, n-1\}$ (tangential directions), let $\zeta \in C_c^\infty(B_1^+)$ be a cutoff function. Define the test function: \begin{align*} \phi := -D_k^{-h}(\zeta^2 D_k^h v). \end{align*} Since $v(y', 0) = 0$ and tangential translation preserves the flat boundary ($v(y' + he_k, y_n)\big|_{y_n=0} = 0$ for $k < n$), the function $\phi$ vanishes on $\{y_n = 0\}$ and belongs to $H^1_0(B_1^+)$. Substituting $\phi$ into the transformed weak formulation and applying discrete integration by parts (shifting $D_k^{-h}$ onto the other side): \begin{align*} \sum_{r,s=1}^n \int_{B_1^+} D_k^h(A_{rs}\, \partial_{y_r} v)\, \partial_{y_s}(\zeta^2 D_k^h v) \, d\mathcal{L}^n(y) = \int_{B_1^+} \tilde{f}\, D_k^{-h}(\zeta^2 D_k^h v) \, d\mathcal{L}^n(y). \end{align*} This has the same structure as in the [interior regularity](/theorems/95) proof. The ellipticity of $A_{rs}$ and Cauchy's inequality with $\varepsilon$ yield: \begin{align*} \int_{B_1^+} \zeta^2\, |D_k^h \nabla v|^2 \, d\mathcal{L}^n(y) \le C\bigl(\|\tilde{f}\|_{L^2(B_1^+)}^2 + \|v\|_{H^1(B_1^+)}^2\bigr), \end{align*} uniformly in $h$. By the difference quotient characterization of Sobolev spaces, this implies $\partial_{y_k} \partial_{y_l} v \in L^2$ for all second derivatives involving at least one tangential index $k < n$. [guided] Why can we use tangential difference quotients at the boundary but not normal ones? The test function $\phi = -D_k^{-h}(\zeta^2 D_k^h v)$ requires that $D_k^h v$ vanishes on $\{y_n = 0\}$ so that $\phi \in H^1_0(B_1^+)$. For tangential directions $k < n$, the translation $y' \mapsto y' + he_k$ stays on the flat boundary, and the Dirichlet condition $v(y', 0) = 0$ is preserved: \begin{align*} D_k^h v(y', 0) = \frac{v(y' + he_k, 0) - v(y', 0)}{h} = \frac{0 - 0}{h} = 0. \end{align*} For the normal direction $k = n$, the translation $y_n \mapsto y_n + h$ moves the evaluation point off the boundary into the interior, where $v(y', h) \ne 0$ in general. Thus $D_n^h v(y', 0) = h^{-1} v(y', h)$, which does not vanish, and the test function $\phi$ would fail to satisfy the homogeneous Dirichlet condition required for membership in $H^1_0(B_1^+)$. After substitution into the weak formulation, the argument proceeds exactly as in the interior case. The product rule for difference quotients gives: \begin{align*} D_k^h(A_{rs}\, \partial_{y_r} v) = A_{rs}^h\, D_k^h \partial_{y_r} v + (D_k^h A_{rs})\, \partial_{y_r} v, \end{align*} which isolates the principal term involving $|D_k^h \nabla v|^2$. Uniform ellipticity of $A_{rs}^h$ (the shifted coefficients inherit the same ellipticity constant) bounds the principal term from below: \begin{align*} \sum_{r,s} \int_{B_1^+} A_{rs}^h\, D_k^h \partial_{y_r} v\, \zeta^2\, D_k^h \partial_{y_s} v \, d\mathcal{L}^n \ge \theta' \int_{B_1^+} \zeta^2\, |D_k^h \nabla v|^2 \, d\mathcal{L}^n. \end{align*} The remainder terms, involving bounded quantities $(D_k^h A_{rs})$ paired with at most one factor of $\zeta D_k^h \nabla v$, are absorbed using Cauchy's inequality with $\varepsilon = \theta'/2$, yielding the uniform bound that implies tangential second derivatives exist in $L^2$. [/guided] [/step] [step:Recover the normal derivative $\partial_{y_n y_n} v$ from the PDE] From the strong form of the transformed equation, isolate the pure normal second derivative: \begin{align*} A_{nn}\, \partial_{y_n y_n} v = \tilde{f} + \sum_{(r,s) \ne (n,n)} \partial_{y_r}(A_{rs}\, \partial_{y_s} v) - \sum_{r=1}^n B_r\, \partial_{y_r} v - C\, v. \end{align*} Dividing by $A_{nn}$ (which satisfies $A_{nn}(y) \ge \theta' > 0$ by ellipticity): \begin{align*} |\partial_{y_n y_n} v| \le C\bigl(|\tilde{f}| + \textstyle\sum_{(r,s) \ne (n,n)} |\partial_{y_r y_s} v| + |\nabla v| + |v|\bigr). \end{align*} The right-hand side involves only tangential/mixed second derivatives (controlled in the previous step), $\tilde{f} \in L^2$, and lower-order terms in $L^2$. Therefore $\partial_{y_n y_n} v \in L^2(B_1^+)$, and consequently $v \in H^2(B_1^+)$. [/step] [step:Transform back and assemble the global result] Since $\Phi$ is a $C^2$-diffeomorphism, the chain rule for second derivatives shows $u \in H^2(U \cap B(x_0, r))$. Covering $\partial U$ with finitely many such balls $\{B(x_0^{(j)}, r_j)\}_{j=1}^N$ and using the [Elliptic Interior Regularity](/theorems/95) theorem for points away from $\partial U$, we obtain $u \in H^2(U)$ with: \begin{align*} \|u\|_{H^2(U)} \le C\bigl(\|f\|_{L^2(U)} + \|u\|_{L^2(U)}\bigr). \end{align*} [/step]