Approximation of $C^1$-smooth equivariant maps by infinitely smooth ones Let $M,N$ be smooth closed manifolds acted by a finite group $G$. Let $f\colon M\to N$ be a $C^1$-smooth $G$-equivariant map.
Is it true that for any $\varepsilon>0$ there exists a $C^\infty$-smooth $G$-equivariant map
$$f_\varepsilon\colon M\to N$$ such that $\|f-f_{\varepsilon}\|_{C^1}<\varepsilon$, where the $C^1$-norm is taken with respect to some  Riemannian metrics on $M$ and $N$?
A reference would be most helpful.
 A: One option is to use the harmonic map flow developed by Eels and Sampson [1]. In a certain sense this is a (non-linear) analog of the heat equation for maps $M \to N$.
Endow the manifolds $M$ and $N$ with two smooth Riemannian metrics $g$ and $h$, which additionally we may assume invariant under the action of $G$. Then you can evolve $u_0 := f$ along a family of mappings $(u(t,\cdot) \mid 0 \leq t < \epsilon)$ which solve the so-called harmonic map flow for some short time.
Specifically
\begin{equation}
\begin{cases}
\partial_t u &=& \tau_g(u) \\
u(0,\cdot) &=& f
\end{cases}
\end{equation}
where $\tau_g(u) = \mathrm{tr} \nabla \mathrm{d} u$ is known as the tension field.
This is the gradient flow for the energy functional \begin{equation}E(u) = \int_M \lvert \mathrm{d} u \rvert^2 \, \mathrm{d vol}_g . \end{equation} The (short-time) existence, uniqueness and smoothness of the solution is guaranteed by virtue of the PDE being strictly parabolic.
In particular, for all $t \in (0,\epsilon)$, $u(t,\cdot): M \to N$ is smooth, and as $t \to 0$,
\begin{equation}
\lvert u(t,\cdot) - f \rvert_{C^1} \to 0.
\end{equation}
The equivariance of $u$ under $G$—given that of the initial datum $f$—follows from the uniqueness of the solution. This is because the PDE is also invariant under the action of the group.
For an explicit confirmation, let $\gamma \in G$ be arbitrary and define $u_1: x \mapsto \gamma u_0(\gamma^{-1} x)$. Then on the one hand the family $(\gamma u(t,\gamma^{-1} \cdot) \mid 0 \leq t < \epsilon)$ is the harmonic map flow from $u_1$. On the other hand $u_1 = u_0$ by assumption, and therefore
\begin{equation}
\gamma u(t,\gamma^{-1} \cdot) = u(t,\cdot)  \text{ for all $0 \leq t < \epsilon$.}
\end{equation}
[1] J. Eels and J.H. Sampson. Harmonic mappings of Riemannian manifolds. Amer. J. Math. 86 (1964) 109-160.
