Computing harmonic maps between Riemannian manifolds

Abstract

In our previous paper (Gaster et al., 2018, arXiv:1810.11932), we showed that the theory of harmonic maps between Riemannian manifolds, especially hyperbolic surfaces, may be discretized by introducing a triangulation of the domain manifold with independent vertex and edge weights. In the present paper, we study convergence of the discrete theory back to the smooth theory when taking finer and finer triangulations, in the general Riemannian setting. We present suitable conditions on the weighted triangulations that ensure convergence of discrete harmonic maps to smooth harmonic maps, introducing the notion of (almost) asymptotically Laplacian weights, and we offer a systematic method to construct such weighted triangulations in the two-dimensional case. Our computer software Harmony successfully implements these methods to compute equivariant harmonic maps in the hyperbolic plane.

1 Introduction

Let M and N be Riemannian manifolds, let us assume M compact and N complete. A harmonic map $f \colon M \to N$ is a critical point of the energy functional

$$ \begin{align*} E(f) = \frac{1}{2} \int_M \Vert \mathop{}\!\mathrm{d} f\Vert^2 \mathop{}\!\mathrm{d} v. \end{align*} $$

Equivalently, f has vanishing tension field $\tau (f) = 0$ , a nonlinear generalization of the Laplace operator that can be defined as the trace of the Riemannian Hessian: $ \tau (f) = \nabla (\mathop {}\!\mathrm {d} f)$ . When N is compact and has negative sectional curvature, there exists a harmonic map $M \to N$ in any homotopy class of smooth maps, and it is unique unless it is constant or maps to a geodesic. This foundational result due to Eells–Sampson [7] and Hartman [9] can be understood in terms of the convexity properties of the energy. Essentially, the curvature assumption on N implies that the energy functional is convex on any component of the space of smooth maps $\mathcal {C}^{\infty }(M,N)$ , which guarantees convergence of the gradient flow—also called heat flow in this setting—from any initial smooth map to the energy minimizer.

In our previous work [8], which mostly specialized to surfaces, we showed that the theory can be appropriately discretized by meshing the domain manifold with a triangulation and assigning two independent systems of weights, on the set of vertices and edges, respectively. One of the main results is the strong convexity of the discrete energy functional, from which we derive convergence of the discrete heat flow to the unique discrete harmonic map. (The second focus of [8] is on center of mass methods, which we do not discuss in the present paper.) While that paper was concerned with a fixed discretization, the purpose of the present paper is to study the convergence of the discrete theory back to the smooth theory when one takes finer and finer meshes.

After introducing the discretization setup in Section 2, in Section 3 we discuss special conditions on weighted triangulations in order to adequately capture the local geometry of the domain manifold. We define Laplacian systems of weights, which aim to produce a good approximation of the Laplacian (i.e., tension field) by the discrete Laplacian. As a fundamental example, we introduce our favorite volume vertex weights and cotangent edge weights.

In Section 4, we study fine sequences of meshes (with maximum edge length converging to zero), and the approximation of the relevant smooth objects by their discrete counterparts. A key requirement for the sequence is to be crystalline, meaning that all angles of the triangulation stay bounded away from zero. We also strategically weaken the notion of Laplacian weights to (almost) asymptotically Laplacian weights. We show that for such sequences of weighted meshes, there is convergence of the discrete volume form, tension field, energy density, and energy to their smooth counterparts.

In Section 5, we study the convergence of discrete maps to smooth harmonic maps. If the discrete energy is sufficiently convex, and the sequence of meshes is almost asymptotically Laplacian, we prove that (the center of mass interpolations of) the discrete harmonic maps converge to the unique smooth harmonic map in $\operatorname {\mathrm {L}}^2$ . We expect the strong convexity assumption to hold in a very broad setting, and have proved it in the two-dimensional case in [8]. Pending stronger assumptions, we also show convergence in $\operatorname {\mathrm {L}}^{\infty }$ , and in energy. Furthermore, we show that the discrete heat flow starting from any discretized map converges to the smooth harmonic map when both the time index and the space index run to $+\infty $ , provided a CFL-type condition is satisfied. This theorem may be seen as a constructive implementation of the theorem of Eells–Sampson and Hartman.

Section 6 describes how to systematically construct almost asymptotically Laplacian sequences of meshes, so that our previous theorems can apply, at least in the two-dimensional case. These are quite simply constructed by iterated midpoint geodesic subdivision from an initial triangulation of the domain manifold, and taking the volume weights on vertices and cotangent weights on edges. Proving the required Laplacian qualities requires some delicate Riemannian geometry estimates, naturally building on the Euclidean case; we largely relegate these to Appendix A to avoid burdening our exposition. It is quite remarkable how the conditions for our constructed sequences to be almost asymptotically Laplacian are barely met, and in turn how these conditions are barely sufficient for our main convergence theorem (Theorem 5.1) to hold.

Putting together the main theorems in Sections 5 and 6 (Theorems 5.1, 5.23, and 6.7), we obtain explicit constructions of sequences of discretizations that ensure convergence to the desired harmonic map. Here is a sample theorem summarizing our main results for surfaces:

Theorem Let M and N be compact Riemannian $2$ -manifolds of negative Euler characteristics, and assume N has negative sectional curvature. Consider a sequence of meshes on M obtained by iterated midpoint subdivision with all angles bounded away from $0$ and from $\pi /2$ , and equip it with the area vertex weights and cotangent edge weights. Let $\mathcal {C}$ be a component of $\mathcal {C}^{\infty }(M,N)$ of nonzero degree, and let $v_n$ be the unique discrete harmonic map in the corresponding discrete homotopy class. Then $v_n$ converges to the unique harmonic map $w \in \mathcal {C}$ in the $\operatorname {\mathrm {L}}^2$ topology.▪

