Let \(e\)  be a smooth, divergence-free vector field, defined in
a neighbourhood of \(a_j\) . Let
\(I_\eta (u_*)(a_j):= \int _{\partial B_\eta (a_j)}(D_{\nu } u_*, \, D_{e}u_*)_g \, \mathrm {d}{H}^1- \frac{1}{2} \int _{\partial B_\eta (a_j)}\left| D u_* \right|^2_g (\nu , \, e)_g\, \mathrm {d}{H}^1\) 
The definition (REF ) of \(j(u_*)\) 
implies, via standard computations, that
\(\left| D u_* \right|_g^2= \left| j(u_*) \right|_g^2\) , \(j(u_*) = \mathrm {d}^*\Psi \)  and
\(D_v u_* = j(u_*)(v) \, iu_*\)  for any smooth vector field \(v\) .
Moreover, \(j(u_*) = \mathrm {d}\Psi (\mathbf {a}, \, \mathbf {d})\)  by
Definition REF .
For simplicity, we write \(\Psi \)  instead of \(\Psi (\mathbf {a}, \, \mathbf {d})\)  from now on.
Then, we have
\( \begin{split}I_\eta (u_*)(a_j)= \int _{\partial B_\eta (a_j)}\mathrm {d}^*\Psi (\nu ) \, \mathrm {d}^*\Psi (e) \, \mathrm {d}{H}^1- \frac{1}{2}\int _{\partial B_\eta (a_j)}\left| \mathrm {d}^*\Psi  \right|^2_g (\nu , \, e)_g \, \mathrm {d}{H}^1\end{split}\) 
Now,
\(\mathrm {d}^*\Psi = -\star \mathrm {d}(\star \Psi )\stackrel{(\ref {S_j})}{=} -\star \mathrm {d}(S_j + 2\pi d_j G(\cdot , a_j))\) 
and thus
\( \begin{split}\vert \mathrm {d}^*\Psi \vert ^2_g& = \vert \mathrm {d}(S_j + 2\pi d_j G(\cdot ,a_j))\vert ^2_g= \vert \nabla (S_j + 2\pi d_j G(\cdot ,a_j))\vert ^2_g \\& = \vert \nabla S_j\vert ^2_g+ 4\pi ^2 d_j^2 \, \vert \nabla G(\cdot ,a_j)\vert ^2_g+ 4\pi d_j \, (\nabla S_j, \, \nabla G(\cdot ,a_j))_g\end{split}\) 
On the other hand, for any vector field \(X\)  we have
\(\begin{split}\mathrm {d}^*\Psi (X) = -\star \mathrm {d}(S_j + 2\pi d_j G(\cdot , a_j))\,(X)&\stackrel{(\ref {eq:grad_ruotato})}{=}-(i\nabla (S_j + 2\pi d_j G(\cdot , a_j)), \, X) \\&\hspace{7.11317pt}= (\nabla (S_j + 2\pi d_j G(\cdot , a_j)), \, iX)\end{split}\) 
because in the tangent plane at each point, \(i\)  is a rotation
and its adjoint is equal to \(i^{-1} = -i\) . As a consequence,
\(\begin{split}\mathrm {d}^*\Psi (\nu ) \, \mathrm {d}^*\Psi (e)&= \nabla _{i\nu }(S_j + 2\pi d_j G(x,a_j))\, \nabla _{i e}(S_j + 2\pi d_j G(x,a_j)) \\&= \nabla _{i\nu }S_j(\cdot ) \, \nabla _{ie}S_j(\cdot )+ 2\pi d_j \, \nabla _{i\nu }S_j(\cdot ) \, \nabla _{ie} G(\cdot ,a_j) \\&\qquad + 2\pi d_j \, \nabla _{i\nu } G(\cdot ,a_j) \, \nabla _{i e}S_j(\cdot )+ 4\pi ^2 d_j^2 \, \nabla _{i\nu } G(\cdot ,a_j) \, \nabla _{ie} G(\cdot ,a_j)\end{split}\) 
We write the Green function \(G\)  in the form (REF ).
For ease of notation, we let \(L_j := \log \operatorname{dist}_g(\cdot , \, a_j)\) 
and \(H_j := H(\cdot , \, a_j)\) , so that
\(2\pi d_j G(\cdot , \, a_j) = -d_j L_j + 2\pi d_j H_j\) 
The equations (REF ) and (REF )
imply respectively
\( \begin{split}\vert \mathrm {d}^*\Psi \vert ^2_g&= \left| \nabla S_j \right|^2_g + d_j^2 \left| \nabla L_j \right|^2_g+ 4\pi ^2 d_j^2 \left| \nabla H_j \right|^2_g \\&\qquad - 4\pi d_j^2 \, (\nabla L_j, \, \nabla H_j)_g- 2 d_j \, (\nabla S_j, \, \nabla L_j)_g+ 4\pi d_j \, (\nabla S_j, \, \nabla H_j)_g\end{split}\) 
and
\( \begin{split}\mathrm {d}^*\Psi (\nu ) \, \mathrm {d}^*\Psi (e)& = \nabla _{i\nu }S_j \, \nabla _{ie}S_j- d_j \, \nabla _{i\nu }S_j \, \nabla _{ie} L_j+ 2\pi d_j \, \nabla _{i\nu } S_j \, \nabla _{ie} H_j \\&\qquad - d_j \, \nabla _{i\nu } L_j \, \nabla _{ie} S_j+2\pi d_j \, \nabla _{i\nu } H_j \, \nabla _{ie} S_j+ d_j^2 \, \nabla _{i\nu } L_j \, \nabla _{ie} L_j \\&\qquad -2\pi d_j^2 \, \nabla _{i\nu } L_j \, \nabla _{ie} H_j-2\pi d_j^2 \, \nabla _{i\nu } H_j \, \nabla _{ie} L_j+ 4\pi ^2 d_j^2 \, \nabla _{i\nu } H_j \, \nabla _{ie} H_j\end{split}\) 
The function \(H\)  is of class \(C^1\)  in \(M\times M\) 
and, by Equation (REF ), \(S_j\)  is smooth
in a neighbourhood of \(a_j\) . In particular, both \(\nabla H_j\) 
and \(\nabla S_j\)  are locally bounded around \(a_j\) .
Therefore, we can write Equations (REF )
and (REF ) as
\( \left| \mathrm {d}^*\Psi  \right|^2_g = d_j^2 \left| \nabla L_j \right|^2_g- 2d_j(\nabla (S_j+2\pi d_j H_j), \, \nabla L_j)_g+ \mathrm {O}(1)\) 
and
\( \begin{split}\mathrm {d}^*\Psi (\nu ) \, \mathrm {d}^*\Psi (e)= d_j^2 \, \nabla _{i\nu } L_j \, \nabla _{ie} L_j&- d_j \, \nabla _{i\nu }(S_j + 2\pi d_j H_j) \, \nabla _{ie} L_j \\&- d_j \, \nabla _{i\nu } L_j \, \nabla _{ie} (S_j + 2\pi d_j H_j)+ \mathrm {O}(1)\end{split}\) 
respectively. Gauss' lemma implies that
\( \nabla \left(\log \operatorname{dist}_g(\cdot , \, a_j)\right)= \frac{\nabla \operatorname{dist}_g(\cdot , \, a_j)}{\operatorname{dist}_g(\cdot , \, a_j)}= \frac{\nu }{\operatorname{dist}_g(\cdot , \, a_j)} \qquad \textrm {in } B_\eta (a_j)\setminus \lbrace a_j\rbrace \) 
so on the circle \(\partial B_\eta (a_j)\)  we have
\(\nabla L_j = \nu /\eta \) , \(\nabla _{i\nu } L_j = 0\)  and
\(\left| \mathrm {d}^*\Psi  \right|^2 &= \frac{d_j^2}{\eta ^2}- \frac{2d_j}{\eta } \nabla _\nu (S_j+2\pi d_j H_j)+ \mathrm {O}(1)  \\\mathrm {d}^*\Psi (\nu ) \, \mathrm {d}^*\Psi (e)&= - \frac{d_j}{\eta } \nabla _{i\nu }(S_j + 2\pi d_j H_j) (\nu , \, ie)+ \mathrm {O}(1) \) 

