2.1 Discretization setup

Our discretization setup is the following. (We also refer to [Reference Gaster, Loustau and Monsaingeon8, Section 2] for more details, although it focuses on the equivariant setting and $\mathbb {H}^2$ .) A mesh on M is any topological triangulation; we denote by $\mathcal {G}$ the embedded graph that is the $1$ -skeleton. A mesh (or its underlying graph) is called geodesic if all edges are embedded geodesic segments.

Denote $\mathcal {V} = \mathcal {G}^{(0)}$ and $\mathcal {E} = \mathcal {G}^{(1)}$ the set of vertices and (unoriented) edges of $\mathcal {G}$ . We shall equip $\mathcal {G}$ with vertex weights $(\mu _x)_{x \in \mathcal {V}}$ and edge weights $(\omega _{xy})_{\{x,y\} \in \mathcal {E}}$ . For now, these weights are two arbitrary and independent collections of positive numbers. Such a biweighted graph allows one to develop a discrete theory of harmonic maps $M \to N$ as follows:

• • The system of vertex weights defines a measure $\mu _{\mathcal {G}} = (\mu _x)_{x \in \mathcal {V}}$ on $\mathcal {V}$ . Since $\mathcal {G}$ is embedded in M, $\mu _{\mathcal {G}} $ can also be seen as a discrete measure on M supported by the set of vertices.

• • A discrete map from M to N along $\mathcal {G}$ is a map $\mathcal {V} \to N$ . The space $\operatorname {\mathrm {Map}}_{\mathcal {G}} (M, N)$ of such maps is a smooth finite-dimensional manifold with tangent space

  $$ \begin{align*} \operatorname{\mathrm{T}}_f \operatorname{\mathrm{Map}}_{\mathcal{G}} (M, N) = \Gamma(f^* \operatorname{\mathrm{T}} N) := \bigoplus_{x \in \mathcal{V}} \operatorname{\mathrm{T}}_{f(x)} N. \end{align*} $$
  It carries a smooth $\operatorname {\mathrm {L}}^2$ -Riemannian metric given by
  $$ \begin{align*} \langle V, W \rangle := \int_M \langle V_x, W_x\rangle \, \mathop{}\!\mathrm{d} \mu_{\mathcal{G}}(x) = \sum_{x\in \mathcal{V}} \mu_x \langle V_x, W_x\rangle, \end{align*} $$
  and an associated $\operatorname {\mathrm {L}}^2$ distance given by
  $$ \begin{align*} d(f,g)^2 := \int_M d(f(x), g(x))^2 \, \mathop{}\!\mathrm{d} \mu_{\mathcal{G}}(x) = \sum_{x \in \mathcal{V}} \mu_x d(f(x), g(x))^2, \end{align*} $$
  where $d(f(x), g(x))$ denotes the Riemannian distance in N.

• • The discrete energy density of a discrete map $f \in \operatorname {\mathrm {Map}}_{\mathcal {G}} (M, N)$ is the discrete nonnegative function $e_{\mathcal {G}} (f) \in \operatorname {\mathrm {Map}}_{\mathcal {G}} (M, \mathbb {R})$ defined by

  $$ \begin{align*} e_{\mathcal{G}} (f)_{x} = \frac{1}{4\mu_x} \sum_{y\sim x} \omega_{xy} \, d(f(x),f(y))^2. \end{align*} $$

• • The discrete energy functional on $\operatorname {\mathrm {Map}}_{\mathcal {G}} (M, N)$ is the map $E_{\mathcal {G}} \colon \operatorname {\mathrm {Map}}_{\mathcal {G}} (M, N) \to \mathbb {R}$ given by

  (2.1) $$ \begin{align} \begin{aligned} E_{\mathcal{G}} (f) &= \int_M e_{\mathcal{G}} (f) \mathop{}\!\mathrm{d} \mu_{\mathcal{G}} \\[3pt] &= \frac12 \sum_{x\sim y} \omega_{xy} \, d(f(x),f(y))^2. \end{aligned} \end{align} $$
  A discrete harmonic map is a critical point of $E_{\mathcal {G}} $ .

  Remark 2.1 The discrete energy functional does not depend on the choice of vertex weights, and neither does the harmonicity of a discrete map. When M is two-dimensional, this reflects the fact that the energy functional $E \colon \mathcal {C}^{\infty }(M,N) \to \mathbb {R}$ only depends on the conformal structure on S.

• • The discrete tension field of $f \in \operatorname {\mathrm {Map}}_{\mathcal {G}} (M, N)$ is $\tau _{\mathcal {G}} (f) \in \Gamma (f^* \operatorname {\mathrm {T}} N)$ defined by

  $$ \begin{align*} \tau_{\mathcal{G}} (f)_{x} = \frac{1}{\mu(x)} \sum_{y \sim x} \omega_{xy} \overrightarrow{f(x) f(y)}. \end{align*} $$

  Notation 2.2 Throughout the paper, we abusively denote $\overrightarrow {xy} := \exp _{x}^{-1}(y)$ (whenever well-defined), where $\exp _x$ is the Riemannian exponential map.

  In [Reference Gaster, Loustau and Monsaingeon8, Proposition 2.21], we show the discrete first variational formula:
  $$ \begin{align*} \tau_{\mathcal{G}} (f) = - \operatorname{\mathrm{grad}} E_{\mathcal{G}} (f). \end{align*} $$
  In particular, f is harmonic if and only if $\tau _{\mathcal {G}} (f) = 0$ . This is equivalent to the property that for all $x\in \mathcal {V}$ , $f(x)$ is the center of mass of its neighbors’ values (more precisely of the system $\{f(y), \omega _{xy}\}$ for y adjacent to x [Reference Gaster, Loustau and Monsaingeon8, Proposition 2.22]).

• • Given $u_0 \in \operatorname {\mathrm {Map}}_{\mathcal {G}} (M, N)$ and $t>0$ , the discrete heat flow with fixed stepsize t is the sequence $(u_k)_{k \geqslant 0}$ defined by

  $$ \begin{align*} u_{k+1} = \exp(t \, \tau_{\mathcal{G}} (u_k)). \end{align*} $$
  The discrete heat flow is precisely the fixed stepsize gradient descent method for the discrete energy functional $E_{\mathcal {G}} $ .

One of the main theorems of [Reference Gaster, Loustau and Monsaingeon8] is that if $S = M$ and N are closed oriented surfaces of negative Euler characteristics and $u_0$ has nonzero degree, then the discrete heat flow converges as $k \to +\infty $ to the unique minimizer of $E_{\mathcal {G}} $ in the same homotopy class with exponential convergence rate. See [Reference Gaster, Loustau and Monsaingeon8, Theorem 4.5] for more details.

2.2 Midpoint subdivision of a mesh

Assume $(M,g)$ is equipped with a geodesic mesh and denote by $\mathcal {G}$ the associated graph. One can define a new mesh called the midpoint subdivision (or refinement) as follows. For comfort, let us assume $M = S$ is two-dimensional; the definition is easily generalized. Define a new geodesic graph $\mathcal {G}'$ by adding to the vertex set of $\mathcal {G}$ all the midpoints of edges of $\mathcal {G}$ , and adding new edges so that every triangle in $\mathcal {G}$ is subdivided as four triangles in $\mathcal {G}'$ (see [Reference Gaster, Loustau and Monsaingeon8, Definition 2.2]). This clearly defines a new geodesic triangulation of S whose $1$ -skeleton is $\mathcal {G}'$ . See Figure 1 for an illustration of an invariant mesh in $\mathbb {H}^2$ and its refinement generated by the software Harmony.

Figure 1:

A mesh of the Poincaré disk model of $\mathbb {H}^2$ on the left and its midpoint refinement on the right. Both are invariant under the action of a Fuchsian group $\Gamma $ , yielding meshes on a closed hyperbolic surface S of genus $2$ . The highlighted central region is a fundamental domain. The blue circle arcs are the axes of generators of $\Gamma \approx \pi _1 S$ .

Evidently, this subdivision process may be iterated, thus one can define the refinement of order n of a geodesic mesh. Meshes obtained by successive midpoint refinements will be our standard support for approximating a smooth manifold by discrete data. Properties of such meshes will be further discussed in Section 6.

2.3 Interpolation

2.3.1 Generalities

Assume $(M,g)$ is equipped with a geodesic mesh and denote by $\mathcal {G}$ the associated graph. A continuous map $f \colon M \to N$ is piecewise smooth along $\mathcal {G}$ if f is smooth in restriction to any simplex of the mesh.

Note that there is a forgetful (restriction) map

$$ \begin{align*} \pi_{\mathcal{G}} \colon \mathcal{C}(M,N) \to \operatorname{\mathrm{Map}}_{\mathcal{G}} (M, N), \end{align*} $$

which assigns to any continuous map $f \colon M \to N$ its restriction to the vertex set of $\mathcal {G}$ . A first definition of an interpolation scheme would be a right inverse $\iota _{\mathcal {G}} $ of the map $\pi _{\mathcal {G}} $ .

Of course, a natural requirement to add is that $\iota _{\mathcal {G}} $ is a continuous map whose image is contained in the subspace of piecewise smooth maps along $\mathcal {G}$ . In the Euclidean setting, there is one canonical choice for interpolation, namely linear interpolation. In the general Riemannian setting, there is no such obvious choice. For our purposes, we will view center of mass interpolation as the preferred interpolation, although there are other natural options (e.g., harmonic interpolation), which we will not discuss.

There is a subtle deficiency in the above definition of interpolation scheme when N is not simply connected: one would like to require that $\iota _{\mathcal {G}} \circ \pi _{\mathcal {G}} $ preserves homotopy classes of maps, but that is not possible. This problem can be solved by defining an interpolation scheme as attached to the choice of a homotopy class:

Definition 2.1 Let $\mathcal {C}\subset \mathcal {C}(M,N)$ be a connected component. An interpolation scheme $\iota _{\mathcal {G}} $ is a continuous right inverse of $\pi _{\mathcal {G}} \big |_{\mathcal {C}}$ , whose image consists of piecewise smooth maps along $\mathcal {G}$ .

Note that this definition still does not allow one to define the homotopy class of a discrete map. A more elegant way to deal with deficiency, which we favored in [Reference Gaster, Loustau and Monsaingeon8], is to work equivariantly in the universal covers.

2.3.3 Center of mass interpolation

We refer to [Reference Gaster, Loustau and Monsaingeon8, Section 5.1] for generalities on centers of mass, also called barycenters, in metric spaces and Riemannian manifolds.

For comfort, let us assume that $S = M$ is two-dimensional; it is quite straightforward to generalize what follows to higher dimensions. First, we describe interpolation between triples of points. Let $A, B$ , and C be three points on the surface $(S,g)$ . We assume that these three points are sufficiently close, more precisely that they lie in a strongly convex geodesic ball B, i.e., any two points of B are joined by a unique minimal geodesic segment in S and this segment is contained in B. In particular, there is a uniquely defined triangle $T \subseteq S$ with vertices A, B, and C, and with geodesic boundary. Any point $P \in T$ can uniquely be written as the center of mass of $\{(A, \alpha ), (B, \beta ), (C, \gamma )\}$ , where $\alpha , \beta , \gamma \in [0,1]$ and $\alpha + \beta + \gamma = 1$ . Let similarly $A'$ , $B'$ , and $C'$ be three sufficiently close points in the Riemannian manifold $(N, h)$ . Then there is a unique center of mass interpolation map $f \colon A B C \to N$ such that for any point $P \in T$ as above, $f(P)$ is the center of mass of $\{(A', \alpha ), (B', \beta ), (C', \gamma )\}$ . In other words, f is the identity map in barycentric coordinates.

Clearly, given a discrete map $f \in \operatorname {\mathrm {Map}}_{\mathcal {G}} (S, N)$ , one can define its center of mass interpolation triangle by triangle following the procedure above. Although there seems to be a restriction on the size of the triangles in S and their images by f in N for the interpolation to be well-defined, one can work equivariantly in the universal covers as explained in Section 2.3.2 and the restriction disappears as long as S has nonpositive sectional curvature, or $\mathcal {G}$ is sufficiently fine, i.e., has small maximum edge length, and N has nonpositive sectional curvature.

Definition 2.3 Assume $(M,g)$ has nonpositive sectional curvature, or $\mathcal {G}$ is sufficiently fine, and N has nonpositive sectional curvature. The discussion above yields a center of mass interpolation scheme

$$ \begin{align*} \iota_{\mathcal{G}} \colon \operatorname{\mathrm{Map}}_{\mathcal{G}} (M, N) \to \mathcal{C}(M,N).\\[-20pt] \end{align*} $$

