I'm too lazy to type-up the proof myself, so I'll send you to a reference.
Chang, S.-Y. A., Wang, L. and Yang, P. C. (1999), "Regularity of harmonic maps". CPAM has the proof in Section 3. Once you get $C^{1,\gamma}$ you immediately get RHS is in $C^\gamma$ and the rest follow by standard elliptic regularity.
Note that the structure of the equation (RHS being of the form $d(u\cdot du)$) is only used for Wente's lemma. For the upgrade of regularity one uses a Caccioppoli type inequality.
(BTW, the Chang-Wang-Yang result bypasses the Hardy space estimates. For that the result can be found in the original paper of Helein, though I'd guess the material is also in his book if you don't read French.)