I'm currently trying to understand Dennis Deturck's and Jerry Kazdan's proof of the existence of isothermal coordinates for $2$-dimensional Riemannian manifolds in Some regularity theorems in Riemannian geometry, p. 251, but I'm kinda stuck. On a Riemannian manifold $(M, g)$, coordinate functions $(x_1, \dots, x_n)$ are called isothermal if the metric is locally conformal to the Euclidean metric in these coordinates, that is, if $g$ is locally of the shape $d s^2 = \lambda(x_1, \dots, x_n) (d x_1^2 + \dots + d x_n^2)$ for some $\lambda(x_1, \dots, x_n) > 0$.
Firstly, Deturck and Kazdan prove that on any Riemannian manifold, harmonic coordinate functions (where the Lacplace operator of the function is $= 0$) exist everywhere locally. Now to prove the existence of isothermal coordinates in the $2$-dimensional case, this allows to choose a harmonic function $u$ with non-vanishing derivative as the first component. So far I can follow.
Now in the next step, they choose the second coordinate function $v$ using the Hodge star operator $\star$. Since $d u$ is a $1$-form and the Hodge star sends $1$-forms to $1$-forms, $\star d u $ is also a $1$-form. Now because $\star du$ is closed, the Poincaré lemma ensures the existence of a coordinate function $v$ satisfying $d v = \star du$. Finally, they show that $v$ is also harmonic through a simple calculation.
Now my question is: how can I see that these coordinate functions $(u, v)$ constructed this way are actually isothermal? Deturck and Kazdan don't seem to elaborate on that unless I have missed it somewhere, so maybe it's very easy and just I'm missing something trivial, but I've been trying to wrap my head around why the off-diagonal entries of the metric become $= 0$ when representing it with respect to these coordinate functions. I'm not sure how exactly one can deduce that from the fact that $u$ and $v$ are harmonic and that $d v = \star d u$.
The relation $dv = \star du$ implies that $du, dv$ are orthogonal, since (letting $\omega$ denote the Riemannian area form) the definition of the Hodge star gives $$\langle du, dv \rangle \omega= \langle du, \star du\rangle \omega = du \wedge du = 0.$$
In other words, the matrix $g^{ij}$ representing the inverse metric in the coordinates $(u,v)$ is diagonal. Since the inverse of a diagonal matrix is diagonal, this means that $g_{ij}$ is diagonal too.