[proofplan]
The proof is a direct application of differentiation under the integral sign, justified by the [Dominated Convergence Theorem](/theorems/511). We form the difference quotient of the [convolution](/page/Convolution) $f_\varepsilon$, pass the difference quotient to the smooth kernel $\eta_\varepsilon$ using the [Mean Value Theorem](/theorems/186), and obtain an integrable dominator from the compact support and $L^\infty$ bound of the kernel derivatives. Iterating over all multi-indices yields $f_\varepsilon \in C^\infty(U_\varepsilon)$ with the stated derivative formula.
[/proofplan]
[step:Write the difference quotient and pass it to the kernel]
Fix $x \in U_\varepsilon$ and a coordinate direction $e_j$ for $1 \le j \le n$. For $|h|$ small enough that $x + he_j \in U_\varepsilon$, the difference quotient of $f_\varepsilon$ is:
\begin{align*}
\frac{f_\varepsilon(x + he_j) - f_\varepsilon(x)}{h} = \int_{B(0,\varepsilon)} \frac{\eta_\varepsilon(y + he_j) - \eta_\varepsilon(y)}{h} \, f(x - y) \, d\mathcal{L}^n(y),
\end{align*}
where we used the definition $f_\varepsilon(x) = \int_{B(0,\varepsilon)} \eta_\varepsilon(y) f(x - y) \, d\mathcal{L}^n(y)$ and the fact that for $y \in B(0, \varepsilon)$ and $x \in U_\varepsilon$, we have $x - y \in U$.
[guided]
We want to show that $f_\varepsilon$ is differentiable and that the derivative passes to the kernel. Starting from the definition of $f_\varepsilon$ as a [convolution](/page/Convolution):
\begin{align*}
f_\varepsilon(x) = \int_{B(0,\varepsilon)} \eta_\varepsilon(y) \, f(x - y) \, d\mathcal{L}^n(y),
\end{align*}
the difference quotient becomes:
\begin{align*}
\frac{f_\varepsilon(x + he_j) - f_\varepsilon(x)}{h} = \int_{B(0,\varepsilon)} \eta_\varepsilon(y) \, \frac{f(x + he_j - y) - f(x - y)}{h} \, d\mathcal{L}^n(y).
\end{align*}
We cannot pass the limit inside directly because $f \in L^1_{\mathrm{loc}}(U)$ may not be differentiable. Instead, we use a different form of the convolution. Writing $f_\varepsilon(x) = \int_{B(0,\varepsilon)} \eta_\varepsilon(x - z) \, f(z) \, d\mathcal{L}^n(z)$ and substituting $y = x - z$, we can equivalently pass the difference quotient to the kernel:
\begin{align*}
\frac{f_\varepsilon(x + he_j) - f_\varepsilon(x)}{h} = \int_{B(0,\varepsilon)} \frac{\eta_\varepsilon(y + he_j) - \eta_\varepsilon(y)}{h} \, f(x - y) \, d\mathcal{L}^n(y).
\end{align*}
The strategy is to pass the difference quotient to $\eta_\varepsilon$ rather than $f$, because $\eta_\varepsilon \in C_c^\infty(\mathbb{R}^n)$ is smooth and we can control its derivatives.
[/guided]
[/step]
[step:Apply the Mean Value Theorem and dominated convergence to pass to the limit]
By the [Mean Value Theorem](/theorems/186) applied to the $C^\infty$ function $\eta_\varepsilon$, for each $y \in B(0,\varepsilon)$ there exists $\xi_h$ between $y$ and $y + he_j$ such that:
\begin{align*}
\frac{\eta_\varepsilon(y + he_j) - \eta_\varepsilon(y)}{h} = \partial_{y_j}\eta_\varepsilon(\xi_h).
\end{align*}
In particular, $\left|\frac{\eta_\varepsilon(y + he_j) - \eta_\varepsilon(y)}{h}\right| \le \|\partial_j \eta_\varepsilon\|_{L^\infty(\mathbb{R}^n)}$ uniformly in $h$ and $y$. The pointwise limit as $h \to 0$ is $\partial_{y_j}\eta_\varepsilon(y)$.
The dominating function is:
\begin{align*}
g: B(0,\varepsilon) &\to [0, \infty) \\
y &\mapsto \|\partial_j \eta_\varepsilon\|_{L^\infty(\mathbb{R}^n)} \cdot |f(x - y)|.
\end{align*}
This is integrable over $B(0,\varepsilon)$: the constant $\|\partial_j \eta_\varepsilon\|_{L^\infty(\mathbb{R}^n)} < \infty$ because $\eta_\varepsilon \in C_c^\infty(\mathbb{R}^n)$, and $\int_{B(0,\varepsilon)} |f(x-y)| \, d\mathcal{L}^n(y) < \infty$ because $f \in L^1_{\mathrm{loc}}(U)$ and the set $\{x - y : y \in B(0,\varepsilon)\} = B(x,\varepsilon)$ is a compact subset of $U$ (since $x \in U_\varepsilon$).
By the [Dominated Convergence Theorem](/theorems/511):
\begin{align*}
\partial_{x_j} f_\varepsilon(x) = \int_{B(0,\varepsilon)} \partial_{y_j} \eta_\varepsilon(y) \, f(x - y) \, d\mathcal{L}^n(y).
\end{align*}
[guided]
We need to verify the hypotheses of the [Dominated Convergence Theorem](/theorems/511) for the sequence of integrands $h \mapsto \frac{\eta_\varepsilon(y + he_j) - \eta_\varepsilon(y)}{h} \cdot f(x-y)$.
**Pointwise convergence:** For each fixed $y$, since $\eta_\varepsilon \in C^\infty$:
\begin{align*}
\lim_{h \to 0} \frac{\eta_\varepsilon(y + he_j) - \eta_\varepsilon(y)}{h} = \partial_{y_j}\eta_\varepsilon(y).
\end{align*}
**Integrable dominator:** By the [Mean Value Theorem](/theorems/186), for each $h \neq 0$ and each $y$, there exists $\xi_h$ on the segment between $y$ and $y + he_j$ such that:
\begin{align*}
\left|\frac{\eta_\varepsilon(y + he_j) - \eta_\varepsilon(y)}{h}\right| = |\partial_{y_j}\eta_\varepsilon(\xi_h)| \le \|\partial_j \eta_\varepsilon\|_{L^\infty(\mathbb{R}^n)}.
\end{align*}
Therefore the full integrand is bounded by $\|\partial_j \eta_\varepsilon\|_{L^\infty} \cdot |f(x-y)|$. We verify integrability: $\|\partial_j \eta_\varepsilon\|_{L^\infty} < \infty$ because $\eta_\varepsilon \in C_c^\infty(\mathbb{R}^n)$, and the integral $\int_{B(0,\varepsilon)} |f(x-y)| \, d\mathcal{L}^n(y) < \infty$ because $x \in U_\varepsilon$ implies $B(x,\varepsilon) \subseteq U$, so $f$ is integrable over this compact subset by the hypothesis $f \in L^1_{\mathrm{loc}}(U)$.
Both hypotheses are satisfied, so the [Dominated Convergence Theorem](/theorems/511) gives:
\begin{align*}
\partial_{x_j} f_\varepsilon(x) = \lim_{h \to 0} \int_{B(0,\varepsilon)} \frac{\eta_\varepsilon(y + he_j) - \eta_\varepsilon(y)}{h} \, f(x-y) \, d\mathcal{L}^n(y) = \int_{B(0,\varepsilon)} \partial_{y_j}\eta_\varepsilon(y) \, f(x-y) \, d\mathcal{L}^n(y).
\end{align*}
Note: the sign is consistent because in the convolution form $f_\varepsilon(x) = \int \eta_\varepsilon(x-y) f(y) \, d\mathcal{L}^n(y)$, we have $\partial_{x_j}[\eta_\varepsilon(x-y)] = (\partial_j \eta_\varepsilon)(x-y)$. Under the substitution $y \mapsto x - y$, this becomes $\partial_{y_j}\eta_\varepsilon(y)$ with no sign change.
[/guided]
[/step]
[step:Iterate to all multi-indices and conclude $f_\varepsilon \in C^\infty(U_\varepsilon)$]
The right-hand side $x \mapsto \int_{B(0,\varepsilon)} \partial_{y_j}\eta_\varepsilon(y) \, f(x-y) \, d\mathcal{L}^n(y)$ is itself a [convolution](/page/Convolution) of $\partial_j \eta_\varepsilon \in C_c^\infty(\mathbb{R}^n)$ with $f \in L^1_{\mathrm{loc}}(U)$, so the same argument applies to this function. Iterating over all multi-indices $\alpha$, with the dominating function $\|D^\alpha \eta_\varepsilon\|_{L^\infty} \cdot |f(x-y)| \cdot \mathbb{1}_{B(0,\varepsilon)}(y)$ at each stage:
\begin{align*}
D^\alpha f_\varepsilon(x) = \int_{B(0,\varepsilon)} D^\alpha \eta_\varepsilon(y) \, f(x - y) \, d\mathcal{L}^n(y) = (D^\alpha \eta_\varepsilon * f)(x)
\end{align*}
for every multi-index $\alpha$ and every $x \in U_\varepsilon$. Since each $D^\alpha f_\varepsilon$ is continuous (being a convolution of a $C_c^\infty$ function with a locally integrable function, with the same dominated convergence argument verifying continuity), we conclude $f_\varepsilon \in C^\infty(U_\varepsilon)$.
[guided]
Why can we iterate indefinitely? At each stage, we have established that $D^\alpha f_\varepsilon(x) = (D^\alpha \eta_\varepsilon * f)(x)$. This function has the same structure as the original $f_\varepsilon$: it is a [convolution](/page/Convolution) of a $C_c^\infty$ kernel (namely $D^\alpha \eta_\varepsilon$) with the same $f \in L^1_{\mathrm{loc}}(U)$. Therefore the identical dominated convergence argument applies to differentiate once more, replacing $\eta_\varepsilon$ with $D^\alpha \eta_\varepsilon$ and using $\|D^{\alpha + e_j}\eta_\varepsilon\|_{L^\infty} \cdot |f(x-y)|$ as the new dominator.
Since $\eta_\varepsilon \in C_c^\infty(\mathbb{R}^n)$, the bound $\|D^\alpha \eta_\varepsilon\|_{L^\infty} < \infty$ holds for every multi-index $\alpha$, so the induction never breaks down. The result is that $D^\alpha f_\varepsilon$ exists and is continuous on $U_\varepsilon$ for every $\alpha$, giving $f_\varepsilon \in C^\infty(U_\varepsilon)$ with the formula:
\begin{align*}
D^\alpha f_\varepsilon(x) = \int_{B(0,\varepsilon)} D^\alpha \eta_\varepsilon(y) \, f(x - y) \, d\mathcal{L}^n(y) = (D^\alpha \eta_\varepsilon * f)(x).
\end{align*}
This completes the proof that mollification produces smooth functions with derivatives given by convolving against derivatives of the kernel.
[/guided]
[/step]