I posted this question in MSE but got no response (even after giving a bounty), so I am trying here.
Let $M,N$ be smooth $d$-dimensional Riemannian manifolds.
Suppose $f:M \to N$ is a differentiable isometry ($df_p$ is an isometry at every $p \in M$). I do not assume $f$ is $C^1$.
Is it true that $f$ must be smooth?
This is false when assuming $f$ is only differentiable almost everywhere.
One counter-example is $f(x)=c(x)+x$ where $c$ is the Cantor function. In this case $M=[0,1],N=[0,2]$, $f'=1$ a.e; $M,N$ are not isometric, but both are flat.
Gromov showed similar examples exist for "many" non-flat $M$ and $N=\mathbb{R}^d$.
Some partial results:
By the inverse function theorem for everywhere differentiable maps, $f$ is a local homeomorphism.
We would like to prove it's a local isometry w.r.t the intrinsic distances, then use the Myers-Steenrod theorem (or use the fact geodesics locally minimize length to deduce $f$ maps geodesics to geodesics, hence it factors through its differential via the exponential maps. this works if $f \in C^1$, see here for details).
The devil is in the details: We need to choose suitable length structures on $M,N$ such that $f$ will become an arcwise isometry. We cannot use the class of $C^1$ paths since $f$ does not necessarily map $C^1$ maps to $C^1$ maps.
Also, nothing promises us that for a differentiable path $\gamma$ such that $\|\dot \gamma(t)\|$ is integrable, $\|\dot {f \circ \gamma}(t)\|$ will be integrable.