This construction and the discrete heat flow is implemented in our freely available computer software Harmony, which is presented in our previous paper [8]. Harmony computes the unique harmonic map from the hyperbolic plane to itself that is equivariant with respect to the actions of two Fuchsian groups, which can be selected by the user via Fenchel–Nielsen coordinates.

Much of the theory and techniques that we develop are well-known in the Euclidean setting, such as the discrete heat flow method or the cotangent weights popularized by Pinkall–Polthier [13]. This paper builds upon the Euclidean theory using fine meshes on Riemannian manifolds. However, there are notable differences from the Euclidean setting: First, the Laplace equation is linear in the Euclidean setting, allowing finite element methods. Second, we restrict to compact manifolds without boundary, in contrast to Euclidean domains where boundary conditions are prescribed. Finally, there are important consequences of negative curvature, including the strong convexity of the energy functional and the uniqueness of harmonic maps, that we exploit in the present project.

The program to discretize the theory of harmonic maps between Riemannian manifolds, and to obtain convergence back to the smooth theory, remains unfinished. Celebrated work on the discretized theory includes [2, 6, 11], while convergence to the smooth harmonic map has been analyzed for submanifolds of $\mathbb {R}^n$ notably by Bartels [1]. The present paper seems to have some overlap with Bartels’ work, though our setting is more intrinsic and geometric in nature. A perhaps more powerful approach than ours to prove convergence of discrete harmonic maps to smooth harmonic maps would consist in finding a discrete version of Bochner’s formula and possibly Moser’s Harnack inequality: see Remark 5.13.

A note to the reader: Although this paper is the sequel of [8], the two papers can be read independently. We also point out that Sections 5 and 6 in this paper can be read independently.

2 Setup

Throughout the paper, let $(M,g)$ and $(N,h)$ be smooth connected complete Riemannian manifolds. These will be our domain and target respectively. We will typically assume that M is compact and oriented, and that N is Hadamard (complete, simply connected, with nonpositive sectional curvature). Although most of the paper holds in this generality, we are especially interested in the case where $S = M$ is two-dimensional. For background on the smooth theory of harmonic maps $M \to N$ , please refer to [8, Section 1].

2.1 Discretization setup