We denote $\widehat {f} := \iota _{\mathcal {G}} (f)$ the center of mass interpolation of a discrete map $f \in \operatorname {\mathrm {Map}}_{\mathcal {G}} (M, N)$ .

Theorem 2.3 Assume M has nonpositive sectional curvature, or $\mathcal {G}$ is sufficiently fine, and N has nonpositive sectional curvature. Then

1. (i) For any $f \in \operatorname {\mathrm {Map}}_{\mathcal {G}} (M, N)$ , the interpolation $\widehat {f}$ maps each edge of $\mathcal {G}$ to a geodesic segment in N (and does so with constant speed).

2. (ii) For any $f \in \operatorname {\mathrm {Map}}_{\mathcal {G}} (M, N)$ , the interpolation $\widehat {f}$ is piecewise smooth along $\mathcal {G}$ .

3. (iii) The map $\iota _{\mathcal {G}} \colon \operatorname {\mathrm {Map}}_{\mathcal {G}} (M, N) \to \mathcal {C}(M,N)$ is $1$ -Lipschitz for the $\operatorname {\mathrm {L}}^{\infty }$ distance on both spaces.

Proof For comfort, let us write the proof when $M = S$ is two-dimensional. The proof of (i) is immediate. For (ii), recall that the center of mass P as above is characterized by

$$ \begin{align*} \alpha \overrightarrow{PA} + \beta \overrightarrow{PB} + \gamma \overrightarrow{PC} = \vec{0} \end{align*} $$

