Let \(B\subset \!\subset M\lbrace a_1, \, \ldots , \, a_n\rbrace \) 
be a geodesic ball, of radius smaller than the
injectivity radius of \(M\) . By Lemma REF 
and Fubini's theorem, we can make sure that
\( \begin{split}\int _{B} \left(\left| D u_\varepsilon  \right|^2_g + \frac{1}{2\varepsilon ^2}(1 - \left| u_\varepsilon  \right|_g^2)^2\right)\mathrm {vol}_g+ \int _{\partial B} \left(\left| D u_\varepsilon  \right|^2_g + \frac{1}{2\varepsilon }(1 - \left| u_\varepsilon  \right|_g^2)^2\right) \mathrm {d}{H}^1 \le C\end{split}\) 
The constant \(C\)  does depend on the ball \(B\)  (in particular,
on the distance of \(B\)  from the singular points \(a_1\) , ..., \(a_n\) ),
but not on \(\varepsilon \) . Moreover, Lemma REF  implies
\(\left| u_\varepsilon  \right|_g \rightarrow 1 \qquad \textrm {uniformly in } \overline{B}.\) 
Step 1 We let \(\rho _\varepsilon :=\vert u_\varepsilon \vert _g\)  and \(v_\varepsilon :=\frac{u_\varepsilon }{\vert u_\varepsilon \vert _g}\) .
The field \(v_\varepsilon \)  is well-defined and belongs to \(H^1_{\operatorname{tan}}(B)\) ,
thanks to (REF ).
By lifting results (see e.g. [1]}), there exists a scalar
function \(\alpha _\varepsilon \in H^1(B)\)  such that
\(v_\varepsilon = \cos (\alpha _\varepsilon ) \hat{e}_1 + i\sin (\alpha _\varepsilon ) \hat{e}_1,\) 
where \(\left\lbrace \hat{e}_1,\hat{e}_2\right\rbrace =\left\lbrace \hat{e}_1,i\hat{e}_1\right\rbrace \) 
is an orthonormal moving frame in \(B\)  such that \(\mathrm {div}\,\mathbb {A}=0\) 
(see (REF )).
We rewrite (REF ) in terms of \(\rho _\varepsilon \)  and \(\alpha _\varepsilon \) .
First of all, note that
\(\Delta _g u_\varepsilon = \Delta _g (\rho _\varepsilon v_\varepsilon ) = \rho _\varepsilon \Delta _g v_\varepsilon + v_\varepsilon \Delta _g \rho _\varepsilon + 2 g^{ij}\partial _i \rho _\varepsilon D_j v_\varepsilon .\) 
Moreover (see [1]})
\(D v_\varepsilon = (\mathrm {d}\alpha _\varepsilon -\mathcal {A})iv_\varepsilon ,\) 
and (\(\mathbb {A}:=\mathcal {A}^{\sharp }\) )
\(\Delta _g v_\varepsilon = (\Delta \alpha _\varepsilon ) iv_\varepsilon -\vert \nabla \alpha -\mathbb {A}\vert ^2_g v_\varepsilon .\) 
Therefore,
\(\Delta _g u_\varepsilon = \Big (\rho _\varepsilon \Delta \alpha _\varepsilon + 2(\nabla \alpha _\varepsilon -\mathbb {A},\nabla \rho _\varepsilon )_g\Big )iv_\varepsilon + \Big (\Delta _g \rho _\varepsilon -\vert \nabla \alpha _\varepsilon -\mathbb {A}\vert ^2_g \Big )v_\varepsilon .\) 
We decompose \(f_\varepsilon = (f_\varepsilon , \, v_\varepsilon )_g\,v_\varepsilon + (f_\varepsilon , \, iv_\varepsilon )_gi v_\varepsilon =: f_\varepsilon ^1 \, v_\varepsilon + f_\varepsilon ^2 \, i v_\varepsilon \) . Since we have assumed that
\(f_\varepsilon \)  is bounded in \(L^2_{\operatorname{tan}}(M)\) , we have
\(\left\Vert  f^1_\varepsilon  \right\Vert _{L^2(B)} + \left\Vert  f^2_\varepsilon  \right\Vert _{L^2(B)} \le C.\) 
We can thus rewrite (REF ) as the following system
\({\left\lbrace \begin{array}{ll}\displaystyle \mathrm {div}\,(\rho _\varepsilon ^2(\nabla \alpha _\varepsilon -\mathbb {A})) = f^1_\varepsilon ,\\\\\displaystyle \Delta _g \rho _\varepsilon + \frac{\rho _\varepsilon }{\varepsilon ^2}(1-\rho _\varepsilon ^2) -\rho _\varepsilon \vert \nabla \alpha _\varepsilon \vert ^2_g = f^2_\varepsilon .\end{array}\right.}\) 

