Convergence of discrete period matrices and discrete holomorphic integrals for ramified coverings of the Riemann sphere

We consider the class of compact Riemann surfaces which are ramified coverings of the Riemann sphere $\hat{\mathbb{C}}$. Based on a triangulation of this covering of the sphere $\mathbb{S}^2\cong \hat{\mathbb{C}}$ and its stereographic projection, we define discrete (multi-valued) harmonic and holomorphic functions. We prove that the corresponding discrete period matrices converge to their continuous counterparts. In order to achieve an error estimate, which is linear in the maximal edge length of the triangles, we suitably adapt the triangulations in a neighborhood of every branch point. Finally, we also prove a convergence result for discrete holomorphic integrals for our adapted triangulations of the ramified covering.


Introduction
Smooth holomorphic functions can be characterized in different ways. In particular, the real and imaginary part of any holomorphic function is harmonic and both are related by the Cauchy-Riemann equations. This perspective naturally led to linear discretizations of harmonic and holomorphic functions, starting with results for square grids, see [CFL28,Isa41,Duf53]. Lelong-Ferand further developed this theory of discrete harmonic and holomorphic functions in [Fer44,LF55]. MacNeal and Duffin generalized these notions in [Mac49,Duf56,Duf59,Duf68]. In particular, they considered arbitrary triangulations in the plane and discovered the cotan-weights. The cotan-Laplacian is also considered for triangle meshes, for example for surfaces in discrete differential geometry, see [PP93], or for applications in computer graphics, see for example [MDSB03]. Further properties and theorems of the smooth theory of holomorphic functions have found recently discrete analogues in the discrete linear theory in [BG16,BG17].
Mercat generalized in [Mer01] the discrete linear theory from planar subsets to discrete Riemann surfaces and introduced in [Mer02,Mer07] discrete period matrices. In [BMS11] numerical experiments are considered to compute discrete period matrices for polyhedral surfaces explicitly and compare them to known period matrices for the corresponding smooth surfaces. A convergence proof for the class of polyhedral surfaces was obtained in [BS16].
The interest in numerical computation of period matrices is for example motivated by the computation of finite-genus solutions of integrable differential equations. As Riemann surfaces may be represented as algebraic curves, this is often taken as a starting point for computing discrete period matrices. Recent results in this context include [GSST98,DvH01,FK15,FK17,MN17].
In this article, we take a different approach and consider Riemann surfaces which are ramified coverings of the Riemann sphereĈ. Based on a triangulation of this covering of the sphere S 2 ∼ =Ĉ or its stereographic projection, discrete period matrices can be obtained from this discrete data. Furthermore, we prove convergence of the discrete period matrices to their continuous counterparts (Theorem 2.5). In particular, we obtain an error estimate, which is linear in the maximal edge length of the triangles if we adapt the triangulations in a neighborhood of every branch point. The details of our 'adapted triangulations' will be explained in Section 2.3.
The convergence of discrete analytic functions to their continuous counterparts remains an important issue, although several results have been proved by now. In particular, for the linear theory, convergence was first shown for the square lattice [CFL28,LF55] and recently for more general quadrilateral lattices [CS11,Sko13,BS16]. In this article, we prove the convergence of discrete holomorphic integrals (Abelian integrals of first kind) obtained from suitable triangulations of the ramified covering to their continuous counterparts (Theorem 2.6).
Our main results are stated in Section 2 and proved in Section 3. The proof is inspired by [BS16] and uses energy estimates which allow to prove the convergence of the discrete period matrices directly. Our results are also applied to improve the convergence results of [BS16] in Section 5. Finally, in Section 6, we present some numerical experiments.

Convergence results for discrete period matrices and discrete holomorphic integrals for ramified coverings ofĈ
In the following, we consider any compact Riemann surface R of genus g ≥ 1 which allows a branched covering map f : R →Ĉ. Using this covering map as a local chart, we always locally identify points in R with their images inĈ. Then for points inĈ \ {∞} we apply the standard stereographic projection to the complex plane C. This map from R to C is denoted by P r R and gives a local chart in a neighborhood about every point, except at branch points. For further use, we fix a radius > 1 such that the images of all branch points, except possibly ∞, have a distance at most /2 from the origin. Let T = T R be a triangulation of R such that all branch points are vertices. We assume that every triangle is contained in only one sheet of the covering. We will mostly consider this triangulation via its (local) image under the chart P r R . In this sense, without further mention, we always identify this triangulation with the corresponding (multi-sheeted) triangulation onĈ (which is the image f (T ) under the covering map) and with the (multi-sheeted) image of this triangulation of C under the map P r R , excluding the vertex at infinity. We assume that this triangulation is a locally planar embedding in the complex plane C or equivalently in the Riemann sphereĈ, except at the branch points and possibly at ∞. From now on, we consider the vertices of the triangulation as points of C, except ∞, that is, we always apply P r R , i.e. the covering map and the standard stereographic projection. The edges connecting incident vertices will be straight line segments or circular arcs in C, depending on the following distinction.
(i) All triangles with at least two vertices in the open disk B (0) of radius about the origin are geodesic, that is Euclidean triangles. We always consider these triangles to be embedded in C.
(ii) All triangles whose vertices are all contained in the complement C \ B (0) are preimages of a geodesic triangulation with Euclidean triangles under the map z → 1/z. Therefore, these triangles are in general bounded by circular arcs.
Note that we mostly consider the images of these triangles under the map z → 1/z which are Euclidean triangles embedded into C in a neighborhood of the origin. (iii) The remaining triangles in the 'boundary region' are consequently in general bounded by two straight lines and one circular arc. These triangles will be called boundary triangles and denoted by F . Finally, we assume that the edge lengths of all boundary triangles are strictly smaller than max{ /2, 1}. As in the first case, we always consider these triangles to be embedded in C.
We denote by V, E, E, F the sets of vertices, edges, oriented edges, and faces of T R , respectively, and identify them locally with their images under the map P r R .

