[proofplan] We prove $u \in H^2(V)$ by showing that the difference quotients $D_k^h \nabla u$ are uniformly bounded in $L^2(V)$, which by the characterization of Sobolev spaces via difference quotients implies that the second weak derivatives exist and lie in $L^2(V)$. The key tool is the test function $v = -D_k^{-h}(\zeta^2 D_k^h u)$, where $\zeta$ is a smooth cutoff supported in $U$ with $\zeta \equiv 1$ on $V$. Substituting into the weak formulation and using discrete integration by parts transfers the difference operator onto the coefficients. Uniform ellipticity bounds the principal term from below, and Cauchy's inequality with $\varepsilon$ absorbs the remainder terms. [/proofplan] [step:Set up the cutoff function and difference quotient test function] Fix $V \subset\subset U$ and choose a cutoff function $\zeta \in C_c^\infty(U)$ satisfying $\zeta \equiv 1$ on $V$, $0 \le \zeta \le 1$, and $\zeta \equiv 0$ near $\partial U$. Since $u \in H^1(U)$ and we do not yet know whether second derivatives exist, we use difference quotients. For $k \in \{1, \dots, n\}$ and $0 < |h| < \operatorname{dist}(V, \partial U)$, define: \begin{align*} D_k^h u(x) := \frac{u(x + h e_k) - u(x)}{h}, \end{align*} where $e_k$ is the $k$-th standard basis vector. Define the test function: \begin{align*} v := -D_k^{-h}(\zeta^2 D_k^h u). \end{align*} Since $u \in H^1(U)$ and $\zeta$ has compact support in $U$, the function $v$ belongs to $H^1_0(U)$ and is a valid test function in the weak formulation. [guided] Why use difference quotients instead of directly differentiating? We only know $u \in H^1(U)$, so $\nabla u \in L^2(U)$ but we have no information about second derivatives. Difference quotients $D_k^h u$ approximate $\partial_{x_k} u$ and are well-defined for $H^1$ functions. Concretely, the difference quotient is: \begin{align*} D_k^h u(x) = \frac{u(x + h e_k) - u(x)}{h}, \end{align*} and the key fact is that if $\|D_k^h g\|_{L^2(V)}$ remains bounded uniformly in $h$ for some $g \in L^2$, then $\partial_{x_k} g$ exists as a weak derivative in $L^2(V)$. Applying this characterization to $g = \partial_{x_i} u$ for each $i$: if we can show that $\|D_k^h \nabla u\|_{L^2(V)}$ remains bounded as $h \to 0$, the characterization of Sobolev spaces via difference quotients guarantees that $\partial_{x_k} \nabla u \in L^2(V)$, which is exactly the second-derivative regularity we seek. Why the test function $v = -D_k^{-h}(\zeta^2 D_k^h u)$? The cutoff $\zeta$ localizes the argument to $V$, while the factor $\zeta^2$ (rather than $\zeta$) ensures that the absorption step works — when the product rule is applied to $\partial_{x_j}(\zeta^2 D_k^h u)$, the cross term produces $\zeta |\nabla \zeta|$ rather than $|\nabla \zeta|^2$, which is critical for the Cauchy inequality to close. The backward difference $D_k^{-h}$ is placed on the outside so that the discrete integration by parts identity \begin{align*} \int_U \phi \, D_k^{-h} \psi \, d\mathcal{L}^n = -\int_U (D_k^h \phi) \, \psi \, d\mathcal{L}^n \end{align*} transfers the difference operator onto the coefficients in the bilinear form, producing terms involving $D_k^h(a_{ij} \partial_{x_i} u)$. This is the discrete analogue of integrating by parts twice in the classical regularity argument. [/guided] [/step] [step:Substitute into the weak formulation and apply discrete integration by parts] The weak formulation states $B[u, v] = (f, v)_{L^2(U)}$. Substituting the test function $v = -D_k^{-h}(\zeta^2 D_k^h u)$, we use the discrete integration by parts identity: \begin{align*} \int_U \phi\, D_k^{-h} \psi \, d\mathcal{L}^n = -\int_U (D_k^h \phi)\, \psi \, d\mathcal{L}^n, \end{align*} which holds whenever one factor has compact support. This shifts $D_k^{-h}$ from the test function onto the coefficients, yielding: \begin{align*} \sum_{i,j=1}^n \int_U D_k^h(a_{ij}\, \partial_{x_i} u)\, \partial_{x_j}(\zeta^2 D_k^h u) \, d\mathcal{L}^n = \int_U f\, D_k^{-h}(\zeta^2 D_k^h u) \, d\mathcal{L}^n + \text{lower-order terms}. \end{align*} [/step] [step:Extract the principal term using the difference quotient product rule] Apply the product rule for difference quotients, $D_k^h(fg) = f^h D_k^h g + (D_k^h f) g$, where $f^h(x) := f(x + he_k)$. The principal contribution on the left-hand side becomes: \begin{align*} \sum_{i,j=1}^n \int_U a_{ij}^h\, (D_k^h \partial_{x_i} u)\, (D_k^h \partial_{x_j} u)\, \zeta^2 \, d\mathcal{L}^n + \text{Rem}, \end{align*} where $\text{Rem}$ contains terms involving $D_k^h a_{ij}$ (which are bounded since $a_{ij} \in C^1(U)$) multiplied by at most one factor of $D_k^h \nabla u$. By uniform ellipticity applied to $a_{ij}^h$ (which satisfies the same ellipticity bound $\theta$ since the shifted coefficients inherit the pointwise condition): \begin{align*} \sum_{i,j=1}^n \int_U a_{ij}^h\, (D_k^h \partial_{x_i} u)\, (D_k^h \partial_{x_j} u)\, \zeta^2 \, d\mathcal{L}^n \ge \theta \int_U \zeta^2 |D_k^h \nabla u|^2 \, d\mathcal{L}^n. \end{align*} [/step] [step:Absorb remainder terms via Cauchy's inequality with $\varepsilon$] All remainder and right-hand side terms contain at most one factor of $\zeta D_k^h \nabla u$ paired with controlled quantities ($\nabla u$, $u$, $f$, bounded coefficient derivatives). Applying Cauchy's inequality with $\varepsilon$: \begin{align*} \theta \int_U \zeta^2 |D_k^h \nabla u|^2 \, d\mathcal{L}^n \le \varepsilon \int_U \zeta^2 |D_k^h \nabla u|^2 \, d\mathcal{L}^n + C(\varepsilon)\bigl(\|\nabla u\|_{L^2(U)}^2 + \|f\|_{L^2(U)}^2\bigr). \end{align*} Choosing $\varepsilon = \theta/2$ and absorbing: \begin{align*} \frac{\theta}{2} \int_U \zeta^2 |D_k^h \nabla u|^2 \, d\mathcal{L}^n \le C\bigl(\|u\|_{H^1(U)}^2 + \|f\|_{L^2(U)}^2\bigr). \end{align*} [guided] Why can the remainder terms be absorbed? Each remainder term has the schematic form $\int (\text{bounded}) \cdot |\zeta D_k^h \nabla u| \cdot |\text{something in } L^2|$. Cauchy's inequality $ab \le \varepsilon a^2 + \frac{1}{4\varepsilon} b^2$ puts $\varepsilon$ in front of the $|\zeta D_k^h \nabla u|^2$ term — which is exactly the term we want to bound — and $\frac{1}{4\varepsilon}$ in front of a controllable $L^2$ quantity. Choosing $\varepsilon = \theta/2$ ensures that the $\varepsilon$-contribution is strictly less than the $\theta$ on the left-hand side, allowing absorption. After absorption: \begin{align*} \frac{\theta}{2} \int_U \zeta^2 |D_k^h \nabla u|^2 \, d\mathcal{L}^n \le C\bigl(\|u\|_{H^1(U)}^2 + \|f\|_{L^2(U)}^2\bigr). \end{align*} The constant $C$ depends on the $C^1$ norms of the coefficients, $n$, and $\theta$, but not on $h$. [/guided] [/step] [step:Conclude $u \in H^2(V)$ via the difference quotient characterization] Since $\zeta \equiv 1$ on $V$: \begin{align*} \int_V |D_k^h \nabla u|^2 \, d\mathcal{L}^n \le \frac{2C}{\theta}\bigl(\|u\|_{H^1(U)}^2 + \|f\|_{L^2(U)}^2\bigr). \end{align*} This bound is uniform in $h$ for all $0 < |h| < \operatorname{dist}(V, \partial U)$ and all $k \in \{1, \dots, n\}$. By the characterization of [Sobolev spaces](/page/Sobolev%20Space) via difference quotients, the weak derivatives $\partial_{x_k} \nabla u$ exist and belong to $L^2(V)$. Since this holds for every $k$, all second-order weak derivatives $\partial_{x_k} \partial_{x_i} u \in L^2(V)$, so $u \in H^2(V)$. The estimate $\|u\|_{H^2(V)} \le C(\|f\|_{L^2(U)} + \|u\|_{L^2(U)})$ follows from the uniform bound above together with the interior $H^1$ estimate (which bounds $\|u\|_{H^1(U)}$ by $C(\|f\|_{L^2(U)} + \|u\|_{L^2(U)})$). [/step]