Step 2 Since \(\rho _\varepsilon \ge 1/2\)  and \(f^1_\varepsilon \in L^2(B)\) , it turns out that \(\alpha _\varepsilon \in H^2(B)\)  and thus \(h_\varepsilon :=f^2_\varepsilon +\rho _\varepsilon \vert \nabla \alpha _\varepsilon \vert ^2_g\in L^2\) .
Then, we multiply the equation for \(\rho _\varepsilon \)  with \(1-\rho _\varepsilon \)  and we integrate on \(B\) . We get
\(\int _{B}\vert \nabla \rho _\varepsilon \vert ^2_g \mathrm {vol}_g+ \frac{1}{\varepsilon ^2}\int _{B}\rho _\varepsilon (1-\rho _\varepsilon )^2(1+\rho _\varepsilon )\mathrm {vol}_g\nonumber \\= \int _{B}h_\varepsilon (1-\rho _\varepsilon ) \mathrm {vol}_g-\int _{\partial B}\partial _\nu \rho _\varepsilon (1-\rho _\varepsilon )\, \mathrm {d}{H}^1.\) 

Since \(\rho _\varepsilon \ge 1/2\)
\(\frac{1}{\varepsilon ^2}\int _{B}\rho _\varepsilon (1-\rho _\varepsilon )^2(1+\rho _\varepsilon )\mathrm {vol}_g\ge \frac{3}{4\varepsilon ^2}\int _{B}(1-\rho _\varepsilon )^2\mathrm {vol}_g.\) 
On the other hand, since \(h_\varepsilon \)  is in \(L^2\)  we get
\(\Big \vert \int _{B}h_\varepsilon (1-\rho _\varepsilon )\mathrm {vol}_g\Big \vert \le \frac{\delta }{2\varepsilon ^2}\int _{B}(1-\rho _\varepsilon )^2\mathrm {vol}_g+ \frac{\varepsilon ^2}{2\delta }\int _{B} h_\varepsilon ^2 \mathrm {vol}_g\) 
and (thanks to (REF ) and (REF ))
\(\int _{\partial B}\partial _\nu \rho _\varepsilon (1-\rho _\varepsilon )\, \mathrm {d}{H}^1 \le C\Vert \rho _\varepsilon -1\Vert _{L^\infty (\partial B)}.\) 
Choosing \(\delta =1/2\) , we conclude that
\(\int _{B}\vert \nabla \rho _\varepsilon \vert ^2_g \mathrm {vol}_g+ \frac{1}{4\varepsilon ^2}\int _{B}(1-\rho _\varepsilon )^2\mathrm {vol}_g\le C\big (\varepsilon ^2 \Vert h_\varepsilon \Vert _{L^2}^2 + \Vert \rho _\varepsilon -1\Vert _{L^\infty (\partial B)}\big ).\) 
Therefore,
\(\rho _\varepsilon \rightarrow 1 \qquad \textrm {strongly in } H^1(B)\) 
Moreover, since
\((1-x^2)^2 = (1-x)^2(1+x)^2\le 4(1-x)^2,\) 
for \(x\in (1/2,1)\) , we conclude that
\(\lim _{\varepsilon \rightarrow 0}\frac{1}{\varepsilon ^2}\int _{B}(1-\vert u_\varepsilon \vert _g^2)^2 \, \mathrm {vol}_g=0\) 

Step 3 The first equation of (REF ) can be
equivalently written as
\(\Delta _g\alpha _\varepsilon = (1 - \rho _\varepsilon ^2)\Delta _g\alpha _\varepsilon + \underbrace{(\nabla (1 - \rho _\varepsilon ^2), \, \nabla \alpha _\varepsilon )_g- \mathrm {div}\,(\rho _\varepsilon ^2\mathbb {A}) + f_\varepsilon ^2}_{=: g_\varepsilon }\) 
Both \(\rho _\varepsilon \)  and \(\alpha _\varepsilon \)  are bounded in \(H^1(B)\) ,
due to the energy estimate (REF ),
while \(f_\varepsilon ^1\)  is bounded in \(L^2(B)\)  by (REF ).
Therefore, \(g_\varepsilon \)  is bounded in \(L^2(B)\) .
Let \(B^\prime \subset \!\subset B\)  be a slightly smaller ball.
Thanks to (REF ), we can apply
elliptic regularity estimates and conclude that \(\alpha _\varepsilon \) 
is bounded in \(H^2(B^\prime )\) . By compact Sobolev embedding,
it follows that \(\alpha _\varepsilon \)  converges strongly in \(H^1(B^\prime )\) 
(possibly up to extraction of a subsequence).
Therefore, since
\(D_{\hat{e}_k}v_\varepsilon = (\mathrm {d}\alpha _\varepsilon (\hat{e}_k) -\mathcal {A}(\hat{e}_k))i v_\varepsilon ,\) 
for \(k=1, \, 2\) , we get that \(D_{\hat{e}_k}v_\varepsilon \)  is strongly convergent in \(L^2(B^\prime )\) .
As
\(D_{\hat{e}_k} u_\varepsilon =\rho _\varepsilon D_{\hat{e}_k}v_\varepsilon + v_\varepsilon (\nabla \rho _\varepsilon , \, \hat{e}_k)_g,\) 
thanks to (REF ) we conclude that
\(D u_\varepsilon \)  converges strongly in \(L^2(B^\prime )\) .

