If $\rho:\mathbb{R}^n\to\mathbb{R}$ is a defining function for a 2-dimensional surface $M$ in $\mathbb{R}^n$,i.e., $M=\{\rho=0\}$, how can one construct a local parametrization of $M$ using this defining function?
For simplicity, assume $M$ is a regular surface, $\rho$ is smooth, $|\nabla\rho|\neq0$ along $\rho=0$, etc.
I understand that one needs some kind of implicit mapping theorem but am not able to make my arguments precise.
It might be difficult to getting an explicit formula in the general case, but here's how it goes: assume that $\rho\colon \Bbb R^n \to \Bbb R$ is smooth, has $0$ as a regular value, and $\rho^{-1}[\{0\}]=M$. Let $p \in M$. Since $\nabla\rho(x) \neq 0$, one of the partial derivatives of $\rho$ does not vanish at $p$. Assume that $$\frac{\partial \rho}{\partial x^n}(p) \neq 0.$$This will give us that $M$ is locally a graph with domain in the $(x^1,\ldots, x^{n-1})$-hyperplane. The argument for the other cases is the same. Define $\varphi\colon \Bbb R^{n} \to \Bbb R$ by $$\varphi(x^1,\ldots, x^{n}) = (x^1,\ldots, x^{n-1}, \rho(x^1,\ldots, x^n)).$$We have $$D\varphi(x^1,\ldots, x^n) = \begin{pmatrix} {\rm Id}_{n-1} \\ \hline \nabla\rho(x_1,\ldots,x^n), \end{pmatrix}$$so that $$\det D\varphi(p) = \frac{\partial \rho}{\partial x^n}(p) \neq 0,$$and the Inverse Function Theorem applies: we get $U \ni p$ open such that $\varphi\colon U \to \varphi[U]$ is a diffeomorphism. The inverse necessarily has the form $$\varphi^{-1}(x^1,\ldots, x^n) = (x^1,\ldots, x^{n-1}, g(x^1,\ldots,x^n))$$for some smooth function $g$. Now, if $\pi\colon \Bbb R^n \to \Bbb R^{n-1}$ forgets the last entry of a vector, you may define ${\bf X}\colon \pi[U \cap M]\to U \cap M$ by $${\bf X}(x^1,\ldots, x^{n-1}) = (x^1,\ldots, x^{n-1}, g(x^1,\ldots, x^{n-1},0)),$$which is a parametrization for $M$ (this indeed takes values in $M$ since $\rho(x^1,\ldots, x^{n-1},g(x^1,\ldots, x^{n-1},0))=0$).
Edit: $D\varphi(p)$ is upper triangular with $1$'s in the diagonal except for the last entry, which is $(\partial \rho/\partial x^n)(p)$. When you wrote $\nabla \rho$ instead of $D\rho$ I assumed that you were using the word "surface" loosely here meaning codimension $1$. If the codimension of $M$ is $n-2$, you have $\rho \colon \Bbb R^n \to \Bbb R^{n-2}$ instead, and you can assume WLOG that the last order $n-2$ submatrix of $D\rho(p)$ is non-singular. Define $$\varphi(x^1,\ldots, x^n) = (x^1,x^2, \rho^1(x_1,\ldots, x^n),\ldots, \rho^{n-2}(x_1,\ldots, x^n))$$instead -- $\det D\varphi(p)$ will be equal to the $(n-2)\times (n-2)$ non-zero subdeterminant of $D\rho(p)$ we chose, and the IVT applies. From here on the proof goes just like in the original answer.