# Points where harmonic functions fail to give a coordinates system

Let $$\Omega$$ be a bounded domain in $$\mathbb R^n$$ with a smooth boundary and let $$g$$ be a smooth Riemannian metric on $$\Omega$$. Let $$f_1,f_2,\ldots,f_n$$ be non-constant smooth functions on $$\partial \Omega$$ that are linearly independent from each other. Given any $$j=1,2,\ldots,n$$, we denote by $$u_j$$, the unique harmonic function in $$\Omega$$ (i.e $$\Delta_gu=0$$ on $$\Omega$$) with Dirichlet data $$f_j$$.

Let us consider the set of points $$p$$ where the set $$u_1,\ldots,u_n$$ fails to give a coordinate system near $$p$$ (that is to say, $$\det J=0$$ where $$J_{jk}=\partial_j u_k(p)$$). Can we say that there is only a finite number of such points in $$\Omega$$? If not, can we say that the measure of such points is zero?

To rule out the trivial case that one of the $$f_j$$ is a nonzero constant, I interpret the linear independence condition as saying that no nontrivial linear combination of the $$f_j$$ is constant.
The set $$\{\det J = 0\}$$ is not necessarily finite. Consider for example $$u_1 = r\sin(\theta)$$ and $$u_2 = r^2\sin(2\theta)$$ in $$B_1 \subset \mathbb{R}^2$$, which both vanish on the $$x$$ axis (hence have parallel gradient there), and have linearly independent boundary data.
In dimension $$n = 2$$, the set $$\{\det J = 0\}$$ has measure zero. One can argue as follows: assume not, and let $$u,\,v$$ be the pair of harmonic functions. Since $$u_y + iu_x$$ and $$v_y + iv_x$$ are holomorphic, the sets $$\{|\nabla u| = 0\}$$ and $$\{|\nabla v| = 0\}$$ are locally finite, and furthermore, the imaginary part of their ratio, $$\frac{\nabla^{\perp}u \cdot \nabla v}{|\nabla v|^2},$$ is harmonic. By hypothesis, this function vanishes on a set of positive measure, and by its harmonicity it vanishes identically in $$\{|\nabla v| > 0\}$$. We conclude in this set that $$\nabla u = \lambda(x)\nabla v$$. Harmonicity of $$u$$ and $$v$$ implies that $$\nabla \lambda \cdot \nabla v = 0$$, and differentiating the relation $$u_x = \lambda v_x$$ in $$y$$, the relation $$u_y = \lambda v_y$$ in $$x$$, and subtracting gives $$\nabla \lambda \cdot \nabla^{\perp}v = 0$$. It follows that $$\lambda$$ is constant in $$\{|\nabla v| > 0\}$$ and hence that $$u - \lambda v$$ is constant, contradicting linear independence of the boundary data.