Our discretization setup is the following. (We also refer to [8, 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 [8, 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 [8, 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 [8] 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 [8, 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 [8, 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 [8], is to work equivariantly in the universal covers.

2.3.2 Working equivariantly

Fix a homotopy class $\mathcal {C}$ of a continuous map $M \to N$ , which induces a group homomorphism $\rho \colon \pi _1 M \to \pi _1 N$ . Recall that any $f \in \mathcal {C}$ admits a $\rho $ -equivariant lift between universal covers $\tilde {f} \colon \tilde {M} \to \tilde {N}$ . The mesh $\mathcal {M}$ on M also lifts to a $\pi _1 M$ -invariant geodesic mesh $\tilde {\mathcal {M}}$ of $\tilde {M}$ . As usual, one has to take more care with basepoints on M and N—and use more notation—to make this story complete.

Definition 2.2 The discrete homotopy class $\mathcal {C}_{\mathcal {G}} := \operatorname {\mathrm {Map}}_{\tilde {\mathcal {G}}, \rho }(\tilde {M}, \tilde {N})$ is defined as the space of $\rho $ -equivariant discrete maps $\tilde {M} \to \tilde {N}$ along $\tilde {\mathcal {G}}$ .

One can then define an interpolation scheme as a continuous right inverse of $\pi _{\mathcal {G}} $ on $\mathcal {C}_{\mathcal {G}} $ . For the purposes of this paper, however, all of the convergence analysis can be performed on the quotient manifolds. The presentation is chosen with ease in mind, and so we overlook the subtlety above. Nevertheless, we point out that there are other benefits to the equivariant setting:

• • It allows one to consider equivariance with respect to group homomorphisms $\rho \colon \pi _1 M \to \operatorname {\mathrm {Isom}}(\tilde {N})$ that are not necessarily induced by continuous maps from M to a quotient of $\tilde {N}$ , e.g., nondiscrete representations $\rho $ .

• • Computationally, it is easier to work in the universal covers. This is the point of view that we chose when coding the software Harmony.

This explains our present change in perspective from the equivariance throughout [8].

2.3.3 Center of mass interpolation

We refer to [8, 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 [8, 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 [8, Lemma 5.3].▪

3 Systems of weights

We follow the discretization setup of Section 2 and seek systems of vertex and edge weights on $\mathcal {G}$ that adequately capture the local geometry of M, in the sense that they ensure a good approximation of the theory of smooth harmonic maps from M to any other Riemannian manifold.

Throughout this section, $(M,g)$ is any Riemannian manifold equipped with a geodesic mesh $\mathcal {M}$ . We denote as usual $\mathcal {G}$ the associated graph.

3.1 Laplacian weights

Definition 3.1 A system of vertex weights $(\mu _x)_{x \in \mathcal {V}}$ and edge weights $(\omega _{xy})_{\{x,y\} \in \mathcal {E}}$ on the graph $\mathcal {G}$ is called Laplacian (to third order) at a vertex $x\in \mathcal {V}$ if, for any linear form $L \in \operatorname {\mathrm {T}}_x^*M$ :

1. (1) (First-order condition)

   $$ \begin{align*} \frac{1}{\mu_x}\sum_{y \sim x} \omega_{xy} \, \overrightarrow{xy} = 0. \end{align*} $$

2. (2) (Second-order condition)

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

3. (3) (Third-order condition)

   $$ \begin{align*} \frac{1}{\mu_x}\sum_{y \sim x} \omega_{xy} \; L(\overrightarrow{xy})^3 = 0. \end{align*} $$

The biweighted graph $(\mathcal {G}, (\mu _x), (\omega _{xy}))$ is called Laplacian if it is Laplacian at any vertex.

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

Remark 3.1 As we shall see, the defining properties of Laplacian weights (or their characterization Proposition 3.4) are remarkably versatile. Perhaps the most obvious motivation for their definition is Theorem 4.14, but we will also use it in different ways, e.g., for Lemma 4.10 or Theorem 4.17.

Remark 3.2 A biweighted graph being Laplacian to first order, i.e., satisfying condition (1), is equivalent to the the fact that each vertex of $\mathcal {G}$ is the weighted barycenter of its neighbors. Theorem 3.3 provides many examples of Laplacian graphs to first order.

Theorem 3.3 Assume $M = S$ is two-dimensional and has nonpositive curvature. Any biweighted graph $\mathcal {G}$ underlying a topological triangulation of S admits a unique map to S that is Laplacian to first order, i.e., whose image graph equipped with the same weights is Laplacian to first order.

Proof Note that a map $f \colon \mathcal {G} \to S$ being Laplacian to first order is equivalent to f having zero discrete tension field, i.e., f being discrete harmonic. By [8, Theorem 3.20], the discrete energy functional in this setting is strongly convex, in particular it has a unique critical point.▪

The following seemingly stronger characterization of Laplacian weights is immediate:

Proposition 3.4 A system of weights on $\mathcal {G}$ is Laplacian at $x \in \mathcal {V}$ if and only if for any finite-dimensional vector space W.

1. (1) For any linear map $L \colon \operatorname {\mathrm {T}}_x M \to W$ :

   $$ \begin{align*} \sum_{y \sim x} \omega_{xy} \; L(\overrightarrow{xy}) = 0. \end{align*} $$

2. (2) For any quadratic form q on $\operatorname {\mathrm {T}}_x M$ with values in W:

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

3. (3) For any cubic form $\sigma $ on $\operatorname {\mathrm {T}}_xM$ with values in W:

   $$ \begin{align*} \sum_{y \sim x} \omega_{xy} \;\sigma(\overrightarrow{xy}) = 0. \end{align*} $$

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).

Theorem 3.6 Let N be any Riemannian manifold of nonpositive sectional curvature. Equip the space of discrete maps $\operatorname {\mathrm {Map}}_{\mathcal {G}} (M, N)$ with the $\operatorname {\mathrm {L}}^2$ distance associated to the volume weights. Then the center of mass interpolation map $\iota _{\mathcal {G}} \colon \operatorname {\mathrm {Map}}_{\mathcal {G}} (M, N) \to \mathcal {C}(M,N)$ is L-Lipschitz with respect to the $\operatorname {\mathrm {L}}^2$ distance on both spaces, with $L = \sqrt {1+\dim M}$ . When M is Euclidean (flat), the Lipschitz constant can be upgraded to $L =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 [13] 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 [4].

The following result noticed by Pinkall–Polthier [13] 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.

Definition 3.2 Let $(S,g)$ be a Riemannian surface equipped with a geodesic mesh with underlying graph $\mathcal {G}$ . The edge weights on $\mathcal {G}$ defined as in (3.4) are the system of cotangent weights.

As a direct application of Lemma 3.7, we obtain:

Proposition 3.8 Let $(S,g)$ be a flat surface with a geodesic mesh. Let $\mathcal {G}$ be the underlying graph equipped with the cotangent edge weights. For any piecewise affine map $f \colon S \to \mathbb {R}^n$ , the smooth energy $E(f) := \frac {1}{2}\int _S \Vert \mathop {}\!\mathrm {d} f \Vert ^2 \mathop {}\!\mathrm {d} v$ coincides with the discrete energy $E_{\mathcal {G}} (f)$ defined in (2.1).

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 [2, 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.

Proposition 3.9 Suppose that $(S,g)$ is a flat surface. Then the cotangent weights associated to any triangulation of S are Laplacian to first order.

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:

Proposition 3.10 Let $(S,g)$ be a flat surface. Let $(\mathcal {G}_n)_{n\in \mathbb {N}}$ be a sequence of graphs obtained by iterated midpoint subdivision from a given initial triangulation. Equip $\mathcal {G}_n$ with the area vertex weights and cotangent edge weights. Then $\mathcal {G}_n$ satisfies the second-order Laplacian condition at any vertex except maybe at the vertices of $\mathcal {G}_0$ .

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:

Definition 3.3 Consider a vertex x with valence six in a Euclidean graph.

• • We say that x has hexaparallel symmetry if the set of vectors $\{\vec {xy} ~\colon ~ y \sim x\}$ is in the $\operatorname {\mathrm {GL}}(2,\mathbb {R})$ -orbit of $\{\pm (1,0),\pm (1,1),\pm (0,1)\}$ . Equivalently, the neighbors of x are the vertices of a hexagon whose opposite sides are pairwise parallel and of the same length. See Figure 5a.

  

  Figure 5: Hexaparallel and semi-hexaparallel symmetry.

• • We say that x has semi-hexaparallel symmetry if the neighbors of x may be cyclically labeled $\{y_1, \dots , y_6\}$ and divided into two overlapping sets $\{y_1, y_2, y_3, y_4\}$ and $\{y_4, y_5, y_6, y_1\}$ , each being part of a potential hexaparallel configuration. See Figure 5b.

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:

Lemma 3.11 Any geodesic graph $\mathcal {G}$ in $\mathbb {R}^2$ equipped with the area vertex weights and cotangent edge weights satisfies the second-order Laplacian condition at any (semi-)hexaparallel vertex x.

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.▪

Corollary 3.12 Suppose that $(S,g)$ is a flat surface. Let $(\mathcal {G}_n)_{n\in \mathbb {N}}$ be a sequence of graphs obtained by iterated midpoint subdivision of an initial triangulation $\mathcal {G}_0$ . Equip $\mathcal {G}_n$ with the area vertex weights and the cotangent edge weights. Then $\mathcal {G}_n$ is Laplacian at any interior vertex.

Proof The first- and third-order conditions are trivial due to central symmetry of the neighbors around the vertex x and the fact that linear and cubic functions are odd. (Alternatively, the first-order condition holds by Proposition 3.9.) The second-order condition holds by Lemma 3.11.▪

Remark 3.13 We shall see in Section 6 that in the general Riemannian setting, the cotangent weights will satisfy similar Laplacian properties asymptotically for very fine meshes.

Remark 3.14 While being the best choice of edge weights, the cotangent weights generally do not satisfy the second-order Laplacian condition at vertices with no (semi-)hexaparallel symmetry. Taking finer and finer triangulations will not help with this defect. At such vertices, which generically exist for topological reasons, the discrete Laplacian of a smooth function can not be expected to approximate its Laplacian. This is somewhat unsettling, but it is an intrinsic difficulty to the discretization of the Laplacian. Providing suitable assumptions that neverthless guarantee convergence of discrete harmonic maps to smooth harmonic maps is the central aim of this paper.

4 Sequences of meshes

In this section, we enhance the previous section by considering sequences of meshes on a Riemannian manifold $(M,g)$ . The idea is to capture the local geometry of M sufficiently well provided the mesh is sufficiently fine. This allows a relaxation of the Laplacian weights conditions, which are too stringent for a fixed mesh of an arbitrary Riemannian manifold. We introduce the notions of asymptotically Laplacian and almost asymptotically Laplacian systems of weights, with the aim that these weakened conditions can still be used to demonstrate the convergence theorems we are after.

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 [5], we define:

Definition 4.1 The sequence $(\mathcal {M}_n)_{n \in \mathbb {N}}$ is called fine provided $\lim \limits _{n \to +\infty } r_n = 0$ .

Notation 4.1 For the remainder of the paper, we drop the subscript $r := r_n$ for ease in notation.

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*} $$

Definition 4.2 The sequence $(\mathcal {M}_n)_{n \in \mathbb {N}}$ is called crystalline if there exists a uniform lower bound for the thickness of simplices in $\mathcal {M}_n$ .

Proposition 4.2 Let $(\mathcal {M}_n)_{n \in \mathbb {N}}$ be a fine sequence of meshes. The following are equivalent:

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

2. (ii) There exists a uniform positive lower bound for all angles between adjacent edges in $\mathcal {M}_n$ .

Proof For brevity, we only sketch the proof; the detailed proof would include proper Riemannian estimates: see Appendix A.

First one checks that (i) $\Leftrightarrow $ (ii) in the Euclidean setting. This is an elementary calculation: for a single triangle (or n-simplex), one can bound its radius in terms of its smallest angle. One then generalizes to an arbitrary Riemannian manifold M by arguing that a very small triangle (or n-simplex) in M has almost the same radius and angles as its Euclidean counterpart in a normal chart. The fact that we only consider fine sequences of meshes means that we can assume that all simplices are arbitrarily small, making the previous argument conclusive.▪

Remark 4.3 In [3, Theorem A], it is shown that the sequence of meshes obtained by midpoint subdivision of a triangle of constant sectional curvature is fine and crystalline, and one concludes immediately the same conclusion for any initial triangulation of a compact hyperbolic surface. An initial preprint of this paper contained a faulty proof of the same conclusion for a triangulation of an arbitrary compact Riemannian surface. While we still expect this conclusion to hold, this subtle question is beyond the scope of this paper. See Section 6.2 for further discussion.

Theorem 4.4 Assume that M is compact and the sequence of meshes $(\mathcal {M}_n)_{n \in \mathbb {N}}$ on M is fine and crystalline. Denote by $\mathcal {G}_n$ the graph underlying $\mathcal {M}_n$ and $r = r_n$ its maximum edge length.

1. (i) The volume vertex weights $\mu _{x, n}$ of $\mathcal {G}_n$ are $\Theta \left (r^{\dim M}\right )$ (uniformly in x).

2. (ii) The number of vertices of $\mathcal {G}_n$ is $\big | \mathcal {V}_n \big | = \Theta \left (r^{-\dim M}\right )$ . More generally, the number of k-simplices of $\mathcal {G}_n$ is $\Theta (r^{-\dim M})$ .

3. (iii) The combinatorial diameter $\operatorname {\mathrm {diam}} \mathcal {G}_n$ of the graph $\mathcal {G}_n$ is $\Theta \left (r^{-1}\right )$ .

4. (iv) The combinatorial surjectivity radius $\operatorname {\mathrm {surj}} \operatorname {\mathrm {rad}} \mathcal {G}_n$ (see below) of the graph $\mathcal {G}_n$ is $\Theta \left (r^{-1}\right )$ .

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.

Notation 4.5 In this paper, we use the notation $f = O(g)$ and $f = o(g)$ in the usual sense, we use the notation $f = \Omega (g)$ for $g = O(f)$ , and $f = \Theta (g)$ for [ $f = O(g)$ and $f = \Omega (g)$ ].

Proof For (i), recall that the volume vertex weight at x is the sum of the volumes of the simplices adjacent to x (divided by $\dim M$ ). Since the sequence is fine, the diameter of all simplices is going to $0$ uniformly in x. On first approximation, the volume of any such vertex is approximately equal to its Euclidean counterpart (say, in a normal chart). Since the lengths of all edges are within $[\alpha r, r]$ for some constant $\alpha>0$ and all angles are bounded below by Proposition 4.2, this volume is $\Theta (r^{\dim M})$ .

For (ii), simply notice that $\sum _{x \in \mathcal {V}_n} \mu _{x, n} = \operatorname {\mathrm {Vol}}(M)$ by Proposition 3.5 and use (i). The generalization to k-simplices is immediate since the total number of k-simplices is clearly $\Theta \left (\big |\mathcal {V}_n\big |\right )$ .

For (iii), let us first show that $\operatorname {\mathrm {diam}} \mathcal {G}_n = \Omega \left (r^{-1}\right )$ . Let x and y be two fixed points in M and denote L the distance between them. For all $n\in \mathbb {N}$ , there exists vertices $x_n$ and $y_n$ in $\mathcal {V}_n$ that are within distance r of x and y, respectively, so their distance in M is $d(x_n, y_n) \geqslant L - 2r$ . Denoting $k_n$ the combinatorial distance between $x_n$ and $y_n$ , one has $d(x_n, y_n) \leqslant k_n r$ by the triangle inequality. We thus find that $k_n r \geqslant L - 2 r$ , hence $\operatorname {\mathrm {diam}} \mathcal {G}_n \geqslant k_n \geqslant L r^{-1} - 2$ so that $\operatorname {\mathrm {diam}} \mathcal {G}_n = \Omega \left (r^{-1}\right )$ . Finally, let us show that $\operatorname {\mathrm {diam}} \mathcal {G}_n = O\left (r^{-1}\right )$ . Let $x_n$ and $y_n$ be two vertices that achieve $\operatorname {\mathrm {diam}} \mathcal {G}_n$ . Let $\gamma _n$ be a length-minimizing geodesic from $x_n$ to $y_n$ . Of course, the length of $\gamma _n$ is bounded above by the diameter of M. There is a sequence of simplices $\Delta _1, \dots , \Delta _{k_n}$ such that $x\in \Delta _1$ , $y\in \Delta _{k_n}$ , and any two consecutive simplices are adjacent. Since the valence of any vertex is uniformly bounded (because of a lower bound on all angles), the number of simplices within a distance $\leqslant r_{\min }$ of any point of M is bounded above by a constant C. This implies $k_n \leqslant C L(\gamma )/r_{\min }$ , so that $k_n \leqslant C (\operatorname {\mathrm {diam}} M) \alpha \, r^{-1}$ . Following edges along the simplices $\Delta _i$ , one finds a path of length $(\dim M - 2) k_n$ from x to y, therefore, $\operatorname {\mathrm {diam}} \mathcal {G}_n \leqslant (\dim M -2) C (\operatorname {\mathrm {diam}} M) \alpha \, r^{-1}$ .

For the proof of (iv), the injectivity radius of M provides a lower bound for $\operatorname {\mathrm {surj}}\operatorname {\mathrm {rad}} \mathcal {G}_n$ of the form $\Omega (r^{-1})$ , and $\operatorname {\mathrm {diam}} \mathcal {G}_n$ provides an upper bound. The details are left to the reader.▪

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 [5] we have:

Lemma 4.6 If $(\mathcal {M}_n)_{n\in \mathbb {N}}$ is a sequence of meshes that is fine and crystalline, then for any piecewise smooth function $f \colon M \to \mathbb {R}$ , the center of mass interpolation $\widehat {f_n}$ converges to f for the piecewise $\mathcal {C}^1$ topology.

Proof As for Proposition 4.2, the proof can be conducted in two steps: First in the Euclidean setting, where the center of mass interpolation $\widehat {f_n}$ is just the piecewise linear approximation of $f_n$ . This proof is done in, e.g., [5]. One then generalizes to an arbitrary Riemannian manifold M by arguing that for very fine triangulations, the center of mass interpolation $\widehat {f_n}$ is very close to the piecewise linear approximation of $f_n$ in a normal chart.▪

Remark 4.7 Any interpolation scheme satisfying the conclusion of Lemma 4.6, as well as Theorems 2.3 and 3.6 (or asymptotic versions thereof), would make the machinery work to prove our upcoming main theorems. One could therefore enforce these properties as the definition of a good sequence of interpolation schemes.

Corollary 4.8 Let $f \colon M \to N$ be a $\mathcal {C}^1$ map between Riemannian manifolds. Assume that M is compact and equipped with a fine and crystalline sequence of meshes $(\mathcal {M}_n)_{n\in \mathbb {N}}$ . The center of mass interpolation $\widehat {f_n}$ converges to f in $\operatorname {\mathrm {L}}^{\infty }(M,N)$ and $E(f) = \lim _{n\to +\infty } E(\widehat {f_n})$ .

Remark 4.9 One would like to say that $\widehat {f_n}$ converges to f in the Sobolev space $\operatorname {\mathrm {H}}^1(M,N)$ , but this space is not well-defined. Actually, $\operatorname {\mathrm {H}}^1(M,N)$ may be defined as the subspace of $\operatorname {\mathrm {L}}^2(M,N)$ consisting of $\operatorname {\mathrm {L}}^2$ maps with finite energy, but it is unclear how to define the $\operatorname {\mathrm {H}}^1$ topology. Nevertheless we can say something in that direction: $\widehat {f_n} \to f$ in $\operatorname {\mathrm {L}}^2(M,N)$ and $E(\widehat {f_n}) \to E(f)$ . One should think of the energy as the $\operatorname {\mathrm {L}}^2$ norm of the derivative, but this “norm” does not induce a distance.

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*} $$

Remark 4.11 For ease of notation, we drop the dependence in n when writing $\mu _x$ and $\omega _{xy}$ above.

Remark 4.12 Before writing the proof, let us clarify the quantifiers in Lemma 4.10 (as well as in Theorems 4.14 and 4.17): The statement is that there exists a constant $M>0$ independent of n such that at any vertex x of $\mathcal {G}_n$ where the system of weights is Laplacian, $\frac {1}{\mu _x}\sum _{y\sim x}\omega _{xy} \leqslant M r^{-2}$ .

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} $$

Notation 4.15 We denote $f_n := \pi _{\mathcal {G}_n}(f)$ , the discretization of f along $\mathcal {G}_n$ (i.e., restriction to $\mathcal {G}_n^{(0)}$ ).

Remark 4.16 The proof below shows that in (4.5), the $O(r^2)$ function depends on f, but may be chosen independent of x if M is compact.

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:

Proposition 4.21 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 any system of asymptotic volume vertex weights converge weakly to the volume density on M.

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.

Theorem 4.22 Let M be a compact Riemannian manifold and let $(\mathcal {M}_n)_{n \in \mathbb {N}}$ be a fine sequence of meshes equipped with a system of asymptotic volume vertex weights. For any complete Riemannian manifold N of nonpositive sectional curvature, the center of mass interpolation map $\iota _n \colon \operatorname {\mathrm {Map}}_{\mathcal {G}_n}(M, N) \to \mathcal {C}(M,N)$ is $L_n$ -Lipschitz with respect to the $\operatorname {\mathrm {L}}^2$ distance on both spaces, with $L_n = 1 + o(1)$ .

Definition 4.5 Let M be a Riemannian manifold. Consider a sequence of geodesic meshes $(\mathcal {M}_n)_{n\in \mathbb {N}}$ with mesh size $r = r_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}}$ asymptotically Laplacian provided that

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

2. (ii) The vertex weights are asymptotic volume weights (see Definition 4.4).

3. (iii) The system of vertex and edge weights on $\mathcal {G}_n$ is Laplacian up to $O\left (r^2\right )$ at all vertices.

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).

Definition 4.6 Let M be a compact Riemannian manifold of dimension m. We say that the sequence of biweighted graphs $(\mathcal {G}_n)_{n\in \mathbb {N}}$ is almost asymptotically Laplacian if it satisfies conditions (i) and (ii) of Definition 4.5, and the modified version of (iii):

1. (iii’) There is a decomposition $\mathcal {V}_n = \bigsqcup _{k=0}^2 \mathcal {V}_n^{(k)}$ , with $\mu _n\left (\mathcal {V}_n^{(k)}\right ) = O(r^{k})$ , so that the system of vertex and edge weights on $\mathcal {G}_n$ is Laplacian up to $O\left (r^{2-k}\right )$ on $\mathcal {V}_n^{(k)}$ .

Remark 4.23 Any asymptotically Laplacian sequence of meshes is almost asymptotically Laplacian: take $\mathcal {V}_n^{(0)} = \mathcal {V}_n$ and $\mathcal {V}_n^{(1)} = \mathcal {V}_n^{(2)} = \emptyset $ .

Remark 4.24 In application, the set $\mathcal {V}_n^{(k)}$ will be the vertices contained in the codimension k-skeleton of a fixed triangulation of M (and not contained in $\mathcal {V}_n^{(k+1)}$ ).

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

Theorem 4.25 Let M be a compact Riemannian manifold. Consider a sequence of geodesic meshes $(\mathcal {M}_n)_{n\in \mathbb {N}}$ , with mesh sizes $r=r_n$ , and equip the underlying graphs $\mathcal {G}_n$ with a system of vertex and edge weights. Let $f \colon M \to N$ be any smooth map to another Riemannian manifold.

1. (1) If $(\mathcal {G}_n)_{n\in \mathbb {N}}$ is asymptotically Laplacian, then $\left \Vert \; \tau (f) - \tau _{\mathcal {G}_n}(f_n) \; \right \Vert {}_{\infty } = O\left (r^2 \right )$ . A fortiori,

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

2. (2) If $(\mathcal {G}_n)_{n\in \mathbb {N}}$ is almost asymptotically Laplacian, then

   (4.7) $$ \begin{align} \left\Vert \tau(f) - \tau_{\mathcal{G}_n}(f_n) \right\Vert{}_2 = O\left(r\right). \end{align} $$
   Furthermore, if $\vec {V} \in \operatorname {\mathrm {T}}_{f_n} \operatorname {\mathrm {Map}}_{\mathcal {G}_n}(M, N)$ is a unit tangent vector such that $\Vert \vec {V} \Vert _{\mathcal {V}_n^{(2)}} = o(1)$ , then
   (4.8) $$ \begin{align} \left \langle \; \tau(f) - \tau_{\mathcal{G}_n}(f_n ) \; , \vec{V} \; \right \rangle = o(r). \end{align} $$

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.26 We keep the setup of Theorem 4.25.

1. (1) If $(\mathcal {G}_n)_{n\in \mathbb {N}}$ is Laplacian or asymptotically Laplacian, then

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

2. (2) If $(\mathcal {G}_n)_{n\in \mathbb {N}}$ is almost asymptotically Laplacian, with decomposition $\mathcal {V}_n = \bigsqcup _{k=0}^2 \mathcal {V}_n^{(k)}$ , then

   $$ \begin{align*} \left| e(f)(x) - e_{\mathcal{G}_n}(f_n)(x) \right|= O\left(r^{2-k}\right) \end{align*} $$
   for every $x\in \mathcal {V}_n^{(k)}$ .

Proof The proof is easily adapted from the proof of Theorem 4.17.▪

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*} $$

Remark 4.28 Of course, Theorem 4.27 also holds for Laplacian and asymptotically Laplacian sequences of meshes, given the hierarchy between these conditions.

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 Convergence to smooth harmonic maps

Let $(M,g)$ be a compact Riemannian manifold and let $(N,h)$ be a Riemannian manifold of nonpositive sectional curvature. which does not contain any flats (totally geodesic flat submanifolds). Consider a connected component $\mathcal {C}$ of the space of smooth maps $\mathcal {C}^{\infty }(M,N)$ that does not contain any map of rank everywhere $\leqslant 1$ . For instance, take any connected component of maps whose topological degree is nonzero when $\dim M = \dim N$ . When N is compact, a celebrated theorem of Eells-Sampson implies that $\mathcal {C}$ contains a harmonic map w [7], and by Hartman [9], the harmonic map w is unique.

In this section, we show that one can obtain the harmonic map $w \in \mathcal {C}$ as the limit of discrete harmonic maps $u_n$ along a sequence of meshes $(\mathcal {M}_n)_{n\in \mathbb {N}}$ , provided that:

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

2. (ii) The discrete energy functional $E_n$ is sufficiently convex on the discrete homotopy class $\mathcal {C}_n$ . We expect that this is the case when N is compact and has negative sectional curvature, and have showed it in the two-dimensional case in our previous work [8] (see Section 5.6).

3. (iii) The sequence of meshes is either Laplacian (Definition 4.3), or one of the weaker versions (Definitions 4.5 and 4.6). In the next and final section Section 6, we systematically construct such sequences.

We then show convergence of the discrete heat flow $u_{kn}$ to the smooth harmonic map w, when the time and space discretization indices k and n simultaneously run to $+\infty $ , provided the adequate CFL condition is satisfied (see Section 5.5).

5.1 Strong convexity of the discrete energy

Please refer to [8, 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 [8, 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 [8, 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.

Remark 5.2 Under the assumptions of Theorem 5.1, w must be the unique smooth harmonic map in $\mathcal {C}$ , the minimizer of the energy functional.

Remark 5.3 The case $c=1$ and $\dim M=2$ is especially salient in light of [8], which guarantees that $c=1$ does hold when $\dim M=2$ in a broad setting: see Section 5.6 for details.

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.

Remark 5.9 We believe that the restriction $\dim M = 2$ when $c = 1$ is superfluous. Indeed, we expect that Lemma 5.8 is true in any dimension. However, proving it requires generalizations of the technical estimates of Appendix B when $\dim M> 2$ . We reserve this (possibly) for a future paper, as well as discussing cotangent weights and the constructions of Section 6 to dimensions $>2$ .

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.

Theorem 5.10 In the setup of Theorem 5.1, if $\dim M=2$ and $c=0$ , then $\widehat {v_n} \to w$ in $\operatorname {\mathrm {L}}^{\infty }(M,N)$ .

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)$ .▪

Remark 5.11 We believe that $c=0$ holds in great generality (see Section 5.1).

Remark 5.12 We believe that the restriction $\dim M =2$ (also possibly $c=0$ ) is superfluous, but are unable to omit it in the current stage of our work. See Remark 5.13 for a related discussion.

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.

Remark 5.13 The obstacle to show that $E_n(v_n)$ and $E(\widehat {v_n})$ are asymptotic would be lifted by showing that the sequence $(\widehat {v_n})_{n\in \mathbb {N}}$ has a uniformly bounded Lipschitz constant, but this would be a very strong result. It would in fact enable us to prove Theorem 5.1 for any asymptotically Laplacian sequence of meshes, with no assumption involving c, with a completely different method involving a Rellich–Kondrachov theorem. In the smooth setting, a uniform Lipschitz bound is achieved by using the Bochner formula and Moser’s Harnack inequality (see, e.g., [10] and [12, Section 2.2]). This is an essential feature of the heat flow and the theory of harmonic maps. While developing a discrete Bochner formula and a discrete Moser’s Harnack inequality is certainly a worthwhile project, it is also beyond the scope of this paper.

Theorem 5.14 In the setup of Theorem 5.1, if $c<2$ , then $E_n(v_n) \to E(w)$ . If moreover $\dim M = 2$ , $c = 0$ , and the sequence of meshes is asymptotically Laplacian, then we also have $E(\widehat {v_n}) \to E(w)$ .

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 [8], 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 [8] 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} $$

Remark 5.16 The assumption on the upper bound of the Hessian is reasonable when compared to the Euclidean setting due to scaling considerations. When N is a hyperbolic surface, we have $\beta _{n,K} = O(r^{-2})$ by [8, Proposition 3.17], which satisfies the assumption but is surely not optimal.

Remark 5.17 The CFL condition (5.6) is most likely far from optimal.

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 [8, 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.

Remark 5.20 We could similarly show $\operatorname {\mathrm {L}}^{\infty }$ convergence (resp. convergence of the energy) of $\widehat {u_{k,n}}$ to w under the assumptions of Theorem 5.10 (resp. Theorem 5.14) and suitable CFL conditions.

5.6 Application to surfaces

When M and N are both two-dimensional, our previous work [8] 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*} $$

Proof The estimate for $\alpha _n$ is an immediate consequence of [8, Theorem 3.20]. The estimate for $\beta _n$ is an immediate consequence of [8, Proposition 3.17]. Note that [8, Proposition 3.17] is only stated for a hyperbolic metric, but it can be extended to any Riemannian metric of curvature bounded below, which is always the case on a compact manifold.▪

Remark 5.22 The estimate $\alpha _n = \Omega (r)$ based on [8, Theorem 3.20] only assumes that N has nonpositive sectional curvature. When N has negative curvature (bounded away from zero by compactness), we expect that a better bound $\alpha _n = \Omega (r^c)$ with $c<1$ is possible to achieve, in fact we conjecture that $\alpha _n = \Omega (1)$ .

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

Theorem 5.23 If the sequence of meshes $(\mathcal {M}_n)_{n\in \mathbb {N}}$ is fine, crystalline, and almost asymptotically Laplacian, then the sequence of interpolations $(\widehat {v_n})_{n\in \mathbb {N}}$ of the discrete harmonic maps $(v_n)$ converges to the unique harmonic map $w \in \mathcal {C}$ in $\operatorname {\mathrm {L}}^2(M,N)$ , and $E(w) = \lim _{n \to +\infty } E_n(v_n)$ .

Furthermore, the discrete heat flow $(\widehat {u_{k,n}})_{k, n \in \mathbb {N}}$ from any initial condition $u \in \mathcal {C}$ converges to w in $\operatorname {\mathrm {L}}^2(M,N)$ when both $k, n \to +\infty $ , provided the CFL condition $k = \Omega \left (\log (r^{-1}) r^{-3} \right )$ holds.

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:

Theorem 5.24 In the setup of Theorem 5.23, assuming $\alpha _n = \Omega (1)$ , the convergence of $\widehat {v_n}$ to w is uniform. If, moreover, the sequence of meshes is asymptotically Laplacian, then we also have $E(w) = \lim _{n \to +\infty } E(\widehat {v_n})$ .

6 Construction of Laplacian sequences

Most of our convergence theorems in Sections 4 and 5 require a Laplacian sequence of meshes (as in Definition 4.3), or one of the weaker variants (Definition 4.5, Definition 4.6). Indeed, one should only expect convergence for weighted graphs that reasonably capture the geometry of M.

In this section, we construct a sequence of weighted meshes on any Riemannian surface and prove that it is always almost asymptotically Laplacian, and discuss cases where more can be said. This construction is very explicit: in fact, it is implemented in our software Harmony in the case of hyperbolic surfaces. The construction can simply be described: take a fine, crystalline sequence of meshes obtained by midpoint subdivision (Section 2.2) and equip it with the volume vertex weights (Section 3.2) and the cotangent weights (Section 3.3).

Remark 6.1 It is possible to generalize this construction to higher-dimensional manifolds, most likely with similar results. We reserve this analysis maybe as part of a future paper. In Euclidean space, the formula for higher-dimensional cotangent weights is given in [4].

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).

Remark 6.2 Finding an initial triangulation of S that is acute is far from an easy task, even for a flat surface. The reader may refer to [15] for more background on this active subject.

Definition 6.1 A $\Delta $ -sequence is a sequence of meshes $(\mathcal {M}_n)_{n\in N}$ with the associated biweighted graphs $(\mathcal {G}_n)_{n\in \mathbb {N}}$ constructed as above.

Remark 6.3 We think of “ $\Delta $ ” here as standing for either “Laplacian” or “simplex.”

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 [3], 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 [8, 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.

Theorem 6.4 Let $(M,g)$ be a compact Riemannian manifold. Let $\delta>0$ . There is an $R\geqslant 0$ so that the following holds: If $\mathcal {M}_0$ is a triangulation of M with largest side length $\leqslant R$ and whose angles are all between $\delta $ and $\leqslant \frac {\pi }{2} - \delta $ , then the sequence $(\mathcal {M}_n)_{n\in \mathbb {N}}$ of iterated subdivisions of $\mathcal {M}$ is fine, crystalline, and with angles bounded uniformly bounded away from $\frac \pi 2$ .

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.

Corollary 6.5 Let $(\mathcal {M}_n)_{n\ \in \mathbb {N}}$ be a strongly acute $\Delta $ -sequence in $(S,g)$ . Then all edge weights of $\mathcal {G}_n$ are $\Theta (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.

Definition 6.2 Recall that $\mathcal {V}_n \subseteq S$ denotes the set of vertices of $\mathcal {G}_n$ . Consider the decomposition $\mathcal {V}_n = \mathcal {V}_n^{(0)} \sqcup \mathcal {V}_n^{(1)} \sqcup \mathcal {V}_n^{(2)}$ , where:

• • $\mathcal {V}_n^{(2)}$ consists of the vertices that are also elements of $\mathcal {V}_0$ , called initial vertices.

• • $\mathcal {V}_n^{(1)}$ consists of the vertices that are located on the edges of the initial triangulation $\mathcal {M}_0$ , and are not elements of $\mathcal {V}_n^{(2)}$ , called boundary vertices.

• • $\mathcal {V}_n^{(0)}$ consists of all other vertices, called interior vertices.

Lemma 6.6 We have $\mu \left (\mathcal {V}_n^{(k)}\right ) = \Theta (r^k)$ for $k \in \{0, 1, 2\}$ .

Proof The cardinal $\big | \mathcal {V}_n^{(2)} \big |$ is clearly constant, while it is easy to show by induction that $\big | \mathcal {V}_n^{(1)} \big | = \Theta (2^n)$ and $\big | \mathcal {V}_n^{(0)} \big | = \Theta (4^n)$ . We also have $r_n = \Theta (2^{-n})$ by Proposition A.15 and $\mu _{x, n} = \Theta \left (r_n^2\right )$ for any $x \in \mathcal {V}_n$ by Theorem 4.4. The desired estimates follow.▪

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:

Theorem 6.7 Any strongly acute $\Delta $ -sequence in a closed Riemannian surface $(S,g)$ is almost asymptotically Laplacian.

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.

Remark 6.14 In retrospect, it is remarkable—almost miraculous—how the conditions for a $\Delta $ -sequence to be almost asymptotically Laplacian are barely met, and in turn how these conditions are barely sufficient for the main convergence theorem (Theorem 5.1) to hold, at least in the $c=1$ case. Seeing how delicate the analysis is, the reader should not be too surprised that it took us many failed attempts until we were able to achieve the right definitions and results.