$\newcommand{\lap}{\Delta}$ $\newcommand{\al}{\alpha}$ $\newcommand{\be}{\beta}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\pl}{\partial}$
Let $(M,g),(N,\eta)$ be smooth Riemannian manifolds. Then the equation for harmonic maps $M \to N$, expressed in coordinates, is $$ \lap_{g} f^{\al}+g^{ij}\pl_if^{\be} \pl_j f^{\ga}(\Ga_{\ga \be}^{\al} \circ f)=0, \tag{1}$$
where $i,j$ are indices for coordinates on $M$, $\alpha,\be,\ga$ are indices on $N$. In particular, the Christoffel symbols $\Ga_{\ga \be}^{\al} $ are those of the Levi-Civita connection of $N$.
Question: I understand that this equation is regarded as "semilinear elliptic of second order". Are there any regularity results for this equation?
A google search with the keywords "semilinear elliptic" + "regularity" found only results for equations of the form of $$ \lap u=F(u,x),$$ where $F$ only depends on $x,u(x)$. The harmonic equation is something like $ \lap u=F(du,u,x).$
I am curious, since I know that regularity theory for harmonic maps is done in a different way: Instead of considering the coordinate equation, you isometrically embed the target in Euclidean space, and use the second fundamental form. I guess the reason not to work with the coordinate version $(1)$ is that this does not give the strong regularity results obtained by the extrinsic approach.
Is it known that such analogous regularity results do not hold for the coordinate version?
(e.g any continuous weakly harmonic map-in the extrinsic sense- is $C^{\infty}$. Does something like that holds for the coordinate equation?)
Edit:
Let me clarify something. The strong versions of intrinsic and extrinsic harmonicity are the same. However, I am not sure that their "weak versions" would be the same; that is, suppose we had a way to define a weak solution of equation $(1)$, which we call "intrinsic weakly harmonic". Are the two notions of weak harmonicity (the standard extrinsic, and the new intrinsic) conicide?
Of course, if every continuous weak solution of $(1)$ is smooth, then both notions of weak harmonicity coincide (at least for continuous maps, since they both turn out to be smooth, hence strongly hamronic).
But it's possible that not every weak solution of the "coordinate version" is smooth, and that the weak notions differ. So, I am asking if we can use standard PDE theory directly on the coordinate version, but I think that such an attempt might lead to a different notion of weak harmonicity. At this moment I am not even sure if there is an accepted way to define a weak version of the coordinate equation.
Is there such a way?