[proofplan]
We first record the global Calderón-Zygmund estimate on $\mathbb{R}^n$: second derivatives are obtained from the Laplacian by second-order Riesz transforms, whose $L^p$ boundedness holds for $1<p<\infty$. We then localize by multiplying $u$ with a smooth cutoff supported in a ball and equal to $1$ on a smaller ball. The cutoff computation produces lower-order terms involving $\nabla u$ and $u$; the gradient term is absorbed by a local interpolation inequality and a standard shrinking-radius argument. Finally, finitely many balls covering $\overline V$ convert the local ball estimate into the stated estimate on $V$.
[/proofplan]
[step:Use the global Calderón-Zygmund estimate on compactly supported functions]
Let $w \in C_c^\infty(\mathbb{R}^n)$ be a real-valued smooth compactly supported function. For $1 \le i,j \le n$, define the Fourier multiplier $T_{ij}$ on Schwartz functions by
\begin{align*}
\widehat{T_{ij}g}(\xi)=m_{ij}(\xi)\hat g(\xi), \qquad m_{ij}(\xi)=\frac{\xi_i\xi_j}{|\xi|^2}\ \text{for}\ \xi \ne 0,\quad m_{ij}(0)=0.
\end{align*}
With the symmetric [Fourier transform](/page/Fourier%20Transform) convention, differentiation gives
\begin{align*}
\widehat{\partial_{x_i}\partial_{x_j}w}(\xi)=-\xi_i\xi_j\hat w(\xi)
\end{align*}
and
\begin{align*}
\widehat{\Delta w}(\xi)=-|\xi|^2\hat w(\xi).
\end{align*}
Hence
\begin{align*}
\widehat{\partial_{x_i}\partial_{x_j}w}(\xi)=m_{ij}(\xi)\widehat{\Delta w}(\xi).
\end{align*}
If $n\ge 2$, we invoke the Calderón-Zygmund theorem for second-order Riesz transforms on $\mathbb{R}^n$: for $1<p<\infty$, each multiplier $m_{ij}$ defines a bounded operator on $L^p(\mathbb{R}^n)$. This is the standard $L^p$ [boundedness theorem](/theorems/181) for the singular integral operators whose first-order multiplier identity is recorded in [Riesz Transform Multiplier](/theorems/3163). The hypotheses are satisfied because $1<p<\infty$ and $\Delta w\in L^p(\mathbb{R}^n)$ for $w\in C_c^\infty(\mathbb{R}^n)$. This external analytic input gives a constant $A_{n,p}>0$ such that
\begin{align*}
\|\partial_{x_i}\partial_{x_j}w\|_{L^p(\mathbb{R}^n)} \le A_{n,p}\|\Delta w\|_{L^p(\mathbb{R}^n)}.
\end{align*}
If $n=1$, then $i=j=1$ and $\Delta w=\partial_{x_1}\partial_{x_1}w$, so the same estimate holds with $A_{1,p}=1$ without invoking the singular-integral theorem.
Summing over $1 \le i,j \le n$ gives a constant $A=A(n,p)>0$ satisfying
\begin{align*}
\|D^2w\|_{L^p(\mathbb{R}^n)} \le A\|\Delta w\|_{L^p(\mathbb{R}^n)}.
\end{align*}
Let $w\in W^{2,p}(\mathbb{R}^n)$ have compact support. Choose $\theta\in C_c^\infty(\mathbb{R}^n)$ with $\theta=1$ on an open neighbourhood of $\operatorname{supp} w$, and let $(\eta_m)_{m=1}^\infty$ be a [standard mollifier](/page/Standard%20Mollifier) sequence on $\mathbb{R}^n$. Define $w_m:=\theta(\eta_m*w)$. Since convolution with a standard mollifier converges to the identity in $W^{2,p}(\mathbb{R}^n)$ and multiplication by the fixed smooth compactly supported function $\theta$ is continuous on $W^{2,p}(\mathbb{R}^n)$, we have $w_m\in C_c^\infty(\mathbb{R}^n)$ and $w_m\to w$ in $W^{2,p}(\mathbb{R}^n)$. In particular $D^2w_m\to D^2w$ and $\Delta w_m\to\Delta w$ in $L^p(\mathbb{R}^n)$, so passing to the limit in the displayed estimate gives the same estimate for every compactly supported $w \in W^{2,p}(\mathbb{R}^n)$.
[/step]
[step:Derive a cutoff estimate on concentric balls]
Let $B(x_0,R) \subset U'$ be an open ball, let $0<r<s<R$, and choose $\eta \in C_c^\infty(B(x_0,s))$ such that $0 \le \eta \le 1$ on $\mathbb{R}^n$, $\eta=1$ on $B(x_0,r)$, and
\begin{align*}
|\nabla \eta| \le \frac{C_1}{s-r}, \qquad |\Delta \eta| \le \frac{C_2}{(s-r)^2},
\end{align*}
where $C_1=C_1(n)>0$ and $C_2=C_2(n)>0$ are constants coming from a fixed radial cutoff profile. Define $w \in W^{2,p}(\mathbb{R}^n)$ by extending $\eta u$ by zero outside $B(x_0,s)$.
Since $\eta=1$ on $B(x_0,r)$, the Hessians agree there, so
\begin{align*}
\|D^2u\|_{L^p(B(x_0,r))} \le \|D^2w\|_{L^p(\mathbb{R}^n)}.
\end{align*}
Applying the global estimate to $w$ gives
\begin{align*}
\|D^2u\|_{L^p(B(x_0,r))} \le A\|\Delta(\eta u)\|_{L^p(B(x_0,s))}.
\end{align*}
The product rule in distributions gives
\begin{align*}
\Delta(\eta u)=\eta\Delta u+2\sum_{k=1}^n(\partial_{x_k}\eta)(\partial_{x_k}u)+u\Delta\eta.
\end{align*}
Because $-\Delta u=f$ in $\mathcal{D}'(U)$ and $B(x_0,s)\subset U$, this becomes
\begin{align*}
\Delta(\eta u)=-\eta f+2\sum_{k=1}^n(\partial_{x_k}\eta)(\partial_{x_k}u)+u\Delta\eta.
\end{align*}
Taking the $L^p(B(x_0,s))$ norm and using the displayed cutoff bounds yields
\begin{align*}
\|D^2u\|_{L^p(B(x_0,r))} \le A\|f\|_{L^p(B(x_0,s))}+\frac{2AC_1}{s-r}\|\nabla u\|_{L^p(B(x_0,s))}+\frac{AC_2}{(s-r)^2}\|u\|_{L^p(B(x_0,s))}.
\end{align*}
[guided]
The purpose of the cutoff is to convert a local problem into a compactly supported global problem, because the global Calderón-Zygmund estimate applies on $\mathbb{R}^n$. We choose radii $0<r<s<R$ and a function $\eta \in C_c^\infty(B(x_0,s))$ with $\eta=1$ on $B(x_0,r)$. The constants in the derivative bounds have the expected scale:
\begin{align*}
|\nabla \eta| \le \frac{C_1}{s-r}, \qquad |\Delta \eta| \le \frac{C_2}{(s-r)^2}.
\end{align*}
Define $w:\mathbb{R}^n \to \mathbb{R}$ by $w=\eta u$ on $B(x_0,s)$ and $w=0$ on $\mathbb{R}^n\setminus B(x_0,s)$. Since $\eta$ is compactly supported in $B(x_0,s)$ and $u \in W^{2,p}_{\mathrm{loc}}(U)$, we have $w \in W^{2,p}(\mathbb{R}^n)$ with compact support.
On $B(x_0,r)$ the cutoff equals $1$, so $w=u$ there and therefore $D^2w=D^2u$ there. This gives
\begin{align*}
\|D^2u\|_{L^p(B(x_0,r))} \le \|D^2w\|_{L^p(\mathbb{R}^n)}.
\end{align*}
The global estimate applies to $w$, so
\begin{align*}
\|D^2u\|_{L^p(B(x_0,r))} \le A\|\Delta w\|_{L^p(\mathbb{R}^n)}.
\end{align*}
Because $w$ is supported in $B(x_0,s)$, the right-hand side is
\begin{align*}
A\|\Delta(\eta u)\|_{L^p(B(x_0,s))}.
\end{align*}
Now we compute the localized Laplacian. The distributional product rule is valid because $\eta$ is smooth and $u \in W^{2,p}_{\mathrm{loc}}(U)$:
\begin{align*}
\Delta(\eta u)=\eta\Delta u+2\sum_{k=1}^n(\partial_{x_k}\eta)(\partial_{x_k}u)+u\Delta\eta.
\end{align*}
The equation $-\Delta u=f$ in $\mathcal{D}'(U)$ says exactly that $\Delta u=-f$ as a distribution on every open subset of $U$, in particular on $B(x_0,s)$. Hence
\begin{align*}
\Delta(\eta u)=-\eta f+2\sum_{k=1}^n(\partial_{x_k}\eta)(\partial_{x_k}u)+u\Delta\eta.
\end{align*}
Taking $L^p(B(x_0,s))$ norms, using $0\le \eta\le 1$, and using the derivative bounds for $\eta$, we obtain
\begin{align*}
\|D^2u\|_{L^p(B(x_0,r))} \le A\|f\|_{L^p(B(x_0,s))}+\frac{2AC_1}{s-r}\|\nabla u\|_{L^p(B(x_0,s))}+\frac{AC_2}{(s-r)^2}\|u\|_{L^p(B(x_0,s))}.
\end{align*}
This is the basic localized estimate. Its only defect is the first-order term $\|\nabla u\|_{L^p}$, which will be removed by interpolation.
[/guided]
[/step]
[step:Absorb the first-order term by a shrinking-radius iteration]
We use the following local first-derivative interpolation inequality for $W^{2,p}$ functions, stated here in the exact form needed. There is a constant $I=I(n,p)>0$ such that for every ball $B(x_0,s)\subset\mathbb{R}^n$, every $0<\varepsilon\le 1$, and every $v\in W^{2,p}(B(x_0,s))$,
\begin{align*}
\|\nabla v\|_{L^p(B(x_0,s))} \le \varepsilon s\|D^2v\|_{L^p(B(x_0,s))}+I\varepsilon^{-1}s^{-1}\|v\|_{L^p(B(x_0,s))}.
\end{align*}
This is the correctly scaled form of the unit-ball interpolation inequality. Indeed, for $V\in W^{2,p}(B(0,1))$ the unit-ball estimate is $\|\nabla V\|_{L^p(B(0,1))}\le \varepsilon\|D^2V\|_{L^p(B(0,1))}+I\varepsilon^{-1}\|V\|_{L^p(B(0,1))}$. Applying this to $V(y)=v(x_0+sy)$ and using the scaling of [Lebesgue measure](/page/Lebesgue%20Measure) gives the displayed estimate on $B(x_0,s)$.
Fix radii $0<\rho<R$ with $B(x_0,R)\subset U'$. Define $\Phi:(0,R)\to[0,\infty)$ by
\begin{align*}
\Phi(t)=\|D^2u\|_{L^p(B(x_0,t))}.
\end{align*}
Let $\rho\le r<s<R$. Applying the cutoff estimate and then the interpolation estimate to $v=u|_{B(x_0,s)}$ with
\begin{align*}
\varepsilon_{r,s}=\min\left\{1,\frac{s-r}{4AC_1s}\right\}
\end{align*}
gives
\begin{align*}
\Phi(r) \le \frac{1}{2}\Phi(s)+A\|f\|_{L^p(B(x_0,R))}+C_3(r,s)\|u\|_{L^p(B(x_0,R))},
\end{align*}
where the corrected scaled interpolation estimate gives
\begin{align*}
C_3(r,s)=\frac{2AC_1}{s-r}I\varepsilon_{r,s}^{-1}s^{-1}+\frac{AC_2}{(s-r)^2}.
\end{align*}
Since $\rho\le s\le R$ and $\varepsilon_{r,s}^{-1}\le 1+4AC_1s(s-r)^{-1}\le (1+4AC_1R)(s-r)^{-1}$ whenever $s-r\le 1$, while the case $s-r>1$ is bounded by the same expression after increasing the constant, every term in $C_3(r,s)$ is bounded by $K(s-r)^{-2}$ for a constant $K=K(n,p,\rho,R)>0$. Therefore
\begin{align*}
\Phi(r) \le \frac{1}{2}\Phi(s)+\frac{K}{(s-r)^2}\left(\|f\|_{L^p(B(x_0,R))}+\|u\|_{L^p(B(x_0,R))}\right)
\end{align*}
for all $\rho\le r<s<R$.
Before applying the iteration lemma, note that $\Phi$ is bounded on $[\rho,R)$. Indeed, for every $t<R$,
\begin{align*}
B(x_0,t)\subset B(x_0,R)\subset U',
\end{align*}
and $u\in W^{2,p}_{\mathrm{loc}}(U)$ with $U'\subset\subset U$ gives $D^2u\in L^p(U')$. Hence
\begin{align*}
\Phi(t)\le \|D^2u\|_{L^p(B(x_0,R))}<\infty
\end{align*}
uniformly for $t\in[\rho,R)$.
We now apply the shrinking-radius lemma, which is another standard iteration input: if a bounded nonnegative function $\Psi$ on $[\rho,R)$ satisfies
\begin{align*}
\Psi(r)\le \theta\Psi(s)+\frac{K_0}{(s-r)^q}M
\end{align*}
for all $\rho\le r<s<R$, with $0<\theta<1$, $K_0>0$, $q>0$, and $M\ge0$, then
\begin{align*}
\Psi(\rho)\le C_4K_0(R-\rho)^{-q}M,
\end{align*}
where $C_4=C_4(\theta,q)>0$. To prove the lemma, choose a geometric ratio $\lambda\in(0,1)$ sufficiently close to $1$ that $\theta\lambda^{-q}<1$, set $r_0=\rho$ and $r_j=\rho+(R-\rho)(1-\lambda^j)$ for $j\geq 1$, iterate the displayed inequality along the sequence $(r_j)$, and sum the resulting convergent geometric series; the residual term tends to zero because $\Psi$ is bounded.
Applying this lemma to $\Psi=\Phi$, $\theta=1/2$, $q=2$, and
\begin{align*}
M=\|f\|_{L^p(B(x_0,R))}+\|u\|_{L^p(B(x_0,R))}
\end{align*}
yields a constant $C_5=C_5(n,p,\rho,R)>0$, obtained from $K=K(n,p,\rho,R)$ and the shrinking-radius constant $C_4(1/2,2)$, such that
\begin{align*}
\|D^2u\|_{L^p(B(x_0,\rho))} \le C_5\left(\|f\|_{L^p(B(x_0,R))}+\|u\|_{L^p(B(x_0,R))}\right).
\end{align*}
[guided]
The obstruction in the cutoff estimate is the first-order term $\|\nabla u\|_{L^p}$. We remove it by trading a small multiple of $\|D^2u\|_{L^p}$ for a larger multiple of $\|u\|_{L^p}$. Precisely, for every ball $B(x_0,s)$, every $0<\varepsilon\le1$, and every $v\in W^{2,p}(B(x_0,s))$, the local interpolation estimate, stated here in the exact form used, gives
\begin{align*}
\|\nabla v\|_{L^p(B(x_0,s))} \le \varepsilon s\|D^2v\|_{L^p(B(x_0,s))}+I\varepsilon^{-1}s^{-1}\|v\|_{L^p(B(x_0,s))}.
\end{align*}
Why does the factor $s$ appear in front of the Hessian term? It is forced by scaling. If $V:B(0,1)\to\mathbb{R}$ is defined by $V(y)=v(x_0+sy)$, then $\|\nabla V\|_{L^p(B(0,1))}=s^{1-n/p}\|\nabla v\|_{L^p(B(x_0,s))}$, $\|D^2V\|_{L^p(B(0,1))}=s^{2-n/p}\|D^2v\|_{L^p(B(x_0,s))}$, and $\|V\|_{L^p(B(0,1))}=s^{-n/p}\|v\|_{L^p(B(x_0,s))}$. Applying the unit-ball interpolation inequality to $V$ and dividing by $s$ gives exactly the displayed estimate. We apply this to $v=u|_{B(x_0,s)}$, which is allowed because $u\in W^{2,p}_{\mathrm{loc}}(U)$ and $B(x_0,s)\subset U$.
Fix $0<\rho<R$ with $B(x_0,R)\subset U'$ and define
\begin{align*}
\Phi(t)=\|D^2u\|_{L^p(B(x_0,t))}, \qquad 0<t<R.
\end{align*}
For $\rho\le r<s<R$, choose
\begin{align*}
\varepsilon_{r,s}=\min\left\{1,\frac{s-r}{4AC_1s}\right\}.
\end{align*}
This choice is made so that the coefficient produced by the gradient term is at most $1/2$: multiplying the Hessian part of the interpolation estimate by $2AC_1(s-r)^{-1}$ gives $2AC_1\varepsilon_{r,s}s(s-r)^{-1}\le 1/2$. Substituting the interpolation estimate into the cutoff estimate gives
\begin{align*}
\Phi(r) \le \frac{1}{2}\Phi(s)+A\|f\|_{L^p(B(x_0,R))}+C_3(r,s)\|u\|_{L^p(B(x_0,R))},
\end{align*}
where the corrected scaled interpolation estimate gives
\begin{align*}
C_3(r,s)=\frac{2AC_1}{s-r}I\varepsilon_{r,s}^{-1}s^{-1}+\frac{AC_2}{(s-r)^2}.
\end{align*}
The constants are now explicit: the factor $(s-r)^{-1}$ comes from $|\nabla\eta|$, the factor $(s-r)^{-2}$ comes from $|\Delta\eta|$, and the term involving $\varepsilon_{r,s}^{-1}s^{-1}$ comes from the correctly scaled interpolation inequality. Since $\rho\le s\le R$, the lower bound on $s$ controls the factor $s^{-1}$ and the upper bound on $s$ controls the factor appearing in $\varepsilon_{r,s}^{-1}$. Thus all these contributions are bounded by
\begin{align*}
\frac{K}{(s-r)^2}
\end{align*}
for a constant $K=K(n,p,\rho,R)>0$. Therefore
\begin{align*}
\Phi(r) \le \frac{1}{2}\Phi(s)+\frac{K}{(s-r)^2}\left(\|f\|_{L^p(B(x_0,R))}+\|u\|_{L^p(B(x_0,R))}\right).
\end{align*}
The iteration lemma also requires boundedness of the function being iterated at all radii in $[\rho,R)$. This holds here because $B(x_0,t)\subset B(x_0,R)\subset U'$ for every $t<R$, while $u\in W^{2,p}_{\mathrm{loc}}(U)$ and $U'\subset\subset U$. Therefore $D^2u\in L^p(U')$ and
\begin{align*}
\Phi(t)\le \|D^2u\|_{L^p(B(x_0,R))}<\infty
\end{align*}
for every $t\in[\rho,R)$.
The final point is important: a finite iteration would leave an uncontrolled term involving $\Phi$ at the last radius. Instead we use the shrinking-radius lemma. The lemma says that if a bounded nonnegative function $\Psi$ on $[\rho,R)$ satisfies
\begin{align*}
\Psi(r)\le \theta\Psi(s)+\frac{K_0}{(s-r)^q}M
\end{align*}
for every $\rho\le r<s<R$, with $0<\theta<1$, then
\begin{align*}
\Psi(\rho)\le C_4K_0(R-\rho)^{-q}M.
\end{align*}
To prove it, choose a geometric ratio $\lambda\in(0,1)$ sufficiently close to $1$ so that $\theta\lambda^{-q}<1$, use the radii $r_0=\rho$ and $r_j=\rho+(R-\rho)(1-\lambda^j)$ for $j\geq 1$, and sum the resulting geometric series. This starts the iteration at the endpoint where the desired estimate is needed, rather than only at radii larger than $\rho$. The careful choice of $\lambda$ is what makes the residual term vanish and avoids any restriction such as $\theta 2^q<1$. Applying this lemma to $\Psi=\Phi$ with $q=2$ gives
\begin{align*}
\|D^2u\|_{L^p(B(x_0,\rho))} \le C_5\left(\|f\|_{L^p(B(x_0,R))}+\|u\|_{L^p(B(x_0,R))}\right),
\end{align*}
with $C_5=C_5(n,p,\rho,R)>0$.
[/guided]
[/step]
[step:Cover $\overline V$ by balls compactly contained in $U'$]
Because $V\subset\subset U'$, the compact set $\overline V$ has positive distance from $\mathbb{R}^n\setminus U'$. Choose finitely many points $x_1,\dots,x_N \in V$ and radii $0<r_\ell<R_\ell$ for $1\le \ell\le N$ such that
\begin{align*}
\overline V \subset \bigcup_{\ell=1}^N B(x_\ell,r_\ell)
\end{align*}
and
\begin{align*}
B(x_\ell,R_\ell)\subset U'
\end{align*}
for every $\ell$. Applying the local ball estimate to each pair $B(x_\ell,r_\ell)\subset B(x_\ell,R_\ell)$, with $\rho=r_\ell$ and $R=R_\ell$, gives constants $C_\ell=C_\ell(n,p,r_\ell,R_\ell)>0$ satisfying
\begin{align*}
\|D^2u\|_{L^p(B(x_\ell,r_\ell))} \le C_\ell\left(\|f\|_{L^p(B(x_\ell,R_\ell))}+\|u\|_{L^p(B(x_\ell,R_\ell))}\right).
\end{align*}
Since each $B(x_\ell,R_\ell)$ is contained in $U'$, we have
\begin{align*}
\|D^2u\|_{L^p(B(x_\ell,r_\ell))} \le C_\ell\left(\|f\|_{L^p(U')}+\|u\|_{L^p(U')}\right).
\end{align*}
Using the finite covering of $\overline V$ and summing the finitely many estimates gives
\begin{align*}
\|D^2u\|_{L^p(V)} \le C\left(\|f\|_{L^p(U')}+\|u\|_{L^p(U')}\right),
\end{align*}
where $C=C(n,p,V,U')>0$ is obtained from the finite family of balls and the constants $C_\ell$. This is exactly the asserted interior Calderón-Zygmund estimate.
[guided]
The local estimate controls $D^2u$ on one smaller ball by the data on a larger ball. To pass from one ball to the arbitrary compactly contained set $V$, we use only compactness. Since $V\subset\subset U'$, the compact set $\overline V$ has positive distance from $\mathbb{R}^n\setminus U'$, so every point of $\overline V$ lies in a small ball whose slightly larger concentric ball still lies inside $U'$.
Choose finitely many points $x_1,\dots,x_N\in U'$ and radii $0<r_\ell<R_\ell$ such that
\begin{align*}
\overline V \subset \bigcup_{\ell=1}^N B(x_\ell,r_\ell)
\end{align*}
and
\begin{align*}
B(x_\ell,R_\ell)\subset U'
\end{align*}
for every $1\le \ell\le N$. Applying the local ball estimate to each pair $B(x_\ell,r_\ell)\subset B(x_\ell,R_\ell)$ gives constants $C_\ell=C_\ell(n,p,r_\ell,R_\ell)>0$ such that
\begin{align*}
\|D^2u\|_{L^p(B(x_\ell,r_\ell))} \le C_\ell\left(\|f\|_{L^p(B(x_\ell,R_\ell))}+\|u\|_{L^p(B(x_\ell,R_\ell))}\right).
\end{align*}
Because each larger ball is contained in $U'$, enlarging the domain in each $L^p$ norm gives
\begin{align*}
\|D^2u\|_{L^p(B(x_\ell,r_\ell))} \le C_\ell\left(\|f\|_{L^p(U')}+\|u\|_{L^p(U')}\right).
\end{align*}
Finally, the finite covering of $\overline V$ implies that the $L^p(V)$ norm is bounded by the sum of the finitely many ball norms. Absorbing the finite number of constants into one constant $C=C(n,p,V,U')>0$ gives
\begin{align*}
\|D^2u\|_{L^p(V)} \le C\left(\|f\|_{L^p(U')}+\|u\|_{L^p(U')}\right),
\end{align*}
which is the desired interior estimate.
[/guided]
[/step]