Discrete harmonic functions
We define weights on the edges E of the triangulation T R essentially by using cotan-weights, but we distinguish three cases for edges e = [x, y] ∈ E corresponding to the different cases above: (i) If both triangles incident to e are contained in the open disk B (0), we use cotan-weights where α e and β e are the angles opposite to the edge e ∈ E in the two adjacent triangles, see Figure 1.
(ii) If both triangles incident to e are contained inĈ\B (0), we consider the image of the triangulation under the map z → 1/z. We define the edge weights by (1) for the image triangles. Note that this amounts to using in (1) the angles opposite to the edge e ∈ E in the two adjacent circular arc triangles.
(iii) If e = [x, y] is incident to a boundary triangle in F , we define the weight similarly as above as a sum c(e) = C 1 + C 2 of two parts corresponding to the two incident triangles ∆ 1 , ∆ 2 . If there is a non-boundary triangle, say ∆ 1 , incident to e, we consider the angle α e in this triangle opposite to e and set C 1 = 1 2 cot α e . The second part C 2 = C [x,y] is defined below in (4) using a suitable interpolation function and the smooth Dirichlet energy. More details and explicit calculations are given in Appendix A.1.
Using our edge weights, we can define discrete harmonicity and a discrete Dirichlet energy for functions u : V → R on the vertices of the triangulation T R . In particular, u is called discrete harmonic if for every vertex x ∈ V there holds y∈V :[x,y]∈E c([x, y])(u(y) − u(x)) = 0. (2) The energy of u is The motivation for our choice of weights, in particular for the choice of weights for boundary triangles, is the following connection of discrete and smooth Dirichlet energies. Recall that for a continuous function on a compact Riemann surface R which is smooth almost everywhere the Dirichlet energy is defined as Then the discrete energy of a function u : V → R is in fact the Dirichlet energy of the continuous interpolation function I T u, defined piecewise on every triangle ∆[x, y, z] as follows: (i) If at least two of the three vertices x, y, z are contained in the open disk B (0), we define I T u| ∆[x,y,z] on this triangle as the linear interpolation of the values of u at the vertices.
(ii) If all vertices x, y, z are contained inĈ \ B (0), we consider the image of the triangle under the map z → 1/z. This is a Euclidean triangle ‹ ∆. Let u = u • 1/z be the corresponding transformed function and denote by fi I T u ∆ the linear interpolation of u on ‹ ∆. We define I T u on the original triangle as the corresponding value of fi  It is easy to see that for every triangulation of a ramified covering R ofĈ as above, I T u is a well-defined continuous function on R. Furthermore, we have Lemma 2.1 (Interpolation lemma).
Proof. We can split the energy according to triangles ∆ f for f ∈ F .
In particular, elementary calculations show that for any triangle ∆[x, y, z] we have where the constants C [x,y] , C [y,z] , C [z,x] only depend on the triangle ∆[x, y, z]. Duffin showed in [Duf59,§ 4] that for Euclidean triangles these constants are one half of the cotan of the opposite angles. Using the conformal invariance of the Dirichlet integral for the triangles inĈ \ B (0) and our choice of weights from (4), we obtain the claim.
Remark 2.2. It will be important to note that the constants C [x,y] , C [y,z] , C [z,x] defined by (4) are only small perturbations of the usual cotan-weights in the following sense. If the maximal edge length h in the boundary triangle ∆[x, y, z] ∈ F is small enough and the angles in ∆[x, y, z] as well as the angles in the Euclidean triangle formed by the vertices x, y, z lie in [δ, π − δ] for some δ > 0, then elementary calculations and estimates show that 1 2 cotα e − C δ, · h ≤ C e ≤ 1 2 cotα e + C δ, · h. for some constant C δ, > 0, whereα e denotes the angle opposite to the edge e in the Euclidean triangle formed by the vertices x, y, z. The details are given in Appendix A.1.
Note that the difference between the the angles in the Euclidean triangle with vertices x, y, z and the corresponding angles in ∆[x, y, z] is of order h. Thus, for h small enough, the corresponding estimates on |C e − 1 2 cot α e | also hold for the actual angles α e in the triangle ∆[x, y, z].

