Solutions to the eikonal equation Theorem. Let V be a $C^\infty$ function on a riemannian manifold $M$ and $p$ be a nondegenerate local minimum with $V(p)=0$. Then there is a unique positive function $\varphi \in C^\infty(U)$ such that $\varphi$ solves the eikonal equation
$$\|\mathrm{grad} \varphi \|^2 = V.$$
Here, $U$ is some open neighborhood of $p$.
I found this statement (at least a quite similar one) in a physics paper without a real proof, just some motivation. It seems highly nontrivial and somehow I am struggling to find a proof anywhere. Can someone give me a good reference?
\Edit: I forgot to write down the (for uniqueness obviously essential) condition that $\varphi$ is positive. Sorry everyone!
 A: Note:  I have realized that, using the Stable Manifold Theorem, one can prove the smoothness of the solution $\phi$ that I describe below.  Thus, I am modifying my answer to incorporate that.  
Local existence and uniqueness of a smooth solution near $p$ satisfying $\phi(p)=0$ and $\phi\ge0$ near $p$ follows from an application of the Stable Manifold Theorem.  Here is the argument.
This is a local question, so we might as well assume that $M=\mathbb{R}^n$, that $p=0$, that $g = g_{ij}dx^idx^j$ satisfies $g_{ij}(0)=\delta_{ij}$, and that the function $V$ has a Taylor expansion $V = h_{ij}x^ix^j + O(|x|^3)$, where $(h_{ij})$ is a symmetric, positive definite matrix.  We are looking for a closed $1$-form $d\phi = f_i\ dx^i$ near $x=0$ so that $\phi$ satisfies the equation $g^{ij}f_if_j = V$ and, at the same time, satisfies $\phi(0)=0$ and $\phi\ge0$ near $x=0$.
Let $p_i$ be the coordinates on $T^*\mathbb{R}^n$ such that the canonical $1$-form has the expression $p_i\ dx^i$.  Then the graph of $d\phi$, described by equations $p_i = f_i(x)$, will be a Lagrangian submanifold for the $2$-form $dp_i\wedge dx^i$ and will lie in the zero locus of the Hamiltonian $H(x,p) = g^{ij}(x)p_ip_j - V(x)$.  Therefore, it will be a union of integral curves of the Hamiltonian vector field
$$
X_H = 2g^{ij}p_i\frac{\partial\ \ }{\partial x^j} 
+ \left(\frac{\partial V}{\partial x^k} 
- \frac{\partial g^{ij}}{\partial x^k}p_ip_j\right)\frac{\partial\ \ }{\partial p_k}.
$$
This graph will have to pass through the unique singular point of $X_H$, i.e., $x = p = 0$ (since $\phi$ clearly must have a critical point at $x=0$ because it vanishes there and is nonnegative nearby), and the linear part of $X_H$ at $x=p=0$ is
$$
Y_H = 2\delta^{ij}p_i\frac{\partial\ \ }{\partial x^j} 
      + 2h_{ij}x^i\frac{\partial\ \ }{\partial p_j}.
$$
The unstable manifold of $Y_H$ is the $n$-dimensional submanifold defined by $p_i = L_{ij} x^j$, where $L$ is the (unique) symmetric positive definite square root of $(h_{ij})$. The stable manifold of $Y_H$ is defined by $p_i = -L_{ij}x^j$.  It follows from the Stable Manifold Theorem that $X_H$ has a smooth $n$-dimensional unstable submanifold $N_+$ given by $p_i = f_i(x) =  L_{ij}x^j + O(|x|^2)$ and a smooth $n$-dimensional stable submanifold $N_-$ given by $p_i = -L_{ij}x^j + O(|x|^2)$.
From the dynamics of $X_H$, it is clear that the only possibility for $d\phi$, when $\phi$ satisfies the above conditions and is at least $C^2$, is to have its graph be $N_+$.  Conversely, taking $d\phi = f_i(x)\ dx^i$ where $p_i = f_i(x)$ is the (necessarily Lagrangian) unstable manifold of $X_H$ and fixing the additive constant by requiring that $\phi(0)=0$ does give a smooth solution to the original equation.  
