[proofplan]
We parametrize the graph by the map $F(x)=(x,u(x))$ and compute its first fundamental form, upward unit normal, and second fundamental form in the coordinate frame induced from $\mathbb{R}^m$. The [mean curvature scalar](/page/Mean%20Curvature) is the contraction of the [second fundamental form](/page/Second%20Fundamental%20Form) with the inverse of the [induced metric](/page/Induced%20Metric). Expanding this contraction and separately differentiating $\nabla u/\sqrt{1+|\nabla u|^2}$ produce the same expression:
\begin{align*}
\frac{1}{W}\sum_{i=1}^m \partial_{x_i}\partial_{x_i}u-\frac{1}{W^3}\sum_{i,j=1}^m \partial_{x_i}u\,\partial_{x_j}u\,\partial_{x_i}\partial_{x_j}u,
\end{align*}
where $W=\sqrt{1+|\nabla u|^2}$.
Since reversing the unit normal only changes the sign of the mean curvature scalar, vanishing of mean curvature is equivalent to the displayed equation.
[/proofplan]
[step:Parametrize the graph and compute the induced metric]
Define the graph parametrization $F: U \to \mathbb{R}^{m+1}$ by
\begin{align*}
F(x)=(x,u(x)).
\end{align*}
For each $i \in \{1,\dots,m\}$, define the coordinate tangent vector field $E_i: U \to \mathbb{R}^{m+1}$ by
\begin{align*}
E_i(x):=\partial_{x_i}F(x)=e_i+\partial_{x_i}u(x)e_{m+1},
\end{align*}
where $e_1,\dots,e_{m+1}$ denote the standard basis vectors of $\mathbb{R}^{m+1}$.
The [induced metric](/page/Induced%20Metric) matrix $g: U \to \mathbb{R}^{m \times m}$ has entries
\begin{align*}
g_{ij}(x) := E_i(x)\cdot E_j(x)
= \delta_{ij}+\partial_{x_i}u(x)\partial_{x_j}u(x).
\end{align*}
Set $W: U \to (0,\infty)$ by
\begin{align*}
W(x)=\sqrt{1+|\nabla u(x)|^2}.
\end{align*}
Then the inverse metric matrix $g^{-1}: U \to \mathbb{R}^{m \times m}$ is given by
\begin{align*}
g^{ij}(x)=\delta_{ij}-\frac{\partial_{x_i}u(x)\partial_{x_j}u(x)}{W(x)^2}.
\end{align*}
Indeed,
\begin{align*}
\sum_{j=1}^m \left(\delta_{ij}+\partial_{x_i}u\,\partial_{x_j}u\right)\left(\delta_{jk}-\frac{\partial_{x_j}u\,\partial_{x_k}u}{W^2}\right)=\delta_{ik}-\frac{\partial_{x_i}u\,\partial_{x_k}u}{W^2}+\partial_{x_i}u\,\partial_{x_k}u-\frac{\partial_{x_i}u\,\partial_{x_k}u|\nabla u|^2}{W^2}.
\end{align*}
Combining the last three terms gives
\begin{align*}
\sum_{j=1}^m \left(\delta_{ij}+\partial_{x_i}u\,\partial_{x_j}u\right)\left(\delta_{jk}-\frac{\partial_{x_j}u\,\partial_{x_k}u}{W^2}\right)=\delta_{ik}+\partial_{x_i}u\,\partial_{x_k}u\left(1-\frac{1+|\nabla u|^2}{W^2}\right).
\end{align*}
Since $W^2=1+|\nabla u|^2$, this equals $\delta_{ik}$.
[guided]
The graph is naturally described by the parametrization $F: U \to \mathbb{R}^{m+1}$ defined by
\begin{align*}
F(x)=(x,u(x)).
\end{align*}
This map is $C^2$ because $u \in C^2(U;\mathbb{R})$. Its coordinate tangent vector fields are obtained by differentiating $F$ in the coordinate directions. For each $i \in \{1,\dots,m\}$, define $E_i: U \to \mathbb{R}^{m+1}$ by
\begin{align*}
E_i(x):=\partial_{x_i}F(x)=e_i+\partial_{x_i}u(x)e_{m+1},
\end{align*}
where $e_1,\dots,e_{m+1}$ are the standard basis vectors of $\mathbb{R}^{m+1}$.
The [induced metric](/page/Induced%20Metric) records inner products of these tangent vectors. Thus, for $i,j \in \{1,\dots,m\}$,
\begin{align*}
g_{ij}(x):=E_i(x)\cdot E_j(x).
\end{align*}
Expanding the Euclidean [inner product](/page/Inner%20Product) gives
\begin{align*}
g_{ij}(x)=\left(e_i+\partial_{x_i}u(x)e_{m+1}\right)\cdot\left(e_j+\partial_{x_j}u(x)e_{m+1}\right).
\end{align*}
Using orthonormality of the standard basis, this becomes
\begin{align*}
g_{ij}(x)=\delta_{ij}+\partial_{x_i}u(x)\partial_{x_j}u(x).
\end{align*}
The matrix form is therefore $g=I+\nabla u\otimes \nabla u$.
Define the positive function $W: U \to (0,\infty)$ by
\begin{align*}
W(x)=\sqrt{1+|\nabla u(x)|^2}.
\end{align*}
The inverse metric is the rank-one inverse of $I+\nabla u\otimes \nabla u$:
\begin{align*}
g^{ij}(x)=\delta_{ij}-\frac{\partial_{x_i}u(x)\partial_{x_j}u(x)}{W(x)^2}.
\end{align*}
To verify this formula, multiply the proposed inverse by $g$. Suppressing the common point $x$ only in this displayed computation, we get
\begin{align*}
\sum_{j=1}^m \left(\delta_{ij}+\partial_{x_i}u\,\partial_{x_j}u\right)\left(\delta_{jk}-\frac{\partial_{x_j}u\,\partial_{x_k}u}{W^2}\right)=\delta_{ik}-\frac{\partial_{x_i}u\,\partial_{x_k}u}{W^2}+\partial_{x_i}u\,\partial_{x_k}u-\frac{\partial_{x_i}u\,\partial_{x_k}u|\nabla u|^2}{W^2}.
\end{align*}
Collecting the terms with the factor $\partial_{x_i}u\,\partial_{x_k}u$ gives
\begin{align*}
\sum_{j=1}^m \left(\delta_{ij}+\partial_{x_i}u\,\partial_{x_j}u\right)\left(\delta_{jk}-\frac{\partial_{x_j}u\,\partial_{x_k}u}{W^2}\right)=\delta_{ik}+\partial_{x_i}u\,\partial_{x_k}u\left(1-\frac{1+|\nabla u|^2}{W^2}\right).
\end{align*}
Since $W^2=1+|\nabla u|^2$, the coefficient in parentheses is zero, and the product is $\delta_{ik}$. Hence the displayed matrix is the inverse metric.
[/guided]
[/step]
[step:Compute the second fundamental form using the upward unit normal]
Define the upward unit normal field $\nu: U \to \mathbb{R}^{m+1}$ by
\begin{align*}
\nu(x)=\frac{1}{W(x)}\left(-\partial_{x_1}u(x),\dots,-\partial_{x_m}u(x),1\right).
\end{align*}
For each $i$,
\begin{align*}
\nu(x)\cdot E_i(x)
=
\frac{-\partial_{x_i}u(x)+\partial_{x_i}u(x)}{W(x)}
=0,
\end{align*}
and
\begin{align*}
|\nu(x)|^2
=
\frac{|\nabla u(x)|^2+1}{W(x)^2}
=1.
\end{align*}
Thus $\nu$ is a unit normal field along the graph.
Define the [second fundamental form](/page/Second%20Fundamental%20Form) coefficients $h_{ij}: U \to \mathbb{R}$ by
\begin{align*}
h_{ij}(x):=\partial_{x_i}\partial_{x_j}F(x)\cdot \nu(x).
\end{align*}
Since
\begin{align*}
\partial_{x_i}\partial_{x_j}F(x)=\partial_{x_i}\partial_{x_j}u(x)e_{m+1},
\end{align*}
we obtain
\begin{align*}
h_{ij}(x)=\frac{\partial_{x_i}\partial_{x_j}u(x)}{W(x)}.
\end{align*}
[/step]
[step:Contract the second fundamental form to identify the mean curvature scalar]
With respect to the upward unit normal $\nu$, we use the convention that the [mean curvature scalar](/page/Mean%20Curvature) is the full metric trace of the [second fundamental form](/page/Second%20Fundamental%20Form), without the normalizing factor $1/m$. Define $H_\nu: U \to \mathbb{R}$ by
\begin{align*}
H_\nu(x):=\sum_{i,j=1}^m g^{ij}(x)h_{ij}(x).
\end{align*}
Changing to the convention with the factor $1/m$, or changing the sign convention for the second fundamental form, multiplies $H_\nu$ by a nonzero constant or by $-1$, and therefore does not change the vanishing condition. Using the formulas for $g^{ij}$ and $h_{ij}$ gives
\begin{align*}
H_\nu=\sum_{i,j=1}^m \left(\delta_{ij}-\frac{\partial_{x_i}u\,\partial_{x_j}u}{W^2}\right)\frac{\partial_{x_i}\partial_{x_j}u}{W}.
\end{align*}
Expanding the contraction gives
\begin{align*}
H_\nu=\frac{1}{W}\sum_{i=1}^m \partial_{x_i}\partial_{x_i}u-\frac{1}{W^3}\sum_{i,j=1}^m \partial_{x_i}u\,\partial_{x_j}u\,\partial_{x_i}\partial_{x_j}u.
\end{align*}
Now compute the divergence of $\nabla u/W$. Since $u \in C^2(U;\mathbb{R})$ and $W \in C^1(U;(0,\infty))$,
\begin{align*}
\operatorname{div}\left(\frac{\nabla u}{W}\right)=\sum_{i=1}^m \partial_{x_i}\left(\frac{\partial_{x_i}u}{W}\right).
\end{align*}
Using the product rule gives
\begin{align*}
\operatorname{div}\left(\frac{\nabla u}{W}\right)=\frac{1}{W}\sum_{i=1}^m \partial_{x_i}\partial_{x_i}u-\frac{1}{W^2}\sum_{i=1}^m \partial_{x_i}u\,\partial_{x_i}W.
\end{align*}
Differentiating $W=(1+|\nabla u|^2)^{1/2}$ yields
\begin{align*}
\partial_{x_i}W
=
\frac{1}{W}\sum_{j=1}^m
\partial_{x_j}u\,\partial_{x_i}\partial_{x_j}u.
\end{align*}
Substituting this identity into the divergence formula gives
\begin{align*}
\operatorname{div}\left(\frac{\nabla u}{W}\right)=\frac{1}{W}\sum_{i=1}^m \partial_{x_i}\partial_{x_i}u-\frac{1}{W^3}\sum_{i,j=1}^m \partial_{x_i}u\,\partial_{x_j}u\,\partial_{x_i}\partial_{x_j}u.
\end{align*}
Comparing this formula with the preceding expression for $H_\nu$, we get
\begin{align*}
\operatorname{div}\left(\frac{\nabla u}{W}\right)=H_\nu.
\end{align*}
[guided]
The [mean curvature scalar](/page/Mean%20Curvature) is obtained by taking the trace of the [second fundamental form](/page/Second%20Fundamental%20Form) with respect to the [induced metric](/page/Induced%20Metric), not with respect to the Euclidean metric on the parameter domain. We use the convention that this scalar is the full trace, without the normalizing factor $1/m$. The convention with the extra factor $1/m$, or the opposite sign convention for the second fundamental form, gives an equation multiplied by a nonzero constant or by $-1$, so the zero set of the mean curvature is unchanged.
First recall why the second fundamental form has the displayed form. For the graph parametrization $F: U \to \mathbb{R}^{m+1}$, the coordinate tangent fields are $E_i: U \to \mathbb{R}^{m+1}$ with
\begin{align*}
E_i(x)=e_i+\partial_{x_i}u(x)e_{m+1}.
\end{align*}
The upward unit normal field is $\nu: U \to \mathbb{R}^{m+1}$ given by
\begin{align*}
\nu(x)=\frac{1}{W(x)}\left(-\partial_{x_1}u(x),\dots,-\partial_{x_m}u(x),1\right).
\end{align*}
This vector is orthogonal to every $E_i(x)$ and has Euclidean norm one because $W(x)^2=1+|\nabla u(x)|^2$. Define $h_{ij}: U \to \mathbb{R}$ by
\begin{align*}
h_{ij}(x):=\partial_{x_i}\partial_{x_j}F(x)\cdot \nu(x).
\end{align*}
Since $\partial_{x_i}\partial_{x_j}F(x)=\partial_{x_i}\partial_{x_j}u(x)e_{m+1}$ and the last component of $\nu(x)$ is $1/W(x)$, we get
\begin{align*}
h_{ij}(x)=\frac{\partial_{x_i}\partial_{x_j}u(x)}{W(x)}.
\end{align*}
The inverse metric coefficients are
\begin{align*}
g^{ij}(x)=\delta_{ij}-\frac{\partial_{x_i}u(x)\partial_{x_j}u(x)}{W(x)^2}.
\end{align*}
Thus, with respect to the upward normal field $\nu$, define $H_\nu: U \to \mathbb{R}$ by
\begin{align*}
H_\nu(x)=\sum_{i,j=1}^m g^{ij}(x)h_{ij}(x).
\end{align*}
Substituting the formulas for $g^{ij}$ and $h_{ij}$ gives
\begin{align*}
H_\nu=\sum_{i,j=1}^m \left(\delta_{ij}-\frac{\partial_{x_i}u\,\partial_{x_j}u}{W^2}\right)\frac{\partial_{x_i}\partial_{x_j}u}{W}.
\end{align*}
Expanding this contraction gives
\begin{align*}
H_\nu=\frac{1}{W}\sum_{i=1}^m \partial_{x_i}\partial_{x_i}u-\frac{1}{W^3}\sum_{i,j=1}^m \partial_{x_i}u\,\partial_{x_j}u\,\partial_{x_i}\partial_{x_j}u.
\end{align*}
Now we compare this expression with the divergence in the theorem statement. Define the vector field $V: U \to \mathbb{R}^m$ by
\begin{align*}
V(x)=\frac{\nabla u(x)}{W(x)}.
\end{align*}
Since $u \in C^2(U;\mathbb{R})$, we have $\nabla u \in C^1(U;\mathbb{R}^m)$, and since $W$ is positive and $C^1$, the vector field $V$ is $C^1$. Therefore its divergence is computed componentwise:
\begin{align*}
\operatorname{div}V
=
\sum_{i=1}^m \partial_{x_i}\left(\frac{\partial_{x_i}u}{W}\right).
\end{align*}
Using the product rule for each summand gives
\begin{align*}
\operatorname{div}V
=
\frac{1}{W}\sum_{i=1}^m \partial_{x_i}\partial_{x_i}u
-
\frac{1}{W^2}\sum_{i=1}^m \partial_{x_i}u\,\partial_{x_i}W.
\end{align*}
It remains to compute $\partial_{x_i}W$. Because
\begin{align*}
W=(1+|\nabla u|^2)^{1/2},
\end{align*}
the chain rule gives
\begin{align*}
\partial_{x_i}W=\frac{1}{2}(1+|\nabla u|^2)^{-1/2}\partial_{x_i}\left(1+\sum_{j=1}^m(\partial_{x_j}u)^2\right).
\end{align*}
Since $W=(1+|\nabla u|^2)^{1/2}$, this becomes
\begin{align*}
\partial_{x_i}W=\frac{1}{2W}\sum_{j=1}^m 2\partial_{x_j}u\,\partial_{x_i}\partial_{x_j}u.
\end{align*}
After cancelling the factor $2$, we obtain
\begin{align*}
\partial_{x_i}W=\frac{1}{W}\sum_{j=1}^m \partial_{x_j}u\,\partial_{x_i}\partial_{x_j}u.
\end{align*}
Substituting this into the divergence computation gives
\begin{align*}
\operatorname{div}\left(\frac{\nabla u}{W}\right)=\frac{1}{W}\sum_{i=1}^m \partial_{x_i}\partial_{x_i}u-\frac{1}{W^2}\sum_{i=1}^m \partial_{x_i}u\left(\frac{1}{W}\sum_{j=1}^m \partial_{x_j}u\,\partial_{x_i}\partial_{x_j}u\right).
\end{align*}
Equivalently,
\begin{align*}
\operatorname{div}\left(\frac{\nabla u}{W}\right)=\frac{1}{W}\sum_{i=1}^m \partial_{x_i}\partial_{x_i}u-\frac{1}{W^3}\sum_{i,j=1}^m \partial_{x_i}u\,\partial_{x_j}u\,\partial_{x_i}\partial_{x_j}u.
\end{align*}
This is exactly the expression obtained for $H_\nu$. Hence
\begin{align*}
H_\nu
=
\operatorname{div}\left(\frac{\nabla u}{\sqrt{1+|\nabla u|^2}}\right).
\end{align*}
[/guided]
[/step]
[step:Conclude that minimality is equivalent to the divergence equation]
Define the graph set $\Gamma_u \subset \mathbb{R}^{m+1}$ by
\begin{align*}
\Gamma_u:=\{(x,u(x)) : x \in U\}.
\end{align*}
The graph $\Gamma_u$ is a [minimal hypersurface](/page/Minimal%20Hypersurface) precisely when its [mean curvature scalar](/page/Mean%20Curvature) vanishes identically for either choice of unit normal. For the upward unit normal $\nu$ computed above,
\begin{align*}
H_\nu
=
\operatorname{div}\left(\frac{\nabla u}{\sqrt{1+|\nabla u|^2}}\right).
\end{align*}
Therefore $H_\nu=0$ on $U$ if and only if
\begin{align*}
\operatorname{div}\left(\frac{\nabla u}{\sqrt{1+|\nabla u|^2}}\right)=0
\end{align*}
on $U$. Replacing $\nu$ by $-\nu$ replaces $H_\nu$ by $-H_\nu$, so the vanishing condition is independent of the choice of orientation. This proves the equivalence and completes the proof.
[guided]
Define the graph set $\Gamma_u \subset \mathbb{R}^{m+1}$ by
\begin{align*}
\Gamma_u:=\{(x,u(x)) : x \in U\}.
\end{align*}
By definition, $\Gamma_u$ is minimal precisely when its [mean curvature scalar](/page/Mean%20Curvature) vanishes identically. The choice of orientation does not affect this vanishing condition: replacing the upward unit normal $\nu$ by the downward unit normal $-\nu$ replaces the scalar mean curvature $H_\nu$ by $-H_\nu$, and a function is identically zero if and only if its negative is identically zero.
From the preceding computation with the upward unit normal, we have the identity
\begin{align*}
H_\nu
=
\operatorname{div}\left(\frac{\nabla u}{\sqrt{1+|\nabla u|^2}}\right)
\end{align*}
on $U$. Therefore $H_\nu=0$ on $U$ if and only if
\begin{align*}
\operatorname{div}\left(\frac{\nabla u}{\sqrt{1+|\nabla u|^2}}\right)=0
\end{align*}
on $U$. This is exactly the equivalence asserted in the theorem statement.
[/guided]
[/step]
