We consider the expression for \(W\)  given by Proposition REF .
We computed already, in Lemma REF , the gradient
of the terms that do not depend on \(\xi \) ; now, we
focus our attention on the term that depends on \(\xi \) .
We take a tangent vector \(v\in T_{a_j} M\) .
We consider the smooth map \(\mathbf {b}\mapsto \Xi (\mathbf {b})\) ,
defined locally in a neighbourhood of \(\mathbf {a}\)  in \(M^n\) ,
such that \(\Xi (\mathbf {a}) = \xi \)  and \(\Xi (\mathbf {b})\in \mathcal {L}(\mathbf {b}, \, \mathbf {d})\) 
for any \(\mathbf {b}\) , as given by Lemma REF .
It suffices to prove that
\( \frac{1}{2}\left(\nabla _{a_j} \int _M \left| \Xi (\mathbf {a}) \right|^2 \, \mathrm {vol}_g, \, v \right)_g= 2\pi d_j \left(i\xi ^{\#}(a_j), \, v\right)_g\) 
Combining (REF ) with Lemma REF ,
the result will follow.
Step 1 Let \(\gamma _1\) , ..., \(\gamma _{2\mathfrak {g}}\)  is a family of
simple, closed, smooth curves whose homology classes generate
the homology group \(H_1(M; \, \mathbb {Z})\) .
As before, we assume (up to a perturbation of the \(\gamma _k\) 's)
that the curves \(\gamma _k\)  do not contain any of the points \(a_j\) .
Let \(\Xi ^\prime := (\nabla _{a_j} \Xi (\mathbf {a}), \, v)_g\)  be the derivative
of \(\Xi \)  with respect to the \(a_j\) -variable, in the direction of \(v\) .
We have seen in Remark REF , Equation (REF ),
that \(\Xi ^\prime \)  is the unique harmonic 1-form such that
\( \int _{\gamma _k} \Xi ^\prime = 2\pi d_j \int _{\gamma _k} \star \mathrm {d}\sigma (\cdot , \, a_j, \, v)\) 
for any \(k\in \lbrace 1, \, \ldots , \, 2\mathfrak {g}\rbrace \) ,
where \(\sigma (\cdot , \, a_j, \, v)\)  is given by (REF ).
From now on, we write \(\sigma \)  instead of \(\sigma (\cdot , \, a, \, v)\) .
Since \(\sigma \)  is smooth in \(M\setminus \lbrace a_j\rbrace \) , the 1-form
\( \beta := \frac{\Xi ^\prime }{2\pi d_j} - \star \mathrm {d}\sigma \) 
is smooth in \(M\setminus \lbrace a_j\rbrace \)  and satisfies
\( \int _{\gamma _k} \beta = 0\qquad \textrm {for any } k\in \lbrace 1, \, \ldots , \, 2\mathfrak {g}\rbrace \) 
due to (REF ). By differentiating (REF ),
and keeping in mind that a harmonic form is closed, we obtain
\(\mathrm {d}\beta = -\mathrm {d}\star \mathrm {d}\sigma = \star \mathrm {d}^*\mathrm {d}\sigma = \star \Delta \sigma = 0\qquad \textrm {in } M\setminus \lbrace a_j\rbrace ,\) 
i.e. \(\beta \)  is closed in \(M\setminus \lbrace a_j\rbrace \) .
Moreover, the homology classes of the curves \(\gamma _1\) ,   ...,   \(\gamma _{2\mathfrak {g}}\) 
are also generators of \(H_1(M\setminus \lbrace a_j\rbrace ; \, \mathbb {Z})\) ,
because the inclusion \(M\setminus \lbrace a_j\rbrace \hookrightarrow M\) 
induces a group isomorphism \(H_1(M\setminus \lbrace a_j\rbrace ; \, \mathbb {Z})\rightarrow H_1(M; \, \mathbb {Z})\) .
This may be seen by writing the surface \(M\)  as a quotient of a polygon
(more precisely, a \((4\mathfrak {g})\) -gon) with suitable identifications
on the boundary; see e.g. [1]}
and Figure REF .
\(M\setminus \lbrace a_j\rbrace \)  retracts by deformation
onto the 1-skeleton of \(M\)  and \(M\)  is obtained by attaching
a 2-cell along a trivial cycle of the 1-skeleton,
so \(H_1(M\setminus \lbrace a_j\rbrace ; \, \mathbb {Z})\simeq H_1(M; \, \mathbb {Z})\) .An
alternative argument: there is an exact sequence
\(H_2(M; \, \mathbb {Z}) \rightarrow H_2(M, \, M\setminus \lbrace a_j\rbrace ; \, \mathbb {Z})\rightarrow H_1(M\setminus \lbrace a_j\rbrace ; \, \mathbb {Z}) \xrightarrow{} H_1(M; \, \mathbb {Z})\rightarrow 0\) 
of relative homology groups with integer coefficients
(see e.g. [2]}). Here \(\chi _*\) 
is the group homomorphism induced by the
inclusion \(\chi \colon M\setminus \lbrace a_j\rbrace \hookrightarrow M\) .
As \(M\)  is closed and orientable, the first arrow is an isomorphism
(see e.g. [2]}); it follows that \(\chi _*\)  is an isomorphism, too.

Then, thanks to (REF ), \(\beta \)  is not only
closed but also exact in \(M\setminus \lbrace a_j\rbrace \) . As a consequence,
there exists a function \(f\in C^\infty (M\setminus \lbrace a_j\rbrace )\)  such that
\( \frac{\Xi ^\prime }{2\pi d_j} = \star \mathrm {d}\sigma + \mathrm {d}f\qquad \textrm {in } M\setminus \lbrace a_j\rbrace \) 

Step 2 Let
\(I := \frac{1}{4\pi d_j}\left(\nabla _{a_j} \int _M \left| \Xi (\mathbf {a}) \right|^2 \, \mathrm {vol}_g, \, v \right)_g\) 
Let \(\eta >0\)  be a small parameter and let \(B_\eta (a_j)\) 
be the geodesic disk with radius \(\eta \)  and center \(a_j\) .
We may evaluate \(I\)  as
\(\begin{split}I = \frac{1}{2\pi d_j} \int _M (\xi , \, \Xi ^\prime )_g \, \mathrm {vol}_g&= \frac{1}{2\pi d_j} \lim _{\eta \rightarrow 0}\int _{M\setminus B_\eta (a_j)}(\xi , \, \Xi ^\prime )_g \, \mathrm {vol}_g\\&\hspace{-5.69046pt} \stackrel{(\ref {gradrenSHg5})}{=}\lim _{\eta \rightarrow 0}\int _{M\setminus B_\eta (a_j)}(\xi , \, \star \mathrm {d}\sigma + \mathrm {d}f)_g \, \mathrm {vol}_g\\&= \lim _{\eta \rightarrow 0}\int _{M\setminus B_\eta (a_j)}(\xi , \, -\mathrm {d}^*(\star \sigma ) + \mathrm {d}f)_g \, \mathrm {vol}_g\end{split}\) 
We integrate by parts in the right-hand side,
by applying (REF ):
\( I = \lim _{\eta \rightarrow 0}\int _{\partial B_\eta (a_j)}\left(-\sigma \xi - f(\star \xi )\right)\) 
The volume integrals are equal to zero because \(\xi \) 
is harmonic, hence \(\mathrm {d}\xi =0\) , \(\mathrm {d}^*\xi =0\) .
The signs account for orientations effects,
as \(\partial B_\eta (a_j)\)  is given the orientation induced by \(B_\eta (a_j)\) ,
not by \(M\setminus B_\eta (a_j)\) . Let \(\nu \)  be the outward-pointing
unit normal to \(\partial B_\eta (a_j)\) . Then, \(i\nu \)  is a unit-norm,
tangent vector field on \(\partial B_\eta (a_j)\) , consistent with the orientation
induced by \(B_\eta (a_j)\) . By definition, we have
\( \int _{\partial B_\eta (a_j)} \left(-\sigma \xi - f(\star \xi )\right)= \int _{\partial B_\eta (a_j)} \left(-\sigma (\xi ^\#, \, i\nu )_g- f((\star \xi )^\#, \, i\nu )_g\right) \mathrm {d}{H}^1\) 
and, by Lemma REF , \((\star \xi )^\# = i\xi ^\#\) .
Therefore, from (REF ) and (REF ) we obtain
\(\begin{split}I &= \lim _{\eta \rightarrow 0} \left(-\int _{\partial B_\eta (a_j)}\sigma \, (\xi ^\#, \, i\nu )_g \, \mathrm {d}{H}^1- \int _{\partial B_\eta (a_j)}f (i\xi ^\#, \, i\nu )_g \, \mathrm {d}{H}^1 \right)\end{split}\) 
or, since \(i\)  is a rotation by \(\pi /2\)  in each tangent plane,
\( \begin{split}I &= \lim _{\eta \rightarrow 0} \left(\int _{\partial B_\eta (a_j)}\sigma \, (i\xi ^\#, \, \nu )_g \, \mathrm {d}{H}^1- \int _{\partial B_\eta (a_j)}f (i\xi ^\#, \, i\nu )_g \, \mathrm {d}{H}^1 \right)\end{split}\) 

<FIGURE>Step 3 In order to proceed further, we need to estimate \(\sigma \) 
and \(f\)  on \(\partial B_\eta (a_j)\) .
We work in normal geodesic coordinates centered at \(a_j\) .
We identify points in \(B_\eta (a_j)\)  with points in \(\mathbb {R}^2\) 
and vector fields with maps \(B_\eta (a_j)\subseteq \mathbb {R}^2\rightarrow \mathbb {R}^2\) 
(in particular, we identify \(v\)  with a constant \(v\in \mathbb {R}^2\) ).
An estimate for \(\sigma \)  is given by Lemma REF :
\( \sigma = \frac{1}{2\pi \eta } (\nu , \, v)_g + \mathrm {O}(1)\qquad \textrm {on } \partial B_\eta (a_j)\) 
As for \(f\) , we combine (REF ) with ()
to obtain
\( \mathrm {d}f = -\star \mathrm {d}\sigma + \mathrm {O}(1)= \frac{1}{2\pi \eta ^2}\left((i\nu , \, v)_g\nu ^\flat + (\nu , \, v)_g\star \nu ^\flat \right)+ \mathrm {O}(\eta ^{-1})\qquad \textrm {on } \partial B_\eta (a_j)\) 
Let \(x_\eta \in \partial B_\eta (a_j)\)  be a reference point, say
\(x_\eta = a_j + (\eta , \, 0)\) . For any \(x\in \partial B_\eta (a_j)\) ,
let \(\gamma _x\)  be the circular arc in \(\partial B_\eta (a_j)\) 
with endpoints \(x_\eta \)  and \(x\) , oriented in the anticlockwise
direction. Then
\(f(x) = f(x_\eta ) + \int _{\gamma _x} \mathrm {d}f\stackrel{(\ref {gradrenSHg9})}{=} f(x_\eta )+ \frac{1}{2\pi \eta ^2}\int _{\gamma _x}(\nu , \, v)_g\star \nu ^\flat + \mathrm {O}(1)\) 
By Gauss' lemma, when we work in geodesic coordinates
we have \(\nu (x) = (x-a_j)/\left| x-a_j \right|\)  and moreover, the metric
tensor satisfies \(g_{ij}(x) = \delta _{ij} + \mathrm {O}(\eta ^2)\) .
Therefore, writing \(x = a_j + \eta (\cos \theta _x, \, \sin \theta _x)\) 
for some \(\theta _x\in (0, \, 2\pi )\) , we may evaluate the integral as
\( \begin{split}f(x) &= f(x_\eta ) + \frac{1}{2\pi \eta } \int _0^{\theta _x}\left(v^1\cos \theta + v^2\sin \theta \right) \mathrm {d}\theta + \mathrm {O}(1) \\&= \frac{1}{2\pi \eta }\left(v^1\sin \theta _x - v^2\cos \theta _x\right)+ C_\eta + \mathrm {O}(1)= -\frac{1}{2\pi \eta }(i\nu (x), \, v)_g + C_\eta + \mathrm {O}(1)\end{split}\) 
for any \(x\in \partial B_\eta (a_j)\) , where \(C_\eta \)  depends on \(\eta \) 
but not on \(x\) .

Step 4 We combine (REF ), (REF )
and (REF ):
\( \begin{split}I = \lim _{\eta \rightarrow 0} \bigg (\frac{1}{2\pi \eta }\int _{\partial B_\eta (a_j)} (i\xi ^\#, \, \nu )_g \, (\nu , \, v)_g \, \mathrm {d}{H}^1&+ \frac{1}{2\pi \eta } \int _{\partial B_\eta (a_j)}(i\xi ^\#, \, i\nu )_g \, (i\nu , \, v)_g \, \mathrm {d}{H}^1 \\+ \frac{C_\eta }{2\pi \eta }\int _{\partial B_\eta (a_j)} (i\xi ^\#, \, i\nu )_g \, \mathrm {d}{H}^1&+ \mathrm {O}(\eta ) \bigg )\end{split}\) 
where \(C_\eta \)  is an (unknown) function of \(\eta \) . However,
Stokes' theorem implies
\(\begin{split}\int _{\partial B_\eta (a_j)} (i\xi ^\#, \, i\nu )_g \, \mathrm {d}{H}^1= \int _{\partial B_\eta (a_j)} \star \xi = \int _{B_\eta (a_j)} \mathrm {d}(\star \xi ) = 0\end{split}\) 
because \(\xi \)  is harmonic.
Therefore, (REF ) simplifies to
\(\begin{split}I &= \lim _{\eta \rightarrow 0} \frac{1}{2\pi \eta }\int _{\partial B_\eta (a_j)}\left(i\xi ^\#, \, (\nu , \, v)_g \, \nu +(i\nu , \, v)_g \, i\nu \right)_g\mathrm {d}{H}^1= \lim _{\eta \rightarrow 0} \frac{1}{2\pi \eta }\int _{\partial B_\eta (a_j)} (i\xi ^\#, \, v)_g \,\mathrm {d}{H}^1\end{split}\) 
As \({H}^1(\partial B_\eta (a_j)) = 2\pi \eta + \mathrm {O}(\eta )\) 
(see e.g. [4]}),
the equality (REF ) follows.

