[proofplan]
We rewrite the minimal graph operator by differentiating the vector field $F: \mathbb{R}^n \to \mathbb{R}^n$ defined by
\begin{align*}
F(p)=\frac{p}{\sqrt{1+|p|^2}}.
\end{align*}
This produces a positive scalar factor times the matrix $a(p)$, so the equation $M[u]=0$ is equivalent to the non-divergence form equation with coefficients $a_{ij}(\nabla u)$. Finally, we compute the quadratic form of $a(p)$ and estimate it directly by separating the component of $\xi$ in the direction of $p$ from the component orthogonal to $p$.
[/proofplan]
[step:Differentiate the minimal surface vector field with respect to its argument]
Define the map $F: \mathbb{R}^n \to \mathbb{R}^n$ by
\begin{align*}
F(p)=\frac{p}{\sqrt{1+|p|^2}}.
\end{align*}
For $1 \leq i,j \leq n$, its $i$-th component is
\begin{align*}
F_i(p) = p_i(1+|p|^2)^{-1/2}.
\end{align*}
Differentiating with respect to $p_j$ gives
\begin{align*}
\partial_{p_j}F_i(p)=\delta_{ij}(1+|p|^2)^{-1/2}-p_i p_j(1+|p|^2)^{-3/2}.
\end{align*}
Factoring out $(1+|p|^2)^{-1/2}$ gives
\begin{align*}
\partial_{p_j}F_i(p)=(1+|p|^2)^{-1/2}\left(\delta_{ij}-\frac{p_i p_j}{1+|p|^2}\right).
\end{align*}
By the definition of $a_{ij}(p)$, this is
\begin{align*}
\partial_{p_j}F_i(p)=(1+|p|^2)^{-1/2}a_{ij}(p).
\end{align*}
[/step]
[step:Apply the chain rule to rewrite $M[u]$ in non-divergence form]
Fix $x \in U$. Since $u \in C^2(U)$, the gradient map $\nabla u: U \to \mathbb{R}^n$ is $C^1$, and the composition $F \circ \nabla u: U \to \mathbb{R}^n$ is $C^1$. By the definition of $M[u]$,
\begin{align*}
M[u](x)=\sum_{i=1}^n \partial_{x_i}\left(F_i(\nabla u(x))\right).
\end{align*}
The chain rule and the preceding derivative formula with $p=\nabla u(x)$ give
\begin{align*}
M[u](x)=\sum_{i=1}^n \sum_{j=1}^n \partial_{p_j}F_i(\nabla u(x))\,\partial_{x_i x_j}u(x).
\end{align*}
Substituting $\partial_{p_j}F_i(\nabla u(x))=(1+|\nabla u(x)|^2)^{-1/2}a_{ij}(\nabla u(x))$ gives
\begin{align*}
M[u](x)=\frac{1}{\sqrt{1+|\nabla u(x)|^2}}\sum_{i,j=1}^n a_{ij}(\nabla u(x))\,\partial_{x_i x_j}u(x).
\end{align*}
The scalar factor $(1+|\nabla u(x)|^2)^{-1/2}$ is strictly positive for every $x \in U$. Therefore $M[u](x)=0$ if and only if
\begin{align*}
\sum_{i,j=1}^n a_{ij}(\nabla u(x))\,\partial_{x_i x_j}u(x)=0.
\end{align*}
Since $x \in U$ was arbitrary, the two equations are equivalent in $U$.
[guided]
The goal is to pass from divergence form to non-divergence form. The divergence form differentiates the composite vector field $F(\nabla u)$, where $F: \mathbb{R}^n \to \mathbb{R}^n$ is defined by
\begin{align*}
F(p)=\frac{p}{\sqrt{1+|p|^2}}.
\end{align*}
Because $u \in C^2(U)$, the map $\nabla u: U \to \mathbb{R}^n$ is $C^1$. Since $F$ is smooth on all of $\mathbb{R}^n$, the composition $F \circ \nabla u: U \to \mathbb{R}^n$ is $C^1$, so differentiating component by component is legitimate.
Fix $x \in U$. By the definition of $M[u]$,
\begin{align*}
M[u](x)=\sum_{i=1}^n \partial_{x_i}\left(F_i(\nabla u(x))\right).
\end{align*}
For each fixed $i$, the function $F_i(\nabla u(\cdot))$ depends on $x_i$ through all components of $\nabla u$. Thus the chain rule gives
\begin{align*}
\partial_{x_i}\left(F_i(\nabla u(x))\right)=\sum_{j=1}^n \partial_{p_j}F_i(\nabla u(x))\,\partial_{x_i x_j}u(x).
\end{align*}
From the previous computation,
\begin{align*}
\partial_{p_j}F_i(\nabla u(x))=(1+|\nabla u(x)|^2)^{-1/2}a_{ij}(\nabla u(x)).
\end{align*}
Substituting this into the sum defining $M[u](x)$ yields
\begin{align*}
M[u](x)=\sum_{i=1}^n \sum_{j=1}^n (1+|\nabla u(x)|^2)^{-1/2}a_{ij}(\nabla u(x))\,\partial_{x_i x_j}u(x).
\end{align*}
Since the scalar factor is independent of the summation indices $i$ and $j$, this becomes
\begin{align*}
M[u](x)=\frac{1}{\sqrt{1+|\nabla u(x)|^2}}\sum_{i,j=1}^n a_{ij}(\nabla u(x))\,\partial_{x_i x_j}u(x).
\end{align*}
The coefficient outside the sum is never zero, because $1+|\nabla u(x)|^2>0$. Therefore multiplying by this positive scalar cannot change the zero set of the expression. Hence $M[u](x)=0$ if and only if
\begin{align*}
\sum_{i,j=1}^n a_{ij}(\nabla u(x))\,\partial_{x_i x_j}u(x)=0.
\end{align*}
Since this equivalence holds for every $x \in U$, the two equations are equivalent on $U$.
[/guided]
[/step]
[step:Compute the quadratic form of the coefficient matrix]
Fix $p \in \mathbb{R}^n$ and $\xi \in \mathbb{R}^n$. By the definition of $a_{ij}(p)$,
\begin{align*}
\sum_{i,j=1}^n a_{ij}(p)\xi_i\xi_j=\sum_{i,j=1}^n \left(\delta_{ij}-\frac{p_i p_j}{1+|p|^2}\right)\xi_i\xi_j.
\end{align*}
Using $\sum_{i,j=1}^n \delta_{ij}\xi_i\xi_j=|\xi|^2$ and collecting the product terms gives
\begin{align*}
\sum_{i,j=1}^n a_{ij}(p)\xi_i\xi_j=|\xi|^2-\frac{1}{1+|p|^2}\left(\sum_{i=1}^n p_i\xi_i\right)\left(\sum_{j=1}^n p_j\xi_j\right).
\end{align*}
Since $p\cdot \xi=\sum_{i=1}^n p_i\xi_i$, this is
\begin{align*}
\sum_{i,j=1}^n a_{ij}(p)\xi_i\xi_j=|\xi|^2-\frac{(p\cdot \xi)^2}{1+|p|^2}.
\end{align*}
Since $(p\cdot \xi)^2 \geq 0$, this identity immediately gives
\begin{align*}
\sum_{i,j=1}^n a_{ij}(p)\xi_i\xi_j \leq |\xi|^2.
\end{align*}
[/step]
[step:Bound the quadratic form from below by separating the parallel direction]
If $p=0$, then $a_{ij}(0)=\delta_{ij}$ for all $1 \leq i,j \leq n$, and therefore
\begin{align*}
\sum_{i,j=1}^n a_{ij}(0)\xi_i\xi_j = |\xi|^2 = \frac{|\xi|^2}{1+|0|^2}.
\end{align*}
Assume now that $p \neq 0$. Define
\begin{align*}
\xi_{\parallel} := \frac{p\cdot \xi}{|p|^2}p, \qquad \xi_{\perp} := \xi-\xi_{\parallel}.
\end{align*}
Then $p\cdot \xi_{\perp}=0$ and $\xi=\xi_{\parallel}+\xi_{\perp}$. The orthogonality relation gives
\begin{align*}
|\xi|^2 = |\xi_{\parallel}|^2 + |\xi_{\perp}|^2, \qquad (p\cdot \xi)^2 = |p|^2|\xi_{\parallel}|^2.
\end{align*}
Substituting these identities into the quadratic-form formula gives
\begin{align*}
\sum_{i,j=1}^n a_{ij}(p)\xi_i\xi_j=|\xi_{\parallel}|^2+|\xi_{\perp}|^2-\frac{|p|^2|\xi_{\parallel}|^2}{1+|p|^2}.
\end{align*}
Combining the two terms involving $|\xi_{\parallel}|^2$ yields
\begin{align*}
\sum_{i,j=1}^n a_{ij}(p)\xi_i\xi_j=\frac{|\xi_{\parallel}|^2}{1+|p|^2}+|\xi_{\perp}|^2.
\end{align*}
Because $1+|p|^2 \geq 1$, we have $|\xi_{\perp}|^2 \geq |\xi_{\perp}|^2/(1+|p|^2)$, and hence
\begin{align*}
\sum_{i,j=1}^n a_{ij}(p)\xi_i\xi_j \geq \frac{|\xi_{\parallel}|^2+|\xi_{\perp}|^2}{1+|p|^2}.
\end{align*}
Using $|\xi|^2 = |\xi_{\parallel}|^2 + |\xi_{\perp}|^2$, we obtain
\begin{align*}
\sum_{i,j=1}^n a_{ij}(p)\xi_i\xi_j \geq \frac{|\xi|^2}{1+|p|^2}.
\end{align*}
Together with the upper bound from the previous step, this proves
\begin{align*}
\frac{|\xi|^2}{1+|p|^2} \leq \sum_{i,j=1}^n a_{ij}(p)\xi_i\xi_j \leq |\xi|^2
\end{align*}
for every $p,\xi \in \mathbb{R}^n$. This completes the proof.
[/step]
