[proofplan]
We prove the estimate after scaling and translating to the unit ball. The graph of $u$ is a minimal hypersurface, and the vertical normal component $v^{-1}$ is a positive Jacobi function, where $v=\sqrt{1+|\nabla u|^2}$. Combining the stability inequality, the Michael-Simon Sobolev inequality on the graph, and Moser iteration gives a local $L^\infty$ bound for $v$ in terms of an intrinsic average of $v$. The remaining point is to bound that average by the height oscillation of the graph, which produces the exponential factor.
[/proofplan]
[step:Normalize the estimate to the unit ball]
Let $\mathcal{L}^n$ denote $n$-dimensional [Lebesgue measure](/page/Lebesgue%20Measure) on $\mathbb{R}^n$. Define the rescaled function $w: B(0,1) \to \mathbb{R}$ by
\begin{align*}
w(y):=R^{-1}\bigl(u(x_0+Ry)-\inf_{B(x_0,R)}u\bigr).
\end{align*}
Then $w\in C^2(B(0,1))$ and $w$ solves
\begin{align*}
\operatorname{div}\left(\frac{\nabla w}{\sqrt{1+|\nabla w|^2}}\right)=0
\end{align*}
on $B(0,1)$, because $\nabla w(y)=\nabla u(x_0+Ry)$ and the divergence transforms by the constant factor $R$. Also
\begin{align*}
\operatorname{osc}_{B(0,1)}w=\frac{\operatorname{osc}_{B(x_0,R)}u}{R}.
\end{align*}
Thus it is enough to prove
\begin{align*}
\sup_{B(0,1/2)}\sqrt{1+|\nabla w|^2}\le C(n)\exp\bigl(C(n)\operatorname{osc}_{B(0,1)}w\bigr).
\end{align*}
After proving this normalized estimate, substituting $w(y)=R^{-1}(u(x_0+Ry)-\inf_{B(x_0,R)}u)$ gives the stated bound.
[/step]
[step:Apply the Jacobi equation to the vertical normal component]
For the normalized function, define $v: B(0,1) \to [1,\infty)$ by
\begin{align*}
v(x):=\sqrt{1+|\nabla w(x)|^2}.
\end{align*}
Let $M\subset \mathbb{R}^{n+1}$ be the graph
\begin{align*}
M:=\{(x,w(x)):x\in B(0,1)\},
\end{align*}
and let $d\mu_M$ denote the induced $n$-dimensional Riemannian volume measure on $M$. Let $\nabla_M$ denote the intrinsic gradient on $M$ with respect to the induced Riemannian metric, and let $\Delta_M$ denote the corresponding Laplace-Beltrami operator. The upward unit normal is the map $\nu: M \to \mathbb{R}^{n+1}$ defined by
\begin{align*}
\nu(x,w(x)):=\frac{(-\nabla w(x),1)}{v(x)}.
\end{align*}
Its vertical component is $\nu_{n+1}=v^{-1}$. Since $w$ solves the minimal surface equation, $M$ is a minimal hypersurface. The [Jacobi equation for minimal hypersurfaces](/page/Jacobi%20Operator) applied to the normal component of the constant vertical ambient vector field gives
\begin{align*}
\Delta_M(v^{-1})+|A|^2v^{-1}=0,
\end{align*}
where $\Delta_M$ is the Laplace-Beltrami operator on $M$ and $A$ is the second fundamental form of $M$. Equivalently, using $v\ge 1$ and differentiating $v^{-1}$,
\begin{align*}
\Delta_M v\ge \frac{|\nabla_M v|^2}{v}.
\end{align*}
Thus $\log v$ is weakly subharmonic on $M$.
[guided]
We introduce the graph because the gradient quantity in the theorem is geometrically natural. Define $v: B(0,1) \to [1,\infty)$ by
\begin{align*}
v(x):=\sqrt{1+|\nabla w(x)|^2}.
\end{align*}
The graph
\begin{align*}
M:=\{(x,w(x)):x\in B(0,1)\}\subset\mathbb{R}^{n+1}
\end{align*}
is a minimal hypersurface precisely because $w$ solves the minimal surface equation. Its upward unit normal is the map $\nu: M \to \mathbb{R}^{n+1}$ defined by
\begin{align*}
\nu(x,w(x)):=\frac{(-\nabla w(x),1)}{v(x)}.
\end{align*}
The last coordinate of this normal is $\nu_{n+1}=v^{-1}$. Let $\nabla_M$ denote the intrinsic gradient on $M$ with respect to the induced Riemannian metric, and let $\Delta_M$ denote the corresponding Laplace-Beltrami operator. The [Jacobi equation for minimal hypersurfaces](/page/Jacobi%20Operator) says that the normal component of every constant ambient vector field satisfies
\begin{align*}
\Delta_M \nu_{n+1}+|A|^2\nu_{n+1}=0,
\end{align*}
where $A$ is the second fundamental form. Substituting $\nu_{n+1}=v^{-1}$ gives
\begin{align*}
\Delta_M(v^{-1})+|A|^2v^{-1}=0.
\end{align*}
Since $v>0$, differentiating $v^{-1}$ on $M$ gives
\begin{align*}
\Delta_M(v^{-1})=-v^{-2}\Delta_Mv+2v^{-3}|\nabla_Mv|^2.
\end{align*}
Multiplying the Jacobi equation by $v^2$ therefore yields
\begin{align*}
\Delta_Mv=2\frac{|\nabla_Mv|^2}{v}+|A|^2v\ge \frac{|\nabla_Mv|^2}{v}.
\end{align*}
This is the differential inequality needed for Moser iteration: it says that $\log v$ is weakly subharmonic, because
\begin{align*}
\Delta_M\log v=\frac{\Delta_Mv}{v}-\frac{|\nabla_Mv|^2}{v^2}\ge 0.
\end{align*}
[/guided]
[/step]
[step:Apply the logarithmic mean-value estimate to $\log v$]
Let $M_r:=M\cap (B(0,r)\times\mathbb{R})$ for $0<r<1$. Since $M$ is a smooth minimal graph, it is stable, so the [stability inequality for minimal hypersurfaces](/page/Stability%20Inequality) holds on every compactly supported [test function](/page/Test%20Function) $\phi\in C_c^1(M_{3/4})$. Since $M$ is minimal, the [Michael-Simon Sobolev inequality](/page/Michael-Simon%20Sobolev%20Inequality) applies with zero mean-curvature term: there is a constant $S(n)>0$ such that every $\phi\in C_c^1(M_{3/4})$ satisfies
\begin{align*}
\left(\int_{M_{3/4}} |\phi|^{n/(n-1)}\,d\mu_M\right)^{(n-1)/n}\le S(n)\int_{M_{3/4}} |\nabla_M\phi|\,d\mu_M.
\end{align*}
We use the following precise local estimate, the logarithmic mean-value inequality obtained from Moser iteration with the stability inequality and the Michael-Simon Sobolev inequality. If $N$ is a smooth stable minimal hypersurface in $\mathbb{R}^{n+1}$, if $N_\rho:=N\cap (B(0,\rho)\times\mathbb{R})$ is compactly contained in the coordinate patch under consideration, if $0<\rho_0<\rho<1$, and if $\psi:N_\rho\to[0,\infty)$ is weakly subharmonic with respect to $\Delta_N$, then
\begin{align*}
\sup_{N_{\rho_0}}\psi\le A(n,\rho_0,\rho)\left(1+\frac{1}{\mu_N(N_\rho)}\int_{N_\rho}\psi\,d\mu_N\right),
\end{align*}
where $A(n,\rho_0,ρ)>0$ depends only on $n$ and on the fixed cutoff geometry determined by $\rho-\rho_0$. This is the exact form of the [Moser iteration](/page/Moser%20Iteration) consequence used here; the volume factor is normalized into the average over $N_\rho$. Apply this statement with $N=M$, $\rho_0=1/2$, $\rho=3/4$, and with the map $\psi:M_{3/4}\to[0,\infty)$ defined by $\psi(x,w(x))=\log v(x)$. The hypotheses hold because $M$ is a smooth stable minimal hypersurface, $M_{3/4}$ lies a positive Euclidean distance from the lateral boundary of the graph over $B(0,1)$, the preceding step proves $\log v$ is weakly subharmonic, and $v\ge 1$ gives $\psi\ge 0$. Hence, for a constant $A_1(n):=A(n,1/2,3/4)>0$,
\begin{align*}
\sup_{M_{1/2}}\log v\le A_1(n)\left(1+\fint_{M_{3/4}}\log v\,d\mu_M\right),
\end{align*}
where
\begin{align*}
\fint_{M_{3/4}}\log v\,d\mu_M:=\frac{1}{\mu_M(M_{3/4})}\int_{M_{3/4}}\log v\,d\mu_M.
\end{align*}
Exponentiating and absorbing $\exp(A_1(n))$ into the multiplicative constant gives a dimensional constant $C_1(n)>0$ such that
\begin{align*}
\sup_{M_{1/2}} v\le C_1(n)\exp\left(C_1(n)\fint_{M_{3/4}}\log v\,d\mu_M\right).
\end{align*}
[guided]
The iteration is applied to $\log v$, not directly to $v$. This distinction matters: the previous step established weak subharmonicity of $\log v$, and an estimate for $\sup \log v$ becomes an exponential estimate for $\sup v$ after exponentiation.
We first record the analytic structure needed for the iteration. Because $M$ is a smooth minimal graph, it is stable; hence the [stability inequality for minimal hypersurfaces](/page/Stability%20Inequality) holds for every compactly supported $\phi\in C_c^1(M_{3/4})$. Because the mean curvature of $M$ is zero, the [Michael-Simon Sobolev inequality](/page/Michael-Simon%20Sobolev%20Inequality) has no mean-curvature error term, so there is a constant $S(n)>0$ such that
\begin{align*}
\left(\int_{M_{3/4}} |\phi|^{n/(n-1)}\,d\mu_M\right)^{(n-1)/n}\le S(n)\int_{M_{3/4}} |\nabla_M\phi|\,d\mu_M
\end{align*}
for every $\phi\in C_c^1(M_{3/4})$.
The precise local consequence of [Moser iteration](/page/Moser%20Iteration) used here is the following logarithmic mean-value estimate. If $N$ is a smooth stable minimal hypersurface in $\mathbb{R}^{n+1}$, if $N_\rho:=N\cap (B(0,\rho)\times\mathbb{R})$ is compactly contained in the coordinate patch under consideration, if $0<\rho_0<\rho<1$, and if $\psi:N_\rho\to[0,\infty)$ is weakly subharmonic with respect to $\Delta_N$, then
\begin{align*}
\sup_{N_{\rho_0}}\psi\le A(n,\rho_0,\rho)\left(1+\frac{1}{\mu_N(N_\rho)}\int_{N_\rho}\psi\,d\mu_N\right).
\end{align*}
Here the constant $A(n,\rho_0,\rho)>0$ is produced by the Michael-Simon Sobolev constant, the stability inequality, and the cutoff functions used between the two radii. For the fixed radii $\rho_0=1/2$ and $\rho=3/4$, this dependence is purely dimensional, and the normalization by $\mu_N(N_\rho)$ is exactly the average appearing in the conclusion.
We apply this result with $N=M$ and with the map $\psi:M_{3/4}\to[0,\infty)$ defined by $\psi(x,w(x))=\log v(x)$. The subharmonicity hypothesis holds by the Jacobi-equation computation in the previous step. The non-negativity hypothesis holds because $v\ge 1$, hence $\log v\ge 0$. The stability and Sobolev hypotheses hold because $M$ is a smooth minimal graph, hence stable, and the Michael-Simon Sobolev inequality was verified above on $M_{3/4}$. The locality hypothesis holds because $M_{3/4}$ is strictly inside the graph over $B(0,1)$, so the cutoff functions used in the iteration can be supported before reaching the lateral boundary. Therefore, with $A_1(n):=A(n,1/2,3/4)$,
\begin{align*}
\sup_{M_{1/2}}\log v\le A_1(n)\left(1+\fint_{M_{3/4}}\log v\,d\mu_M\right),
\end{align*}
where
\begin{align*}
\fint_{M_{3/4}}\log v\,d\mu_M:=\frac{1}{\mu_M(M_{3/4})}\int_{M_{3/4}}\log v\,d\mu_M.
\end{align*}
Exponentiating gives
\begin{align*}
\sup_{M_{1/2}}v\le \exp(A_1(n))\exp\left(A_1(n)\fint_{M_{3/4}}\log v\,d\mu_M\right).
\end{align*}
Absorbing the fixed factor $\exp(A_1(n))$ and the coefficient $A_1(n)$ into a single dimensional constant $C_1(n)>0$ gives
\begin{align*}
\sup_{M_{1/2}} v\le C_1(n)\exp\left(C_1(n)\fint_{M_{3/4}}\log v\,d\mu_M\right).
\end{align*}
[/guided]
[/step]
[step:Bound the logarithmic average by the height oscillation]
Let
\begin{align*}
H:=\operatorname{osc}_{B(0,1)}w.
\end{align*}
Since $0\le w\le H$ on $B(0,1)$ after the normalization, the graph $M_{3/4}$ lies in the cylinder $B(0,3/4)\times[0,H]$. The area element of the graph, measured with respect to the already defined Lebesgue measure $\mathcal{L}^n$ on the base, is
\begin{align*}
d\mu_M(x,w(x))=v(x)\,d\mathcal{L}^n(x).
\end{align*}
We use the following precise radius-dependent form of the [Bombieri-De Giorgi-Miranda logarithmic area estimate](/page/Bombieri-De%20Giorgi-Miranda%20Estimate). If $f\in C^2(B(0,1))$ solves the minimal surface equation, if $0<\rho<1$, if $\Gamma_f:=\{(x,f(x)):x\in B(0,1)\}$, if $\Gamma_{f,\rho}:=\Gamma_f\cap(B(0,\rho)\times\mathbb{R})$, and if $H_0\ge \sup_{B(0,1)}f-\inf_{B(0,1)}f$, then, with the map $v_f:B(0,1)\to[1,\infty)$ defined by $v_f(x)=\sqrt{1+|\nabla f(x)|^2}$,
\begin{align*}
\frac{1}{\mu_{\Gamma_f}(\Gamma_{f,\rho})}\int_{\Gamma_{f,\rho}}\log v_f\,d\mu_{\Gamma_f}\le B(n,\rho)(1+H_0),
\end{align*}
where $B(n,\rho)>0$ depends only on $n$ and the interior radius $\rho$. This is the exact slope-average estimate needed below, with the averaging measure specified as the induced graph measure. Applying this with $f=w$, $\rho=3/4$, and $H_0=H$ gives
\begin{align*}
\fint_{M_{3/4}}\log v\,d\mu_M\le B(n,3/4)(1+H).
\end{align*}
Set $C_2(n):=B(n,3/4)$. The hypotheses are satisfied here because $w\in C^2(B(0,1))$ solves the minimal surface equation, $B(0,3/4)\subset B(0,1)$ is the fixed interior ball, and the normalization gives the cylinder height bound $0\le w\le H$.
[guided]
The remaining input converts geometric height control into slope control on average. Define
\begin{align*}
H:=\operatorname{osc}_{B(0,1)}w.
\end{align*}
Because of the normalization in the first step, $0\le w\le H$ throughout $B(0,1)$. Therefore the part of the graph over $B(0,3/4)$ is contained in the Euclidean cylinder $B(0,3/4)\times[0,H]$.
The graph measure is computed from the parametrization $x\mapsto (x,w(x))$. With respect to $n$-dimensional Lebesgue measure $\mathcal{L}^n$ on the base, the induced volume element is
\begin{align*}
d\mu_M(x,w(x))=v(x)\,d\mathcal{L}^n(x).
\end{align*}
Thus averaging $\log v$ over $M_{3/4}$ is exactly averaging the logarithmic slope with respect to the intrinsic area measure of the graph.
Now apply the radius-dependent form of the [Bombieri-De Giorgi-Miranda logarithmic area estimate](/page/Bombieri-De%20Giorgi-Miranda%20Estimate). The theorem says that if $f\in C^2(B(0,1))$ solves the minimal surface equation, if $0<\rho<1$, if $\Gamma_f:=\{(x,f(x)):x\in B(0,1)\}$, if $\Gamma_{f,\rho}:=\Gamma_f\cap(B(0,\rho)\times\mathbb{R})$, and if $H_0\ge \sup_{B(0,1)}f-\inf_{B(0,1)}f$, then, with the map $v_f:B(0,1)\to[1,\infty)$ defined by $v_f(x)=\sqrt{1+|\nabla f(x)|^2}$,
\begin{align*}
\frac{1}{\mu_{\Gamma_f}(\Gamma_{f,\rho})}\int_{\Gamma_{f,\rho}}\log v_f\,d\mu_{\Gamma_f}\le B(n,\rho)(1+H_0).
\end{align*}
The conclusion is exactly a normalized intrinsic average over the graph, and the constant depends only on $n$ and the fixed radius $\rho$. We verify the hypotheses with $f=w$, $\rho=3/4$, and $H_0=H$. The function $w$ is in $C^2(B(0,1))$ and solves the minimal surface equation. The subdomain $B(0,3/4)$ is a fixed interior ball of $B(0,1)$. The normalization gives $0\le w\le H$, so $\sup_{B(0,1)}w-\inf_{B(0,1)}w=H$. Therefore
\begin{align*}
\fint_{M_{3/4}}\log v\,d\mu_M\le B(n,3/4)(1+H).
\end{align*}
Since the radius $3/4$ is fixed once and for all, we write $C_2(n):=B(n,3/4)$; this is a dimensional constant.
[/guided]
[/step]
[step:Combine the iteration and height estimates]
Combining the two preceding estimates gives
\begin{align*}
\sup_{B(0,1/2)}\sqrt{1+|\nabla w|^2}=\sup_{M_{1/2}}v\le C_1(n)\exp\left(C_1(n)C_2(n)(1+H)\right).
\end{align*}
The preceding two estimates have constants depending only on $n$, because the radii $1/2$ and $3/4$ are fixed. Since $H=\operatorname{osc}_{B(0,1)}w$ and $H\ge 0$, enlarging the dimensional constant gives
\begin{align*}
\sup_{B(0,1/2)}\sqrt{1+|\nabla w|^2}\le C(n)\exp\bigl(C(n)\operatorname{osc}_{B(0,1)}w\bigr).
\end{align*}
Undoing the normalization from the first step yields
\begin{align*}
\sup_{B(x_0,R/2)} \sqrt{1+|\nabla u|^2} \le C(n)\exp\left(C(n)\frac{\operatorname{osc}_{B(x_0,R)}u}{R}\right),
\end{align*}
which is the claimed Moser interior gradient estimate.
[guided]
The previous two steps give complementary bounds: Moser iteration controls the maximum of $v$ by the exponential of the logarithmic average of $v$, and the Bombieri-De Giorgi-Miranda estimate controls that logarithmic average by the height oscillation. Substituting the second estimate into the first gives
\begin{align*}
\sup_{B(0,1/2)}\sqrt{1+|\nabla w|^2}=\sup_{M_{1/2}}v\le C_1(n)\exp\left(C_1(n)C_2(n)(1+H)\right).
\end{align*}
Here $H=\operatorname{osc}_{B(0,1)}w\ge 0$. The factor depending on the additive $1$ in $1+H$ is $\exp(C_1(n)C_2(n))$, which depends only on $n$, so it can be absorbed into the prefactor. Enlarging the dimensional constant gives
\begin{align*}
\sup_{B(0,1/2)}\sqrt{1+|\nabla w|^2}\le C(n)\exp\bigl(C(n)\operatorname{osc}_{B(0,1)}w\bigr).
\end{align*}
Finally, the normalization was
\begin{align*}
w(y)=R^{-1}\bigl(u(x_0+Ry)-\inf_{B(x_0,R)}u\bigr),
\end{align*}
so $\nabla w(y)=\nabla u(x_0+Ry)$ and
\begin{align*}
\operatorname{osc}_{B(0,1)}w=\frac{\operatorname{osc}_{B(x_0,R)}u}{R}.
\end{align*}
The change of variables $x=x_0+Ry$ maps $B(0,1/2)$ onto $B(x_0,R/2)$. Therefore
\begin{align*}
\sup_{B(x_0,R/2)} \sqrt{1+|\nabla u|^2} \le C(n)\exp\left(C(n)\frac{\operatorname{osc}_{B(x_0,R)}u}{R}\right),
\end{align*}
which is the claimed Moser interior gradient estimate.
[/guided]
[/step]