[proofplan]
We regard the graph of $u$ as a complete oriented minimal surface in $\mathbb R^3$ with upward unit normal. The minimal graph equation gives zero mean curvature, and the global conformal parametrization theorem for entire minimal graphs gives conformal coordinates on the whole complex plane. In those coordinates the stereographic Gauss map is a bounded [holomorphic function](/page/Holomorphic%20Function), hence constant by [Liouville's theorem](/page/Liouville's%20Theorem). A constant normal vector forces the graph to lie in an affine plane, and the graph condition then says exactly that $u$ is affine.
[/proofplan]
[step:Convert the analytic equation into a complete oriented minimal graph]
Define the graph parametrization $X_0: \mathbb R^2 \to \mathbb R^3$ by
\begin{align*}
X_0(x_1,x_2) = (x_1,x_2,u(x_1,x_2)).
\end{align*}
Since $u \in C^2(\mathbb R^2)$, the map $X_0$ is a $C^2$ immersion. Define its upward unit normal $N_0: \mathbb R^2 \to S^2$ by
\begin{align*}
N_0(x_1,x_2) = \frac{(-\partial_{x_1}u(x_1,x_2),-\partial_{x_2}u(x_1,x_2),1)}{\sqrt{1+|\nabla u(x_1,x_2)|^2}}.
\end{align*}
The third component satisfies $N_{0,3}(x_1,x_2)>0$ for every $(x_1,x_2) \in \mathbb R^2$. By the minimal graph equation characterization of zero mean curvature, the displayed divergence equation is equivalent to the graph surface $\Sigma := X_0(\mathbb R^2)$ having mean curvature zero with respect to $N_0$.
The induced Riemannian metric $g_0 := X_0^*\langle \cdot,\cdot\rangle_{\mathbb R^3}$ is complete. Let $L_{g_0}$ denote the Riemannian length functional associated to the metric $g_0$ on $\mathbb R^2$, equivalently the length of the corresponding lifted curve in $\Sigma$ under $X_0$. Indeed, if $\gamma:[0,T) \to \mathbb R^2$ is a $C^1$ curve leaving every compact subset of $\mathbb R^2$, then its graph lift $X_0\circ \gamma:[0,T)\to \Sigma$ has length
\begin{align*}
L_{g_0}(X_0\circ \gamma) = \int_0^\top \sqrt{|\gamma'(t)|^2 + |\nabla u(\gamma(t))\cdot \gamma'(t)|^2}\, d\mathcal L^1(t).
\end{align*}
Since the square root is bounded below by $|\gamma'(t)|$ for every $t \in [0,T)$, we obtain
\begin{align*}
L_{g_0}(X_0\circ \gamma) \geq \int_0^\top |\gamma'(t)|\, d\mathcal L^1(t).
\end{align*}
The Euclidean length on the right is infinite for every curve leaving every compact subset of the complete Euclidean plane $\mathbb R^2$. Hence every divergent curve in $(\mathbb R^2,g_0)$ has infinite length, so $g_0$ is complete.
[/step]
[step:Choose global conformal coordinates on the entire graph]
Since $\Sigma$ is an entire minimal graph over $\mathbb R^2$, the projection $\pi:\Sigma\to\mathbb R^2$, defined by $\pi(x_1,x_2,x_3)=(x_1,x_2)$, is a global $C^2$ diffeomorphism with inverse $X_0$. The surface is therefore connected and simply connected. The preceding step shows that the induced metric is complete, and the same step gives zero mean curvature. By the global conformal parametrization theorem for complete entire minimal graphs, the conformal type is the complex plane rather than the disk, so there exists a conformal $C^2$ diffeomorphism $\Phi: \mathbb C \to \Sigma$.
Define the conformal immersion $X: \mathbb C \to \mathbb R^3$ by
\begin{align*}
X(z) = \Phi(z).
\end{align*}
Define the normal map $N: \mathbb C \to S^2$ by
\begin{align*}
N(z) = N_0(\pi(\Phi(z))).
\end{align*}
Because $N_0$ is the upward normal on the graph, $N_3(z)>0$ for every $z \in \mathbb C$.
[/step]
[step:Apply Liouville's theorem to the bounded holomorphic Gauss map]
Define the stereographic Gauss map from the south pole as the function $g: \mathbb C \to \mathbb C$ given by
\begin{align*}
g(z) = \frac{N_1(z)+iN_2(z)}{1+N_3(z)}.
\end{align*}
The denominator never vanishes because $N_3(z)>0$. Since $X:\mathbb C\to\mathbb R^3$ is a conformal minimal immersion, the holomorphicity theorem for the stereographic Gauss map of a conformal minimal immersion implies that $g$ is holomorphic on $\mathbb C$. Moreover, using $N_1(z)^2+N_2(z)^2+N_3(z)^2=1$, we obtain
\begin{align*}
|g(z)|^2 = \frac{N_1(z)^2+N_2(z)^2}{(1+N_3(z))^2}.
\end{align*}
Thus
\begin{align*}
|g(z)|^2 = \frac{1-N_3(z)^2}{(1+N_3(z))^2} = \frac{1-N_3(z)}{1+N_3(z)} < 1.
\end{align*}
Thus $g$ is bounded. By [Liouville's theorem](/theorems/38) for bounded entire holomorphic functions, $g$ is constant on $\mathbb C$.
[guided]
The reason for introducing $g$ is that the minimal surface equation is nonlinear in $u$, while the Gauss map becomes holomorphic after passing to conformal coordinates. We use stereographic projection from the south pole because the upward normal satisfies $N_3>0$. Define $g: \mathbb C \to \mathbb C$ by
\begin{align*}
g(z) = \frac{N_1(z)+iN_2(z)}{1+N_3(z)}.
\end{align*}
This formula is defined everywhere on $\mathbb C$ because $1+N_3(z)>1$ for every $z \in \mathbb C$.
We now verify the hypotheses of the holomorphicity theorem for the Gauss map. The map $X:\mathbb C\to\mathbb R^3$ is a conformal immersion by construction of the global conformal coordinates, and its image has zero mean curvature because the original graph satisfies the minimal graph equation. Therefore the holomorphicity theorem for the stereographic Gauss map of a conformal minimal immersion applies and gives that $g$ is holomorphic on $\mathbb C$.
Next we verify boundedness. Since $N(z)\in S^2$, its components satisfy
\begin{align*}
N_1(z)^2+N_2(z)^2+N_3(z)^2=1.
\end{align*}
Therefore
\begin{align*}
|g(z)|^2 = \left|\frac{N_1(z)+iN_2(z)}{1+N_3(z)}\right|^2.
\end{align*}
Using $|a+ib|^2=a^2+b^2$ for [real numbers](/page/Real%20Numbers) $a$ and $b$, we get
\begin{align*}
|g(z)|^2 = \frac{N_1(z)^2+N_2(z)^2}{(1+N_3(z))^2}.
\end{align*}
Substituting $N_1(z)^2+N_2(z)^2=1-N_3(z)^2$ gives
\begin{align*}
|g(z)|^2 = \frac{1-N_3(z)^2}{(1+N_3(z))^2} = \frac{(1-N_3(z))(1+N_3(z))}{(1+N_3(z))^2} = \frac{1-N_3(z)}{1+N_3(z)}.
\end{align*}
Because $N_3(z)>0$, the last quotient is strictly less than $1$. Hence $|g(z)|<1$ for every $z\in\mathbb C$.
[Liouville's theorem](/theorems/346) for bounded entire holomorphic functions requires a holomorphic function on all of $\mathbb C$ whose absolute value is bounded above by a finite constant. We have just verified holomorphicity and the bound $|g|<1$. Hence Liouville's theorem implies that $g$ is constant on $\mathbb C$.
[/guided]
[/step]
[step:Turn the constant Gauss map into an affine graph]
Stereographic projection from the south pole is injective on $S^2\setminus\{(0,0,-1)\}$. Since $N_3>0$, the image of $N$ lies in this domain, and the constancy of $g$ implies that there exists a fixed vector $a=(a_1,a_2,a_3)\in S^2$ such that $N(z)=a$ for every $z\in\mathbb C$. Transporting this identity through the parametrization $\Phi$ gives $N_0(x_1,x_2)=a$ for every $(x_1,x_2)\in\mathbb R^2$.
From the explicit formula for $N_0$, the third component is
\begin{align*}
a_3 = \frac{1}{\sqrt{1+|\nabla u(x_1,x_2)|^2}}>0,
\end{align*}
so $a_3\neq 0$. The first two components give
\begin{align*}
\partial_{x_1}u(x_1,x_2) = -\frac{a_1}{a_3} \quad \text{and} \quad \partial_{x_2}u(x_1,x_2) = -\frac{a_2}{a_3}
\end{align*}
for every $(x_1,x_2)\in\mathbb R^2$. Define constants $b_1 := -a_1/a_3$, $b_2 := -a_2/a_3$, and $c := u(0,0)$. Then $\nabla u=(b_1,b_2)$ on the connected set $\mathbb R^2$.
Fix $(x_1,x_2)\in\mathbb R^2$, and define the $C^1$ function $h:[0,1]\to\mathbb R$ by $h(t)=u(tx_1,tx_2)$. Applying the [fundamental theorem of calculus](/theorems/632) on $[0,1]$ gives
\begin{align*}
u(x_1,x_2)-u(0,0) = \int_0^1 h'(t)\, d\mathcal L^1(t).
\end{align*}
By the chain rule and the identity $\nabla u=(b_1,b_2)$, we have $h'(t)=b_1x_1+b_2x_2$ for every $t\in[0,1]$. Therefore
\begin{align*}
u(x_1,x_2)-u(0,0) = \int_0^1 \bigl(b_1x_1+b_2x_2\bigr)\, d\mathcal L^1(t) = b_1x_1+b_2x_2.
\end{align*}
Thus
\begin{align*}
u(x_1,x_2)=b_1x_1+b_2x_2+c
\end{align*}
for every $(x_1,x_2)\in\mathbb R^2$. Hence $u$ is affine, and equivalently the entire minimal graph $\Sigma$ is the affine plane $\{(x_1,x_2,b_1x_1+b_2x_2+c):(x_1,x_2)\in\mathbb R^2\}$.
[/step]