Discrete analytic functions, discrete holomorphic integrals and discrete period matrices
In the following, we define discrete analytic functions and discrete holomorphic integrals analogously as in [BS16].
For an oriented edge e ∈ E, we denote by h e ∈ V and t e ∈ V the head and the tail of e, and by l e ∈ F and r e ∈ F the left shore and the right shore of e, respectively, see Figure 1. Two functions u : The pair f = (u : V → R, v : F → R) of two conjugate functions is called a discrete analytic function. We write Ref := u and Imf := v. If both u and v are constant functions, not necessarily equal to each other, we write f = const. A direct checking shows that on simply connected surfaces discrete harmonic functions are precisely real parts of discrete analytic functions. Note that for non-zero weights c(e) = 0 we define the (discrete) energy of a function v : F → R by We will consider multi-valued functions on the vertices and the faces of the triangulation T . Informally, a multi-values function changes its values after performing some nontrivial loop on the surface.
Recall that the Riemann surface R is a branched covering ofĈ with genus g ≥ 1. Denote by p : R → R the universal covering of R and by p : ‹ T → T the induced universal covering of T . Fix a base point z 0 ∈ R and closed paths α 1 , . . . , α g , β , . . . , β g : [0, 1] → R forming a standard basis of the fundamental group π 1 (R, p(z 0 )) such that α 1 β 1 α −1 1 β −1 1 · · · α g β g α −1 g β −1 g is null-homotopic. Each closed path γ : [0, 1] → R with γ(0) = γ(1) = p(z 0 ) determines the deck transformation d γ : R → R, that is, the homeomorphism such that p • d γ = p and d γ (z 0 ) =γ(1), whereγ : R → R is the lift of γ : [0, 1] → S such thatγ(0) = z 0 . The induced deck transformation of ‹ T is also denoted by d γ : ‹ T → ‹ T . A multi-valued function with periods A 1 , . . . , A g , B 1 , . . . , B g ∈ C is a pair of functions f = (Ref : The numbers A 1 , . . . , A g and B 1 , . . . , B g are called the A-periods and the B-periods of the multivalued function f , respectively. Analogously, we can also define multi-valued functions u : ‹ V → R, v : ‹ F → R or u : R → R. Note in particular, that for each multi-valued function u : ‹ V → R and every edge [x, y] ∈ E the difference u(x) − u(y) is well defined. The (discrete) energy of the multi-valued function is Similarly, for each multi-valued function u : R → R, which is smooth on every face of ‹ F , at each point inside a face ∆ ∈ F the gradient ∇u is well defined. The (Dirichlet) energy of the multi-valued function is A multi-valued discrete analytic function is called a discrete holomorphic integral or discrete Abelian integral of the first kind. For each l = 1, . . . , g denote by φ l T = (Reφ l T : ‹ V → R, Imφ l T : ‹ F → R) the unique (up to constant) discrete holomorphic integral with A-periods given by A k = δ kl , where k = 1, . . . , g. The g × g-matrix Π T whose l-th column is formed by the B-periods of φ l T , where l = 1, . . . , g, is called the period matrix of the triangulation T .

Convergence of energy and discrete period matrices
So far, we have defined our notions like discrete energy for a rather general class of triangulations. In view of our convergence results, we now make some additional assumptions. In orde to measure distances and other metric properties, we always consider the images of the triangles in C under the projection P r R . By our assumptions above, these are Euclidean triangles if they are contained in B (0) and we use the standard metric in C. We also apply this metric for boundary triangles in F . For triangles which are mapped toĈ \ B (0), we consider their image under the map z → 1/z and then use again the standard metric. Alternatively, we could work onĈ with the chordal metric.
First we determine the maximal distance between two vertices in a triangle which lies inside B (0) or on the boundary F and the maximal edge length of the triangles outside B (0) after the mapping by 1/z. The maximum of these two numbers is called maximal edge length and denoted by h = h(T ).
Furthermore, we suppose that near all branch points of R the edge lengths are adapted to the singularity which then guarantees an approximation error of order h. In particular, for every branch Furthermore, we assume that all these disks are all contained in we first apply the mapping z → 1/z and assume that r O = 1/(2 ). For all branch points we already have a natural complex structure and charts. In In any case, consider all triangles in C R O which are incident to v O . The aperture of O is the sum of all the face angles at O of the projection of these triangles. Denote by γ O the value 2π divided by the aperture. Note that for branch points we have γ O ∈ {1/n : n = 2, 3, 4, . . . }, so γ O ≤ 1/2.
We demand that the triangles in the neighborhood C O of O have an adapted size: as an additional condition, we demand that • the images under the chart g O of any two incident vertices in C O have maximum distance h.
where |Oz| denotes the distance of z to O in C. Then we deduce from our assumption that the maximal edge length in ∆ is smaller than h · |Oz| 1−γ O .
In Section 5 we explain how our ideas can be used for polyhedral surfaces with more general conical singularities with 0 < γ O ≤ 1/2.
We will always assume that the maximal edge length h is strictly smaller than max{ /4, r O /4, 1}.
A triangulation T which satisfies these additional properties for all its branch points will be called adapted triangulation.
Theorem 2.4 (Energy convergence). For each δ > 0 and each smooth multi-valued harmonic function u : R → R there are two constants Const u,δ,R, , const u,δ,R, > 0 such that for any adapted triangulation T of R with maximal edge length h < const u,δ,R, and minimal face angle > δ we have The assumption on the minimal face angle in the theorem cannot be dropped, see [BS16, Example 4.14].
Based on energy estimates from this theorem, we deduce convergence of discrete period matrices. To this end, recall that R is a Riemann surface which is a branched covering ofĈ. Therefore, a basis of holomorphic integrals φ l R : R → C and the period matrix Π R of R are defined analogously to the discrete case above.
Theorem 2.5 (Convergence of period matrices). For each δ > 0 there are two constants Const δ,R, , const δ,R, > 0 such that for any adapted triangulation T of R with maximal edge length h < const u,δ,R, and minimal face angle > δ we have Both theorems are proved in Section 3.

Convergence of discrete holomorphic integrals
For the next theorem, we need some additional notions similarly as in [BS16]. The discrete holo- Recall that a triangulation T is Delaunay, if for every edge e ∈ E we have α e + β e ≤ π. Let {T n } be a sequence of adapted triangulations of the surface R with maximal edge length h < max{ /4, r 0 /4, 1}. Such a sequence of adapted triangulations is called non-degenerate uniform, if there is a constant Const, not depending on n, such that for each member of the sequence: (A) the angles of each triangle are greater than 1/Const.
. Then we require that after this mapping in each disk of radius equal to the maximal edge length the number of image points of vertices is smaller than Const.
A sequence of functions f n = (Ref n : ‹ V n → R, Imf n : ‹ F n → R) converges to a function f : R → C uniformly on compact subsets if for every compact set K ⊂ R we have Theorem 2.6 (Convergence of holomorphic integrals). Let {T n } be a sequence of non-degenerate uniform adapted triangulations of R with maximal edge length h n → 0 as n → ∞. Let z n ∈ ‹ V n be a sequence of vertices converging to a point z 0 ∈ R. Let ∆ n ∈ ‹ F n be a sequence of faces with its vertices converging to z 0 . Then for each 1 ≤ l ≤ g the discrete holomorphic integrals normalized at z n and ∆ n converge uniformly on each compact set to the holomorphic integral φ l R : R → C normalized at z 0 . This theorem is proved in Section 4.

Proof of convergence of energy and period matrices
In this section, we prove convergence of the discrete energy to the corresponding Dirichlet energy and convergence of discrete period matrices to their continuous counterparts. The main ideas of the proof follow [BS16, Section 4], but we improve the estimates near branch points (which can be considered as special conical singularities) by using the additional properties of the adapted triangulations.
We denote by Const a,b,c a positive constant which only depends on the parameters a, b, c. The symbol Const may denote distinct constants at different places of the text, for example in 2 · Const ≤ Const. Furthermore, we set D k u(z) := max 0≤j≤k In the following, all triangle which are considered are in C after application of P r R .

Energy estimates in a triangle
First we consider only one triangle ∆ of the triangulation T . Let u : ∆ → R be a smooth function which smoothly extends to a neighborhood of ∆. Let I T u : ∆ → R be the corresponding interpolation function defined in Section 2.1. Then we set E ∆ (u) = ∆ |∇u| 2 dxdy and E T ∆ (u) = ∆ |∇I T u| 2 dxdy. Denote by δ the minimal angle of the triangle ∆.
(iii) If the triangle ∆ is a boundary triangle in F then Proof. (i) For triangles contained in B this is Lemma 4.4 in [BS16].
(ii) For triangles contained inĈ \ B , note that the discrete energy is actually defined using the image of ∆ under the map 1/z and the corresponding functionũ = u • 1/z. By conformality, the smooth Dirichlet energy can also be considered on ‹ ∆. Therefore, the same arguments as in (i) apply.
(iii) For boundary triangles we estimate the discrete and the smooth energy separately. For the smooth energy, we have . For the discrete energy, an estimate of E T ∆ (u) is contained in Appendix A.2, see in particular (17).

Energy estimates near a a branch point
Let v O ∈ R be a branch point of R with γ O ≤ 1/2. In this subsection we only consider those triangles of the adapted triangulation T which are completely contained in the neighborhood C R O . Denote by T O the connected component of these triangles which contains v O . For the estimate of the difference of energies for these triangles we consider in particular If O = ∞, we compose as above the map P r R with z → 1/z and consider the corresponding function u • 1/z instead of u. Then the following reasoning also applies to the image branch point at the origin and its 1/(2 )-neighborhood. As Proof. We use Lemma 3.4, sum the inequalities and estimate the integral.
Proof. We deduce from Lemma 3.3 similarly as in the previous lemma that T respectively. We consider the energies on both parts separately.

Convergence of energies
Our assumption on the maximal edge length, the definition of the discrete energy, the compactness of R, and the estimates in Lemma 3.1 imply that For all triangles ∆ ∈ G (i) T we consider the image ‹ ∆ under the map z → 1/z. Using the corresponding map u = u • 1/z we obtain analogously Lemma 3.9. We have E S (u) ≤ Const u, · h, where S denotes the subset of R which is covered by triangles of F .
Proof. This estimate is due to the fact that the derivative of u is bounded away from the branch points. Furthermore, the area of the ring { − h ≤ |z| ≤ + h} is bounded by 4π h and the degree of the branched covering is fixed. Therefore, The proof of this lemma is given in Appendix A.2.
Proof of Theorem 2.4. Summing up the estimates obtained in Lemmas 3.5-3.10 we get the desired result:
Denote u T,P = Reφ T,P , where φ T,P is the discrete holomorphic integral defined in Theorem 3.12 for each vector P ∈ R 2g . Analogously, let φ R,P : R → C be a holomorphic integral whose periods have real parts A 1 , . . . , A g , B 1 , . . . , B g , respectively. Denote u R,P = Reφ R,P .
Lemma 3.13. For every δ > 0 and every vector P ∈ R 2g there are constants Const P,δ,R, , const P,δ,R, > 0 such that for any adapted triangulation T of R with maximal edge length h < const P,δ,R, and minimal face angle δ > 0 we have Proof. From the interpolation lemma 2.1 we know that E T (u T,P ) = E(I T u T,P ) and the interpolation function I T u T,P is continuous and piecewise smooth. Using Lemma 3.11 and its smooth counterpart (Dirichlet's principle) we deduce from Theorem 2.4 that For each l = 1, . . . , g denote by φ l T * = (Reφ l T * : ‹ V → R, Imφ l T * : ‹ F → R) the unique (up to constant) discrete holomorphic integral with A-periods given by A k = iδ kl , where k = 1, . . . , g. The g × g-matrix Π T * whose l-th column is formed by the B-periods of φ l T * divided by i, where l = 1, . . . , g, is called the dual period matrix of the triangulation T .
The following theorem connects the period matrices to the energies.

Lemma 3.14 ([BS16, Lemmas 3.14 & 3.15]). (i) The energy E T (u T,P ) is a quadratic form in the vector P ∈ R 2g with the block matrix
The energy E(u R,P ) is a quadratic form in the vector P ∈ R 2g with the block matrix Combining Lemmas 3.13 and 3.14, we obtain: Corollary 3.15. Let {T n } be a nondegenerate uniform sequence of adapted triangulations of R with maximal edge length tending to zero as n → ∞. Let P n ∈ R 2g be a sequence of 2g-dimensional real vectors converging to a vector P ∈ R 2g . Then E Tn (u Tn,Pn ) → E(u R,P ) as n → ∞.
Proof of Theorem 2.5. Both E Tn (u Tn,P ) and E(u R,P ) are quadratic forms in P ∈ R 2g by Lemma 3.14 with block matrices E T and E R , respectively. Thus by Lemma 3.13 for every δ > 0 there are constants Const δ,R, , const δ,R, > 0 such that for any adapted triangulation T of R with maximal edge length h < const δ,R, and minimal face angle δ > 0 we have E T − E R ≤ Const δ,R, · h. From this inequality we deduce estimates on ReΠ T − ReΠ R and ImΠ T − ImΠ R of the same type, but with different constants which are derived in the following. These estimates complete the proof.

Proof of convergence of discrete holomorphic integrals
The strategy of the proof of Theorem 2.6 follows the corresponding ideas in [BS16, Section 5]. Due to our different setup, we need some modifications.

Equicontinuity
In this section we consider triangulations T of branched coverings with boundary. The main goal is to consider (sufficiently small) intrinsic discs about a branch points or about a regular point and derive an estimate for harmonic functions there. A function u : V → R is discrete harmonic on T if it satisfies (2)  Let T be a non-degenerate uniform adapted triangulation of the branched covering of R. We assume that T is a simply connected part of T . For simplicity, we directly consider the projection of all triangles into C by P r R .

Lemma 4.1 (Equicontinuity lemma). (i) Let T be contained in an open disc
For |z − w| < h < r/3 the same inequality holds with |z − w| replaced by h .
(10) In case (ii), we consider the harmonic function u as defined on the image triangulation T g . The proof only uses the fact that u satisfies a maximum principle which still holds in our case. For the third case, we just work with the triangulation T 1/z and assume without loss of generality that u is defined there.

Lemma 4.2. Let T be a triangulation of a ramified covering with boundary such that all angles are
in [δ, π − δ] for some π/4 > δ > 0. Then there exist constants const δ, , Const δ, > 0 such that for 0 < h < const δ, and every function u : Proof. Let ∆ ∈ F be a triangle with vertices x, y, z ∈ V such that [x, z] is a boundary edge of T . Denote the angle in ∆ at the vertex v ∈ {x, y, z} by α v .
First consider the case that ∆ ∈ F is no boundary triangle. We want to show that holds for some constant Const δ > 0. Thus we only need to consider the case α y > π/2. Take Const δ = 1/(cot 2 δ − 1). As α x , α z > δ and α x + α y + α z = π, elementary calculations imply that This implies (11). If ∆ ∈ F , we know that where α E v denotes the angle at the vertex v in the Euclidean triangle with vertices x, y, z and |r v | ≤ Const δ, , see Appendix A.1. Therefore, there are constants const δ, , ‡ Const δ, > 0 such that for all 0 < h < const δ, we have E T ∆ (u) ≥ ‡ Const δ, |C [x,z] |(u(z) − u(x)) 2 . Take Const δ, = max{ ‡ Const δ, , Const δ }, sum the above inequalities over all such faces and deduce E T (u) − E T (u) ≤ Const δ, · E T (u). Now the claim follows.

Convergence of multi-valued discrete harmonic functions and discrete holomorphic integrals
As a first step, we can deduce that the uniform limit of a sequence of discrete harmonic functions is harmonic. To this end, we say that a sequence of triangulated polygons {T n } approximates a domain Ω ⊂ C, if for n → ∞ the following three quantities tend to zero: the maximal distance from a point of the boundary ∂T n to the set ∂Ω, the maximal distance from a point of ∂Ω to the set ∂T n , and the maximal edge length of the triangulation T n . Theorem 4.4 (Convergence of multi-valued discrete harmonic functions). Let {T n } be a nondegenerate uniform sequence of adapted Delaunay triangulations of R with maximal edge length h n tending to zero as n → ∞. Let z n ∈ ‹ V n be a sequence of vertices converging to a point z 0 ∈ R. Let P n ∈ R 2g be a sequence of vectors converging to a vector P ∈ R 2g . Then the functions u Tn,Pn : ‹ V n → R satisfying u Tn,Pn (z n ) = 0 converge to u R,P : R → R with u R,P (z 0 ) = 0 uniformly on every compact subset.
Proof. We will start with some estimates on compact subsets of R of a special form. Let π = P r R •p : R → C be the local projection map P r R composed with the universal covering p. For v ∈ R denote by ‹ B r (v) ⊂ R the subset which projects for π(v) ∈ C to an open intrinsic disc B r (π(v)) = π( ‹ B r (v)) with radius r about π(v). If π(v) = ∞, we assume that π( ‹ B r (v)) = C \ B 1/r (0). We restrict ourselves to the following cases: • π(v) = O is a branch point and r = r O (v) > 0 its associated radius defined in Section 2.3, • π(v) = ∞ and r = 1/ .
Note that the union of these sets ‹ B r (v) covers R and every compact set K ⊂ R is contained in the union of finitely many of these sets.
Let ‹ B r (v) be one of these sets. Consider those triangles of the given triangulation › T n which are completely contained in ‹ B r (v) and denote by ‹ T n (v, r) the connected component of these triangles which contains v. Choose n 1 such that for all n > n 1 the maximal edge length h n < r min O /200. Consider u Vn(v,r) := u Tn,Pn | Vn(v,r) . By Lemma 4.2 and Corollary 3.15 the sequence of energies E Tn(v,r) (u Vn(v,r) ) is bounded. Thus the Equicontinuity lemma 4.1 implies that the function u Vn(v,r) | Vn∩ B 3 4 r (v) has uniformly bounded differences. That is, there exists a constant Const R,P,δ such that for all n > n 1 and z, w ∈ ‹ V n ∩ ‹ B3 4 r (v) we have |u Tn,Pn (z) − u Tn,Pn (w)| ≤ Const R,P,δ . Lemma 4.1 also implies that the sequence is equicontinuous, that is, there exists a function δ(ε) for ε > 0 such that for each n > n 1 and z, w ∈ ‹ V n ∩ ‹ B 3 4 r (v) with |z − w| < δ(ε) we have |u Tn,Pn (z) − u Tn,Pn (w)| ≤ ε. Now take a sequence of compact sets K 1 ⊂ K 2 ⊂ · · · ⊂ R such that R = ∞ j=1 K j . Assume that K 1 contains all point of the convergent sequence {z n }. Since K 1 is compact, it is contained in the union of finitely many of the sets considered above. Therefore, the sequence u Tn,Pn | Vn∩K 1 is equicontinuous and has uniformly bounded differences (this bound also depends on K 1 ). Furthermore, as all z n ∈ K 1 and u Tn,Pn (z n ) = 0, the sequence u Tn,Pn | Vn∩K 1 is uniformly bounded. We deduce from the Arzelà-Ascoli theorem that there is a continuous function u 1 : K 1 → R and a subsequence {l k } with l 1 = n 1 such that u T l k ,P l k converges to u 1 uniformly on K 1 .
Analogously, we see that there is a continuous function u 1 : K 1 → R and a subsequence {m k } of {l k } with m 1 = n 1 , m 2 = l 2 such that u Tm k ,Pm k converges to u 2 uniformly on K 2 . Clearly, we have u 1 = u 2 on K 1 . This procedure can be continued and eventually we obtain a continuous function u : R → R and a subsequence {n k } of {1, 2, 3, . . . } such that u Tn k ,Pn k converges uniformly to u on each compact subset of R. Also, u has the same periods P as u R,P and u(z 0 ) = 0. Applying Lemma 4.3 to bounded domains not containing any branch point, we see that the limit function u : R → R is harmonic in R except possibly at the branch points. But as u is locally bounded, these singularities can be removed and therefore the continuous function u is in fact harmonic on the whole surface R. Thus u = u R,P by our normalization u(z 0 ) = 0 = u R,P (z 0 ).
Since the limit function u = u R,P is unique, it follows that the whole sequence u Tn,Pn , not just the subsequence u Tn k ,Pn k , converges to u R,P uniformly on every compact subset.
Proof of Theorem 2.6. Let P n , P ∈ R 2g be the periods of the real parts Reφ l Tn : ‹ V → R and Reφ l R : R → R of the discrete and smooth holomorphic integrals, respectively. Then by Theorem 3.12 Reφ l Tn = u Tn,Pn and Reφ l R = u R,P . Theorem 2.5, implies that P n → P as n → ∞. Thus we deduce from Theorem 4.4 that the real parts Reφ l Tn converge to Reφ l R uniformly on every compact subset. Convergence of the imaginary parts is proven analogously due to the following Lemma 4.5. Proof. This follows immediately from (6) together with the definitions of the discrete energies in (3) and (7).

Improved convergences of period matrices and holomorphic integrals for polyhedral surfaces
The techniques applied for adapted triangulations near branch points may also be used to improve the order of convergence of period matrices and holomorphic integrals for polyhedral surfaces compared to the results obtained in [BS16]. A polyhedral surface S is an oriented two-dimensional manifold without boundary which has a piecewise flat metric with isolated conical singularities. An example is the surface of a polyhedron in three-dimensional space. Let T S be a geodesic triangulation of the polyhedral surface S such that all faces are flat triangles. Note in particular, that all singular points of the metric are vertices of T S . On all edges we use cotan weights given by (1). If γ O > 1/2, we do not adapt the triangulation further. But for singularities O with γ O ≤ 1/2 we consider a chart g O , which maps a neighborhood C O of O to a neighborhood of the origin in C. Furthermore, we can introduce as above "polar coordinates" (r, φ) on C O with the origin at the vertex O. We map all vertices in C O to a neighborhood of the origin in C by the chart g O : where |Oz| denotes the distance of z to O in S. As in Section 2.3 we deduce from our assumption that the maximal edge length in ∆ is smaller than h · |Oz| 1−γ O .
Applying the estimates of Sections 3.1 and 3.2, we obtain the following improved versions of Theorems 2.5 and 2.7 of [BS16].
Theorem 5.1 (Energy convergence). For each δ > 0 and each smooth multi-valued harmonic function u : ‹ S → R there are two constants Const u,δ,S , const u,δ,S > 0 such that for any adapted triangulation T of S with maximal edge length h < const u,δ,S and minimal face angle > δ we have Theorem 5.2 (Convergence of period matrices). For each δ > 0 there exist constants Const δ,S , const δ,S > 0 such that for any adapted triangulation T of S with maximal edge length h < const u,δ,S and minimal face angle > δ we have Theorem 5.3 (Convergence of holomorphic integrals). Let {T n } be a sequence of non-degenerate uniform adapted triangulations of S with maximal edge length h n → 0 as n → ∞. Let π : ‹ S → S be the universal covering of S. Denote by ‹ T n the corresponding triangulation of ‹ S such that π( ‹ T n ) = T n . Let z n ∈ ‹ V n be a sequence of vertices converging to a point z 0 ∈ ‹ S. Let ∆ n ∈ ‹ F n be a sequence of faces with its vertices converging to z 0 . Then for each 1 ≤ l ≤ g the discrete holomorphic integrals φ l T = (Reφ l T : ‹ V → R, Imφ l T : ‹ F → R) normalized at z n and w n converge uniformly on each compact set to the holomorphic integral φ l S : ‹ S → C normalized at z 0 .

Numerical experiments
In the following, we present some numerical analysis for our convergence results detailed above. We are very grateful to Stefan Sechelmann for writing software and performing numerical experiments. Mainly, we apply the scheme described in Section 2, but we consider the triangulations on the sphere S 2 ∼ =Ĉ without stereographic projection to C. Furthermore, we use an approximation of the discrete energy E T and of the discrete multi-valued harmonic functions u T,P because we use slightly different weights instead of those given in Section 2.1. In particular, each triangle ∆ of T as an embedded triangle in S 2 before stereographic projection has circular arcs as edges where the circles all pass through the north pole (= ∞) if ∆ is contained in the spherical region corresponding to B (0) or all pass through the south pole (= 0) if ∆ is contained in the spherical region corresponding toĈ \ B (0), respectively. For boundary triangles, there are two types of circular arcs. For practical reasons, we do not work with these triangles in S 2 ⊂ R 3 . Instead, we take the vertices and add straight line segments in R 3 between incident vertices. For every original triangle ∆ on S 2 we obtain a corresponding triangle ∆ S in R 3 , see Figure 2 for some examples of triangulations. As is fixed and the maximum edge length h tends to zero, any angle α in the triangle ∆ and the corresponding angle α S in the triangle ∆ S only differ by an error of order h, in particular |α − α S | ≤ Const δ,ρ · l max (∆). Thus, for uniform Delaunay triangulations which we consider, the weights c S (e) = 1 2 (cot α S e +cot β S e ), using the cotan-formula for the angles of the triangles ∆ S , can be estimated using the original weights c(e) for T . More precisely, there is a constant Const δ, such that c(e) . Thus for 0 < h < 1/Const δ, and some constant Const P,δ,R, > 0 we obtain similarly as in the proof of Lemma 3.13 that As concrete examples we consider two surfaces with known period matrices, namely the torus T of genus 1 with branch points 0.5 + 0.4i, −0.3 + 0.2i, −0.1, 0.1 − 0.2i and Lawson's minimal surface L of genus 2 which corresponds to the hyperelliptic curve µ 2 = λ 6 − 1 with branch points e ikπ/3 , k = 0, 1, . . . , 5. The smooth period matrices are Π T ≈ 0.836 + 0.955i and Π L = i √ 3 Ä 2 −1 −1 2 ä . We compare four different types of triangulations of S 2 which are used as basis for the computations of the discrete period matrices. Figure 2 shows an example for each of these four types. In order to simplify calculations and the construction of cycles, we always use the same triangulation on every sheet of the covering.

Random
We sample points at random on the sphere and then build the corresponding Delaunay triangulation.
Fibonacci The points on a sphere are evenly distributed by means of a Fibonacci spiral. This leads to very 'regular' triangulations with triangles which are almost equilateral, except near branch points, which are in general additional vertices.
These two types of triangulations are directly used for further computations (called Homogeneous Random ♦ and Homogeneous Fibonacci ). In these two cases the estimates of [BS16] apply. Our numerical results are plotted in the lower rows of Figures 3 and 4. The log-log plots show that the error behaves indeed like √ h for 'Homogeneous Random', where h is the maximal edge length. This was also observed in the example studied in Section 7.3 of [BS16]. Nevertheless, for the more regular triangulations using Fibonacci spirals, our numerical evidence indicates a higher order of the error bound, possibly a linear dependence on h.
According to our new idea of adapted triangulations explained in Section 2.3, we refine the examples of the two types of triangulations above in a neighborhood of the branch points by suitably  adding vertices (called Clustering Random • and Clustering Fibonacci ). Our results in Figures 3  and 4 confirm that the error between smooth and discrete period matrices decreases indeed faster in the adapted case. In particular, the log-log plots in the upper row of Figures 3 and 4, respectively, show that for our adapted method the error depends at least linearly on the maximal edge length h (as proven in Theorem 2.5) and is again possibly even of higher order for more regular triangulations using Fibonacci spirals.
Recall that the actual bounds on the approximation error depend on the angles in the triangles which differ for all our examples. In our proof we only use some (rough) estimate of the angles such that our constants in Theorems 3.13 and 2.5 depend only on the minimal angle of the adapted triangulation. We did not study the dependence of the angles in detail, but our numerical results suggest that the order of convergence also depends significantly on the regularity of the triangulations. u v = u(v) for the values of the given smooth function u at the vertices x, y, z. Therefore, we obtain ∆[x,y,z] |∇I T u| 2 = 1 0 1 0 Du ∆ (τ, σ)Dp −1 (p(τ, σ))(Dp −1 (p(τ, σ))) T (Du ∆ (τ, σ)) T | det Dp(τ, σ)|dτ dσ This gives an explicit way to calculate the edge weights. Note that by the same method we can obtain the usual cotan-weights on the Euclidean triangle with vertices x, y, z, if we use s E (t) = y + t(z − y) instead of s(t) for t ∈ [0, 1]. The function s E is the usual linear parametrization of the straight edge from y to z. An important observation is that these seemingly 'complicated' weights are in fact only small perturbations of the usual cotan-weights if the edge length is small enough. To see this, we will estimate the quantities in the above integrals compared to the corresponding quantities for s E .
Proof of the estimate in Remark 2.2. We assume that there is some δ > 0 such that all angles in the triangle ∆[x, y, z] are in [δ, π − δ] and also all angles of the Euclidean triangle with vertices x, y, z are in [δ, π − δ]. In the following, we will always assume that τ ∈ [0, 1] as in the integral terms above. Also, we are not interested in the best possible estimates, any constant, depending only on the indicated parameter, will suffice.
Combining these estimates, we have

A.2 Proof of Lemma 3.10
Note that as u is smooth, there exist constants Const u, , const u, > 0 such that for 0 < h < const u, we have |u(x) − u(y)| 2 /|x − y| 2 ≤ Const δ, · max w∈∆ D 1 u(w) 2 ≤ Const u, for all edges e = [x, y] of boundary triangles in F with edge lengths smaller than h.

A.3 Proof of Equicontinuity Lemma 4.1
First note that condition (D) from Section 2.4 implies that for every path in a non-degenerate uniform adapted triangulation T with consecutive vertices v 0 v 1 . . . v m we have Here e denotes the eccentricity as defined in Section 4.1 and for the last estimate we have used Schwarz's inequality. Now consider a simply connected triangulation T with boundary contained in an open disc B r (v) ⊂ C. This is the assumption for part (i) of the lemma. In the case of part (ii), we consider the image triangulation T g by the chart g O (z) = (z − O) γ O . By abuse of notation, we still denote this image triangulation by T . Also, we denote the vertices of T by Z and W , which are the actual vertices z, w in the first case and the images Z = g O (z), W = g O (w) in the case of a branch point. For simplicity, we assume that the edges between vertices are straight line segments, as we do not need the actual, possibly curved edges.
Let u : V → R be any function which assumes its maximum and its minimum on the boundary for any subgraph of T . Let Z, W be two distinct interior vertices of T . Denote by ZW the straight line segment joining these points. Let dist(ZW, ∂T ) be the Euclidean distance of this straight line segment to the curve of boundary edges. We assume that |Z −W | < r/3 < dist(ZW, ∂T )/3 for some r > 0. Let h denote twice the maximum circumradius of the triangles of T . Let m = r−|Z−W | 2h be the largest integer smaller than r−|Z−W | 2h . We consider auxiliary rectangles R k , k = 1, . . . , m, which are centered at (Z + W )/2 with one pair of sides parallel to ZW with length ZW + 2k · h and other pair of sides orthogonal to ZW of length 2k · h . Then the interior of R k , k = 1, . . . , m, is covered by triangles of T . Denote by V k the set of vertices contained in R k \ R k−1 , where R 0 = ZW . Then any two vertices v A , v B ∈ V k may be connected by a path v 0 v 1 . . . v N with v 0 = v A , v N = v B and all vertices v j ∈ V k as h is larger than any edge length.
Without loss of generality, assume that u(Z) ≥ u(W ). As u assumes its maximum and minimum on the boundary, there exists Z k , W k ∈ V k such that u(Z k ) ≥ u(Z) ≥ u(W ) ≥ u(W k ). The length of the path joining Z k and W k is at most the number of vertices in V k . The set of these vertices can be covered by at most Const · (|Z − W | + 4kh )/h discs of radius h /2. Therefore, by condition (U) of Section 2.4 the number of vertices in V k is less than Const · e · (|Z − W |/h + 4k). Therefore, we can estimate the energy for a path v 0 v 1 . . . v N in V k from v 0 = Z k , v N = W k using (18) This implies the desired inequalities (9) and (10).