(see [Reference Gaster, Loustau and Monsaingeon8, Equation (37)]), where we denote $\overrightarrow {PA} := \exp _{P}^{-1}(A)$ etc. It follows from the implicit function theorem that $(\alpha , \beta , \gamma )$ provide smooth barycentric coordinates on T (resp. $T'$ ). Conclude by observing that $\widehat {f}$ is the identity map in barycentric coordinates.

The proof of (iii) is a little more delicate, and crucially relies on N having nonpositive sectional curvature. Let $f_1, f_2 \in \operatorname {\mathrm {Map}}_{\mathcal {G}} (S, N)$ , we want to show that $d_{\infty }(\widehat {f_1}, \widehat {f_2}) \leqslant d_{\infty }(f_1, f_2)$ . Consider any triangle in $\mathcal {G}$ with vertices $A, B, C \in S$ . Let $p \in S$ be any point inside or on the boundary of the triangle $ABC \subseteq S$ . We denote $A_i = f_i(A)$ , $B_i = f_i(B)$ , $C_i = f_i(C)$ , $P_i = \widehat {f_i}(P)$ for $i\in \{1, 2\}$ . Since p is an arbitrary point on S, we win if we show that $d(P_1, P_2) \leqslant d_{\infty }(f_1, f_2)$ . By definition of the center of mass interpolation, $P_i$ is the center of mass of $\{(A_i, \alpha ), (B_i, \beta ), (C_i, \gamma )\}$ , where $\alpha , \beta , \gamma \in [0,1]$ is some triple with $\alpha + \beta + \gamma = 1$ (namely, the unique triple such that M is the center of mass of $\{(A, \alpha ), (B, \beta ), (C, \gamma )\}$ ). Let $\vec {V}_i = \alpha \overrightarrow {P_i A_i} + \beta \overrightarrow {P_i B_i} + \gamma \overrightarrow {P_i C_i}$ and let $\vec {W} = \alpha \overrightarrow {P_1 A_2} + \beta \overrightarrow {P_1 B_2} + \gamma \overrightarrow {P_1 C_2}$ , where we denote $\overrightarrow {P_i A_i} = \exp _{P_i}^{-1}(A_i)$ , etc. By definition of the center of mass $\vec {V}_i = \vec {0}$ , so we can write $\vec {W} = \vec {W} - \vec {V}_1$ :

(2.2) $$ \begin{align} \vec{W} = \alpha \left( \overrightarrow{P_1 A_2} - \overrightarrow{P_1 A_1} \right) + \beta \left( \overrightarrow{P_1 B_2} - \overrightarrow{P_1 B_1} \right) + \gamma \left( \overrightarrow{P_1 C_2} - \overrightarrow{P_1 C_1} \right). \end{align} $$

Since N has nonpositive sectional curvature, the exponential map $\exp _{P_1} \colon \operatorname {\mathrm {T}}_{P_1} N \to N$ is distance nondecreasing (for this argument to be completely rigorous, we may need to pass to universal covers), so that $\Vert \overrightarrow {P_1 A_2} - \overrightarrow {P_1 A_1} \Vert \leqslant d(A_1, A_2)$ , etc. Using the triangle inequality in (2.2), we find $\Vert \vec {W} \Vert \leqslant d_{\infty }(f_1, f_2)$ . This shows that $d(P_1, P_2) \leqslant d_{\infty }(f_1, f_2)$ by [Reference Gaster, Loustau and Monsaingeon8, Lemma 5.3].▪

3.1 Laplacian weights

Recall that we denote $\overrightarrow {xy} := \exp _x^{-1}y \in \operatorname {\mathrm {T}}_xM$ .

The following seemingly stronger characterization of Laplacian weights is immediate:

Note that we use the metric (inner product) in $\operatorname {\mathrm {T}}_x M$ to define $\operatorname {\mathrm {tr}} q$ . By definition, $\operatorname {\mathrm {tr}} q$ is the trace of the self-adjoint endomorphism associated to q.

3.2 Preferred vertex weights: the volume weights

In this paper, we favor one system of vertex weights associated to any mesh of any Riemannian manifold, the so-called volume weights.

For comfort assume $(M, g) = S$ is two-dimensional, although what follows is evidently generalized to higher dimensions. Let x be a vertex of the triangulation and consider the polygon $P_x \subseteq S$ equal to the union of the triangles adjacent to x. We define the weight of the vertex x by

$$ \begin{align*} \mu_x := \frac{1}{3} \operatorname{\mathrm{Area}}(P_x), \end{align*} $$

where $\operatorname {\mathrm {Area}}(P_x)$ denotes the Riemannian volume (area) of $P_x$ . This clearly defines a system of positive vertex weights $\mu _{\mathcal {G}} := (\mu _x)_{x \in \mathcal {V}}$ . We alternatively see $\mu _{\mathcal {G}} $ as a discrete measure on S supported by the set of vertices, which is meant to approximate the volume density $v_g$ of the Riemannian metric: see Section 4.2. Note that the choice of the constant $\frac {1}{1 +\dim M} = \frac {1}{3}$ in the definition of $\mu _x$ is motivated by the fact that each triangle is counted three times when integrating over S. The next proposition is almost trivial:

Proposition 3.5 Let $(M,g)$ be a closed manifold with an embedded graph $\mathcal {G}$ associated to a geodesic mesh. Let $\mu _{\mathcal {G}} $ be the discrete measure on S defined by the volume weights. Then

$$ \begin{align*} \sum_{x \in \mathcal{V}} \mu_x = \int_M \mathop{}\!\mathrm{d} \mu_{\mathcal{G}} = \int_M \mathop{}\!\mathrm{d} v_g = \operatorname{\mathrm{Vol}}(M,g). \end{align*} $$

Recall that any system of vertex weights endows the space of discrete maps $\operatorname {\mathrm {Map}}_{\mathcal {G}}(M,N)$ with an $\operatorname {\mathrm {L}}^2$ distance (see Section 2.1).

Proof Let us assume $M = S$ is two-dimensional for comfort. Let $f, g \in \operatorname {\mathrm {Map}}_{\mathcal {G}} (M, N)$ , denote by $\widehat {f} := \iota _{\mathcal {G}} (f)$ and $\widehat {g} := \iota _{\mathcal {G}} (g)$ their center of mass interpolations. By definition of the $\operatorname {\mathrm {L}}^2$ distance on $\mathcal {C}(M,N)$ ,

$$ \begin{align*} d(\widehat{f}, \widehat{g})^2 = \int_M d(\widehat{f}(x), \widehat{g}(x))^2 \, \mathop{}\!\mathrm{d} v_g(x)~. \end{align*} $$

Denote by $\mathcal {T}$ the set of triangles in the mesh. The integral is rewritten

(3.1) $$ \begin{align} d(\widehat{f}, \widehat{g})^2 = \sum_{T \in \mathcal{T}} \int_T d(\widehat{f}(x), \widehat{g}(x))^2 \, \mathop{}\!\mathrm{d} v_g(x)~. \end{align} $$

Let $T = ABC$ be any triangle in $\mathcal {T}$ . Following the proof of Theorem 2.3 (iii), for all $x \in T$ there exists $\alpha , \beta , \gamma \in [0,1]$ such that $\alpha + \beta + \gamma = 1$ and

$$ \begin{align*} d(\widehat{f}(x), \widehat{g}(x)) \leqslant \alpha d(f(A), g(A)) + \beta d(f(B), g(B)) + \gamma d(f(C), g(C)). \end{align*} $$

By convexity of the square function, it follows

(3.2) $$ \begin{align} d(\widehat{f}(x), \widehat{g}(x))^2 \leqslant \alpha d(f(A), g(A))^2 + \beta d(f(B), g(B))^2 + \gamma d(f(C), g(C))^2 \end{align} $$

hence

(3.3) $$ \begin{align} d(\widehat{f}(x), \widehat{g}(x))^2 \leqslant d(f(A), g(A))^2 + d(f(B), g(B))^2 + d(f(C), g(C))^2. \end{align} $$

Therefore, we may derive from (3.1)

$$ \begin{align*} \begin{aligned} d(\widehat{f}, \widehat{g})^2 &\leqslant \sum_{T \in \mathcal{T}} \left[d(f(A), g(A))^2 + d(f(B), g(B))^2 + d(f(C), g(C))^2\right] \operatorname{\mathrm{Area}}(T)\\[3pt] &\leqslant \sum_{x \in \mathcal{V}} \sum_{T \in \mathcal{T}_x} d(f(x), g(x))^2 \, \operatorname{\mathrm{Area}}(T_x), \end{aligned} \end{align*} $$

where $\mathcal {T}_x$ denotes the set of triangles adjacent to x. Finally, this is rewritten

$$ \begin{align*} d(\widehat{f}, \widehat{g})^2 \leqslant \sum_{x \in \mathcal{V}} 3 \mu_x \, d(f(x), g(x))^2, \end{align*} $$

where $\mu _x$ is the volume weight at x, i.e., $d(\widehat {f}, \widehat {g})^2 \leqslant 3 d(f, g)^2$ .

If M is Euclidean (flat), the proof can be upgraded to obtain a Lipschitz constant $L = 1$ by keeping the finer estimate (3.2) instead of (3.3), and computing the triangle integral.▪

3.3 Preferred edge weights: the cotangent weights

We also have a favorite system of edge weights, the so-called cotangent weights, although they have the following restrictions:

1. (1) We only define them for two-dimensional Riemannian manifolds, although they have natural higher-dimensional analogs.

2. (2) They are only positive for triangulations having the “Delaunay angle property.” (This includes any acute triangulation.)

These weights have a simple definition in terms of the cotangents of the (Riemannian) angles between edges in the triangulation, and coincide with the weights of Pinkall–Polthier [Reference Pinkall and Polthier13] in the Euclidean case. For more background on the cotangent weights in the Euclidean setting and a formula for their higher-dimensional analogs, please see [Reference Crane4].

The following result noticed by Pinkall–Polthier [Reference Pinkall and Polthier13] is an elementary exercise of plane Euclidean geometry:

Lemma 3.7 Let $T = ABC$ and $T' = A'B'C'$ be triangles in the Euclidean plane. Denote by $f \colon T \to T'$ the unique affine map such that $f(A) = A'$ , etc. Then the energy of f is given by

$$ \begin{align*} \begin{aligned} E(f) &:= \frac{1}{2}\int_T \Vert \mathop{}\!\mathrm{d} f \Vert^2 \mathop{}\!\mathrm{d} v\\[3pt] &= \frac{1}{4} \left(a^{\prime 2} \cot \alpha + b^{\prime 2} \cot \beta + c^{\prime 2} \cot \gamma \right), \end{aligned} \end{align*} $$

where $\alpha $ , $\beta $ , and $\gamma $ denote the unoriented angles of the triangle $ABC$ and $a'$ , $b'$ , and $c'$ denote the side lengths of the triangle $A'B'C'$ as in Figure 2.

Figure 2:

A triangle map in $\mathbb {R}^2$ .

In view of Lemma 3.7, given a surface $(S,g)$ equipped with a geodesic mesh, we define the weight of an edge e by considering the two angles $\alpha $ and $\beta $ opposite to e in the two triangles adjacent to e (see Figure 3), and we put

(3.4) $$ \begin{align} \omega_e := \frac{1}{2}(\cot \alpha + \cot \beta). \end{align} $$

Figure 3:

The weight $\omega _e$ of the edge e is defined in terms of the opposite angles $\alpha $ and $\beta $ .

Note that we use the Riemannian metric g to define the geodesic edges of the graph and the angles between edges.

As a direct application of Lemma 3.7, we obtain:

Note that a priori, the cotangent weights are not necessarily positive. Clearly, they are positive for acute triangulations (all of whose triangles are acute). More generally, the cotangent weights are positive if and only if the triangulation has the property that, for any edge e, the two opposite angles add to less than $\pi $ . This is simply because

$$ \begin{align*} \omega_e = \frac{1}{2}(\cot \alpha + \cot \beta) = \frac{\sin(\alpha + \beta)}{2 \sin \alpha \sin \beta}. \end{align*} $$

We call this the Delaunay angle property. In the Euclidean setting (for a flat surface), this property is equivalent to the triangulation being Delaunay, i.e., the circumcircle of any triangle does not contain any vertex in its interior [Reference Bobenko and Springborn2, Lemma 9 and Proposition 10].

3.4 Laplacian qualities of cotangent weights

In the two-dimensional Euclidean setting, in addition to Proposition 3.8, the cotangent weights enjoy some good—and other not so good—Laplacian properties, although this is much less obvious.

Proof Let x be a vertex and consider the polygon $P = P_x$ equal to the union of the triangles adjacent to x. Since in the flat case the exponential map $\exp _x$ is a local isometry, without loss of generality we can assume that P is contained in the Euclidean plane $\operatorname {\mathrm {T}}_x S \approx \mathbb {R}^2$ and $x = O$ .

Suppose that the vertices of P are given in cyclic order by $(A_i)$ , and that we have angles $\alpha _i$ , $\beta _i$ , and $\gamma _i$ as in Figure 4. By definition, the weight of the edge $O A_i$ is given by $\omega _i := \frac {1}{2}\left (\cot \beta _{i-1} + \cot \gamma _i\right )$ .

Figure 4:

The triangles of P at O.

Now consider the identity map $f \colon P \to \mathbb {R}^2$ . It has constant energy density $e(f) = 2$ , therefore, the total energy of f is $E = 2 \operatorname {\mathrm {Area}}(P)$ . On the other hand, E is the sum of the energies of f in restriction to the triangles forming P. By Lemma 3.7 this is

(3.5) $$ \begin{align} E = \frac{1}{4} \sum_i \left[\cot \alpha_i \Vert \overrightarrow{A_i A_{i+1}} \Vert^2 + \cot \beta_i \Vert \overrightarrow{O A_{i+1}} \Vert^2 + \cot \gamma_i \Vert \overrightarrow{ O A_i} \Vert^2 \right]. \end{align} $$

So far, we assumed that O is the origin in $\mathbb {R}^2$ , but of course the argument is valid if O is any point. In fact, let us see the energy E above as a function of $O \in \mathbb {R}^2$ when all the other points $A_i \in \mathbb {R}^2$ are fixed. We compute the infinitesimal variation of E under a variation O. On the one hand, $\dot {E}(O) = 0$ since $E(O) = 2 \operatorname {\mathrm {Area}}(P)$ is constant. On the other hand, (3.5) yields

(3.6) $$ \begin{align} \begin{aligned} \dot{E}(O) &= -\frac{1}{4} \sum_i \left[\frac{\dot{\alpha}_i}{\sin^2 \alpha_i} \Vert \overrightarrow{A_i A_{i+1}} \Vert^2 + \frac{\dot{\beta}_i}{\sin^2 \beta_i} \Vert \overrightarrow{O A_{i+1}} \Vert^2 + \frac{\dot{\gamma}_i}{\sin^2 \gamma_i} \Vert \overrightarrow{ O A_i} \Vert^2 \right] \\[3pt] &\quad - \frac{1}{2} \sum_i \left\langle \dot{O} \, , \, \cot \beta_i \overrightarrow{O A_{i+1}} + \cot \gamma_i \overrightarrow{O A_{i}} \right\rangle. \end{aligned} \end{align} $$

We claim that the first sum in (3.6) vanishes. Indeed, first observe that the law of sines yields

$$ \begin{align*} \frac{\Vert \overrightarrow{A_i A_{i+1}} \Vert^2 }{\sin^2 \alpha_i} = \frac{\Vert \overrightarrow{O A_{i+1}} \Vert^2}{\sin^2 \beta_i} = \frac{ \Vert \overrightarrow{ O A_i} \Vert^2}{\sin^2 \gamma_i} = D^2, \end{align*} $$

where D is the diameter of the triangle $O A_i A_{i+1}$ ’s circumcircle, so the first sum is rewritten

$$ \begin{align*} \begin{aligned} \sum_i \left[ D^2 \left(\dot{\alpha}_i + \dot{\beta}_i + \dot{\gamma}_i\right) \right], \end{aligned} \end{align*} $$

and $\dot {\alpha }_i + \dot {\beta }_i + \dot {\gamma }_i = 0,$ since $\alpha _i + \beta _i + \gamma _i = \pi $ is constant. Thus (3.6) is rewritten

$$ \begin{align*} \begin{aligned} \dot{E}(O) &= - \frac{1}{2} \sum_i \left\langle \dot{O} \, , \, \cot \beta_i \overrightarrow{O A_{i+1}} + \cot \gamma_i \overrightarrow{O A_{i}} \right\rangle\\[3pt] &=- \left\langle \dot{O} \, , \,\sum_i \omega_i \overrightarrow{O A_{i}} \right\rangle. \end{aligned} \end{align*} $$

In other words: $\operatorname {\mathrm {grad}} E(O) = - \sum _i \omega _i \overrightarrow {O A_{i}}$ . Since this must be zero (recall that $E(O)$ is constant), O is indeed the barycenter of its weighted neighbors $\{A_i, \omega _i\}$ .▪

It is not true in general that cotangent weights are Laplacian to second order. However, for triangulations obtained by midpoint refinement, it is true for almost all vertices:

The proof is based on the observation that any vertex of $\mathcal {G}_n$ is either an initial vertex (vertices of $\mathcal {G}_0$ ), a boundary vertex (vertices that are located on edges of the initial triangulation) or an interior vertex (all other vertices), and that the latter two satisfy a strong symmetry condition, which we call (semi-)hexaparallel symmetry:

It is straightforward to check by induction that a plane Euclidean graph obtained by iterated midpoint subdivision is hexaparallel at any interior vertex and semi-hexaparallel at any boundary vertex. Thus Proposition 3.10 reduces to:

Proof We need to show the second-order condition: for any quadratic form q on $\mathbb {R}^2$ ,

(3.7) $$ \begin{align} \frac{1}{\mu_x} \sum_{y\sim x} \omega_{xy} \ q(y-x) = 2 \ \operatorname{\mathrm{tr}} q. \end{align} $$

First, we argue that the semi-hexaparallel case derives from the hexaparallel case. Note that the left-hand side of (3.7) is invariant by the central symmetry at x, since a quadratic function is even. If x has semi-hexaparallel symmetry, we can create two hexaparallel configurations as in Figure 5b, both satisfying (3.7). Taking the half-sum of the two equations then yields the desired result.

Assume from now on that x has hexaparallel symmetry. Denote $y_1, \dots , y_6$ the neighbors in cyclic order. We may choose a complex coordinate on $\mathbb {R}^2 \approx \mathbb {C}$ so that $x=0$ and $y_1 = 1$ . Denote $z=a+bi$ the coordinate of $y_2$ . The hexaparallel condition implies that $y_3 = z-1$ , $y_4 = -1$ , $y_5 = -z$ , and $y_6 = 1-z$ . Let the oriented angles $\angle (y_1,y_2)$ , $\angle (y_2,y_3)$ , and $\angle (y_3,y_4)$ be denoted by $\alpha $ , $\beta $ , and $\gamma $ , respectively. For any $w \in \mathbb {C}$ , we have $\cot (\arg w) = \frac {\operatorname {Re}(w)}{\operatorname {Im}(w)}$ . Therefore we may compute:

$$ \begin{align*} \begin{aligned} \cot \alpha &= \frac{\operatorname{Re}}{\operatorname{Im}}(z)=\frac a b, \\[3pt] \cot \beta &= \frac{\operatorname{Re}}{\operatorname{Im}}\left(\frac{z-1}{z}\right) = \frac {a^2+b^2-a}{b}, \\[3pt] \cot \gamma &= \frac{\operatorname{Re}}{\operatorname{Im}}\left(\frac{1}{1-z}\right) = \frac{1-a}{b}. \end{aligned} \end{align*} $$

Since $\mu _x =\frac 13 ( 6\cdot b/2)=b$ , we get

$$ \begin{align*} \begin{aligned} \frac {1}{\mu_x} \sum_{y\sim x} \omega_{xy} \ q(y-x) &= \frac 2{b} \left( \cot \alpha \cdot q(z-1) + \cot \beta \cdot q(1) + \cot \gamma \cdot q(z) \right) \\[3pt] & = \frac 2b \left( \frac ab \cdot q(z-1) + \frac{a^2+b^2-a}{b} \cdot q(1) + \frac{1-a}{b} \cdot q(z) \right). \end{aligned} \end{align*} $$

The latter is equal to $2$ , $0$ , and $2$ when $q=\mathop {}\!\mathrm {d} x^2$ , $\mathop {}\!\mathrm {d} x\mathop {}\!\mathrm {d} y$ , or $\mathop {}\!\mathrm {d} y^2$ , respectively, as desired.▪

4.1 Fine and crystalline sequences of meshes

Let $(\mathcal {M}_n)_{n \in \mathbb {N}}$ be a sequence of geodesic meshes of a Riemannian manifold $(M,g)$ . Denote by $r_n$ the “mesh size,” i.e., the longest edge length of $\mathcal {M}_n$ . Following [Reference de Saint-Gervais5], we define:

Given a bounded subset $D \subseteq M$ , one calls:

• • diameter of D the supremum of the distance between two points of D, denoted $\operatorname {\mathrm {diam}}(D)$ .

• • radius of D the distance from the center of mass of D to its boundary, denoted $\operatorname {\mathrm {radius}}(D)$ .

• • thickness of D the ratio of its radius and diameter, denoted $\operatorname {\mathrm {thick}}(D)$ :

  $$ \begin{align*} \operatorname{\mathrm{thick}}(D) := \frac{\operatorname{\mathrm{radius}}(D)}{\operatorname{\mathrm{diam}}(D)}. \end{align*} $$

The surjectivity radius at a vertex x of a graph $\mathcal {G}$ is the smallest integer $k \in \mathbb {N}$ such that there exists a vertex at combinatorial distance k from x all of whose neighbors are at combinatorial distance $\leqslant k$ from x. The surjectivity radius of the graph $\mathcal {G}$ , denoted $\operatorname {\mathrm {surj}}\operatorname {\mathrm {rad}} \mathcal {G}$ , is the minimum of its surjectivity radii over all vertices.

For a continuous map $f \colon M \to \mathbb {R}$ , denote $f_n := \pi _n(f)\in \operatorname {\mathrm {Map}}_{\mathcal {G}_n}(M, N)$ the discretization of f: this is just the restriction of f to the vertex set of $\mathcal {G}_n$ . As in [Reference de Saint-Gervais5] we have:

The following lemma will be useful in Section 4.3 and again in Section 4.4.

Lemma 4.10 Assume that the sequence of meshes $(\mathcal {M}_n)_{n \in \mathbb {N}}$ on M is fine and crystalline. Let $\mathcal {G}_n$ be the graph underlying $\mathcal {M}_n$ and $r = r_n$ its maximum edge length. If $\mathcal {G}_n$ is equipped with a system of vertex and edge weights that is Laplacian at some vertex x, then

$$ \begin{align*} \frac{1}{\mu_x}\sum_{y\sim x}\omega_{xy} = O\left(r^{-2}\right). \end{align*} $$

Proof Apply condition (2) of Proposition 3.4 to the quadratic form $q = \Vert \cdot \Vert ^2$ :

(4.1) $$ \begin{align} \frac{1}{\mu_x}\sum_{y\sim x}\omega_{xy} \, d(x,y)^2 = 2m, \end{align} $$

where $m = \dim M$ . The fact that the sequence of meshes is fine and crystalline implies that there exists a uniform lower bound for the ratio of lengths in the triangulation. Thus there exists a constant $\alpha> 0$ such that for any neighbor vertices x and y in $\mathcal {G}_n$ :

(4.2) $$ \begin{align} \alpha \, r \leqslant d(x,y) \leqslant r. \end{align} $$

It follows from (4.1) and (4.2) that

$$ \begin{align*} \frac{1}{\mu_x}\sum_{y\sim x}\omega_{xy} \leqslant \frac{2m}{\alpha^2 r^2}.\\[-40pt] \end{align*} $$

▪

4.2 Convergence of the volume form

Let $(M,g)$ be a Riemannian manifold, let $(\mathcal {M}_n)_{n \in \mathbb {N}}$ be a sequence of meshes with the underlying graphs $(\mathcal {G}_n)_{n \in \mathbb {N}}$ . We equip $\mathcal {G}_n$ with the volume vertex weights defined in Section 3.2. These define a discrete measure $\mu _n$ on M supported by the set of vertices $\mathcal {V}_n = \mathcal {G}_n^{(0)}$ .

Theorem 4.13 If M is any Riemannian manifold and $(\mathcal {M}_n)_{n \in \mathbb {N}}$ is any fine sequence of meshes, then the measures $(\mu _n)_{n\in \mathbb {N}}$ on M defined by the volume vertex weights converge weakly-* to the volume density on M:

$$ \begin{align*} \int_M f \, \mathop{}\!\mathrm{d} \mu_n \stackrel{n\to +\infty}{\longrightarrow} \int_M f \, \mathop{}\!\mathrm{d} \mu \end{align*} $$

for any $f \in \mathcal {C}_c^0(M, \mathbb {R})$ (continuous function with compact support), where $\mu $ denotes the measure on M induced by the volume form $v_g$ .

Proof Recall that a continuity set $A \subseteq M$ is a Borel set such that $\mu (\partial A) = 0$ . Since any compact set has finite $\mu $ -measure, it is well-known that the weakly-* convergence of $\mu _n$ to $\mu $ is equivalent to

$$ \begin{align*} \mu_n(A) \stackrel{n\to +\infty}{\longrightarrow} \mu(A) \end{align*} $$

for any bounded continuity set A. Let thus A be any bounded continuity set. Denote by $B_n$ the union of all simplices that are entirely contained in A, and by $C_n$ the union of all simplices that have at least one vertex in A. We obviously have $B_n \subseteq A \subseteq C_n$ , and by definition of $\mu _n$ we have:

(4.3) $$ \begin{align} \mu(B_n) \leqslant \mu_n(A) \leqslant \mu(C_n). \end{align} $$

On the other hand, clearly we have $C_n - B_n \subseteq N_{\varepsilon _n}(\partial A)$ , where we have denoted $ N_{\varepsilon _n}(\partial A)$ the $\varepsilon _n$ -neighborhood of $\partial A$ , with $\varepsilon _n = 2r$ here. (As usual we denote $r = r_n$ the maximal edge length in $\mathcal {M}_n$ .) By continuity of the measure $\mu $ , we know that $\lim _{n\to +\infty } \mu (N_{\varepsilon _n}(\partial A)) = \mu (\partial A) = 0$ . Note that we used the boundedness of A, which guarantees that $\mu (N_{\varepsilon _n}(\partial A)) < +\infty $ . It follows:

(4.4) $$ \begin{align} \lim_{n \to +\infty} \mu(C_n - B_n) = 0~. \end{align} $$

Since $B_n \subseteq A \subseteq C_n$ , (4.4) implies that $\lim _{n \to +\infty } \mu (B_n) = \lim _{n \to +\infty } \mu (C_n) = \mu (A)$ , and we conclude with (4.3) that $\lim _{n \to +\infty } \mu _n(A) = \mu (A)$ .▪

4.3 Convergence of the tension field

Now we consider another Riemannian manifold N and a smooth function $f \colon M \to N$ .

Consider a fine and crystalline sequence of meshes $(\mathcal {M}_n)_{n\in \mathbb {N}}$ on M, with mesh size (i.e., maximum edge length) $r = r_n$ , and underlying graph $\mathcal {G}_n$ .

Theorem 4.14 Assume that the sequence of meshes $(\mathcal {M}_n)_{n \in \mathbb {N}}$ on M is fine and crystalline. If $\mathcal {G}_n$ is equipped with a system of vertex and edge weights that is Laplacian at some vertex x, then

(4.5) $$ \begin{align} \tau_{\mathcal{G}_n}(f_n)_x - \tau(f)_x = O\left(r^2\right). \end{align} $$

Proof Consider $F:=\exp _{f(x)}^{-1} \circ f \circ \exp _x : \operatorname {\mathrm {T}}_x M \to \operatorname {\mathrm {T}}_{f(x)}N$ . For $y \sim x$ , denote $v = v_y := \exp _x^{-1}y$ . By Taylor’s theorem, we have

(4.6) $$ \begin{align} \exp_{f(x)}^{-1}f(y) = F(v) = (\mathop{}\!\mathrm{d} F)_{|0}(v) + \frac{1}{2} (\mathop{}\!\mathrm{d}^2 F)_{|0}(v,v) + \frac 16 (\mathop{}\!\mathrm{d}^3 F)_{|0}(v,v,v) + O\left(r^4\right). \end{align} $$

This implies

$$ \begin{align*} \begin{aligned} \tau_{\mathcal{G}} (f)(x) & = \frac 1{\mu_x} \sum_{y\sim x} \omega_{xy} \, \exp_{f(x)}^{-1}f(y) \\[3pt] & = \frac 1{\mu_x} \sum_{y\sim x}\omega_{xy} \, (\mathop{}\!\mathrm{d} F)_{|0}(v) + \frac 1{2\mu_x} \sum_{y\sim x}\omega_{xy}\, (\mathop{}\!\mathrm{d}^2 F)_{|0}(v,v)\\[3pt] & \quad + \frac 1{6 \mu_x} \sum_{y\sim x}\omega_{xy} \, (\mathop{}\!\mathrm{d}^3 F)_{|0}(v,v,v)+\frac{1}{\mu_x}\sum_{y\sim x}\omega_{xy} \, O\left(r^4\right). \end{aligned} \end{align*} $$

By conditions (1) and (3) of Proposition 3.4, the first and third sums above vanish, while the second sum is rewritten with condition (2):

$$ \begin{align*} \tau_{\mathcal{G}} (f_{\mathcal{G}} )(x) = \operatorname{\mathrm{tr}}\left(\mathop{}\!\mathrm{d}^2F_{|0}\right) + \frac{1}{\mu_x}\sum_{y\sim x}\omega_{xy} \, O\left(r^4\right). \end{align*} $$

Note that $\operatorname {\mathrm {tr}}\left (\mathop {}\!\mathrm {d}^2F_{|0}\right ) = \operatorname {\mathrm {tr}}\left (\nabla ^2f_{|x}\right ) =\tau (f)(x)$ , and conclude with Lemma 4.10.▪

4.4 Convergence of the energy

We keep the setting of Section 4.3: $f \colon M \to N$ is a smooth function between Riemannian manifolds, and M is equipped with a sequence of meshes $(\mathcal {M}_n)_{n \in \mathbb {N}}$ that is fine and crystalline.

4.4.1 Convergence of the energy density

Theorem 4.17 Assume that the sequence of meshes $(\mathcal {M}_n)_{n \in \mathbb {N}}$ on M is fine and crystalline. Assume $\mathcal {G}_n$ is equipped with a system of vertex and edge weights. Then

$$ \begin{align*} e_{\mathcal{G}_n}(f_n) = e(f) + O\left(r^2\right) \end{align*} $$

on the set of vertices where $\mathcal {G}_n$ is Laplacian.

Recall that we denote $f_n := \pi _n(f)$ the discretization of f along $\mathcal {G}_n$ .

Remark 4.18 Remark 4.16 holds again for Theorem 4.17.

Proof Assume $\mathcal {G}_n$ is Laplacian at x. Using (4.6) again, denoting $v_y = \exp _x^{-1}y$ , we find that

$$ \begin{align*} e_{\mathcal{G}} (f)_x & = \frac1{4\mu_x} \sum_{y\sim x} \omega_{xy} \, \Vert F(v_y) \Vert ^2 \\[3pt] & = \frac1{4 \mu_x} \sum_{y\sim x} \omega_{xy} \, \left\Vert ( \mathop{}\!\mathrm{d} F)_{|0}(v_y) + \frac12 (\mathop{}\!\mathrm{d}^2 F)_{|0}(v_y) + O\left( r^3\right) \right\Vert{}^2 \nonumber \\[-14pt] \nonumber \\ & = \frac1{4\mu_x} \sum_{y\sim x} \omega_{xy} \,\Vert (\mathop{}\!\mathrm{d} F)_{|0}(v_y) \Vert^2 + \frac1{4\mu_x} \sum_{y\sim x} \omega_{xy}\, \langle (\mathop{}\!\mathrm{d} F)_{|0}(v_y) , (\mathop{}\!\mathrm{d}^2 F)_{|0}(v_y) \rangle \\[3pt] & \quad + \frac1{4\mu_x} \sum_{y\sim x} \omega_{xy} \,O\left(r^4 \right). \end{align*} $$

Condition (3) of Proposition 3.4 implies that the second sum vanishes. Lemma 4.10 implies that the third sum is $O\left (r^2 \right )$ . By condition (2) of Proposition 3.4, the remaining first sum is rewritten

$$ \begin{align*} \frac1{4\mu_x} \sum_{y\sim x} \omega_{xy} \,\Vert (\mathop{}\!\mathrm{d} F)_{|0}(v_y) \Vert^2 = \frac12 \Vert (\mathop{}\!\mathrm{d} F)_{|0}\Vert^2 = e(f)_x, \end{align*} $$

since $\operatorname {\mathrm {tr}}(L^2) = \Vert L \Vert ^2$ for any linear form L. We thus get

$$ \begin{align*} e_{\mathcal{G}} (f)_x = e(f)_x + O\left(r^2\right).\\[-40pt] \end{align*} $$

▪

4.4.2 Convergence of the energy

Recall that the energy is $E(f) := \int _M e(f) \mathop {}\!\mathrm {d} \mu $ . The convergence of the discrete energy is now an easy consequence of the weakly-* convergence of measures $\mu _n \to \mu $ and the uniform convergence of the energy densities $e_{\mathcal {G}_n}(f_n) \to e(f)$ . This is the classical combination of weak convergence and strong convergence.

Definition 4.3 Let $(M,g)$ be a Riemannian manifold. Consider a sequence of geodesic meshes $(\mathcal {M}_n)_{n\in \mathbb {N}}$ , and equip the underlying graphs $\mathcal {G}_n$ with a system of positive vertex and edge weights. We call the sequence of biweighted graphs $(\mathcal {G}_n)_{n\in \mathbb {N}}$ Laplacian provided that:

1. (i) The sequence of meshes $(\mathcal {M}_n)_{n\in \mathbb {N}}$ is fine and crystalline.

2. (ii) For every $n\in \mathbb {N}$ , the vertex weights on $\mathcal {G}_n$ are given by the volume weights (see Section 3.2).

3. (iii) For every $n\in \mathbb {N}$ , the system of vertex and edge weights on $\mathcal {G}_n$ is Laplacian.

Theorem 4.19 Let M be a Riemannian manifold and let $(\mathcal {M}_n)_{n\in \mathbb {N}}$ be a Laplacian sequence of meshes. For any smooth $f \colon M \to N$ with compact support:

$$ \begin{align*} \lim_{n \to +\infty} E_{\mathcal{G}_n}(f_n) = E(f). \end{align*} $$

Recall that we denote $f_n := \pi _{\mathcal {G}_n}(f)$ the discretization of f along $\mathcal {G}_n$ .

Proof By Theorem 4.13,

$$ \begin{align*} E(f) = \lim_{n \to + \infty} \int_M e(f) \mathop{}\!\mathrm{d} \mu_n. \end{align*} $$

By Theorem 4.17, on the support of $\mu _n$ , $e(f) = e_{\mathcal {G}_n}(f_n) + O\left (r^2\right )$ . It follows that

$$ \begin{align*} E(f) = \lim_{n \to + \infty} \int_M e_{\mathcal{G}_n}(f_n) \mathop{}\!\mathrm{d} \mu_n, \end{align*} $$

in other words, $E(f) = \lim _{n \to +\infty } E_{\mathcal {G}_n}(f_n)$ .▪

Remark 4.20 The proof of Theorem 4.19 hints that $E(f) = E_{\mathcal {G}_n}(f_n) + O\left (r^2\right )$ , provided that the convergence of $\mu _n$ to $\mu $ is sufficiently fast. Improvements of this estimate can occur in more restricted situations: for instance, when both the target and the domain are hyperbolic surfaces:

$$ \begin{align*} E(f) = E_{\mathcal{G}_n}(f_n) + O\left(r^4\right). \end{align*} $$

This can be proven by carrying out involved calculations in the hyperbolic plane, which we spare.

4.5 Weak Laplacian conditions

It is clear from the proofs of the main results in the previous subsections that the Laplacian conditions for sequences of meshes can be weakened and still produce the same results, or at least some of them, with minimal changes in the proofs. This is a useful generalization, for it is very stringent to require a sequence of weighted graphs $(\mathcal {G}_n)$ to be Laplacian for all n. Instead we start by asking that the sequence is merely asymptotically Laplacian in the following sense.

Definition 4.4 Let M be a Riemannian manifold. Consider a sequence of geodesic meshes $(\mathcal {M}_n)_{n\in \mathbb {N}}$ , and equip the underlying graphs $\mathcal {G}_n$ with a system of positive vertex weights $\{\mu _{x, n}\}$ . We call the sequence of weight systems $(\{\mu _{x, n}\})_{n \in \mathbb {N}}$ asymptotic volume weights provided that

$$ \begin{align*} \mu_{x, n} = (1 + o(1)) \, \widehat{\mu}_{x, n} \end{align*} $$

for some function $o(1)$ independent of x, where $\widehat {\mu }_{x, n}$ denote the volume weights (see Section 3.2).

The following proposition is an immediate consequence of Theorem 4.13:

It is immediate to show that for asymptotic volume weight, Theorem 3.6 holds with a Lipschitz constant $L_n = \sqrt {1+\dim M} + o(1)$ . Although this is sufficient for the needs of this paper (see Lemma 5.5), let us state in the next theorem that the result can be improved to $L_n = 1 + o(1)$ . The proof follows from Theorem 3.6 by writing an expansion of the volume form in normal coordinates, we skip it for brevity.

Explicitly, (iii) means that for all $x \in \mathcal {V}_n$ and $L \in \operatorname {\mathrm {T}}_x^* M$ :

1. (1)

   $$ \begin{align*} \frac{1}{\mu_x} \sum_{y \sim x} \omega_{xy} \, \overrightarrow{xy} = O\left(r^2\right), \end{align*} $$

2. (2)

   $$ \begin{align*} \frac{1}{\mu_x} \sum_{y\sim x} \omega_{xy} \, L(\overrightarrow{xy})^2 = 2\| L \|^2 \left(1 + O\left(r^2\right)\right), \end{align*} $$

3. (3)

   $$ \begin{align*} \frac{1}{\mu_x} \sum_{y \sim x} \omega_{xy} \, L(\overrightarrow{xy}) ^3 = \| L \|^3 O\left(r^2\right). \end{align*} $$

The $O(r^2)$ functions above should be independent of x and L. Note again that to alleviate notations, we drop the dependence in n when writing r, $\mu _x$ , and $\omega _{xy}$ .

It is immediate to check that the proofs of Theorems 4.14, 4.17, and 4.19 apply to asymptotically Laplacian sequences of graphs. Alas, it is still unreasonable to expect to be able to construct asymptotically Laplacian sequences in general. Fortunately, the notion may be further slightly weakened while keeping the validity of the most important theorems, and allowing the systematic construction of such sequences in Section 6 (at least in the two-dimensional case).

The following theorems are generalized or weakened versions of Theorems 4.14, 4.17, and 4.19.

Note that we use the discrete measure $\mu _n$ on the vertex set of $\mathcal {G}_n$ in order to define the $\operatorname {\mathrm {L}}^2$ -norm on spaces of discrete maps along $\mathcal {G}_n$ .

Proof When $(\mathcal {G}_n)_{n\in \mathbb {N}}$ is Laplacian, (1) is an immediate consequence of Theorem 4.14. When $(\mathcal {G}_n)_{n\in \mathbb {N}}$ is merely asymptotically Laplacian, the proof of Theorem 4.14 is still valid up to $O\left (r^2 \right )$ .

For the proof of (2), let $\mathcal {V}_n^{(k)}$ be the subset of $\mathcal {V}_n$ of mass $O(r^{k})$ where $\mathcal {G}_n$ is Laplacian up to $O\left (r^{2-k}\right )$ . By tracing the proof of Theorem 4.14, one quickly sees that $\tau (f) = \tau _{\mathcal {G}_n}(f_n) + O\left (r^{2-k}\right )$ on $\mathcal {V}_n^{(k)}$ , for each $k \in \{0, 1, 2\}$ . The decomposition $\mathcal {V}_n = \bigsqcup _{k=0}^2 \mathcal {V}_n^{(k)}$ implies

$$ \begin{align*} \begin{aligned} \Vert \tau(f) - \tau_{\mathcal{G}_n}(f_n ) \Vert^2 & = \sum_{k=0}^2 \Vert \tau(f) - \tau_{\mathcal{G}_n}(f_n ) \Vert_{\mathcal{V}_n^{(k)}}^2 \\[3pt] &\leqslant \sum_{k=0}^2 \Vert \tau(f) - \tau_{\mathcal{G}_n}(f_n ) \Vert_{\infty, \mathcal{V}_n^{(k)}}^2 ~\mu(\mathcal{V}_n^{(k)}) \\[3pt] &\leqslant \sum_{k=0}^2 O\left(r^{4-2k}\right) O\left(r^{k}\right) = O(r^2). \end{aligned} \end{align*} $$

For the second estimate, write similarly

$$ \begin{align*} \begin{aligned} \left\langle \tau(f) - \tau_{\mathcal{G}_n}(f_n ) \, , \, \vec{V} \right\rangle & = \sum_{k=0}^2 \left\langle \tau(f) - \tau_{\mathcal{G}_n}(f_n ) \, , \, \vec{V} \right\rangle_{\mathcal{V}_n^{(k)}} \\[3pt] &\leqslant \sum_{k=0}^2 \Vert \tau(f) - \tau_{\mathcal{G}_n}(f_n ) \Vert_{\mathcal{V}_n^{(k)}} ~ \Vert \vec{V} \Vert_{\mathcal{V}_n^{(k)}} \\[3pt] &\leqslant O\left(r^{2}\right) \cdot 1 + O\left(r^{3/2}\right) \cdot 1 + O\left(r\right) \cdot o(1) = o(r). \end{aligned} \\[-38pt] \end{align*} $$

▪

Theorem 4.27 We keep the setup of Theorem 4.25. If $(\mathcal {G}_n)_{n\in \mathbb {N}}$ is almost asymptotically Laplacian,

$$ \begin{align*} \lim_{n \to +\infty} E_{\mathcal{G}_n}(f_n) = E(f)~. \end{align*} $$

Proof By definition of almost asymptotically Laplacian, the sequence of measures $(\mu _n)_{n\in \mathbb {N}}$ converges weakly-* to the measure $\mu $ on M, therefore

(4.9) $$ \begin{align} \begin{aligned} E(f) = \int_M e(f) \mathop{}\!\mathrm{d} \mu = \lim_{n \to +\infty} \int_M e(f) \mathop{}\!\mathrm{d} \mu_n. \end{aligned} \end{align} $$

Let $\mathcal {V}_n = \bigsqcup _{k=0}^2 \mathcal {V}_n^{(k)}$ be the decomposition of the vertices of $\mathcal {G}_n$ granted by Definition 4.6. By Theorem 4.26,

$$ \begin{align*} \begin{aligned} \int_M e(f) \mathop{}\!\mathrm{d} \mu_n = \sum_{k=0}^2 \int_{ \mathcal{V}_n^{(k)} } e(f) \mathop{}\!\mathrm{d} \mu_n = \sum_{k=0}^2 \int_{\mathcal{V}_n^{(k)}} e_{\mathcal{G}_n}(f_n) + O\left( r^{2-k} \right) \mathop{}\!\mathrm{d} \mu_n . \end{aligned} \end{align*} $$

It follows:

$$ \begin{align*} \begin{aligned} \int_M e(f) \mathop{}\!\mathrm{d} \mu_n &= \int_M e_{\mathcal{G}_n}(f_n) \mathop{}\!\mathrm{d} \mu_n + \sum_{k=0}^2 O\left(r^{k}\right) \, O\left(r^{2-k}\right) = E_{\mathcal{G}_n}(f_n) + O\left(r^2\right) . \end{aligned} \end{align*} $$

In particular, we find that $\int _M e(f) \mathop {}\!\mathrm {d} \mu _n = E_{\mathcal {G}_n}(f_n) + o\left (1\right )$ . Injecting this into (4.9) yields the desired result $E(f) = \lim _{n \to +\infty } E_{\mathcal {G}_n}(f_n)$ .▪

5.1 Strong convexity of the discrete energy

Please refer to [Reference Gaster, Loustau and Monsaingeon8, Section 3.1] for the definition of convex, strictly convex, and strongly convex functions on Riemannian manifolds. In a nutshell, these notions are generalized from the one-dimensional case by restricting to geodesics; the convexity [resp. $\alpha $ -strong convexity] of a smooth function is characterized by its Hessian being $\geqslant 0$ [resp. $\geqslant \alpha g$ where g is the Riemannian metric].

Keeping the same setup as above, assume moreover that N is compact and has negative sectional curvature. In this case, we expect that the discrete energy functional $E_{\mathcal {G}} \colon \mathcal {C}_{\mathcal {G}} \to \mathbb {R}$ is $\alpha _{\mathcal {G}} $ -strongly convex for any biweighted graph $\mathcal {G}$ on M underlying a mesh, for some $\alpha _{\mathcal {G}}>0$ . In our previous paper, we proved this statement when M and N are two-dimensional. The estimates we obtained (see [Reference Gaster, Loustau and Monsaingeon8, Theorem 3.20 and Proposition 3.14]) imply that, when $\mathcal {G}$ is equipped with volume weights, $\alpha _{\mathcal {G}} = \Omega \left ( \operatorname {\mathrm {diam}}(\mathcal {G})^{-1} \right )$ . Further, when $(\mathcal {M}_n)_{n \in \mathbb {N}}$ is a fine and crystalline sequence of meshes of M and mesh sizes $r=r_n$ , with underlying graphs $\mathcal {G}_n$ , discrete energy functionals $E_n:=E_{\mathcal {G}_n}$ , and moduli of convexity $\alpha _n:=\alpha _{\mathcal {G}_n}$ , Theorem 4.4 implies that we have the estimate $\alpha _n = \Omega \left ( r \right )$ .

In fact, we conjecture that the smooth energy $E \colon \mathcal {C} \to \mathbb {R}$ is $\alpha $ -strongly convex for some $\alpha> 0$ (see [Reference Gaster, Loustau and Monsaingeon8, Section 3.2] for a discussion), and we expect that $\alpha = \lim _{n \to +\infty }{\alpha _n}$ for any asymptotically Laplacian sequence of meshes $(\mathcal {M}_n)_{n \in \mathbb {N}}$ . In particular, the sequence $(\alpha _n)_{n \in \mathbb {N}}$ should be $\Omega (1)$ in great generality (see Notation 4.5 for the notations $\Omega $ and $\Theta $ ).

5.2 $\operatorname {\mathrm {L}}^2$ convergence

The main theorem of this section is:

Theorem 5.1 Let M and N be any Riemannian manifolds with M compact and N complete with nonpositive sectional curvature. Let $\mathcal {C}$ be a connected component of $\mathcal {C}^{\infty }(M,N)$ containing a harmonic map w. Consider a fine and crystalline sequence of meshes $(\mathcal {M}_n)_{n\in \mathbb {N}}$ of M with mesh size $r=r_n$ and underlying graphs $(\mathcal {G}_n)_{n\in \mathbb {N}}$ that satisfy:

1. (i) The sequence $(\mathcal {G}_n)_{n\in \mathbb {N}}$ is almost asymptotically Laplacian.

2. (ii) The discrete energy $E_n \colon \operatorname {\mathrm {Map}}_{\mathcal {G}_n}(M,N) \to \mathbb {R}$ is $\alpha _n$ -strongly convex on $\mathcal {C}_n$ , with $\alpha _n = \Omega \left ( r^c \right )$ .

Denote $v_n \in \operatorname {\mathrm {Map}}_{\mathcal {G}_n}(M,N)$ , the minimizer of $E_n$ on $\mathcal {C}_n$ and $\widehat {v_n}$ its center of mass interpolation.

If $c<1$ , then

$$ \begin{align*} \widehat{v_n} \xrightarrow[n \to +\infty]{} w \quad \text{in } \operatorname{\mathrm{L}}^2(M, N). \end{align*} $$

Moreover, the conclusion still holds if $c=1$ and $\dim M=2$ , assuming $(\mathcal {G}_n)_{n\in \mathbb {N}}$ has uniformly bounded ratio between edge weights.

Proof The proof is a combination of a few key ideas that we emphasize using in-proof lemmas. The bulk of the hard work has been done in the previous sections, which we will refer to for the proof of these lemmas.

Let $w_n := \pi _n(w)\in \operatorname {\mathrm {Map}}(\mathcal {G}_n, N)$ denote the discretization of w (restriction of w to the vertex set of $\mathcal {G}_n$ ). We also denote $\widehat {w_n}$ the center of mass interpolation of $w_n$ .

Lemma 5.4 We have $\widehat {w_n} \to w$ in $\operatorname {\mathrm {L}}^2(M,N)$ when $n \to +\infty $ , moreover $E(\widehat {w_n}) \to E(w)$ .▪

Proof This is an immediate consequence of Corollary 4.8, which we can invoke since M is compact and the sequence of meshes $(\mathcal {M}_n)_{n\in \mathbb {N}}$ is fine and crystalline.▪

Lemma 5.5 There exists a constant $L>0$ such that

$$ \begin{align*} d(\widehat{w_n}, \widehat{v_n}) \leqslant L\, d(w_n, v_n), \end{align*} $$

where $d(\widehat {w_n}, \widehat {v_n})$ and $d(w_n, v_n)$ indicate the $\operatorname {\mathrm {L}}^2$ distances in $\mathcal {C}(M, N)$ and $\operatorname {\mathrm {Map}}_{\mathcal {G}_n}(M, N)$ .

Proof This follows immediately from Theorem 4.22.▪

Lemma 5.6 Let R be a complete Riemannian manifold and $F \colon R \to \mathbb {R}$ be a $\mathcal {C}^2 \ \alpha $ -strongly convex function. Then F has a unique minimizer $x^*$ , and for all $x\in R$

(5.1) $$ \begin{align} d(x,x^*) \leqslant \frac{ \left| \left \langle \; \operatorname{\mathrm{grad}} F(x) \;, \vec{V} \; \right \rangle \right|}{\alpha}, \end{align} $$

where $\vec {V} $ is a unit tangent vector in the direction $\exp _x^{-1}(x^*)$ , in particular

(5.2) $$ \begin{align} d(x,x^*) \leqslant \frac{\Vert \operatorname{\mathrm{grad}} F(x) \Vert}{\alpha}. \end{align} $$

We also have

(5.3) $$ \begin{align} 0 \leqslant F(x) - F(x^*) \leqslant \frac{\left\Vert\operatorname{\mathrm{grad}} F(x)\right\Vert{}^2}{\alpha}. \end{align} $$

Proof Recall that on a complete Riemannian manifold, there exists a length-minimizing geodesic between any two points. It is not hard to show that a strongly convex function on a complete (finite-dimensional) Riemannian manifold is proper, hence existence of the minimizer, and uniqueness follows from strict convexity.

The first inequality (5.1) is easy to prove for a function $f \colon \mathbb {R} \to \mathbb {R}$ by integrating $f''(x) \geqslant \alpha $ . For the general case, take a length-minimizing unit geodesic $\gamma \colon \mathbb {R} \to R$ with $\gamma (0) = x^*$ and $\gamma (L) = x$ , and apply the previous result to $f = F \circ \gamma $ . The second inequality (5.3) follows with Cauchy–Schwarz. For (5.3), the one-dimensional case is readily obtained via the mean value theorem, and the general case quickly follows.▪

Lemma 5.7 We have

(5.4) $$ \begin{align} d(w_n, v_n) \leqslant \frac{ \left| \left \langle \; \tau_{\mathcal{G}_n}(w_n) \; , \vec{V} \; \right \rangle \right| }{\alpha_n} \quad \text{ where } \quad \vec{V} = \frac{\exp_{w_n}^{-1}v_n}{\Vert \exp_{w_n}^{-1}v_n \Vert} , \end{align} $$

where d denotes the $\operatorname {\mathrm {L}}^2$ distance in $\operatorname {\mathrm {Map}}_{\mathcal {G}_n}(M, N)$ . In particular,

(5.5) $$ \begin{align} d(w_n, v_n) \leqslant \frac{ \Vert \tau_{\mathcal{G}_n}(w_n) \Vert }{\alpha_n}. \end{align} $$

Proof Apply Lemma 5.6 (5.1) and (5.2) to $R = \operatorname {\mathrm {Map}}_{\mathcal {G}_n}(M,N)$ and $F = E_n$ .▪

At this point, we would like to apply Lemma 5.7 and Theorem 4.25 to conclude that

$$ \begin{align*} d(w_n, v_n) \to 0. \end{align*} $$

Indeed, (5.5) together with (4.7) imply that $d(w_n, v_n) = O(r^{1-c})$ . If $c<1$ , we thus clearly have $d(w_n, v_n) \to 0$ . The equality case $c=1$ is much more subtle. In theory, we can still conclude that $d(w_n, v_n) \to 0$ with (5.4) and (4.8), which together yield $d(w_n, v_n) = o(1)$ . However, to apply (4.8), we need to know that $\Vert \vec {V} \Vert _{\mathcal {V}_n^{(2)}} = o(1)$ . Although we believe this is always true, we only show it when $\dim M = 2$ in this paper.

Lemma 5.8 Assume $\dim M = 2$ . We have $\Vert \vec {V} \Vert _{\mathcal {V}_n^{(2)}} = o(1)$ .

Proof Clearly, $\Vert \vec {V} \Vert _{\mathcal {V}_n^{(2)}}^2 \leqslant \Vert \vec {V} \Vert _{\infty }^2 ~ \mu (\mathcal {V}_n^{(2)})$ , that is

$$ \begin{align*} \Vert \vec{V} \Vert_{\mathcal{V}_n^{(2)}}^2 \leqslant \frac{d_{\infty}(w_n, v_n)^2}{d(w_n, v_n)^2} ~ O(r^2). \end{align*} $$

It appears that we win if we can show that $\frac {d_{\infty }(w_n, v_n)}{d(w_n, v_n)} = o(r^{-1})$ . Unfortunately, the comparison between the $\operatorname {\mathrm {L}}^{\infty }$ distance and the $\operatorname {\mathrm {L}}^2$ distance on $\operatorname {\mathrm {Map}}_{\mathcal {G}_n}(M, N)$ only satisfies $\frac {d_{\infty }(u, v)}{d(u, v)} = O(r^{-1})$ in general. However, this inequality may be slightly improved when v is the discrete energy minimizer. In order to avoid burdening our exposition, we relegate this technical estimate to Appendix B. The desired comparison is given in Corollary B.4 (which requires the uniform bound assumption on ratios of edge weights).▪

We can now smoothly wrap up the proof of Theorem 5.1: write

$$ \begin{align*} d(\widehat{v_n}, w) & \leqslant d(\widehat{v_n}, \widehat{w_n}) + d(\widehat{w_n}, w), && \text{(triangle inequality)}\\[3pt] &\leqslant L\, d(w_n, v_n) + o(1). && \text{(by Lemmas 5.4 and 5.5)} \end{align*} $$

We proved that $d(w_n, v_n) \to 0$ if $c<1$ or $c=1$ and $\dim M = 2$ , so we are done.

5.3 $\operatorname {\mathrm {L}}^{\infty }$ convergence

Under stronger assumptions, we are able to prove uniform convergence in the two-dimensional case by comparing the $\operatorname {\mathrm {L}}^2$ and $\operatorname {\mathrm {L}}^{\infty }$ distances on the space of discrete maps $\operatorname {\mathrm {Map}}_{\mathcal {G}_n}(M,N)$ (also using Corollary 4.8). See Appendix B for details about this comparison.

Proof Write

$$ \begin{align*} d_{\infty}(\widehat{v_n}, w) \leqslant d_{\infty}(\widehat{v_n}, \widehat{w_n}) + d_{\infty}(\widehat{w_n}, w). \end{align*} $$

The second term $d_{\infty }(\widehat {w_n}, w)$ converges to zero by Corollary 4.8. It remains to show that $d_{\infty }(\widehat {v_n}, \widehat {w_n}) \to 0$ . By Theorem 2.3 (iii), $d_{\infty }(\widehat {v_n}, \widehat {w_n}) \leqslant d_{\infty }(v_n, w_n)$ . Using Corollary B.4, we find that $d_{\infty }(v_n,w_n) = o\left ( r^{2-\dim M} \right )$ , and we conclude that $d_{\infty }(v_n, w_n) = o(1)$ .▪

5.4 Convergence of the energy

One would like to discuss convergence of the discrete minimizer $\widehat {v_n}$ to the smooth harmonic map w in the Sobolev space $\operatorname {\mathrm {H}}^1(M,N)$ , say, under the assumptions of Theorem 5.1, but this function space (or rather its topology) is not well-defined, see Remark 4.9. It is however still reasonable to ask whether the energy of $v_n$ converges to the energy of w.

We shall see that it does not cost much to prove that the discrete energy $E_n(v_n)$ converges to $E(w)$ , however, it is much more difficult to show that the energy of the interpolation $E(\widehat {v_n})$ also converges to $E(w)$ . While we believe that $E_n(v_n)$ and $E(\widehat {v_n})$ are asymptotic, proving it is too hard in the current state of our work. We will thus be content with stating the desired convergence result under very restrictive assumptions.

Proof First write that $E(w) = \lim _{n \to +\infty } E_n(w_n)$ by Theorem 4.27. Thus it is sufficient to show that $E_n(w_n)$ and $E_n(v_n)$ are asymptotic. By Lemma 5.6 (5.3) applied to $F = E_n$ , we find that

$$ \begin{align*} 0 \leqslant E_n(w_n) - E_n(v_n) \leqslant \frac{\left\Vert \tau_{\mathcal{G}_n}(w_n) \right\Vert{}^2 }{\alpha_n}, \end{align*} $$

so with (4.7) we find that $\big |E_n(w_n) - E_n(v_n)\big | = O\left (r^{2-c}\right )$ and the claim follows.

For the second claim, first write that $E(w) = \lim _{n \to +\infty } E(\widehat {w_n})$ by Corollary 4.8. Thus it is sufficient to show that $E(\widehat {w_n})$ and $E(\widehat {v_n})$ are asymptotic. One can derive from Theorem 2.3 (iii) and Proposition 4.2 (ii) that for a fine and crystalline sequence of meshes,

$$ \begin{align*} \left| \left\Vert \mathop{}\!\mathrm{d} \widehat{f} (x) \right\Vert - \left\Vert \mathop{}\!\mathrm{d} \widehat{g} (x) \right\Vert \right| = O\left(\frac{d_{\infty}(f,g)}{r}\right)~, \end{align*} $$

uniformly in $f, g \in \operatorname {\mathrm {Map}}_{\mathcal {G}_n}(M,N)$ and in $x \in M$ in the interior of the triangulation, from which it follows $\big | E( \widehat {f}) - E( \widehat {g}) \big | = O\left (\frac {d_{\infty }(f,g)}{r}\right )$ . In our case this gives $\big | E( \widehat {w_n}) - E( \widehat {v_n}) \big | = O\left (\frac {d_{\infty }(w_n,v_n)}{r}\right )$ . By Lemma 5.7, Theorem 4.25 (1), and Corollary B.4, we have $d_{\infty }(w_n, v_n) = o\left (r^{2-c-\frac {\dim M}{2}}\right )$ , so we find $\big | E( \widehat {w_n}) - E( \widehat {v_n}) \big | = o\left (r^{1-c-\frac {\dim M}{2}}\right )$ hence $\big | E( \widehat {w_n}) - E( \widehat {v_n}) \big | = o(1)$ when $c = 0$ and $\dim M = 2$ .▪

5.5 Convergence in time and space of the discrete heat flow

We turn to more practical considerations about how to compute harmonic maps. In the previous subsections, we established that, under suitable assumptions, the discrete harmonic map $v_n$ converges to the smooth harmonic map w. In our previous work [Reference Gaster, Loustau and Monsaingeon8], we showed that for each fixed $n \in \mathbb {N}$ , $v_n$ may be computed as the limit of the discrete heat flow $u_{k, n}$ when $k \to +\infty $ . While this is relatively satisfactory, in practice one cannot wait for the discrete heat flow to converge for each n. Hence, it is preferable to let both indices k and n run to $+\infty $ simultaneously. In the theory of PDEs, this situation with a double discretization in time and space is typical—they call it full discretization, and one expects convergence to the solution provided that the time step and the space step satisfy a constraint, called a CFL condition. We are happy to report a similar result.

We keep the same setup as in the beginning of the section. Let $u \in \mathcal {C}$ be a smooth map, denote by $u_n \in \operatorname {\mathrm {Map}}_{\mathcal {G}_n}(M,N)$ its discretization. For each $n \in \mathbb {N}$ , denote by $(u_{k,n})_{k\in \mathbb {N}}$ the sequence in $\operatorname {\mathrm {Map}}_{\mathcal {G}_n}(M,N)$ obtained by iterating the discrete heat flow from the initial map $u_{0,n} = u_n$ . We recall that the discrete heat flow is defined by

$$ \begin{align*} u_{k+1,n} = u_{k, n} + t_n \tau_{\mathcal{G}_n}(u_{k, n}), \end{align*} $$

where $t_n$ is a suitably chosen time step and we use the notation $x+ v$ for the Riemannian exponential map $\exp _x (v)$ in N. We recall that the discrete heat flow is just a fixed stepsize gradient descent method for the discrete energy functional $E_n$ on the Riemannian manifold $\operatorname {\mathrm {Map}}_{\mathcal {G}_n}(M,N)$ . In particular, strong convexity of the $E_n$ implies convergence of the discrete heat flow to the unique discrete harmonic map $v_n$ with exponential convergence rate. We refer to [Reference Gaster, Loustau and Monsaingeon8] for more details.

Theorem 5.15 Consider the same setup and assumptions as in Theorem 5.1. Also assume that for any constant $K>0$ , the discrete energy $E_n$ has Hessian bounded above by $\beta _{n,K} = O(r^{-d})$ on its sublevel set $\{E_n \leqslant K\}$ , for some $d \geqslant 0$ independent of K. Then

$$ \begin{align*} \widehat{u_{k,n}} \xrightarrow[k,n \to +\infty]{} w \quad \text{in } \operatorname{\mathrm{L}}^2(M, N), \end{align*} $$

provided the CFL condition:

(5.6) $$ \begin{align} k = \Omega\left( \frac{\log(r^{-1})}{r^{c+d}} \right). \end{align} $$

Proof Let us break the proof into a few key steps.

Lemma 5.18 There exists a constant $K>0$ such that

$$ \begin{align*} E_n(u_{k,n}) \leqslant K \end{align*} $$

for all $k,n \in \mathbb {N}$ .▪

Proof For each fixed $n \in \mathbb {N}$ , the discrete energy $E_{n}(u_{k,n})$ is nonincreasing with k, since the discrete heat flow is a gradient descent for the discrete energy. In particular, $E_n(u_{k,n}) \leqslant E_n(u_{0,n})$ . To conclude, we must argue that the sequence $(E_n(u_n))_{n\in \mathbb {N}}$ is bounded. This is true since it converges to $E(u)$ by Theorem 4.27.▪

Lemma 5.19 For every $k, n$ , we have

$$ \begin{align*} d(u_{k,n}, v_n) \leqslant c_n q_n^k, \end{align*} $$

where $c_n = O\left (r^{-c/2}\right )$ and $q_n = 1 - C r^{c + d} + o\left (r^{c+d}\right )$ with $C>0$ .

Proof This is an immediate consequence of [Reference Gaster, Loustau and Monsaingeon8, Theorem 4.1]. Note that for the estimate of $c_n$ , we need to use the fact that $E_n(u_{0,n}) = O(1)$ , which we showed in Lemma 5.18.▪

We now finish the proof of Theorem 5.15. For every $k, n \in \mathbb {N}$ , we have

$$ \begin{align*} d(\widehat{u_{k,n}}, w) \leqslant d(\widehat{u_{k,n}}, \widehat{v_n}) + d(\widehat{v_n}, w). \end{align*} $$

The second term $d(\widehat {v}_n, w)$ converges to zero by Theorem 5.1. As for the first term, we have $d(\widehat {u_{k,n}}, \widehat {v_n}) \leqslant L\, d(u_{k,n}, v_n)$ for some constant $L>0$ by Theorem 4.22. Thus it is enough to show that $d(u_{k,n}, v_n) \to 0$ under the appropriate CFL condition.

Let $(\varepsilon _n)_{n\in N}$ be a sequence of positive real numbers converging to zero to be chosen later. Since $(u_{k,n})$ converges to $v_n$ when $k \to +\infty $ , there exists $k_0(n)$ such that $d(u_{k,n}, v_n) \leqslant \varepsilon _n$ for all $k \geqslant k_0(n)$ . Note that the inequality $k \geqslant k_0(n)$ is the CFL condition that we are after, for a/any choice of $(\varepsilon _n)$ . It is possible to compute $k_0(n)$ explicitly with Lemma 5.19; one finds that

$$ \begin{align*} k_0(n) = \frac{\log(c_n) + \log(\varepsilon_n^{-1})}{\log(q_n^{-1})} \end{align*} $$

is sufficient. With our estimates we get $\log (c_n) = \Theta (\log (r^{-1}))$ and $\log (q_n^{-1}) \sim C r^{c+d}$ . It is easy to choose $\varepsilon _n$ so that $\log (\varepsilon _n^{-1})$ is negligible compared to $\log (r^{-1})$ , e.g., $\varepsilon _n = r_n$ . We thus find $ k_0(n) = \Theta \left ( \frac {\log (r^{-1})}{r^{c+d}} \right )$ as desired.

5.6 Application to surfaces

When M and N are both two-dimensional, our previous work [Reference Gaster, Loustau and Monsaingeon8] gives estimates for the strong convexity of the discrete energy. More precisely, consider the following setup:

Let $S = M$ and N be closed Riemannian surfaces of negative Euler characteristic. Assume N has negative sectional curvature. Assume that S is equipped with a fine and crystalline sequence of meshes $(\mathcal {M}_n)_{n\in \mathbb {N}}$ , equipped with asymptotic volume weights and positive edge weights such that the ratio of any two edge weights is uniformly bounded. Consider a homotopy class of maps $\mathcal {C} \subset \mathcal {C}^{\infty }(M,N)$ of nonzero degree, and its discretization $\mathcal {C}_n$ along each mesh.

Lemma 5.21 The discrete energy functional $E_n \colon \mathcal {C}_n \to \mathbb {R}$ has Hessian bounded below by $\alpha _n$ and above by $\beta _{n,K}$ on any sublevel set $\{E_n \leqslant K\}$ , with

$$ \begin{align*} \begin{aligned} \alpha_n &= \Omega(r), \\ \beta_{n,K} &= O(r^{-2}). \end{aligned} \end{align*} $$

As a consequence of Lemma 5.21 and the previous theorems of this section, we obtain the following theorem for surfaces.

The previous theorems of this section (Theorems 5.10 and 5.14) also show that under the stronger assumption $\alpha _n = \Omega (1)$ (which we believe holds in a very general setting), the conclusions of the previous theorem may be strengthened:

6.1 Description

Let $S = M$ be a two-dimensional compact Riemannian manifold. One could consider complete metrics with punctures and/or geodesic boundary, but for simplicity we assume S is closed.

Consider a sequence of meshes $(\mathcal {M}_n)_{n\in \mathbb {N}}$ with underlying graphs $(\mathcal {G}_n)_{n\in \mathbb {N}}$ defined by:

• • $\mathcal {M}_0$ is any acute triangulation.

• • $\mathcal {M}_{n+1}$ is obtained from $\mathcal {M}_n$ by midpoint subdivision (see Section 2.2).

Furthermore, equip $\mathcal {G}_n$ with the volume vertex weights (Section 3.2) and the cotangent weights (Section 3.3).

6.2 Angle properties

In order for $\Delta $ -sequences to be crystalline and have reasonable edge weights systems, we need to address some questions about the behavior of angles when iterating midpoint subdivision:

1. (1) Do all the angles of the triangulation remain bounded away from zero?

2. (2) Do all angles remain acute?

3. (3) Do all angles remain bounded away from $\frac \pi {2}$ ?

These questions are surprisingly subtle. In the context of applying Theorem 5.23, they are quite natural since: (1) is necessary and sufficient for the sequence of meshes to be crystalline (see Proposition 4.2), (2) is sufficient for the edge weights to remain positive, and (3) is necessary for the ratio of any two edge weights to remain uniformly bounded, a requirement of Theorem 5.23.

As remarked above, (1) has a positive answer for surfaces of constant curvature by the main theorem of [Reference Brunck3], and we expect a positive answer to follow from compactness, without any curvature assumption. Questions (2) and (3) are more subtle, and the answer seems likely to depend on the initial triangulation chosen, even for hyperbolic surfaces. In practice, the triangulations of hyperbolic surfaces produced in our software Harmony (roughly speaking, chosen to maximize the smallest angles, see [Reference Gaster, Loustau and Monsaingeon8, Section 6.2]) seem to always yield positive answers to both (2) and (3). It seems possible, however, to produce contrived initial triangulations whose answers will be negative to both (2) and (3). Nonetheless, the Riemannian estimates we carry out to produce Laplacian qualities of mesh sequences can be used to verify positive answers to (1)–(3) in great generality.

To avoid burdening our presentation with technical Riemannian geometry estimates, we postpone this proof to the Appendix: see Propositions A.13 and A.15 in Section A.2.

We say a sequence of acute triangulations is strongly acute if the angles remain uniformly bounded away from $0$ and from $\frac {\pi }{2}$ . Thus, provided the initial triangulation is sufficiently fine, any sequence of triangulations obtained from iterated refinement as in Theorem 6.4 will be strongly acute.

We record the following easy consequence of Definition 6.1.

6.3 Laplacian qualities

Let $(\mathcal {M}_n)_{n\ \in \mathbb {N}}$ be a $\Delta $ -sequence in $(S,g)$ , denote $(\mathcal {G}_n)_{n\in \mathbb {N}}$ the underlying graphs.

The decomposition $\mathcal {V}_n = \mathcal {V}_n^{(0)} \sqcup \mathcal {V}_n^{(1)} \sqcup \mathcal {V}_n^{(2)}$ thus makes any $\Delta $ -sequence a candidate to be almost asymptotically Laplacian: see Definition 4.6. The main theorem of this section provides a positive answer:

Proof There are several conditions to check: refer to Definition 4.6. Conditions (i) and (ii) are trivially satisfied by definition of a $\Delta $ -sequence.

It remains to check the Laplacian qualities stated in (iii’), namely that $\mathcal {G}_n$ is Laplacian up to $O\left (r^{2-k}\right )$ on $\mathcal {V}_n^{(k)}$ for $k\in \{0,1,2\}$ . For each k, there are three conditions to check: the first-order, second-order, and third-order Laplacian conditions, up to $O\left (r^{2-k}\right )$ (see Definition 4.5 (iii)). There are thus nine conditions to check, some of which can be grouped together.

This first lemma almost comes “for free”:

Lemma 6.8 At any vertex $x \in \mathcal {V}_n$ , the j-th order Laplacian condition (for $j \in \{1, 2, 3\}$ ) holds up to $O(r^{j-2})$ .▪

Proof We have $\mu _x = \Theta \left (r_n^2\right )$ (Theorem 4.4), $\omega _{xy} = \Theta (1)$ (Corollary 6.5) and $\overrightarrow {xy} = O(r)$ for any $y \sim x$ , therefore

$$ \begin{align*} \frac{1}{\mu_x} \sum_{y \sim x} \omega_{xy} L(\overrightarrow{xy})^j = \Vert L \Vert^j O(r^{j-2}). \end{align*} $$

The conclusion easily follows for each $j \in \{1, 2, 3\}$ .▪

In what follows, we will frequently need to compare our present Riemannian setting to its “Euclidean counterpart.” Let us clarify what we typically mean by that. Consider a vertex $x \in \mathcal {V}_n \subseteq S$ and its neighbors $\{y_i\} \subseteq S$ . By working in the normal chart at x, we can imagine that x and $\{y_i\}$ live in the Euclidean plane $\operatorname {\mathrm {T}}_x S$ . In this plane, each edge of the triangulation, a Riemannian geodesic, may be replaced by a Euclidean straight segment, yielding a Euclidean triangulation. One can then define, for instance, the Euclidean cotangent weights associated to this Euclidean triangulation. We shall call the Euclidean cotangent weights $\omega _{xy}^{\mathrm {E}}$ the Euclidean counterparts of the cotangent weights $\omega _{xy}$ .

Lemma 6.9 The cotangent weights $\omega _{xy}$ are within $O(r^2)$ of their Euclidean counterparts $\omega _{xy}^{\mathrm {E}}$ .

Proof This follows directly from the first-order expansion of the cotangent given in Proposition A.5. Note that we need to know that all angles are bounded away from $0$ and $\frac {\pi }{2}$ , which is guaranteed by definition of a strongly acute $\Delta $ -sequence.▪

The fact that the cotangent weights are exactly Laplacian to first order in the Euclidean setting (Proposition 3.9) and the previous lemma allow us to upgrade the $j=1$ case of Lemma 6.8:

Lemma 6.10 At any vertex $x \in \mathcal {V}_n$ , the first-order Laplacian condition holds up to $O(r)$ .

Proof Write

$$ \begin{align*} \sum_{y \sim x} \omega_{xy} \overrightarrow{xy} = \sum_{y \sim x} \left(\omega_{xy} - \omega_{xy}^{\mathrm{E}}\right) \overrightarrow{xy} ~+ ~\sum_{y \sim x} \omega_{xy}^{\mathrm{E}} \overrightarrow{xy}. \end{align*} $$

The first sum is $O(r^3)$ by Lemma 6.9 and the second vanishes by Proposition 3.9 (note that $\overrightarrow {xy} := \exp _x^{-1}(y)$ is equal to its “Euclidean counterpart” $\overrightarrow {xy}^{\mathrm {E}}$ , since we are looking at the normal chart at x). Since $\mu _x = O(r^2)$ , conclude that $\frac {1}{\mu _x} \sum _{y \sim x} \omega _{xy} \overrightarrow {xy} = O(r)$ .▪

As far as the first-order Laplacian condition is concerned, Lemma 6.10 is good enough for vertices $x \in \mathcal {V}_n^{(2)}$ and $x \in \mathcal {V}_n^{(1)}$ . However, for $x \in \mathcal {V}_n^{(0)}$ , we need to upgrade the estimate to $O(r^2)$ . Essentially, this follows from the fact that interior vertices have “almost central symmetry,” and second-order Riemannian estimates. The computations are tedious but fairly straightforward, we condensed them in the proof of the next lemma:

Lemma 6.11 At any interior vertex $x \in \mathcal {V}_n^{(0)}$ , the first-order Laplacian condition holds up to $O(r^2)$ .

Proof We need to push one step further the asymptotic expansion of the cotangent weights mentioned in Lemma 6.9. Order the neighbors of x cyclically, and given a neighbor y, denote $y'$ and $y''$ the previous and the next neighbors. By Proposition A.12, we have

$$ \begin{align*} \omega_{xy} = \omega_{xy}^{\mathrm{E}} + \lambda_{xy} + O(r^3) \quad \text{ with } \quad \lambda_{xy} = \frac{1}{2} \left(\varepsilon_{x y' y} + \varepsilon_{x y'' y} \right), \end{align*} $$

where the notation $\varepsilon _{O A B}$ is defined in Proposition A.12. It follows that

$$ \begin{align*} \sum_{y \sim x} \omega_{xy} \overrightarrow{xy} = \sum_{y \sim x} \omega_{xy}^{\mathrm{E}} \overrightarrow{xy} ~+ ~\sum_{y \sim x} \lambda_{xy} \overrightarrow{xy} ~+ ~O(r^4). \end{align*} $$

The first sum vanishes as in Lemma 6.10. Since $\mu _x = O(r^2)$ , we need to show that $\sum _{y \sim x} \omega _{xy} \overrightarrow {xy} = O(r^4)$ . Hence we win if we show that $\sum _{y \sim x} \lambda _{xy} \overrightarrow {xy} = O(r^4)$ .

We note that any interior vertex $x \in \mathcal {V}_n^{(0)}$ has “almost central symmetry” up to $O(r^3)$ , meaning that its set of neighbors may be divided into pairs $\{y_+, y_-\}$ such that $\overrightarrow {x y_+} + \overrightarrow {x y_-} = O(r^3)$ (equivalently, the central symmetry at x preserves the set of neighbors up to $O(r^3)$ ). This immediately follows from the fact that $x \in \mathcal {V}_n^{(0)}$ has in fact “almost hexaparallel symmetry,” as we shall see in (Lemma 6.12).

Now write

$$ \begin{align*} \begin{aligned} \sum_{y \sim x} \lambda_{xy} \overrightarrow{xy} &= \sum_{\{y_+, y_-\}} \lambda_{xy_+} \overrightarrow{xy_+} + \lambda_{xy_-} \overrightarrow{xy_-}\\ &= \sum_{\{y_+, y_-\}} \left(\lambda_{xy_+} - \lambda_{xy_-}\right) \overrightarrow{xy_+} + \lambda_{xy_-} \left(\overrightarrow{xy_-} + \overrightarrow{xy_+} \right) . \end{aligned} \end{align*} $$

It is not hard to see from the expression of $\lambda _{xy}$ that $\lambda _{xy} = O(r^2)$ , and, due to the almost central symmetry, $\lambda _{xy_+} - \lambda _{xy_-} = O(r^4)$ . (To be fair, it is a few lines of calculations, but let us skip the unnecessary details.) We also have $\overrightarrow {xy_-} + \overrightarrow {xy_+} = O(r^3)$ , we thus derive from the previous identity that $\sum _{y \sim x} \lambda _{xy} \overrightarrow {xy} = O(r^5)$ , which is better than the $O(r^4)$ desired result.▪

At this point, it is good to pause and see that we have proved that the first-order Laplacian condition holds up to $O(r^{2-k})$ on $\mathcal {V}_n^{(k)}$ for all $k\in \{0,1,2\}$ , as required. Let us now turn to the second-order condition. On $\mathcal {V}_n^{(2)}$ , we have already proved that it holds up to $O(1)$ as required: see Lemma 6.8. Let us now show that it holds up to $O(r^2)$ on $\mathcal {V}_n^{(1)}$ (better than the required $O(r)$ ) and on $\mathcal {V}_n^{(0)}$ (as required). Along with Lemma 6.10, this is the most difficult part of the proof.

Lemma 6.12 At any interior vertex $x \in \mathcal {V}_n^{(0)}$ or boundary vertex $x \in \mathcal {V}_n^{(1)}$ , the second-order Laplacian condition holds up to $O(r^2)$ .

Proof Let $x \in \mathcal {V}_n^{(0)}$ be an interior vertex. Using Riemannian estimates, we shall prove that the second-order Laplacian condition holding up to $O(r^2)$ is a consequence of the fact that x has “almost hexaparallel symmetry.” We defined hexaparallel symmetry in the Euclidean setting: see Definition 3.3. This definition naturally extends to the Riemannian setting, using the normal chart at x to bring x and its neighbors back to the Euclidean setting. We further say that x has almost hexaparallel symmetry (up to $O(r^3)$ ) provided that the neighbors of x are within $O(r^3)$ of a hexaparallel configuration. Using Proposition A.8, one quickly shows that any interior vertex has almost hexaparallel symmetry.

Denote $\hat {y}_i$ the hexaparallel configuration around x such that $\hat {y}_i - y_i = O(r^3)$ , and denote $\omega _{x\hat {y}}^{\mathrm {E}}$ the Euclidean counterparts of the cotangent weights $\omega _{x\hat {y}}$ . As in Lemma 6.10, one shows that $\omega _{xy}$ , $\omega _{x\hat {y}}$ , and $\omega _{x\hat {y}}^{\mathrm {E}}$ are all within $O(r^2)$ . Now write

$$ \begin{align*} \frac{1}{\mu_x} \sum_{y\sim x} \omega_{xy} L\left(\overrightarrow{xy}\right)^2 &= \frac{1}{\mu_x} \sum_{y\sim x} (\omega_{xy} - \omega_{x\hat{y}}^{\mathrm{E}}) L\left(\overrightarrow{xy}\right)^2 \\ &\quad + \frac{1}{\mu_x} \sum_{y\sim x} \omega_{x\hat{y}}^{\mathrm{E}} L\left(\overrightarrow{xy} - \overrightarrow{x \hat{y}}\right) L\left(\overrightarrow{xy} + \overrightarrow{x \hat{y}}\right) + \frac{1}{\mu_x} \sum_{y\sim x} \omega_{x\hat{y}}^{\mathrm{E}} L\left(\overrightarrow{x \hat{y}}\right)^2. \end{align*} $$

One quickly sees that the first and second sums are $\Vert L \Vert ^2 O\left (r^2\right )$ . As for the third sum, first note that denoting $\mu _x^{\mathrm {E}}$ the Euclidean area weight at x, we have

$$ \begin{align*} \frac{1}{\mu_x^{\mathrm{E}}} \sum_{y\sim x} \omega_{x\hat{y}}^{\mathrm{E}} L\left(\overrightarrow{x \hat{y}}\right)^2 = 2 \Vert L \Vert^2, \end{align*} $$

by Lemma 3.11. Since $\mu _x = \mu _x^{\mathrm {E}} \left (1 + O\left (r^2\right )\right )$ by Proposition A.11, we find

$$ \begin{align*} \frac{1}{\mu_x} \sum_{y\sim x} \omega_{x\hat{y}}^{\mathrm{E}} L\left(\overrightarrow{x \hat{y}}\right)^2 = 2 \Vert L \Vert^2 \left(1 + O\left(r^2\right)\right) . \end{align*} $$

Gathering all three sums, we find $\frac {1}{\mu _x} \sum _{y\sim x} \omega _{xy} L\left (\overrightarrow {xy}\right )^2 = 2 \Vert L \Vert ^2 \left (1 + O\left (r^2\right )\right )$ as desired.

One conducts a similar proof when x is a boundary vertex: in that case, it has almost semi-hexaparallel symmetry up to $O(r^3)$ , and the proof is similarly derived from the Euclidean case.▪

This concludes the proof that the second-order Laplacian condition holds up to $O(r^{2-k})$ on $\mathcal {V}_n^{(k)}$ for all $k\in \{0, 1, 2\}$ . Let us finally examine the third-order condition. We already proved in Lemma 6.8 that it holds up to $O(r)$ at any vertex, which is good enough for $\mathcal {V}_n^{(2)}$ and $\mathcal {V}_n^{(1)}$ . It remains to prove that it holds up to $O(r^2)$ on $\mathcal {V}_n^{(0)}$ . It actually holds up to $O(r^3)$ :

Lemma 6.13 At any $x \in \mathcal {V}_n^{(0)}$ , the third-order Laplacian condition holds up to $O(r^3)$ .

Proof This is an easy consequence of the almost central symmetry: write

$$ \begin{align*} \begin{aligned} \sum_{y \sim x} \omega_{xy} L(\overrightarrow{xy})^3 &= \sum_{\{y_+, y_-\}} \omega_{xy_+} L(\overrightarrow{xy_+})^3 + \omega_{xy_-} L(\overrightarrow{xy_-})^3\\[3pt] &= \sum_{\{y_+, y_-\}} \left(\omega_{xy_+} - \omega_{xy_-}\right) L(\overrightarrow{xy_+})^3 + \omega_{xy_-} \left(L(\overrightarrow{xy_-})^3 + L(\overrightarrow{xy_+})^3 \right) . \end{aligned} \end{align*} $$

By almost central symmetry, we have $\omega _{xy_+} - \omega _{xy_-} = O(r^2)$ and $\overrightarrow {xy_-} + \overrightarrow {xy_+} = O(r^3)$ . It follows that the first term is $\Vert L \Vert ^3 O(r^5)$ , as is the second term. (For the second term, write $L(\overrightarrow {xy_-}) = L(\overrightarrow {xy_+}) + \Vert L \Vert O(r^3)$ and expand the third power of this identity.) Thus we find that $\sum _{y \sim x} \omega _{xy} L(\overrightarrow {xy})^3 = \Vert L \Vert ^3 O(r^5)$ , therefore, $\frac {1}{\mu _x}\sum _{y \sim x} \omega _{xy} L(\overrightarrow {xy})^3 = \Vert L \Vert ^3 O(r^3)$ as required.▪

This concludes the proof that the third-order Laplacian condition holds up to $O(r^{2-k})$ on $\mathcal {V}_n^{(k)}$ for all $k\in \{0,1,2\}$ . The proof of Theorem 6.7 is now complete.