By substituting (REF ), ()
in (REF ), we obtain
\( \begin{split}I_\eta (u_*)(a_j)= &-\frac{d_j}{\eta } \int _{\partial B_\eta (a_j)}\nabla _{i\nu } (S_j+2\pi d_j H_j) \, (\nu , \, ie)_g \, \mathrm {d}{H}^1 \\&+ \frac{d_j}{\eta } \int _{\partial B_\eta (a_j)}\nabla _\nu (S_j+2\pi d_j H_j) \, (\nu , \, e)_g \, \mathrm {d}{H}^1 \\&\qquad \qquad \qquad - \frac{d_j^2}{2\eta ^2}\int _{\partial B_\eta (a_j)} (\nu , \, e)_g \, \mathrm {d}{H}^1+ \mathrm {O}(\eta )\end{split}\) 
The properties of \(i\)  imply \((\nu , \, ie)_g = -(i\nu , \, e)_g\)  and thus
\( \begin{split}&- \nabla _{i\nu }(S_j +2\pi d_j H_j) \, (\nu , \, ie)_g+ \nabla _\nu (S_j + 2\pi d_j H_j) \, (\nu , \, e)_g \\&\qquad = \big (\nabla _{\nu }(S_j + 2\pi d_j H_j)\nu , \, e\big )_g+ \big (\nabla _{i\nu }(S_j + 2\pi d_j H_j) i\nu , \, e\big ) _g= (\nabla (S_j + 2\pi d_j H_j), \, e)_g.\end{split}\) 
Moreover, by integrating by parts (see (REF )),
we obtain
\( \begin{split}\int _{\partial B_\eta (a_j)} (\nu , \, e)_g \, \mathrm {d}{H}^1= \int _{B_\eta (a_j)} (\mathrm {div}\,e) \, \mathrm {vol}_g= 0\end{split}\) 
because we have assumed that \(\mathrm {div}\,e = 0\) .
Combining (REF ), (REF )
and (REF ), we deduce
\( \begin{split}I_\eta (u_*)(a_j)&= \frac{d_j}{\eta } \int _{\partial B_\eta (a_j)}\nabla _{e} (S_j+2\pi d_j H_j) \, \mathrm {d}{H}^1 + \mathrm {O}(\eta )\end{split}\) 
and, since \(\mathcal {H}^1(\partial B_\eta (a_j))= 2\pi \eta + \mathrm {O}(\eta ^2)\) 
(see e.g. [1]}), we conclude that
\(\lim _{\eta \rightarrow 0} I_{\eta }(u_*)(a_j)= 2\pi d_j\left(\nabla S_j(a_j)+ 2\pi d_j \nabla H(a_j, \, a_j), \, e(a_j)\right)_g\) 
The result follows by Lemma REF .
