[proofplan]
We first record a local $H^1$ estimate, obtained by testing the weak equation with a cutoff multiple of $u$. We then fix a coordinate direction and test the equation with a translated difference-quotient expression, which turns the leading part into a coercive integral involving $\nabla D_{k,h}u$. The coefficient-difference, cutoff, and forcing terms are estimated by Cauchy-Schwarz and the elementary Young inequality, with the local $H^1$ estimate controlling $\nabla u$. The resulting bound is uniform in the difference-quotient scale $h$, so the difference-quotient criterion gives all second weak derivatives on $V$.
[/proofplan]
[step:Extend the weak formulation to compactly supported $H^1$ test functions]
Let $\phi\in H^1_0(U)$ have compact support in $U$. Choose an [open set](/page/Open%20Set) $W\subset\mathbb{R}^n$ such that $\operatorname{supp}\phi\subset W\subset\subset U$. Since $a_{ij}\in L^\infty(W)$, $u\in H^1(W)$, and $f\in L^2(W)$, the functional $\Lambda_W:H^1_0(W)\to\mathbb{R}$ defined by
\begin{align*}
\Lambda_W(\psi):=\sum_{i,j=1}^n \int_W a_{ij}(x)\,\partial_{x_j}u(x)\,\partial_{x_i}\psi(x)\,d\mathcal{L}^n(x)-\int_W f(x)\psi(x)\,d\mathcal{L}^n(x)
\end{align*}
is continuous. Indeed, Cauchy-Schwarz gives continuity of the first term from the $L^\infty(W)$ coefficient bounds and $\nabla u\in L^2(W)$, and continuity of the second term from $f\in L^2(W)$. By the density of smooth compactly supported functions in $H^1_0(W)$, $C_c^\infty(W)$ is dense in $H^1_0(W)$. The weak equation says $\Lambda_W(\psi)=0$ for every $\psi\in C_c^\infty(W)$, hence by continuity $\Lambda_W(\phi)=0$. Therefore
\begin{align*}
\sum_{i,j=1}^n \int_U a_{ij}(x)\,\partial_{x_j}u(x)\,\partial_{x_i}\phi(x)\,d\mathcal{L}^n(x)=\int_U f(x)\phi(x)\,d\mathcal{L}^n(x)
\end{align*}
holds for every $\phi\in H^1_0(U)$ with compact support in $U$.
[/step]
[step:Derive a local $H^1$ estimate from a cutoff test function]
Let $W_0,W\subset\mathbb{R}^n$ be open with $W_0\subset\subset W\subset\subset U$. Choose $\eta\in C_c^\infty(W)$ with $0\le \eta\le 1$ and $\eta=1$ on a neighbourhood of $W_0$. Define
\begin{align*}
A_W := \left(\sum_{i,j=1}^n \|a_{ij}\|_{L^\infty(W)}^2\right)^{1/2}
\end{align*}
and
\begin{align*}
M_\eta := \|\nabla \eta\|_{L^\infty(W)}.
\end{align*}
Using $\phi=\eta^2u$ in the weak formulation gives
\begin{align*}
\sum_{i,j=1}^n \int_W a_{ij}(x)\,\partial_{x_j}u(x)\,\partial_{x_i}(\eta^2u)(x)\,d\mathcal{L}^n(x)=\int_W f(x)\eta(x)^2u(x)\,d\mathcal{L}^n(x).
\end{align*}
Since $\partial_{x_i}(\eta^2u)=\eta^2\partial_{x_i}u+2\eta(\partial_{x_i}\eta)u$, ellipticity gives
\begin{align*}
\theta\int_W \eta(x)^2|\nabla u(x)|^2\,d\mathcal{L}^n(x)\le \left|\int_W f(x)\eta(x)^2u(x)\,d\mathcal{L}^n(x)\right|+2A_WM_\eta\left(\int_W \eta(x)^2|\nabla u(x)|^2\,d\mathcal{L}^n(x)\right)^{1/2}\|u\|_{L^2(W)}.
\end{align*}
Applying Cauchy-Schwarz to the forcing term and the elementary inequality $2ab\le \varepsilon a^2+\varepsilon^{-1}b^2$ with $\varepsilon=\theta/2$, we obtain
\begin{align*}
\|\nabla u\|_{L^2(W_0)}^2 \le C_1\left(\|f\|_{L^2(W)}^2+\|u\|_{L^2(W)}^2\right),
\end{align*}
where $C_1$ depends only on $\theta$, $A_W$, and $M_\eta$.
[/step]
[step:Choose a cutoff and define the translated difference quotients]
Fix open sets $V\subset\subset U'\subset\subset U$. Choose $\zeta\in C_c^\infty(U')$ such that $0\le \zeta\le 1$ and $\zeta=1$ on a neighbourhood of $V$. Let $K=\operatorname{supp}\zeta$, and choose $\rho>0$ such that
\begin{align*}
K+\{te_k:|t|<\rho,\ 1\le k\le n\}\subset U'.
\end{align*}
For $1\le k\le n$ and $0<|h|<\rho$, define the translated-safe open set
\begin{align*}
U'_{k,h}:=\{x\in U':x+he_k\in U'\}.
\end{align*}
Define the translation operator $\tau_{k,h}:H^1(U')\to H^1_{\mathrm{loc}}(U'_{k,h})$ by
\begin{align*}
(\tau_{k,h}v)(x):=v(x+he_k),\qquad x\in U'_{k,h},
\end{align*}
and define the difference-quotient operator $D_{k,h}:H^1(U')\to L^2_{\mathrm{loc}}(U'_{k,h})$ by
\begin{align*}
D_{k,h}v(x):=\frac{v(x+he_k)-v(x)}{h},\qquad x\in U'_{k,h}.
\end{align*}
Set
\begin{align*}
w_h:=D_{k,h}u.
\end{align*}
This function is used only on sets where the translate $x+he_k$ lies in $U'$. More precisely, on the compact set $K=\operatorname{supp}\zeta$ and on its $h$-translate, the choice of $\rho$ ensures that $w_h\in H^1$ locally. The product $\zeta^2w_h$ is first defined on $K$ and then extended by zero to $U'$; with this convention $\zeta^2w_h\in H^1_0(U')$.
[/step]
[step:Test with the backward difference quotient and isolate the coercive term]
Use the admissible [test function](/page/Test%20Function)
\begin{align*}
\phi_h:=-D_{k,-h}(\zeta^2w_h).
\end{align*}
Indeed, $\operatorname{supp}(\zeta^2w_h)\subset K$. For $0<|h|<\rho$, the definition of $\rho$ gives $K-he_k\subset U'$ and $K+he_k\subset U'$. With the convention
\begin{align*}
D_{k,-h}\psi(x)=\frac{\psi(x-he_k)-\psi(x)}{-h},
\end{align*}
the support of $D_{k,-h}(\zeta^2w_h)$ is contained in $K\cup(K+he_k)$, because $D_{k,-h}(\zeta^2w_h)(x)$ can be nonzero only when $x\in K$ or $x-he_k\in K$. Thus $\phi_h\in H^1_0(U')$. The same support statement holds after applying $\partial_{x_i}$, since distributional differentiation cannot enlarge support.
The discrete integration-by-parts identity
\begin{align*}
\int_{U'} g(x)D_{k,-h}\psi(x)\,d\mathcal{L}^n(x)=-\int_{U'} D_{k,h}g(x)\psi(x)\,d\mathcal{L}^n(x)
\end{align*}
holds whenever both $\psi$ and its relevant translate have compact support in $U'$. Here this applies with $\psi=\partial_{x_i}(\zeta^2w_h)$, because its support is contained in $K$ and its translate is contained in $K+he_k\subset U'$. Applying it to the left-hand side of the weak equation gives
\begin{align*}
\sum_{i,j=1}^n \int_{U'} D_{k,h}\!\left(a_{ij}\partial_{x_j}u\right)(x)\,\partial_{x_i}(\zeta^2w_h)(x)\,d\mathcal{L}^n(x)=-\int_{U'} f(x)D_{k,-h}(\zeta^2w_h)(x)\,d\mathcal{L}^n(x).
\end{align*}
Using the product rule for difference quotients,
\begin{align*}
D_{k,h}(a_{ij}\partial_{x_j}u)=(\tau_{k,h}a_{ij})\,\partial_{x_j}w_h+(D_{k,h}a_{ij})\,\partial_{x_j}u.
\end{align*}
The principal part therefore contains
\begin{align*}
\sum_{i,j=1}^n \int_{U'} \zeta(x)^2(\tau_{k,h}a_{ij})(x)\,\partial_{x_j}w_h(x)\,\partial_{x_i}w_h(x)\,d\mathcal{L}^n(x).
\end{align*}
By ellipticity applied at the point $x+he_k$, this term is bounded below by
\begin{align*}
\theta\int_{U'} \zeta(x)^2|\nabla w_h(x)|^2\,d\mathcal{L}^n(x).
\end{align*}
[guided]
The point of this test function is to move one discrete derivative from the test function onto the coefficient-gradient product $a_{ij}\partial_{x_j}u$. Fix a coordinate direction $e_k$ and a scale $h$ with $0<|h|<\rho$. We use
\begin{align*}
\phi_h:=-D_{k,-h}(\zeta^2D_{k,h}u).
\end{align*}
The product $\zeta^2D_{k,h}u$ is defined on $K$, where the translate $x+he_k$ lies in $U'$, and then extended by zero to $U'$. Since $\operatorname{supp}(\zeta^2D_{k,h}u)\subset K$ and, for $0<|h|<\rho$, the sets $K-he_k$ and $K+he_k$ are contained in $U'$, this zero extension belongs to $H^1_0(U')$. With the convention $D_{k,-h}\psi(x)=(\psi(x-he_k)-\psi(x))/(-h)$, the backward difference $D_{k,-h}(\zeta^2D_{k,h}u)$ can be nonzero only when $x\in K$ or $x-he_k\in K$; hence its support is contained in $K\cup(K+he_k)$, a compact subset of $U'$. Differentiating does not enlarge support, so the same translated-support condition holds for $\partial_{x_i}D_{k,-h}(\zeta^2D_{k,h}u)$.
The weak equation gives
\begin{align*}
\sum_{i,j=1}^n \int_{U'} a_{ij}(x)\partial_{x_j}u(x)\partial_{x_i}\phi_h(x)\,d\mathcal{L}^n(x)=\int_{U'} f(x)\phi_h(x)\,d\mathcal{L}^n(x).
\end{align*}
Since weak derivatives commute with difference quotients on these compactly supported translated functions, we have
\begin{align*}
\partial_{x_i}\phi_h=-D_{k,-h}\partial_{x_i}(\zeta^2D_{k,h}u).
\end{align*}
The discrete integration-by-parts identity applies because $\partial_{x_i}(\zeta^2D_{k,h}u)$ and its required translate have compact support in $U'$. It says that a backward difference can be moved from one factor to the other as a forward difference, with a minus sign. The minus sign in the definition of $\phi_h$ was chosen so that the leading term has the positive coercive sign. Thus
\begin{align*}
\sum_{i,j=1}^n \int_{U'} D_{k,h}\!\left(a_{ij}\partial_{x_j}u\right)(x)\partial_{x_i}(\zeta^2D_{k,h}u)(x)\,d\mathcal{L}^n(x)=-\int_{U'} f(x)D_{k,-h}(\zeta^2D_{k,h}u)(x)\,d\mathcal{L}^n(x).
\end{align*}
Now expand the discrete derivative of the product. With $w_h=D_{k,h}u$,
\begin{align*}
D_{k,h}(a_{ij}\partial_{x_j}u)=(\tau_{k,h}a_{ij})\partial_{x_j}w_h+(D_{k,h}a_{ij})\partial_{x_j}u.
\end{align*}
The first summand is the desired one: it contains $\partial_{x_j}w_h$, hence second-order information about $u$. Since
\begin{align*}
\partial_{x_i}(\zeta^2w_h)=\zeta^2\partial_{x_i}w_h+2\zeta(\partial_{x_i}\zeta)w_h,
\end{align*}
the highest-order part is
\begin{align*}
\sum_{i,j=1}^n \int_{U'} \zeta(x)^2(\tau_{k,h}a_{ij})(x)\partial_{x_j}w_h(x)\partial_{x_i}w_h(x)\,d\mathcal{L}^n(x).
\end{align*}
For each $x$ in the integration region, the coefficient matrix is evaluated at $x+he_k$. The ellipticity hypothesis is valid at that point, so with $\xi=\nabla w_h(x)$ it yields
\begin{align*}
\sum_{i,j=1}^n(\tau_{k,h}a_{ij})(x)\partial_{x_j}w_h(x)\partial_{x_i}w_h(x)\ge \theta|\nabla w_h(x)|^2.
\end{align*}
Multiplying by $\zeta(x)^2$ and integrating gives the coercive lower bound
\begin{align*}
\theta\int_{U'} \zeta(x)^2|\nabla w_h(x)|^2\,d\mathcal{L}^n(x).
\end{align*}
[/guided]
[/step]
[step:Bound the cutoff, coefficient, and forcing errors uniformly in $h$]
Define
\begin{align*}
A:=\left(\sum_{i,j=1}^n \|a_{ij}\|_{L^\infty(U')}^2\right)^{1/2}
\end{align*}
and
\begin{align*}
B:=\left(\sum_{i,j,k=1}^n \|\partial_{x_k}a_{ij}\|_{L^\infty(U')}^2\right)^{1/2}.
\end{align*}
Also set $M:=\|\nabla\zeta\|_{L^\infty(U')}$. Define the compact translated cutoff neighbourhood
\begin{align*}
K_\rho:=\{x+te_k:x\in K,\ |t|\le \rho,\ 1\le k\le n\}\subset U'.
\end{align*}
Since $a_{ij}\in C^1(U')$, the [mean value theorem](/theorems/186) along the line segment $x+the_k$, which is contained in $K_\rho$ whenever $x\in K$ and $|h|<\rho$, gives
\begin{align*}
|D_{k,h}a_{ij}(x)|\le \|\partial_{x_k}a_{ij}\|_{L^\infty(U')}.
\end{align*}
The cutoff error from differentiating $\zeta^2$ is supported in $K$ and is bounded by
\begin{align*}
2AM\|\zeta\nabla w_h\|_{L^2(U')}\|w_h\|_{L^2(K)}.
\end{align*}
The coefficient-difference terms are supported in $K$ and are bounded by
\begin{align*}
B\|\nabla u\|_{L^2(K_\rho)}\left(\|\zeta\nabla w_h\|_{L^2(U')}+2M\|w_h\|_{L^2(K)}\right).
\end{align*}
For the forcing term, the standard difference-quotient estimate $\|D_{k,-h}\psi\|_{L^2(U')}\le \|\partial_{x_k}\psi\|_{L^2(U')}$ for $\psi\in H^1_0(U')$ gives
\begin{align*}
\left|\int_{U'} f(x)D_{k,-h}(\zeta^2w_h)(x)\,d\mathcal{L}^n(x)\right|\le \|f\|_{L^2(U')}\left(\|\zeta\nabla w_h\|_{L^2(U')}+2M\|w_h\|_{L^2(K)}\right).
\end{align*}
Finally, the one-dimensional difference-quotient estimate applied on the line segments $\{x+te_k:0\le t\le h\}$ for $x\in K$ gives
\begin{align*}
\|w_h\|_{L^2(K)}\le \|\partial_{x_k}u\|_{L^2(K_\rho)}.
\end{align*}
Combining these bounds with the coercive estimate and applying $2ab\le \varepsilon a^2+\varepsilon^{-1}b^2$ with $\varepsilon>0$ chosen small relative to $\theta$, we obtain
\begin{align*}
\|\zeta\nabla D_{k,h}u\|_{L^2(U')}^2\le C_2\left(\|f\|_{L^2(U')}^2+\|\nabla u\|_{L^2(K_\rho)}^2\right)
\end{align*}
for every $1\le k\le n$ and every $0<|h|<\rho$, where $C_2$ depends only on $n$, $\theta$, $A$, $B$, and $M$.
[/step]
[step:Use the local $H^1$ estimate to remove the gradient of $u$ from the right-hand side]
Apply the local $H^1$ estimate from the earlier step with $W_0$ chosen as an open neighbourhood of $K_\rho$ satisfying $K_\rho\subset W_0\subset\subset U'$ and with $W=U'$, using a cutoff whose geometry depends only on the containment $K_\rho\subset\subset U'$. This gives
\begin{align*}
\|\nabla u\|_{L^2(K_\rho)}^2\le C_3\left(\|f\|_{L^2(U')}^2+\|u\|_{L^2(U')}^2\right).
\end{align*}
Substituting this into the uniform difference-quotient estimate yields
\begin{align*}
\|\zeta\nabla D_{k,h}u\|_{L^2(U')}^2\le C_4\left(\|f\|_{L^2(U')}^2+\|u\|_{L^2(U')}^2\right)
\end{align*}
for all $1\le k\le n$ and all $0<|h|<\rho$. Since $\zeta=1$ on a neighbourhood of $V$, this implies
\begin{align*}
\|\nabla D_{k,h}u\|_{L^2(V)}^2\le C_4\left(\|f\|_{L^2(U')}^2+\|u\|_{L^2(U')}^2\right).
\end{align*}
[/step]
[step:Pass from uniform difference quotients to second weak derivatives]
Fix $1\le i,k\le n$. The preceding estimate gives a constant $C_5>0$ such that
\begin{align*}
\|D_{k,h}(\partial_{x_i}u)\|_{L^2(V)}\le C_5\left(\|f\|_{L^2(U')}+\|u\|_{L^2(U')}\right)
\end{align*}
for every sufficiently small nonzero $h$. Let $\varphi\in C_c^\infty(V)$. Choose $\delta_\varphi>0$ such that $\operatorname{supp}\varphi+te_k\subset V$ whenever $|t|<\delta_\varphi$. For $0<|h|<\min\{\rho,\delta_\varphi\}$, the translated supports remain inside $V$, so the discrete integration-by-parts identity gives
\begin{align*}
\int_V D_{k,h}(\partial_{x_i}u)(x)\varphi(x)\,d\mathcal{L}^n(x)=-\int_V \partial_{x_i}u(x)D_{k,-h}\varphi(x)\,d\mathcal{L}^n(x).
\end{align*}
As $h\to 0$, $D_{k,-h}\varphi\to \partial_{x_k}\varphi$ uniformly on $V$. By weak compactness in reflexive Banach spaces, applied to the reflexive [Hilbert space](/page/Hilbert%20Space) $L^2(V)$, the bounded family $D_{k,h}(\partial_{x_i}u)$ has a weakly convergent subsequence in $L^2(V)$; write its weak limit as $g_{ik}\in L^2(V)$. Passing to the limit along that subsequence yields
\begin{align*}
\int_V g_{ik}(x)\varphi(x)\,d\mathcal{L}^n(x)=-\int_V \partial_{x_i}u(x)\partial_{x_k}\varphi(x)\,d\mathcal{L}^n(x).
\end{align*}
Thus $g_{ik}=\partial_{x_k}\partial_{x_i}u$ in the weak sense on $V$. By weak [lower semicontinuity of the norm](/theorems/215) in a Hilbert space,
\begin{align*}
\|\partial_{x_k}\partial_{x_i}u\|_{L^2(V)}=\|g_{ik}\|_{L^2(V)}\le C_5\left(\|f\|_{L^2(U')}+\|u\|_{L^2(U')}\right).
\end{align*}
This is the standard difference-quotient characterization of Sobolev weak derivatives: a uniform $L^2$ bound on directional difference quotients of $\partial_{x_i}u$ identifies an $L^2$ [weak derivative](/page/Weak%20Derivative) in the $x_k$ direction. Since $i$ and $k$ were arbitrary, all second weak derivatives of $u$ belong to $L^2(V)$.
[/step]
[step:Assemble the $H^2$ estimate on $V$]
The local $H^1$ estimate gives
\begin{align*}
\|u\|_{H^1(V)}\le C_6\left(\|f\|_{L^2(U')}+\|u\|_{L^2(U')}\right).
\end{align*}
The previous step gives the same type of bound for every second weak derivative $\partial_{x_k}\partial_{x_i}u$ on $V$. Summing the finitely many first- and second-order bounds and absorbing the dimension-dependent finite sum into the constant gives
\begin{align*}
\|u\|_{H^2(V)}\le C\left(\|f\|_{L^2(U')}+\|u\|_{L^2(U')}\right).
\end{align*}
The constant $C$ may be chosen in terms of the ellipticity constant, the dimension, the cutoff geometry of the fixed containment $V\subset\subset U'$, and the $L^\infty$ and $C^1$ bounds of the coefficients on the translated cutoff neighbourhood $K_\rho\subset\subset U'$. Since $K_\rho$ is determined by the chosen cutoff $\zeta$, the translation radius $\rho$, and the fixed containment $V\subset\subset U'$, this is exactly the stated dependence on the cutoff geometry and on coefficient bounds on the neighbourhood where the cutoff and all sufficiently small translated difference quotients are supported. This proves $u\in H^2(V)$ and the stated interior estimate.
[/step]