[proofplan]
We prove the estimate by induction on the derivative order $k = |\alpha|$. The base case $k = 0$ is the mean-value inequality. For $k = 1$, we use that each partial derivative $\partial_{x_i} u$ is harmonic, apply the mean-value property followed by the divergence theorem to convert the volume average of $\partial_{x_i} u$ to a boundary integral of $u$, then control the boundary values using the $k = 0$ estimate on a concentric ball. The inductive step repeats this argument: we average a $(k-1)$th derivative on a ball of radius $r/k$, convert to a boundary integral, and apply the induction hypothesis on balls of radius $(k-1)r/k$ centred at boundary points.
[/proofplan]
[step:Establish the base case $k = 0$ via the mean-value inequality]
By the solid mean-value property for harmonic functions and the triangle inequality,
\begin{align*}
|u(x_0)| &\le \frac{1}{\mathcal{L}^n(B(x_0,r))} \int_{B(x_0,r)} |u(y)|\, d\mathcal{L}^n(y) = \frac{1}{\omega_n\, r^n}\, \|u\|_{L^1(B(x_0,r))},
\end{align*}
where $\omega_n = \mathcal{L}^n(B(0,1)) = \pi^{n/2}/\Gamma(1 + n/2)$. This gives the estimate with $C_0 = 1/\omega_n$.
[/step]
[step:Prove the first-order estimate by converting the volume average to a boundary integral]
Fix $i \in \{1, \dots, n\}$. Since $u$ is harmonic, $\partial_{x_i} u$ is also harmonic in $U$. Set $a := r/2$; then $B(x_0, a) \subset B(x_0, r) \subset U$. Applying the mean-value property to $\partial_{x_i} u$ on $B(x_0, a)$ and then the divergence theorem:
\begin{align*}
|\partial_{x_i} u(x_0)| &= \left|\frac{1}{\mathcal{L}^n(B(x_0,a))} \int_{B(x_0,a)} \partial_{x_i} u(y)\, d\mathcal{L}^n(y)\right| \\
&= \left|\frac{1}{\mathcal{L}^n(B(x_0,a))} \int_{\partial B(x_0,a)} u(y)\, \nu_i(y)\, d\mathcal{H}^{n-1}(y)\right|,
\end{align*}
where $\nu_i$ is the $i$th component of the outward unit normal on $\partial B(x_0,a)$. Since $|\nu_i| \le 1$,
\begin{align*}
|\partial_{x_i} u(x_0)| &\le \frac{\mathcal{H}^{n-1}(\partial B(x_0,a))}{\mathcal{L}^n(B(x_0,a))}\, \|u\|_{L^\infty(\partial B(x_0,a))} = \frac{n}{a}\, \|u\|_{L^\infty(\partial B(x_0,r/2))},
\end{align*}
using $\mathcal{H}^{n-1}(\partial B(0,a))/\mathcal{L}^n(B(0,a)) = n/a$. Substituting $a = r/2$ gives $|\partial_{x_i} u(x_0)| \le \frac{2n}{r}\, \|u\|_{L^\infty(\partial B(x_0, r/2))}$.
[guided]
Why do we work on the ball of radius $r/2$ rather than $r$? Because the boundary values of $u$ on $\partial B(x_0, r/2)$ can be controlled by the $L^1$ norm on $B(x_0, r)$: for any $x \in \partial B(x_0, r/2)$, the ball $B(x, r/2)$ is contained in $B(x_0, r)$ (by the triangle inequality), so the $k = 0$ estimate applies on $B(x, r/2)$.
Fix $i \in \{1, \dots, n\}$. Since $u$ is harmonic in $U$, its partial [derivative](/page/Derivative) $\partial_{x_i} u$ is also harmonic in $U$ (the Laplacian commutes with partial derivatives for smooth functions). Set $a := r/2$; then $B(x_0, a) \subset B(x_0, r) \subset U$.
The mean-value property for the harmonic function $\partial_{x_i} u$ on $B(x_0, a)$ gives
\begin{align*}
\partial_{x_i} u(x_0) = \frac{1}{\mathcal{L}^n(B(x_0,a))} \int_{B(x_0,a)} \partial_{x_i} u(y)\, d\mathcal{L}^n(y).
\end{align*}
We now convert the volume integral of a derivative into a boundary integral using the divergence theorem. Applying it to the vector field $y \mapsto u(y)\, e_i$ on $B(x_0, a)$:
\begin{align*}
\int_{B(x_0,a)} \partial_{x_i} u(y)\, d\mathcal{L}^n(y) = \int_{\partial B(x_0,a)} u(y)\, \nu_i(y)\, d\mathcal{H}^{n-1}(y),
\end{align*}
where $\nu_i(y)$ is the $i$th component of the outward unit normal $\nu(y) = (y - x_0)/|y - x_0|$ on $\partial B(x_0,a)$. Taking absolute values and using $|\nu_i| \le 1$:
\begin{align*}
|\partial_{x_i} u(x_0)| &\le \frac{1}{\mathcal{L}^n(B(x_0,a))} \int_{\partial B(x_0,a)} |u(y)|\, d\mathcal{H}^{n-1}(y) \\
&\le \frac{\mathcal{H}^{n-1}(\partial B(x_0,a))}{\mathcal{L}^n(B(x_0,a))}\, \|u\|_{L^\infty(\partial B(x_0,a))}.
\end{align*}
The ratio $\mathcal{H}^{n-1}(\partial B(0,a))/\mathcal{L}^n(B(0,a)) = n\omega_n a^{n-1}/(\omega_n a^n) = n/a$. With $a = r/2$:
\begin{align*}
|\partial_{x_i} u(x_0)| \le \frac{2n}{r}\, \|u\|_{L^\infty(\partial B(x_0, r/2))}.
\end{align*}
[/guided]
[/step]
[step:Control the boundary values using the base case to complete the $k = 1$ estimate]
For any $x \in \partial B(x_0, r/2)$, the triangle inequality gives $B(x, r/2) \subset B(x_0, r)$. Applying the base case on $B(x, r/2)$:
\begin{align*}
|u(x)| \le \frac{1}{\omega_n (r/2)^n}\, \|u\|_{L^1(B(x, r/2))} \le \frac{2^n}{\omega_n\, r^n}\, \|u\|_{L^1(B(x_0, r))},
\end{align*}
where the second inequality uses $B(x, r/2) \subset B(x_0, r)$. Taking the supremum over $x \in \partial B(x_0, r/2)$ and substituting into the previous estimate:
\begin{align*}
|\nabla u(x_0)| \le \frac{2n}{r} \cdot \frac{2^n}{\omega_n\, r^n}\, \|u\|_{L^1(B(x_0, r))} = \frac{C_1}{r^{n+1}}\, \|u\|_{L^1(B(x_0, r))},
\end{align*}
with $C_1 = 2^{n+1} n / \omega_n$.
[/step]
[step:Execute the inductive step from order $k - 1$ to order $k$]
Assume the estimate holds for all multi-indices of order $\le k - 1$. Fix $\alpha$ with $|\alpha| = k$ and write $\alpha = \beta + e_i$ where $|\beta| = k - 1$ and $e_i$ is the $i$th standard basis multi-index. Define $v := D^\beta u$, which is harmonic in $U$ since the Laplacian commutes with partial derivatives for smooth functions. Set $a := r/k$.
Repeating the argument from the $k = 1$ step with $v$ in place of $u$ and averaging radius $a = r/k$:
\begin{align*}
|D^\alpha u(x_0)| = |\partial_{x_i} v(x_0)| \le \frac{n}{a}\, \|v\|_{L^\infty(\partial B(x_0, r/k))} = \frac{nk}{r}\, \|D^\beta u\|_{L^\infty(\partial B(x_0, r/k))}.
\end{align*}
For any $x \in \partial B(x_0, r/k)$, the triangle inequality gives $B(x, (k-1)r/k) \subset B(x_0, r)$. Applying the induction hypothesis at order $k - 1$ on this ball:
\begin{align*}
|D^\beta u(x)| \le \frac{C_{k-1}}{\left(\frac{(k-1)r}{k}\right)^{n+k-1}}\, \|u\|_{L^1(B(x, (k-1)r/k))} \le C_{k-1} \left(\frac{k}{k-1}\right)^{n+k-1} \frac{1}{r^{n+k-1}}\, \|u\|_{L^1(B(x_0, r))}.
\end{align*}
Substituting:
\begin{align*}
|D^\alpha u(x_0)| \le \frac{nk}{r} \cdot C_{k-1} \left(\frac{k}{k-1}\right)^{n+k-1} \frac{1}{r^{n+k-1}}\, \|u\|_{L^1(B(x_0, r))} = \frac{C_k}{r^{n+k}}\, \|u\|_{L^1(B(x_0, r))},
\end{align*}
with the recursive constant
\begin{align*}
C_k := nk \left(\frac{k}{k-1}\right)^{n+k-1} C_{k-1}.
\end{align*}
[guided]
The inductive step follows the same two-part structure as the $k = 1$ case: (i) convert a volume average of a derivative into a boundary integral, then (ii) control the boundary values using the estimate at the previous order.
Assume the estimate holds for all multi-indices of order $\le k - 1$ with the stated constants. Fix $\alpha$ with $|\alpha| = k$. Write $\alpha = \beta + e_i$ where $|\beta| = k - 1$. Define $v := D^\beta u$. Since $u$ is harmonic and smooth, $v$ is also harmonic in $U$: $\Delta v = \Delta D^\beta u = D^\beta \Delta u = 0$.
Why choose averaging radius $a = r/k$? We need the boundary points of $B(x_0, r/k)$ to admit balls of radius $(k-1)r/k$ that fit inside $B(x_0, r)$. For $x \in \partial B(x_0, r/k)$, we have $|x - x_0| = r/k$, so $B(x, (k-1)r/k) \subset B(x_0, r/k + (k-1)r/k) = B(x_0, r)$. This is the geometric constraint that determines the radius.
Applying the mean-value property to the harmonic function $\partial_{x_i} v$ on $B(x_0, a)$ and then the divergence theorem (as in the $k = 1$ step):
\begin{align*}
|D^\alpha u(x_0)| = |\partial_{x_i} v(x_0)| \le \frac{n}{a}\, \|v\|_{L^\infty(\partial B(x_0, a))} = \frac{nk}{r}\, \|D^\beta u\|_{L^\infty(\partial B(x_0, r/k))}.
\end{align*}
Now we control $|D^\beta u(x)|$ for $x \in \partial B(x_0, r/k)$. Since $B(x, (k-1)r/k) \subset B(x_0, r)$, the induction hypothesis at order $|\beta| = k - 1$ on the ball $B(x, (k-1)r/k)$ gives
\begin{align*}
|D^\beta u(x)| &\le \frac{C_{k-1}}{\left(\frac{(k-1)r}{k}\right)^{n+k-1}}\, \|u\|_{L^1(B(x, (k-1)r/k))} \\
&\le C_{k-1} \left(\frac{k}{k-1}\right)^{n+k-1} \frac{1}{r^{n+k-1}}\, \|u\|_{L^1(B(x_0, r))},
\end{align*}
where the second inequality uses $B(x, (k-1)r/k) \subset B(x_0, r)$ to enlarge the domain of integration. Taking the supremum over $x \in \partial B(x_0, r/k)$ and combining:
\begin{align*}
|D^\alpha u(x_0)| \le \frac{nk}{r} \cdot C_{k-1} \left(\frac{k}{k-1}\right)^{n+k-1} \frac{1}{r^{n+k-1}}\, \|u\|_{L^1(B(x_0, r))} = \frac{C_k}{r^{n+k}}\, \|u\|_{L^1(B(x_0, r))},
\end{align*}
with $C_k := nk \left(\frac{k}{k-1}\right)^{n+k-1} C_{k-1}$. This constant depends only on $n$ and $k$, completing the induction.
[/guided]
[/step]
