A conformal map whose Jacobian vanishes at a point is constant? Let $f:M \to N$ be a smooth weakly conformal map between connected $d$-dimensional Riemannian manifolds, i.e. $f$ satisfies $df^Tdf =(\det df)^{\frac{2}{d}} \, \text{Id}_{TM}$. 

Assume $d \ge 3$ and that $df=0$ at some point.  Is it true that $f$ is constant?

A proof for the Euclidean case, can be found in "Geometric Function Theory and Non-linear Analysis", by Iwaniec and Martin. Their proof uses the fact both $M,N$ are Euclidean.
If I am not mistaken, a conformal map (in $d\ge3$) is determined by its 2-jet at a point, so it suffices to prove all the second derivatives at the point where $df$ vanishes are zero. (But maybe this wrong; it's possible that the "2-jet determination" only holds for conformal maps whose differentials are everywhere non-zero. I am not sure.)
Edit:
David E Speyer proved below that if $f$ is conformal and $df_{x_0}=0$ then all the derivatives of $f$ vanish at $x_0$ (for all orders). Thus, if the above statement about the two-jet-determination is true, the only solutions are indeed constants.
Does anybody know whether or not this "two-jet-determination" hold in the case where we allow the Jacobian to vanish?
 A: I am no expert in this field, but I just looked into Harmonic Morphisms Between Riemannian Manifolds by Paul Baird and John C. Wood, specifically into chapter 14. Here they give more general definition of weakly conformal mapping (which it seems is the standard definition). Namely, $f: M \to N$ is weakly conformal at $x\in M$ if
$$
\langle df(X), df(Y) \rangle_{f(x)}= \Lambda(x) \langle X, Y \rangle_x
$$
for some number $\Lambda(x)$.
The Proposition 14.4.3 states that $f$ is weakly conformal at $x$ if one of the three mutually exclusive conditions hold, one of them being $df_x = 0$ and it falls into the non-degenerate cases. Few paragraphs later, on page 442, the authors say that non-degenerate weakly conformal maps behave in a similar way to the Riemannian case, namely any non-degenerate mapping is constant provided $\dim M > \dim N.$ I do not see a proof nearby, but then again, I just skimmed the text.
edit:
I don't know why but in the first time I looked at this question I got the impression that it's about semi-Riemannian manifolds. The aforementioned book actually state what the OP is asking for. Unless I misread the question again.

Theorem 11.4.6: Let $f \colon M \to N$ be a non-constant weakly conformal map between manifolds of the same dimension $n \geq 3.$ Then $f$ has no critical points, i.e., $f$ is a conformal local diffeomorphism.

Let me also add that $2$-jet determinacy holds in quite general contexts, see e.g. Jet-determination of symmetries of parabolic geometries and Another proof of the Liouville theorem, but it seems that the nondegeneracy of $\mathrm{d} f$ is always assumed.
A: I am not at all an expert in this field. However, I believe I can show that there are no examples where $Df$ vanishes to finite order. Thus, in particular there are no analytic solutions and, if the OP is right about two-jet-determination, there are no solutions.
Let $f: M \to N$ be a weakly conformal map between Riemannian manifolds of dimension $d>2$. We write $g_M$ and $g_N$ for the metrics on $M$ and $N$, so the conformal condition is that $(Df)^T g^N (Df) = \lambda g^M$ for some scalar valued function $\lambda$ on $M$. If $f$ is singular at $x_0 \in M$, then we must have $\lambda(x_0)=0$, so we must have $Df=0$. We say that $f$ vanishes to order $k$ at $x_0$ if all $r$-fold derivatives of $f$ vanish at $x_0$ for $r<k$ from some $k$-fold derivative doesn't. So we have shown that, if $f$ is singular at $x_0$, then it vanishes to some order $k>1$. We are going to show that there are no solutions which vanish to finite order $>1$. 
The problem is local, so we can assume that $M$ and $N$ are open subsets of $\mathbb{R}^d$ and we may take $x_0=0$ and $f(x_0)=0$. Furthermore, we may take bases such that $g^M(0)$ and $g^N(0)$ are the standard inner product. These choices still allow us the ability to independently make further orthogonal changes of basis on source and target. We expand $f$ in a Taylor series as $f(x) = f_k(x) + O(|x|^{k+1})$, expand the metrics as $g^M(x) = \mathrm{Id}+ O(|x|)$ and $g^N(y) = \mathrm{Id}+O(|y|)$, and expand the scalar $\lambda$ as $\lambda(x)=\lambda_m(x) +O(|x|^{m+1})$. Here $f_k$ and $\lambda_m$ are nonzero polynomials of degree $k$ and $m$.
The conformal condition is that
$$(Df_k + \cdots)^T (\mathrm{Id}+ \cdots) (Df_k + \cdots) = (\lambda_m + \cdots ) (\mathrm{Id}+ \cdots).$$
Since $f_k \neq 0$, the matrix $Df_k$ is nonzero and thus $(Df_k)^T (Df_k)$ is nonzero. So the leading terms $(Df_k)^T (Df_k)$ and $\lambda_m \mathrm{Id}$ must have matching degrees. We deduce that $m = 2k-2$ and $(Df_k)^T (Df_k) = \lambda_{m} \mathrm{Id}$. 
Thus, the polynomial map $f_k$ is conformal between $\mathbb{R}^d$ and itself, with the standard Euclidean structure. This contradicts the result of Iwaniec and Martin cited by the OP. 

Also, in the case $k=2$, we can verify that there is no solution by a direct computation. I will rename my variables -- my goal is to show that there is no conformal map $f: \mathbb{R}^d \to \mathbb{R}^d$ given by quadratic polynomials such that $(Df)(x)^T (Df)(x) = \lambda(x) \mathrm{Id}$ where $\lambda$ is a polynomial of degree $2$. Let the components of $f$ be $f(x) = \tfrac{1}{2} (x^T P_1 x, x^T P_2 x, \cdots, x^T P_n x)$ and let $\lambda(x) = \tfrac{1}{2} x^T Q x$, where the $P_i$ and $Q$ are symmetric. 
We can break the matrix $Df$ up into rows:
$$Df = \begin{bmatrix} x^T P_1 \\ x^T P_2 \\ \vdots \\ x^T P_n \end{bmatrix}.$$
Our given condition is that $(Df)(x)^T (Df)(x) = \lambda(x) \mathrm{Id}$; since $\lambda(x)$ is nonzero, we can work in the field of rational functions of $x$ and deduce that we also have $(Df)(x) (Df)(x)^T = \lambda(x) \mathrm{Id}$, giving
$$\begin{bmatrix} x^T P_1 \\ x^T P_2 \\ \vdots \\ x^T P_n \end{bmatrix}
\begin{bmatrix} P_1 x & P_2 x & \cdots & P_n x \end{bmatrix} =
\tfrac{1}{2} (x^T Q x) \mathrm{Id}.$$
So, for all $x \in \mathbb{R}^d$, we have
$$x^T P_i P_j x = \begin{cases} \tfrac{1}{2} x^T Q x & i = j \\ 0 & i \neq j \end{cases}.$$
We can't deduce that $P_i P_j = \tfrac{1}{2} Q$ from this, since $P_i P_j$ may not be symmetric. But we do deduce
$$P_i P_j + P_j P_i = \begin{cases} Q & i=j \\ 0 & i \neq j \end{cases}. \quad (\ast)$$
Making an orthogonal change of basis, we may assume that $Q$ is diagonal. Let the eigenvalues of $Q$ be $q_1$, $q_2$, ... $q_k$ with multiplicities $d_1$, $d_2$, ..., $d_k$. Since $Q = 2 P_i^2$, all the $P_i$ commute with $Q$ and thus they are block diagonal with blocks of size $d_1$, $d_2$, ..., $d_k$. Thus, we may solve the equation $(\ast)$ separately in each block, and may thus assume that $Q$ is diagonal; say $Q = q \mathrm{Id}$. Since $Q = 2 P_i^2$, we have $q>0$; rescaling each $P$ by $\tfrac{1}{\sqrt{2 q}}$, we may assume that $Q = \mathrm{Id}$. 
Thus, $(\ast)$ becomes 
$$P_i P_j + P_j P_i = \begin{cases} \mathrm{Id} & i=j \\ 0 & i \neq j \end{cases}. \quad (\ast)$$
In other words, $\mathbb{R}^{d_i}$ is a representation of the Clifford algebra on $d$ generators with $d_i \leq d$. For $d>2$, this is impossible.
(We are talking about the algebra labeled $C\ell_{d,0}$ in Wikipedia's table.) 
A: In a different direction, there are solutions where $f$ is $C^1$ and the metric $g_M$ is $C^0$. Let $\phi : \mathbb{R}_{>0} \to \mathbb{R}_{>0}$ be a $C^1$ function such that $\lim_{s \to 0^+} \phi(s) = \lim_{s \to 0^+} s \tfrac{\phi'(s)}{\phi(s)} = 0$. An example is $\phi(s) = \tfrac{1}{-\log s}$. (Or rather, that function near $0$ smoothed out somehow to not blow up at $s=1$.)
Write $s$ for the function $\sum x_i^2$ on $\mathbb{R}^n$. Let $M=N=\mathbb{R}^n$ and define $f:M \to N$ by $f(\vec{x}) = \phi(s)\ \vec{x}$. Let $g_N$ be the standard metric and define $g_M = \tfrac{1}{\phi(s)^2} f^{\ast} g_N$ on $M \setminus \{ 0 \}$. I claim that $g_M$ extends continuously to $0$. 
We compute: The Jacobian of $f$ at the point $\vec{x}$ is
$$J = \phi(s) \mathrm{Id} + \phi'(s)\  \vec{x} \ \vec{x}^T.$$
Using the identity $\vec{x}^T \vec{x}=s$, we compute:
$$g_M = \tfrac{1}{\phi(s)^2} J J^T = \mathrm{Id} + \left(2 \tfrac{\phi'(s)}{\phi(s)} + \tfrac{s \phi'(s)^2}{\phi(s)^2} \right) \vec{x} \ \vec{x}^T .$$
We will show that the second term zpproaches $0$ as $\vec{x} \to 0$. Note that each entry of $\vec{x}^T \vec{x}$ is $O(|\vec{x}|^2) = O(s)$. So each entry in the second term is bounded by 
$$2 \tfrac{s \phi'(s)}{\phi(s)} + \left( \tfrac{s \phi'(s)}{\phi(s)} \right)^2$$
and we assumed $\tfrac{s \phi'(s)}{\phi(s)} \to 0$.
Unfortunately, the hypothesis that $\tfrac{s \phi'(s)}{\phi(s)} \to 0$ gives $\tfrac{\phi'}{\phi} = o(s^{-1})$ so $\log |\phi| = o(\log s)$ and $\phi$ must approach $0$ more slowly than any $s^{\epsilon}$. So $f(x,0,0,\ldots,0)$ has derivative $0$ at $0$ but is larger than any $x^{1+\epsilon}$. So this method can't be pushed to $C^2$ functions.
