Idea of the proof: We prove the estimate by induction on $k$. The case $k=0$ is the mean-value formula. For $k=1$ we use that $u_{x_i}$ is harmonic, average it on the ball of radius $r/2$, convert the average to a boundary integral, and control the boundary values of $u$ by the $L^1$-average over $B(x_0,r)$. For the inductive step we repeat the same argument with $v := \partial^\beta u$ ($|\beta|=k-1$), but we average on the ball of radius $r/k$ so that $B\!\left(x,\tfrac{k-1}{k}r\right) \subset B(x_0,r)$ for $x \in \partial B(x_0,r/k)$.

**Step 0 (base, $k=0$).** \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) && \text{(mean-value for $|u|$)} \\ &= \frac{1}{\frac{\pi^{n/2}}{\Gamma(1+\frac{n}{2})} r^n} \, \|u\|_{L^1(B(x_0,r))} && \text{(volume of $B(x_0,r)$)} \\ &= \frac{C_0}{r^n} \, \|u\|_{L^1(B(x_0,r))}. \end{align*}

**Step 1 (case $k=1$).** Fix $i \in \{1,\dots,n\}$ and set $a := r/2$. Since $u_{x_i}$ is harmonic, \begin{align*} |u_{x_i}(x_0)| &= \left| \frac{1}{\mathcal{L}^n(B(x_0,a))} \int_{B(x_0,a)} u_{x_i}(y) \, d\mathcal{L}^n(y) \right| && \text{(mean-value for $u_{x_i}$)} \\ &= \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| && \text{(divergence theorem)} \\ &\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))} && \text{($|\nu_i|\le 1$)} \\ &= \frac{n}{a} \, \|u\|_{L^\infty(\partial B(x_0,a))} && \text{($\frac{\mathcal{H}^{n-1}(\partial B(0,a))}{\mathcal{L}^n(B(0,a))}=\frac{n}{a}$)} \\ &= \frac{2n}{r} \, \|u\|_{L^\infty(\partial B(x_0,r/2))}. \end{align*} For any $x \in \partial B(x_0,r/2)$ we have $B\!\left(x,\tfrac{r}{2}\right) \subset B(x_0,r)$, hence by the $k=0$ case, \begin{align*} |u(x)| &\le \frac{1}{\mathcal{L}^n(B(x,\tfrac{r}{2}))} \int_{B(x,\tfrac{r}{2})} |u(y)| \, d\mathcal{L}^n(y) \\ &\le \frac{1}{\frac{\pi^{n/2}}{\Gamma(1+\frac{n}{2})} (\tfrac{r}{2})^n} \, \|u\|_{L^1(B(x_0,r))} = \frac{2^n C_0}{r^n} \, \|u\|_{L^1(B(x_0,r))}. \end{align*} Combining, \begin{align*} |\nabla u(x_0)| \le \frac{2n}{r} \cdot \frac{2^n C_0}{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\, C_0$.

**Step 2 (induction step).** Assume the estimate holds for orders $\le k-1$ with the constants stated. Fix a multiindex $\alpha$ with $|\alpha|=k$ and write $\alpha = \beta + e_i$, $|\beta|=k-1$. Put $a := r/k$ and $v := \partial^\beta u$ (harmonic in $U$). Then \begin{align*} |\partial^\alpha u(x_0)| &= |(v)_{x_i}(x_0)| = \left| \frac{1}{\mathcal{L}^n(B(x_0,a))} \int_{\partial B(x_0,a)} v(y)\, \nu_i(y) \, d\mathcal{H}^{n-1}(y) \right| && \text{(mean-value for $v_{x_i}$ + divergence)} \\ &\le \frac{\mathcal{H}^{n-1}(\partial B(x_0,a))}{\mathcal{L}^n(B(x_0,a))} \, \|v\|_{L^\infty(\partial B(x_0,a))} = \frac{n}{a} \, \| \partial^\beta u \|_{L^\infty(\partial B(x_0,r/k))} \\ &= \frac{n k}{r} \, \| \partial^\beta u \|_{L^\infty(\partial B(x_0,r/k))}. \end{align*} For every $x \in \partial B(x_0,r/k)$ we have the inclusion \begin{align*} B\!\left(x, \frac{k-1}{k} r \right) \subset B(x_0,r), \end{align*} so by the induction hypothesis (order $k-1$) applied on $B\!\left(x, \tfrac{k-1}{k} r\right)$, \begin{align*} |\partial^\beta u(x)| &\le \frac{C_{k-1}}{\left( \frac{k-1}{k} r \right)^{n+k-1}} \, \|u\|_{L^1(B(x_0,r))} = 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*} Taking the supremum over $x \in \partial B(x_0,r/k)$ and inserting into (*) gives \begin{align*} \max_{|\alpha|=k} |\partial^\alpha u(x_0)| &\le \frac{n k}{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 \begin{align*} C_k := n k \left( \frac{k}{k-1} \right)^{n+k-1} C_{k-1}. \end{align*} This completes the induction and the proof.