Harmonic functions without critical points (global isothermal coordinates) Let $(M,g)$ be a compact, orientable Riemannian surface with non-empty boundary $\partial M$.
Question. Does there always exist a smooth function $f:M\rightarrow \mathbb{R}$ with $\Delta_gf=0$ in the interior of $M$ and $d_pf\neq 0$ for all $p\in M$?

*

*For $M\subset \mathbb{R}^2$, equipped with the Euclidean metric, this is clear (just take an affine linear map).

*Given a finite subset $P\subset M$, one can always find a smooth function $f:M\rightarrow \mathbb{R}$ such that, for all $p\in P$ we have $\Delta f = 0$ near $p$ and $d_pf\neq 0$. (This follows from a well known trick used to construct isothermal coordinates, see below.)

*If $M$ is contractible, this is equivalent to the existence of global isothermal coordinates on $(M,g)$.


Proof of 2)

  Extend $M$ to a closed surface $(N,g)$ and write $R:H_\perp^3(N)\rightarrow \mathbb{R}^{\vert P \vert}$ for the map that sends $f$ to the list $(\vert d_p f\vert: p\in P)\in \mathbb{R}^{\vert P \vert}$. Here $H^s_\perp(N)$ consists of Sobolev-functions of regularity $s$ with zero mean (i.e. $\perp\{\mathrm{constants}\}$). Further write $F:H_\perp^1(N)\rightarrow H_\perp^3(N)$ for the map that sends $h$ to the solution $f$ of $\Delta f = h$. Now the point is that the set $D_P=\{h\in H_\perp^1(N)\cap C^\infty(N):h \text{ vanishes near } P\}$ is dense in $H_\perp^1(N)$, which implies that $RF(D_P)\subset \mathbb{R}^{\vert P\vert}$ is dense (as both $R$ and $F$ are continuous and surjective). In particular there is a vector in $RF(D_P)\subset\mathbb{R}^{\vert P\vert}$ with all coordinates non-zero and an $R$-preimage $f\in F(D_P)$ satisfies the desired requirements.

 A: Let $k\ge 1$ be the number of boundary components of $M$ and $\gamma$ the genus of the closed surface $\Sigma =\Sigma_M$ obtained by gluing a disk to each boundary component.
Lemma.
Suppose $\gamma =0$ and $k\ge 1$. Then there is a conformal diffeomorphism $F:M\rightarrow F(M)\subset \mathbb{R}^2$ onto a a smooth domain $F(M)$. For $k=1$ this can be chosen to be the unit-disk.
Proof. Extend $g$ arbitrarily to all of $\Sigma$ and denote the extension also with $g$. We now use the fact that all Riemannian metrics on $S^2$ are conformally equivalent (or equivalently there is only one complex structure) to obtain a conformal diffeomorphism $f:(\Sigma,g)\rightarrow (S^2,g_0)$ (with $g_0$  the standard round metric). Without loss of generality we may assume that the north-pole $p$ does not lie in the image $f(M)$ and we consider the stereographic projection $P:(S^2\backslash p,g_0)\xrightarrow{\sim} \mathbb{C}$, which is well known to be conformal. In particular, $F=P\circ f\vert_{M}:(M,g)\rightarrow \mathbb{C}$ is a conformal diffeomorphism onto the image $F(M)\subset \mathbb{C}\equiv \mathbb{R}^2$. The assertion about $k=0$ is just the Riemann mapping theorem. q.e.d.
This answeres the question for $\gamma =0$: Just take any harmonic map $f_0:F(M)\rightarrow \mathbb{R}$ without critical points and pull it back via $F$. As $F$ is a conformal map between surfaces, $f=F^*f_0$ is harmonic as well (Section 2.3 in these notes)
and it clearly has no critical points.
For $\gamma >0$ I actually suspect that the question has to be answered negatively, but I don't know how to prove it yet. I think that in that case $(M,g)$ necessarily has a closed geodesic $c:S^1\rightarrow M$ and one might get a contradiction from looking at $f\circ c$, where $f$ is the supposed harmonic map.
