Harmonic coordinate functions in minimal surfaces 
Let $\Sigma\subset\mathbb{R}^3$ be a regular surface. Prove that $\Sigma$ is minimal if and only if its coordinate functions are harmonic.

I know that every regular surface admits an isothermal parametrization. Assuming that a parametrization $X:U\subset\mathbb{R}^2\to\Sigma$ is isothermal (i.e., $<X_u,X_u>=<X_v,X_v>$ and $<X_u,X_v>=0$), I know how to prove $X_{uu}+X_{vv}=2<X_u,X_u>^2HN$, where $H$ is the mean curvature and $N=\frac{X_u\wedge X_v}{||X_u\wedge X_v||}$. So obviously $H\equiv 0 \iff X_{uu}+X_{vv}\equiv 0$.
The problem doesn't say anything about $X$ being isothermal, so I'm stuck with the general case. I tried to use brute force, but it was such a ridiculous amount of work that I eventuay gave up. 
Is there some clever way around it?
 A: The correct result is "if $X$ is isothermal, then $\Sigma$ is minimal if and only if $X$ is harmonic". 
Consider the following counter-example: take $w_1,w_2 \in \Bbb R^3$ linearly independent such that $\|w_1\| = 1$ and $\|w_2\| = 2$, and put $X\colon \Bbb R^2 \to \Sigma \doteq {\rm span}(w_1,w_2)$ given by $X(u,v) = uw_1+vw_2$. Clearly $\Sigma$ is minimal but $\langle X_u,X_u\rangle \neq \langle X_v,X_v\rangle$.
The bottom line is: you already "won" at $H \equiv 0  \iff X_{uu}+X_{vv} = 0$.
A: I am not sure if you have background of Riemannian geometry, following is what I do.
Suppose this surface $(\Sigma^2,g)$ which is parametrized by $$F=(f_1,f_2): \ \ (x^1,x^2)\in D\subset{R}^2\rightarrow R^3$$
If we can prove the fact that $\Delta_\Sigma F=(\Delta_{\Sigma}f_1,\Delta_\Sigma f_2)= 2\vec{H}$, then the answer follows directly. In fact, it is not hard to see that $\Delta_{\Sigma}F$ is normal to surface since
$$\begin{eqnarray*}<\Delta_{\Sigma}F,\partial_lF>&=&<g^{ij}\frac{\partial^2F}{\partial x^i\partial x^j}-g^{ij}\Gamma_{ij}^k\partial_kF,\partial_l F>\\
&=&g^{ij}<\partial_i\partial_jF,\partial_lF>-g^{ij}\Gamma_{ij}^k<\partial_kF,\partial_lF>\\
&=&g^{ij}g(\nabla_{\partial_i}\partial_j,\partial_l)-g^{ij}\Gamma_{ij}^k<\partial_kF,\partial_lF>\\
&=&g^{ij}\Gamma_{ij}^kg_{kl}-g^{ij}\Gamma_{ij}^kg_{kl}\\
&=&0
\end{eqnarray*}$$
That means we just choose the normal part of $\Delta_\Sigma F$, so
$$
\begin{eqnarray*}\Delta_{\Sigma}F&=&(\Delta_{\Sigma}F)^N\\
&=& (\frac{1}{\sqrt{|g|}}\partial_i(\sqrt{|g|}g^{ij}\partial_jF) )^N\\
&=& (g^{ij}\partial_i\partial_j F)^N+(\frac{1}{\sqrt{|g|}}\partial_i(\sqrt{|g|}g^{ij})\partial_jF)^N\\
&=&g^{ij}(\partial_i\partial_jF)^N\\
&=& 2\vec{H}
\end{eqnarray*}$$
