Asymptotics of the Determinant of Discrete Laplacians on Triangulated and Quadrangulated Surfaces

Consider a surface \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document}Ω with a boundary obtained by gluing together a finite number of equilateral triangles, or squares, along their boundaries, equipped with a vector bundle with a flat unitary connection. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega ^{\delta }$$\end{document}Ωδ be a discretization of this surface, in which each triangle or square is discretized by a bi-periodic lattice of mesh size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta $$\end{document}δ, possessing enough symmetries so that these discretizations can be glued together seamlessly. We show that the logarithm of the product of non-zero eigenvalues of the discrete Laplacian acting on the sections of the bundle is asymptotic to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} A|\Omega ^{\delta }|+B|\partial \Omega ^{\delta }|+C\log \delta +D+o(1). \end{aligned}$$\end{document}A|Ωδ|+B|∂Ωδ|+Clogδ+D+o(1).Here A and B are constants that depend only on the lattice, C is an explicit constant depending on the bundle, the angles at conical singularities and at corners of the boundary, and D is a sum of lattice-dependent contributions from singularities and a universal term that can be interpreted as a zeta-regularization of the determinant of the continuum Laplacian acting on the sections of the bundle. We allow for Dirichlet or Neumann boundary conditions, or mixtures thereof. Our proof is based on an integral formula for the determinant in terms of theta function, and the functional Central limit theorem.


Introduction
Let be a connected surface, possibly with boundary, obtained by gluing finitely many equal equilateral triangles, or squares, along their boundaries. Thus, may have conical singularities and piece-wise straight boundary with corners; the cone and wedge angles either all belong to π 3 k, or to π 2 k, k ∈ N; we refer to the former situation as a triangulation and the latter one as a quadrangulation. The boundary of will be decomposed into two parts, ∂ D and ∂ N , that will carry Dirichlet and Neumann boundary conditions respectively; we assume that there are finitely many points separating the two. We assume that is equipped with a finite rank vector bundle with a unitary flat connection ϕ. We Fig. 1. An example of a surface , glued of five equilateral triangles, and its discretization δ by triangular lattice. In this case, has one conical singularity of angle 5π 3 and five boundary corners of angle 2π 3 . Note that the discrete triangles are glued so that the local graph structure at the edges is no different from that in the bulk of the triangles furthermore allow to have finitely many punctures, distinct from the set of conical singularities, and allow ϕ to have monodromy around those punctures. Our main results and techniques are new already in the case of a trivial line bundle, and the reader interested in the simplest situation may think of this case.
We summarize our setup briefly here, referring to Sect. 2 for the detailed definitions. By a lattice, we mean an (infinite) undirected planar graph with non-negative weights on edges embedded bi-periodically in the plane C. Given a lattice which is symmetric under reflections z → z and under rotation by π 2 or π 3 (such as, e.g., the square lattice in the case of quadrangulations and the triangular lattice in the case of triangulations), we can discretize by this lattice scaled to have small mesh δ. We denote the discretized surface by δ , see Fig. 1. In what follows, by a lattice-dependent constant, we mean a quantity that does not depend on δ, but may depend on the underlying lattice and weights (and, in the case of D p , on local geometric data); a lattice-independent, or universal quantity, is allowed to depend only on , the boundary conditions, and the connection ϕ. The parallel transport of ϕ along the edges gives a discrete connection, and we can consider the corresponding discrete Laplacian δ ,ϕ , with Dirichlet/Neumann boundary conditions approximating those in . We will denote the rank of ϕ by d. The subject of the present paper is the asymptotics of det δ ,ϕ as δ → 0, where det stands for the product of all non-zero eigenvalues. We make the following conventions regarding the parameters, see Sect. 2 for details. The mesh size δ is defined so that the lattice has δ −2 vertices per unit area. Thus, for a given lattice, δ = δ 0 N for a constant δ 0 and an integer N . We denote by | δ | the size of δ , measured in the number of fundamental "plaquettes" (triangles or squares) in their discretization; thus as δ → 0, | δ | is of order δ −2 . Similarly, |∂ D δ |, |∂ N δ | denote the size of Dirichlet and Neumann parts of the boundary, measured in the number of plaquette sides, thus |∂ D δ |, |∂ N δ | have order δ −1 . We will also assume that the weights on the graph are normalized so that the random walk on δ converges to the standard Brownian motion.
• C is a lattice-independent constant given by where ,ϕ stands for the Friedrichs extension the Laplacian on acting on sections of ϕ, see Sect. 3; the sum is over the set of conical singularities, corners, and punctures of , and the values C p for a cone C α of angle α, corners ϒ α D , ϒ α N , ϒ N D of angle α with Dirichlet, Neumann, changing Neumann-to-Dirichlet boundary conditions, and a puncture P M with monodromy operator M are given by • the constant D has the form where D p are lattice-dependent constants (entirely determined by the lattice and angle, boundary conditions, or the monodromy at a puncture p); and det ζ ,ϕ is the zeta-regularized determinant of ,ϕ , see Sect. 3.
The constants A, B D , B N depend linearly on the rank d of the bundle, i.e., when divided by d, they only depend on the underlying lattice and the weights.
There is a lot of previous work on the subject. Discretizations of Laplacian were studied by Dodziuk [10] on a class of discretizations of Riemannian manifolds in arbitrary dimension. He established [10, Section 5] the convergence of the eigenvalues and the spectral zeta function of a combinatorial Laplacian with weights inherited from the Riemannian metric on the underlying manifold. Cardy and Peschel [5] conjectured that the asymptotic of the partition function of any critical 2D lattice model on a Riemann surface should take a form similar to (1.1). Since the determinants of the discrete Laplacian and its vector bundle versions are partition functions of a number of lattice models, such as dimers and double dimers, discrete GFF, spanning trees and cycle-rooted spanning forests [3,11,24,25,27,28,30,31,37], our results can be viewed as a rigorous proof of a particular case of the Cardy-Peschel conjecture. Duplantier and David [13] computed the asymptotics of the determinant of the discrete square lattice Laplacian on a torus and a rectangle; their results were extended to cylinder, Möbius strip and Klein bottle by Brankov-Priezzhev and Izmailyan-Oganesyan-Hu [4,23]. The approach in these papers is based on the fact that these geometries have large groups of symmetries acting on them, and hence the discrete Laplacian can be diagonalized explicitly. The analyze separately, see the key formula (4.1). The main idea is to regularize the discrete heat kernel by subtracting the heat kernel on one of a discretized model surface-the full plane, the half-plane, a punctured plane, a cone, or a corner-that matches the geometry of δ locally. This immediately isolates the volume term A · | δ | in (1.1), and also the boundary, cone and puncture contributions, whose asymptotics can be analyzed separately by studying heat kernels on model surfaces. On the other hand, the regularized trace of the heat kernel receives only Brownian scale contributions, and after rescaling converges to its continuous counterpart, more or less, by the local Central limit theorem. We derive this convergence from the functional CLT in the plane and the parabolic Harnack inequality of Delmotte. To complete the asymptotic analysis, we need to pass to the limit under the integral; to this end, we employ a large deviation estimate for small t and a uniform spectral gap bound for large t.
Thus, essentially, we only use three ingredients: the functional Central limit theorem, parabolic regularity, and the fact that microscopically, our lattice approximation "looks the same at all places". The only reason we do not consider more general surfaces (e.g., allowing for conical singularities with arbitrary angles, or for genuinely curved surfaces) is that those do not admit nice discretizations with this last property; see Remark 9.2. To see the difficulty it entails, note that the constant A, being essentially the free energy per lattice site in the infinite volume limit, is lattice-dependent; if the lattice is different in different places of δ , then A will also fluctuate, and the volume term in the asymptotics may be hard to control with meaningful precision. This is why we believe that approximation schemes such as Dodziuk's are too general to admit asymptotic formulae like (1.1). On the positive side, apart from the volume term (4.3), the rest of our analysis, at least in the absence of the boundary, does not use regularity of the lattice in any essential way. This leaves hope that the method can be applied to more general surfaces and approximations in sufficiently integrable case, e.g. on isoradial graphs.
Note that if we only want to study the the difference log det δ ,ϕ 1 −log det δ ,ϕ 2 , the aforementioned volume term cancels out, as do most of other terms in the asymptotics, and and our method yields a simple proof the asymptotics 2(dim ker ϕ 1 −dim ker ϕ 2 ) log δ+ log det ζ ,ϕ 1 − log det ζ ,ϕ 2 , using only the functional CLT and the parabolic regularity. Such differences are of interest since they often compute interesting topological observables in the models, see e.g. [11,12,[24][25][26], and there are a number of results in this direction. Dubédat and Gheissari [12] proved convergence for tori in a different way under even weaker assumptions and Kassel-Kenyon [25] proved a similar result on general Riemann surfaces, without identifying the limit. A remark on cancellation of the singular terms in this setting is also made by Kenyon [28], for trivial line bundles with different simply-connected domains, and Finski [18], for general bundles.
Some of the ingredients of our proof could be alternatively established by the methods [17,18]; e.g., the convergence of theta function in the Brownian scale follows easily from the convergence of rescaled eigenvalues and the uniform Weyl law, as does the spectral gap estimate of Lemma 5.5. Arguably, our proofs via the functional CLT are more streamlined; for example, cones and conical singularities, that seem to be a source of technical difficulties in [17], do not enter this part of our proof at all, since the Brownian motion avoids them almost surely. On the other hand, some intermediate results and techniques in [17,18] are of independent interest.
Our method gives the explicit values of A, B D = −B N and D p as sums of integrals of discrete (continuous time) heat kernels on the discretizations of the plane, half-plane, corners, or cones. Note, however, that to determine those constant, it suffices to compute the asymptotics (1.1) in some simple geometry; for example, a torus for A, or a cylinder or a square (with corresponding boundary conditions) for B D , B N . Thus, on the square lattice, it follows from [13,Eq. (4.24) and (5.14)] that where G is the Catalan's constant 1 − 1 3 2 + 1 5 2 − . . . Note that [13] computes B D and B + N on two different discretizations of a square (see Fig. 4), shifted with respect to each other by δ 2 + i δ 2 , and B D is related to B + N by the Uniform Spanning tree duality. By contrast, our identity B D = −B N holds for a fixed discretization. An expression for A in terms of polylogarithms is known, by a different method, for isoradial graphs with critical weights, including the hexagonal and the triangular lattices [29, Theorem 1.1]. In Sect. 8, we compute closed-form expressions for B D = −B N for triangular lattice (rotated in two different ways), for the square lattice rotated by 45 • , and give an alternative computation in the case of a non-rotated square lattice, recovering the result of [13]. Note that the specific values of the constants are sensitive to conventions such as the definition of | δ | and the normalization of weights, hence e.g. the extra − log 2 in A above compared to [13,18,28].
In our proof, we employ the functional Central limit theorem with an error bound that we derive from Einmahl's multidimensional KMT coupling [14]. The only place the error bound is used in earnest is when we sharpen the asymptotics of the contribution of cones and corners to (1.1), from C p log δ + o(log δ) in [18,28] to C p log δ + D p + o (1). (We also use the error bound it in the computation of the boundary contribution; however, there we do not need the functional CLT, so e.g. Berry-Esseen bounds would suffice.) Since any power law error bound suffices for that, one could use e.g. the Skorokhod embedding instead, cf. [34,Theorem 3.4.2]; see also [34,Section 7.2] for simple proofs of the KMT coupling on the square and the triangular lattices. The same error bound allows one to track the rate of convergence in Lemma 5.6, Corollary 6.4, and (6.2), (6.6), improving the o(1) in (1.1) to O(δ ρ ) with ρ > 0. The recent independent work [21] gives explicitly the error bound on the square lattice; it also shows that similar methods allow one to treat polygonal boundaries with any rational slopes.
The assumption that the weights are symmetric is only used in the proofs of technical Lemmas 5.1-5.6, where we found it convenient to use the parabolic Harnack inequality of Delmotte [9]; we note that we use the continuous time version that is significantly simpler than the discrete time one. With some work, these lemmas, modified accordingly, can be given alternative proofs, allowing one to lift the symmetry assumption. The same applies to the unitarity of ϕ, which can probably be relaxed under the assumption that the real parts of the eigenvalues of the Laplacian remain positive.

Laplacians, Heat Kernels and Zeta Functions
Although we are mainly interested in the discretizations of the Riemann surfaces, we start from a more general setup. Let G = (V, E) be a connected finite undirected graph with weights w e > 0 assigned to its edges. A vector bundle of rank d over G is a collection of d-dimensional complex vector spaces V x attached to its vertices x ∈ V. A connection on a vector bundle is a collection of linear isomorphisms ϕ xy = ϕ −1 yx : V x → V y for each pair x ∼ y of adjacent vertices. From now on, we will only consider unitary connections, that is, each V x is equipped with an inner product and the maps ϕ xy are unitary. A section f of a vector bundle is a choice of an element f (x) ∈ V x for each x ∈ V; we will take the liberty to refer to it as a section of ϕ in order to lighten the notation. When the vector bundle is unitary, the linear space of it sections comes with the inner product f ; g = x∈V f (x); g(x) .
The Laplace operator G,ϕ acts on sections of the vector bundles by the formula Note that with respect to the inner product as above, we have which shows that G,ϕ is non-negative and self-adjoint and therefore diagonalizable with real non-negative eigenvalues. We define the heat operator associated to G,ϕ by This is again a linear operator acting on the linear space of sections of the vector bundle. By linearity, we can write (P The trace of P G,ϕ t is called the theta function. By computing the trace first as the sum of eigenvalues, and then as the sum of the diagonal elements, we get the theta inversion identity (2.1) The zeta function associated to G,ϕ is the Mellin transofrm of G,ϕ (t), defined for s ∈ C as Let k = dim ker G,ϕ . By subtracting k from both sides of the theta inversion formula, multiplying by t s−1 and integrating, one gets note that the integral converges at infinity because of the positivity of the eigenvalues. Our main motivation for studying the zeta function is the identity We will need a probabilistic interpretation of the heat kernel. Let γ t denote the continuous time random walk on the weighted graph (G; w), that, being at x ∈ V, moves following exponential clock to an adjacent vertex y with intensity w xy . We denote ϕ γ [0;t] := ϕ x n−1 x n • · · · • ϕ x 0 x 1 , where γ 0 = x 0 ∼ x 1 ∼ · · · ∼ x n = γ t are the vertices visited consecutively by γ up to time t. We denote by P x and E x the probability and the expectation with respect to this random walk started at x. We have the following lemma: Proof. Denote the right-hand side byP G,ϕ (x, y, t). We observe that ∂ tP G,ϕ (x, y, t) + G,ϕ xP G,ϕ (x, y, t) = 0, and the same system of ODEs with the same initial conditions holds for P G,ϕ (x, y, t), and the solution is unique.
We use (2.4) to extend the definition of the heat kernel to infinite graphs.
In what follows, the graph G will be a discretization of a triangulated or a quadrangulated surface , as in the introduction. To discretize , choose an infinite locally finite weighted connected graph C δ 0 with vertices embedded in the plane; we assume that the embedded weighted graph is bi-periodic with periods either 1 and 1 2 + √ 3 2 i (for triangulation case), or 1 and i (for quadrangulation case) and has δ −2 0 vertices per unit area. We moreover assume that C δ 0 is preserved under rotations by π/3, respectively, π/2, around the origin, and under reflections with respect to the real line. We denote 2 i} (respectively, the unit square {0, 1, 1 + i, i}). We denote by δ the discrete surface obtained by discretizing each triangle/square in with T δ = C δ ∩ T ; since C δ has all the symmetries of C δ 0 , these discretizations can be naturally glued together. We will interchangeably use δ and N = δ 0 δ −1 as mesh parameters of the discretization. We do not require C δ 0 to be properly embedded, i.e., edges are allowed to intersect; however, we do assume that the graph obtained by removing edges connecting a vertex strictly inside T with one strictly outside T is still connected, so that δ is connected.
As an example, Z + iZ and 1 2 + i 2 + Z + iZ can serve as C δ 0 (with δ 0 = 1 in this case) in the quadrangulated case, leading to two different families of discretizations (one of them will have vertices at cone tips). The square lattice rotated by 45 • yields another discretization. The triangular and the hexagonal lattices can serve as C δ 0 in the triangulation case, with δ 0 = It is easy to show (see Lemma 10.2 below) that the continuous time random walk on C δ satisfies the Central limit theorem, i.e., converges, as δ → 0, to a Brownian motion, whose covariance matrix must be scalar, because of the symmetries. We will assume that the weights w xy are chosen so that this matrix is the identity, i.e., 2t )dy for any disc B. This can always be achieved by simultaneously multiplying all the weights by a common factor, which only affects the values of the lattice-dependent constants in (1.1).
Let us comment on the boundary conditions. We will assume that if a point p ∈ ∂ of Dirichlet-to-Neumann change is at a corner, then, at δ , the boundary conditions also change at the corresponding corner, and if p is an inner point of a side of one of the triangles/squares, then it is approximated by sequence of points p δ → p at distance m δ /N , m δ ∈ N, from the corner of that triangle (or square). To define an action of δ ,ϕ on a section f with Dirichlet (respectively, Neumann) boundary condition at a boundary segment l, we extend f across l by f (z ) = ∓ f (z), where z → z is the reflection with respect to l; if there are vertices on the Dirichlet boundary, we only consider sections f that are zero at those vertices. This procedure may lead to non-symmetric weights that are gauge equivalent to symmetric ones; the above formulae are not affected. We adopt the convention that when writing the sum over x ∈ δ (as e.g. in (2.1)) we do not include vertices lying on the Dirichlet part of the boundary, but do include those on the Neumann part.
An approximation p δ to a puncture p will be realized as a point of disjoint from any edge of δ ; we moreover insist that it is an image under scaling of a fixed point in a fundamental domain of C δ 0 .

Heat Kernels and Theta Functions in the Continuum
We recall a continuous version of the above theory, in which the random walk γ δ on δ is replaced by the Brownian motion γ t on , reflected at ∂ N and absorbed at ∂ D . This is only needed for the interpretation of the constant term in the asymptotics (1.1); a reader willing to accept (3.3-3.5) as the definition of the zeta-reguralized determinant of the Laplacian may skip most of the rest. Recall [38] that the Laplacian on surfaces with conical singularities is not essentially self-adjoint; hence some care is needed when specifying its self-adjoint extension, see also [17,Section 2.3].
We start by constructing the Brownian motion on , which can be done by elementary means. Letˆ be two copies of glued along the boundary with conical singularities removed; we can define the Brownian motionγ onˆ by coupling it to the Brownian motionγ in the plane, liftingγ toˆ by local isometries. The lifting is well defined at least up to the first timeγ hits a conical singularity ofˆ , but this can only happen wheň γ hits a point of a triangular or a square lattice, i.e., with probability 0. Hence, almost surely for Leb(ˆ ) a.e. starting point,γ t is defined for all t. We then define γ t to beγ t reflected to and stopped upon hitting ∂ D .
Lemma 3.1. The Markov process γ t is symmetric with respect to the Lebesgue measure on . Moreover, for each fixed t > 0, its transition kernel P (·, ·, t) is bounded, and it is smooth away from ∂ and the cone tips.
The following construction is general for symmetric Markov processes, see [45,Section 1.4] or [20,Section 1.4]. Let L 2,ϕ ( ) be the set of L 2 sections of ϕ. Define the Markov semi-group P ,ϕ t From the symmetry of γ t and the unitarity of ϕ, we see that P ,ϕ t are self-adjoint; also, P ,ϕ ≤ P pointwise. Hence, by Cauchy-Schwarz, i.e., for each t > 0, P ,ϕ t is a self-adjoint contraction on L 2,ϕ ( ). The semi-group P ,ϕ t is strongly continuous: it is enough to check that P ,ϕ t f t→0 −→ f in L 2 for a continuous f , in which case the convergence clearly holds uniformly. Put D The form E ,ϕ is non-negative and closed [45,Lemma 4.2], and we denote by ,ϕ the unique non-negative self-adjoint operator associated to E ,ϕ , called the generator of {P ,ϕ t } t>0 . In terms of the common spectral projections E λ for P is Hilbert-Schmidt, and hence P ,ϕ 2t is trace class for all t > 0. Therefore, we can define the spectral theta function by and, integrating, we arrive at the continuous theta inversion identity, We define the spectral zeta function by postponing the proof of convergence to Sect. 9. Here k = dim ker ,ϕ ; we can also describe k explicitly (for a connected ) as follows: if ∂ D = ∅, then k = 0, otherwise, k is the dimension of the maximal trivial sub-bundle. This description follows from Lemma 5.5, which uses the latter definition k, and whose conclusion, together with (5.4), implies that both δ (t) − k and (t) − k tend to zero as t → ∞. In particular, k = dim ker ,ϕ = dim ker δ ,ϕ = k δ . Alternatively, one could note from Lemma A.1 that, if ψ ∈ ker ,ϕ , then ∇ ϕ ψ, ∇ ϕ ψ = 0, that is, ψ is a covariant constant, also known as a flat section, cf. [17, Corollary 2.6].
As in the discrete case, we have ζ ,ϕ (s) = i λ −s i whenever either the series or the integral (3.4) converges absolutely. Moreover, in in Sect. 9, we analytically continue x is a plane, half-plane, a infinite cone or an infinite wedge, depending on the local geometryof δ near x ζ ,ϕ into a neighborhood of the origin, cf. [6]. In analogy with (2.3), this allows one to define the zeta-regularized determinant of ,ϕ (cf. [41,42]) by It is possible to describe the form E ,ϕ and its domain D(E ,ϕ ) (or, equivalently, the generator ,ϕ ) more explicitly, which we postpone to the Appendix.

The Key Formula for the Determinant of the Discrete Laplacian
For notational simplicity, we first assume there are no punctures, and also that the lattice is such that there are no vertices at conical singularities. We will then discuss the necessary modifications in the general case. Given x ∈ , we define x to be one of the model surfaces, namely, a plane, a halfplane, an infinite cone, or an infinite wedge, that agrees with locally near x, see Fig. 2. More precisely, fix a small r > 0 in such a way that the 2r -neighborhoods of the tips of the conical singularities and the boundary components do not overlap or self-overlap. We define r α := r/ sin(α/2) if α < π and r α = r otherwise. We then define x to be the wedge of angle α if x is at distance at most r α from the tip of a corner (see Fig. 3) with angle α, else, if x at distance ≤ r from ∂ (or a conical singularity), we define x to be the half-plane (respectively, the cone); else, x is a plane.
In a similar way (with the same r independent of δ), we define the infinite graphs δ x discretizing each respective x . These graphs, when they have a boundary, come equipped with boundary conditions inherited from δ . We define the heat kernel P in each of the discrete model domains δ x by (2.4), with ϕ the trivial connection on the rank one bundle, and the random walk being stopped at the Dirichlet boundary and reflected at the Neumann one.
We start from (2.2) and rewrite it as where w x := y∼x w xy . This identity is valid for 0 < e s < 1, since each integral defines an analytic function in that region (see Lemma 6.2 below for large t bounds). Moreover, all the integrals but the last two are in fact analytic at least in −1 < e s < 1; indeed it follows from (2. , where e −w x t is the probability that the random walk does not move at all by time t, and O(t 2 ) is the probability that it makes at least two steps, which is bounded from above by The last two terms are equal to This is our key formula; analyzing it term by term will lead to (1.1). The last two terms are already explicit, and we see the −2k = −2 dim ker ,ϕ contribution to the logarithmic term in (1.1). For the first three integrals, going back to the probabilistic interpretation of the heat kernel, we observe that only the walks with δ −1 steps contribute, hence, these terms converge to their continuous counterparts by Central limit theorem, see Sect. 5 and Corollary 6.4 for details. As for the fourth term, we note that the summands are zero unless x is r -close to a conical singularity or to the boundary. Thus, the whole sum only depends on the number of conical singularities and their angles, and on the geometry of the boundary. We treat it in Sect. 6.
Turning to the fifth term in (4.1), denote and observe that this quantity only depends on the vertex of T δ 0 corresponding to x in the discretization procedure, in particular in the vertex-transitive case this is just a constant. Assume first that there are no vertices of C δ 0 on ∂ T δ 0 . Then, subdividing the vertices of δ into scaled copies of T δ 0 , the fifth term above gives the leading term of the asymptotics (1.1): If there are vertices of T δ 0 on ∂ T but not in its corners, then the contribution of those should be included in the definition of A with weights 1 2 . This leads to a miscount for the contribution of the vertices on ∂ , which we absorb into the |∂ δ | term by re-defining the constants B N , B D , cf. the second term in (6.8) below. Similarly, if there are vertices at the corners of T δ 0 , they should be counted with weight 1 6 or 1 4 , which leads to a miscount for boundary corners and cones which we absorb into D p .
In the case there are punctures, in an r -neighborhood of a puncture p, we define x to be the punctured plane C δ \{ p}, equipped with the connection ϕ p obtained by first restricting ϕ to the neighborhood of p and then extending it to a flat connection on the whole C δ \ {p}. We then simply use P C δ \{ p},ϕ p instead of d · P δ x in the above formulae. If there are vertices at the corners of T (and thus at the conical singularities), then the asymptotics of the heat kernel at a conical singularity p of angle α reads P This results in additional constant contributions to the asymptotics that can be absorbed into D p .

Contributions from the CLT Part
The goal for this section is to prove convergence of the first two terms in the key formula. We start with five fairly standard Lemmas, whose proof is deferred to Sect. 10. We denote † := ∂ ∪ C ∪ P

Lemma 5.2 (Short time large diameter bound).
For every ε > 0, there are constants C, c > 0 such that, for all x ∈ δ , all t > 0 and all δ < c,

Lemma 5.4 (Hölder regularity of heat kernels).
There exists a number θ > 0 such that, for any η > 0, there exists a constant C η with the following property: if dist(x, y) < 1 2 dist(x, † ) < 1 2 η and t > η 2 , then where, to make sense of the left-hand side, we identify the vector spaces V y , y ∈ B(x, η 2 ), using a trivialization of ϕ over B(x, η 2 ). Lemma 5.5 (Spectral gap). There are constants C, c > 0 independent of δ, such that for all δ small enough, one has We are in the position to prove convergence of the first two terms in the key formula: Lemma 5.6. We have the following convergence results: Proof. In view of Lemma 5.5, for (5.4), it suffices to prove the convergence for any fixed T > 0. Let 0 < η < η 0 , and let {ψ j } be a partition of unity for such that diam(supp ψ j ) < η for all j. We write and split the sum according to whether j ∈ J : For j ∈ J , we may apply Lemma 5.4 to any x, y ∈ supp ψ j . Averaging (5.2) with weights , multiplying by t −1 , and integrating yields Summing these bounds over j ∈ J yields the upper bound C η 0 | δ |δ 2 η θ , which goes to zero uniformly in δ as η → 0. Finally, with the initial point X chosen at random with P(X = x) = ψ j (x)/S j , and we pick a trivialization of ϕ over supp ψ j . The expression inside the expectation is continuous with respect to the path γ δ , therefore, by Lemma 5.1, (5.6) converges to its continuous counterpart In view of the bounds we have collected, we have that Since the continuous heat kernels satisfy the suitable counterparts of Lemmas 5.3-5.4 (for instance, as a consequence of the discrete bounds and the convergence), an argument as above gives that the latter quantity is equal to as required. For (5.5), the same argument as above, applied to δ and each of δ x , gives the convergence of the integral from t 0 to 1, for any fixed t 0 > 0. Hence, it suffices to show that the integral from 0 to t 0 converges to 0 as t 0 → 0 uniformly in δ. We can write where γ andγ are random walks on δ and δ x , respectively, coupled in such a way that they coincide up until τ r := min{s : dist(γ δ −2 s , x) ≥ r }. On the event τ r ≥ t, the expression in the expectation is zero. Hence, we can write where we have used Lemma 5.3 and then Lemma 5.2. Summing over x gives a bound on the integrand in the left-hand side of (5.5) that is independent of δ and integrable at 0. This concludes the proof.

Local Contributions
In this section, we compute the asymptotics of the local term Each of the model surfaces = H, ϒ α , C α has scaling acting on it, and I . Thus, decreasing δ by going from N to N + 1 is tantamount to adding new terms to the sum, corresponding to those x whose distance to a conical singularity, the boundary, or a puncture is between r N N +1 and r . The asymptotics of those new terms is governed by Central limit theorem. We postpone the proof of the following Lemmas to Sect. 10: Lemma 6.1 (Local CLT with error bound). If is one of the model surfaces, δ its discretization, and ε > 0, then there exist q > 0 and C > 0 such that for all δ, all t ∈ (δ q , δ −q ) and x, y at distance at least ε from the tip (if is a wedge or a cone). Lemma 6.2 (Uniform tail bound for the heat kernel). If δ is one of the model surfaces, then there exists C > 0 such that, for any δ, t > 0 and x, y ∈ δ , Let 1 x 0 be a continuous model surface equipped with boundary conditions. Let 2 be another model surface that contains an isometric copy B (x 0 , η) of the ball B(x 0 , η) ⊂ 1 not containing tips of a wedge or a cone, with corresponding parts of the boundary having the same boundary conditions, and let δ 1,2 be their discretizations that respect the isometry. We will denote, for x ∈ B(x 0 , η), where we identify the points in B(x 0 , η) with their isomorphic copies. We use a similar notation for discretizations δ 1,2 of 1,2 We have the following Lemma: Corollary 6.3. In the above setup, there exist ρ > 0 and C > 0 such that for all δ and all x ∈ B(x 0 , η 2 ).
Proof. Let q be as in Lemma 6.1. At small times t ≤ δ q/2 , we repeat the argument in the end of the proof of Lemma 5.6 for d = 1 and δ 1,2 instead of δ , δ x ; this gives and, integrating, At large times t > δ −q , we use Lemma 6.2 to get Clearly, similar estimates, with right-hand side divided by δ 2 , hold for continuous heat kernels, e.g., as a consequence of convergence. At intermediate times δ q 2 < t < δ −q , we apply Lemma 6.1 to each of 1,2 separately to get Combining all the estimates above yields the result.
Proof. As in the proof of Lemma 5.6, we change the variable t → δ −2 t. We then use Lemma 5.3 ensure that the sum and the integral over the η 0 -neighborhood of the boundary and singularities tends to zero with η 0 , uniformly in δ. This allows to consider the sum over the complement of that neighborhood only. In that region, we argue as in the proof of Corollary 6.3 that the summand converges to the integrand uniformly: apply Lemma 6.1 on the integral from 1 to δ −q and Lemma 6.2 to the integral from δ −q to infinity. The only difference with the proof of Corollary 6.3 is that now we do not need to deal with small t.

Conical singularities.
Let us compute the contribution of an r -neighborhood of the tip of a conical singularity with angle α to (6.1). Changing the scale to δ 0 , we see that x∈C α,δ :|x|≤r where | · | denotes the distance to the tip; that is, decreasing δ for a fixed r simply results in adding new terms to the sum. The asymptotics of those terms as |x| → ∞ can be read off Corollary 6.3: if r N ≤ |x| < r (N + 1), then where 1 is any point at distance 1 from the tip, and we used rotational invariance of , and the relation δ N = δ 0 . The error term O(|x| −2−ρ ) sums to a constant over C α,δ 0 , and recalling that C δ 0 has δ −2 0 vertices per unit area, we have Taking into account that log(r N) = log r + log δ 0 − log δ, we conclude whereD α is a (lattice-dependent) constant.

Boundary segments.
Let l ⊂ ∂ be a side of a triangle or a square comprising ; we introduce a local coordinate in which l is identified with (0; 1) ⊂ ∂H. Let l δ be the corresponding segment of ∂ δ . Let α 0,1 be the angles of the wedges at its endpoints 0 and 1, and denoteα 0,1 := min{α 0,1 /2; π/2}. We consider the contribution to (6.1) of the points that are at distance at most r from l δ , but at the distance greater than r α 0,1 = r/ sin(α 0,1 ) from its endpoints 0 and 1, respectively. This contribution reads Fig. 3; the boundary conditions in H δ above are inherited from l δ . We first treat the sum over R δ r , which we can split into N = δ 0 δ −1 strips R δ r (k) := { m x ≤ r, k N ≤ Fig. 3. A decomposition of a neighborhood of a boundary segment. The rectangle R δ r includes two sectors Y δ 0,1 and a triangle δ 0 ; in this case, δ 1 = ∅ since the corresponding angle is greater than π . The shaded sectors of radii r α = r/ sin α 2 for α = 2π/3 and r α = r for α = 4π/3 > π are the regions for which δ x is a wedge; in the white part of R δ r , δ x = H δ e x < k+1 N } that all give equal contributions. As in the cone case, we see that decreasing δ is tantamount to adding new terms to the sum over R δ r (0), i.e., x∈R δ r (0) where we can compute, with the sign s = ±1 depending on the boundary conditions as s N = +1 and s D = −1, The contributions of Y δ 0,1 to (6.3) will cancel the corresponding contributions from corners, thus we will leave them as they are for a while. The contribution of δ 0,1 is computed as in the cone case, applying Corollary 6.3 and then using Brownian scaling and shift invariance of I H C (i) : and similarly y∈ δ 6.3. Boundary corners. We parameterize a boundary corner ϒ α by a local coordinate z so that ϒ α = {z ∈ C : 0 < argz < α}, and denote we denote by Y δ left (respectively, Y δ right ) the set Y δ 1 (respectively, Y δ 0 ) corresponding to the boundary segment adjacent to ϒ α on the left (respectively, on the right). We also denote Y δ which is non-empty if and only if α > π. The contribution of ϒ α to (6.1) can be written as whereĤ δ stands for the upper-half plane H δ rotated counterclockwise by α − π around π (so that its boundary coincides with the left boundary of the corner), and the boundary conditions ϒ α,δ , H δ ,Ĥ δ are inherited from those in δ . The first three terms yield, similarly to the computations above, Observe that when collecting the contributions to (6.1) along ∂ δ , the last two terms in (6.4) cancel out the corresponding terms in (6.3). The total contribution of the rneighborhood of ∂ δ to (6.1) is thereforê where C p depends only on the angle and the boundary conditions b,b ∈ {D, N } on the segments adjacent to the p ϒ α bb as .
Taking into account the discussion at the end of Sect. 4, we have where the sum is over the vertices in one boundary edge of the fundamental domain T δ 0 .

Punctures.
Using a suitable modification of Lemma 6.1, similarly to the conical singularity case, we have, in the local coordinate where p = 0,

Explicit Computations for the Logarithmic Term
In this section, we compute the integrals involving heat kernels that contribute to the logarithmic term of the asymptotics. The results are not new. Namely, as pointed out in [18,21], the constant C in (1.1) is related to the spectral zeta-function by C = −2ζ (0), see Remark 9.1. The value of ζ (0) has been computed in a much greater generality by Cheeger, see [6,Theorem 4.4] and the discussion thereafter.
Here, we propose an alternative computation based on the following identity for the heat kernel on the universal cover of a punctured plane: Lemma 7.1. We have, for the heat kernelP := P C\{0} , Proof. We use the Brownian loop measure of Lawler and Werner, see [35] or [33, Section 5.6], defined as a σ -finite measure on the space of unrooted closed loops in a Riemann surface by where P is the heat kernel in with Dirichlet boundary conditions, and μ ,z,t is the Brownian probability measure on the paths from z to z of duration t. We will only need the conformal invariance of this measure, see [35,Proposition 6] which is stated for planar domains, but the proof, being a local computation, extends verbatim to Riemann surfaces. Consider the annular region A r,α = {r ≤ |z| ≤ 1}/{z ∼ e iα z} in the cone of angle α,Ã r,α its universal cover, and let E denote the set of loops in A r,α that wind around the annulus once counterclockwise. The map φ : z → −i log z maps A r,α onto the cylinder O α,r = S r /{z ∼ z + α}, where S r = {0 ≤ m z ≤ − log r }, so, by the conformal invariance, we have Note that by scaling invariance, the total μ C α measure of the loops that wind around C α and intersect a given circle |z| = r does not depend on r . We claim that it is also finite. Indeed, by conformal invariance, we can pass to the cylinder O α = C/{z ∼ z + α}, when the circle is mapped to l h := {w : m w = h}, and then use that for some c, C > 0, These two bounds imply that the contribution to (7.1) from the points z with | m z − h| > 1 is finite. Since the probability that a bridge with a small t winds around O α is exponentially small, the region | m z − h| ≤ 1 also gives a finite contribution.
Hence, up to O(1) as r → 0, the left-hand side of (7.2) equals We conclude by comparing this to the right-hand side of (7.2), which is, up to O(1), Since C \ {0} also covers each of the cones C α C/{z ∼ e iα z}, we have We now can compute, using that ∞ k=1 1/k 2 = π 2 6 , By the reflection principle, we have P ϒ α (x, y, t) = P C 2α (x, y, t) ± P C 2α (x,ȳ, t) for ϒ α = ϒ α N and ϒ α = ϒ α D respectively, and hence, using that k∈Z (x + k) −2 = π 2 sin −2 π x, we get Using that, by reflection principle applied to the line argz = α, we have x, e 2iαȳ , t ; (7.5) a similar straightforward but tedious computation yields Now we are ready to collect the values of C ϒ α for α ≤ π :Ĉ DD α andĈ N N α consist of two equal terms given by (7.4), with the cotangent terms canceling out the corresponding terms in (6.7), whileĈ DN α is given by the sum of (7.7) and (7.8) above. If α > 2π , then we need, in addition, to compute the contribution of the third term in (6.5) which is done similarly, and change the integration limits in (7.4), (7.7-7.8) to π 2 ; we leave it to the reader to check that the answer is (unsurprisingly) given by the same analytic expression in α.
In the case of a puncture, we can compute

Explicit Computation of the Constant B for Some Lattices
In this section, we show how the exact values of B N = −B D can be computed for several nice lattices, see Figure 4. We assume d = 1, since, as remarked before, the constants A, B D , B N depend linearly on d. We note that the values B + N and B D are related to each other by planar UST duality, see [13,Section 5], therefore, in fact, each of B + N , B N , B + D , B D can be deduced from the results in [13]. Similarly, it can be deduced from duality that B ♦ N = B ♦ D , and hence B ♦ N = B ♦ D = 0. Here, we give a direct selfcontained computation.
For each of those lattices, the constantB N = −B D can be expressed in the form Let γ Z t be the continuous time symmetric random walk on Z that moves with intensity 1. Note that in the cases * = +, , ♦, , we have X t . In the cases * = + and * = , X t and Y t are independent. We compute P(γ Z t is even) = 1 2 (1 + e −2t ) =: K Z 2 (t) and P(γ Z t = 0) = e −t I 0 (t) =: K Z (t) where I 0 (t) = 1 π π 0 e t cos φ dφ is the modified Bessel function of the first kind (see [7] or the computation in the end of this section). Recall from Sect. 4 that in all cases, (8.2) which follows from the formula for the Laplace transform of the Bessel function, , by taking the derivative of (8.2) with respect to a and matching the behavior as a → ∞.
With these observations, we are ready to evaluate the right-hand side of (8.1) for * = +, , ♦, : and plugging into (6.8) gives In the case = , X t is equal in law to a scaled copy of a continuous time walk on Z jumping by −2, −1, 1, 2 with intensities To identify the universal constant term in (1.1), we derive an analog of the key formula in the continuum. This is slightly more delicate than in the discrete because the continuous heat kernel is more singular at t = 0, in particular, there is no s for which the integral ∞ 0 P C (x, x, t)t s−1 dt converges. Therefore, we perform the analytic continuation in two steps. Starting with (3.4), we write The first two terms converge for all s and rest converge when e s > 1. The last two terms evaluate to d| | 2π(s−1) − k s . Let's have a closer look at the third term. To construct its analytic continuation to s = 0, we first remark that it is, in fact, already analytic for e s > 1 2 . Indeed, for x ∈ H, let ρ := m x; by Brownian scaling, Hence, the integral in the third term in (9.1) over x : x H converges absolutely for all s with e s > 1 2 . For x a cone, the same scaling argument leads to the same bound with ρ the distance to the tip, thus, the contribution of those x converges absolutely for e s > 0. Finally, near a corner, we break the integral down as in (6.4) and treat the first three terms as in the cone case and the last two as in the half-plane case. Now, for 1 2 < e s < 1, we can write , x, t) t s−1 dt, cf. the notation in Sect. 6. The last two terms give a contribution that is analytic over e s < 1, hence our task is to analytically continue I x C (x, s). We split it into contributions of neighborhoods of cones, boundary segments, corners and punctures, and leverage the fact that for the scaling on each of x , one has I x C (x, s) = a 2−2s I x  C (ax, s). For x C α C/{z ∼ e iα z}, we thus get x∈C α :|x|<r This is analytic for e s ≤ 1 when divided by (s), which is the only thing we care about; when we eventually evaluate the derivative at 0, we get some value that can be absorbed into the constant D C α . For the contributions of boundary segments and boundary corners, we split their neighborhoods as in Sects. 6.2 and 6.3. The integral over R r yields m x≤r,0≤ e x≤1 For the other contributions, note that those of Y 0,1 cancel out as in Sects. 6.2-6.3, and other contributions can be treated as in the cone case, eventually contributing a constant that can be absorbed into D ϒ α ; same applies to punctures. Dividing (9.1) by (s) and differentiating at s = 0, we get Remark 9.1. The equation (9.1) also allows one to see that in fact, C = −2ζ (0). Since 1 (s) ∼ s near the origin, ζ (0) only receives the contributions from those terms in the right-hand side of (9.1) that have a pole at the origin, that is, the third term and − k s . The third term is split into the contributions from near boundary segments, i.e., where x = H, and the contributions from neighborhoods of the punctures, the cone tips and the corners. The former is evaluated at (9.3) and has no pole at the origin, and the contribution from a cone tip is computed in (9.2), with the residue at s = 0 matching the coefficient found in (6.2). A similar result holds for boundary corners and punctures.
We are in the position to put everything together and prove Theorem 1.1: Proof of Theorem 1.1. We look at the key formula (4.1) term by term. The first three terms converge to the first three terms of (9.4) by Lemma 5.6 and Corollary 6.4. The fourth term can be broken into the contributions of neighborhoods of conical singularities, punctures, and the boundary, whose asymptotics is given by (6.2), (6.9), (6.6), and the constants C p are made explicit in Sect. 7. The fifth term gives the "volume" contribution that is discussed in the end of Sect. 4.
Remark 9.2. The underlying triangulation or quadrangulation structure of was only used in the discretization procedure, but otherwise it plays no role in the proof. While we found no elegant way to state Theorem 1.1 in a more general form that would account for that, we give an example: let T ω δ 1 ,ω δ 2 = C δ /(ω δ 1 Z + ω δ 2 Z) be a sequence of discretized tori, whose periods ω δ 1,2 converge as δ → 0 to ω 1,2 ∈ C \ {0} with ω 2 /ω 1 / ∈ R. Then, the asymptotics (1.1) holds, and since there are no boundary or corners/cones, it takes the form This extends the results in [7,8,13,19]. The above proof applies verbatim; also note that the symmetries of the lattice (other than the double periodicity) are not needed here, as long as the embedding is such that the random walk converges to the Brownian motion.

Proof of the Lemmas
In preparation for the proof of Lemmas 3.1 and A.1, we prove the following bounds on the transition kernel Pˆ (x, y, t) of for the processγ defined in Sect. 3, in particular, establishing the existence of this transition kernel.
Proof. Ifγ is the Brownian motion in C started at the origin and τ = min{t : |γ (t)| ≥ r }, then, by symmetry, P(|γ t | ≥ r |τ < t) ≥ 1 2 . Therefore, Define the sequence of stopping times . Conditionally on F(γ [0,τ 2k−1 ] ) and on the event t ∈ [τ 2k−1 , τ 2k ], the distribution of γ t is that of the (time-shifted) Brownian motion started at a point on ∂ B 3R 4 (x 0 ) and conditioned to stay in B R (x 0 ). IfP(x, y, t) denotes the heat kernel of this conditioned Brownian motion, then we have, for some C > 0, Indeed, by Brownian scaling, we may assume R = 1. Write P D (x, y, t) = i e −λ i t ψ i (x)ψ i (y) and Q(x, t) = D P D (x, y, t)dy, where ψ i and λ i are normalized eigenfunctions and eigenvalues of the Laplacian in the unit disc D with zero boundary conditions. We haveP(x, y, t) = P D (x, y, t)/Q(x, t), and P D (x, y, t) ≤ P C (x, y, t) = 1 2π t e − |x−y| 2 2t . From this, it follows thatP(x, y, t) is bounded for small t andP(x, y, t) t→∞ −→ ψ 1 (y) for large t, and (10.4) follows. We arrive at which proves (10.1). To prove (10.2), let τ R = min{s : |γ s − x 0 | = R} and note that the Brownian motionsγ ,γ are coupled to coincide up to τ, so that and (10.2) follows by (10.1) and (10.3).
Proof of Lemma 3.1. We first check thatγ t is a Leb(ˆ )-symmetric Markov process. For x 0 ,ŷ 0 ∈ˆ and small ε x,y > 0, fix a ball B ε (x 0 ) ⊂ C and an isometryσ of . If its liftγ ends up in B ε y (ŷ 0 ), we denote by σ y (γ ,σ ) the isometry of B ε y (y 0 ) to B ε y (ŷ 0 ) obtained by extendingσ alongγ (otherwise, σ y (γ ,σ ) is undefined.) We denote by Y the set of alľ y 0 obtained in this way, and by y the set of all isometries σ y (γ ,σ ) modulo shifts. Ifγ ends up in B ε x (x 0 ), we define σ x (γ ,σ ), X and x similarly. Note that X , Y are discrete sets in the plane, and x , y are finite with | x | = | y |. We now write, for the coupled Brownian motionsγ ,γ , where ∼ stands for equality of isometries modulo shifts. By the reversibility of the planar Brownian notion, denoting byγ −1 the time-reversal ofγ , Also, P C (x,y, t) = P C (y,x, t), and if we shift every point in Y to some fixed point in this set, then the pointx 0 gets shifted to every point in X . Collecting all these observations together, we end up with which implies thatγ t is Leb(ˆ )-symmetric. The transition kernel of the reflected process can be written asP(x, y, t) +P(x, y, t), implying symmetry, and killing upon hitting a closed subset preserves the class of symmetric processes. For the boundedness, the bound (10.1) and the symmetry imply that Pˆ (x, y, t) ≤ C R −2 if either x or y are at distance at least R from any conical tip ofˆ . Now, fix a small R > 0, and let x ∈ B R (x 0 ), where x 0 is a conical tip. If y / ∈ B R (x 0 ), then Pˆ (x, y, t) ≤ C R −2 by the strong Markov property with respect to exit time from B R (x 0 ). Else, coupling the Brownian motions inˆ and in the infinite cone C α x 0 up to exiting B R (x 0 ), we similarly obtain Pˆ (x, y, t) ≤ C R −2 + P C x 0 (x, y, t).
Coupling to the Brownian motion in C, we have P C x 0 (x, y, t) ≤ y∈Y P C (x,y, t) ≤ const(t, α) where Y is the finite set of endpoints of paths in C starting atx that lift to a path from x to y.
For the smoothness, if B R (x 0 ) is Euclidean, y ∈ B R 2 (x 0 ) and x / ∈ B R (x 0 ), we can write as in the proof of (10.1): and note that all the derivatives ofP(x, y, t) are uniformly bounded over Since R is at our disposal, this proves smoothness for x = y; for x = y, one can use a similar decomposition, on the event that t < τ 2 , we can use that the heat kerned in the disc is smooth.
Proof of Lemma 5.2. It suffices to prove the bound for the walk on C δ (by passing first toˆ δ as in the proof of Lemma 5.1 and then assuming by Markov property that x is at distance at least ε/10 from conical singularities). If t > δ 3 2 , we use Lemma 10.2 and bound P(diam γ δ . The first term converges to 0 super-polynomially in t while the second one bounded from above by C · exp(−(tδ −2 ) ) ≤ C · exp(−t − 3 ) provided that t ≤ 1 and δ Pickĉ > 0 in such a way that for all δ small enough, the random walk on C δ needs at least K := ĉδ −1 steps to reach diameter ε. The probability of this is bounded by P(X ≥ K ), where X is a Poisson random variable with mean M K , and M = t , with c a constant depending on C δ 0 and ε. If M < 1 2 , then, for any α > 0, we have by Stirling bound Since 1 − M + log M < 0 for M < 1, we can pick α > 0 such that the exponential is bounded by 1 for all M < 1 2 , i.e., P(X ≥ K ) ≤ 2(c t 1 3 ) α ĉδ −1 . For δ small enough, the exponent is at least 20, and so we have P(diam γ δ [0,tδ −2 ] > ε) ≤ P(X ≥ K ) ≤ 2t 20 6 provided that t < (c ) −6 .
For x ∈ δ and r > 10N −1 > 0, denote Q(x, r ) := [0; r 2 ] × B(x, r ). The parabolic Harnack inequality (PHI) asserts that there exists a constant C H such that, for any δ, any x, r such that B(x, r ) ∩ ∂ D δ = ∅, and any u positive and satisfying . In our setting, PHI follows from [9, Theorem 1.7]. Delmotte uses normalized Laplcian in which the random walk jumps at rate one; however since he allows for jumps from a vertex to itself, the two setups are equivalent; note that we rescale the graph distance but don't rescale time, hence the additional factor of δ 2 in (10.7).
Of the three conditions of [9, Theorem 1.7], the volume doubling condition DV (C 1 ) and uniform ellipticity conditions (c) are obvious in our setting: they state that |B(x, r )| ≤ C 1 |B(x, 2r )| for any x, r , and min y∼x w yx ≥ cw x for any x, respectively. This third one, the Poincaré inequality P(C 2 ), asserts that Recall (see e.g. [15]) the classical proof of the Poincaré inequality for a ball B: by from which the Poincaré inequality follows by integrating over x, y ∈ B. This proof extends to the discrete settings of balls in C δ , simply by replacing integration with summation, in particular, integration over a segment [x y] with summation over H δneighborhood of [x y] for a large enough fixed H. In the final step, we sum with weights w x w y ; note that w xy in the right-hand side can be ignored since they are uniformly bounded away from 0. For balls in an infinite cone or an infinite wedge, it suffices to map the cone or the wedge by a bi-Lipschitz map to the plane or the half-plane and apply the same proof, increasing H if necessary. Since we only apply PHI to balls in model surfaces, these cases are all we need. A typical example of a function u that PHI is applied to is P δ (x 0 , x, δ −2 t).
We recall the standard argument that PHI implies Hölder regularity of solutions to (10.7), namely, there exist θ > 0 and C Höl > 0 such that any δ > 0 and for Q(x, r ) as above, one has To prove (10.9), note that ifû is u normalized so that inf Q(x,r )û = 0 and sup Q(x,r )û = 1, and sup Q − (x,r )û ≥ 1 2 , then Using PHI as in the proof of Lemma 5.3, we see that this implies P δ (x, y, δ −2 t) < Cδ 2 (1 − η) t for any x, y. Summing this bound over x = y ∈ δ yields the desired result. Case 2. Suppose that has no Dirichlet boundary and ϕ is the trivial line bundle, thus k = 1. We claim that the heat kernel P δ (x, y, δ −2 t) at time t = 1 is uniformly contracting in the total variation distance, i.e., there exists η > 0 such that for any x 1 , x 2 , and any δ small enough, one has and hence, by PHI, inf x inf y∈B(y 0 ,r ) P δ (x, y, δ −2 ) ≥ C H c ψ δ 2 . This gives the desired improvement on the trivial bound of 1 on the LHS of (10.10). Now, iterating (10.10), we see that 1 for all δ small enough, independently of t, as required. Case 3. Suppose that has no Dirichlet boundary and k = 0. Pick any x 0 ∈ δ , and (non-contractible) loops β (1) , . . . β (n) rooted at x 0 such that {ϕ(γ j )} n j=0 do not have a common eigenvector of eigenvalue 1; if such loops did not exist, then the translations of the common eigenvector would form a covariant constant. Pick a a small r such that B(x 0 , r ) is contractible, and, for x, y ∈ B(x 0 , r ), denote, for i = 0, 1, . . . , n, P where β (0) is a contractible loop, and we identify the points of B(x 0 , r ) in order to compute the homotopy type of a non-closed path. There exists a constant c > 0 such that cδ 2 ≤ P δ (i) (x, y, δ −2 ) ≤ c −1 δ 2 for each i and each x, y ∈ B(x 0 , r ) and all δ small enough; the upper bound follows from Lemma 5.3 and the lower one is done exactly as in Case 2. We can write where in the last term, we used that ||ϕ(γ δ [0,δ −2 t] )|| ≤ 1 as ϕ is unitary. We claim that Indeed, the fraction is strictly smaller than 1 unless all ϕ(β (i) )v are non-negative multiples of each other, but since ϕ(β 0 ) = Id, this would mean that v is a common eigenvector of eigenvalue 1. Applying (10.12) to (10.11), we get, for any x ∈ B(x 0 , r ) and all δ small enough, ]}, and we have used that P(E) ≥ y∈B(x 0 ,r ) P δ (0) (x, y, δ −2 ) is uniformly bounded from below. Since the bounds obtained, in fact, did not depend on x 0 , the proof is now completed as in Case 1. Case 4. Suppose that has no Dirichlet boundary. Let ϕ 0 be the trivial sub-bundle of the maximal dimension of ϕ, which is a direct sum of trivial line bundles. Since ϕ is unitary, we have ϕ = ϕ 0 ⊕ ϕ ⊥ 0 for ϕ ⊥ 0 the (point-wise) orthogonal complement to ϕ 0 ; moreover, ϕ ⊥ 0 has no trivial line sub-bundles. Applying Case 2 to ϕ 0 and Case 3 to ϕ 1 concludes the proof.
Using this to bound the oscillation in the Hölder bound (10.9), we get, for t ≥ δ 2κ , provided that κ is small enough (in which case also ε δ ν so that the Hölder bound applies). A similar estimate holds for the continuous heat kernel. Combining this with (10.15) gives the claim. fails after the moment the former hits the tip, in particular, the distributions of the first time they leave the tip will be different. However, forcingγ δ δ −2 t to leave the tip simultaneously with γ δ δ −2 t yields a coupling of γ δ δ −2 t toγ δ δ −2 t+τ such that still γ δ δ −2 t ∈ ∪ k σ k (B(y, ε)) if and only ifγ δ δ −2 t+τ ∈ ∪ kσ k (B(y, ε)) . Here τ a random variable given by a sum of N i. i. d. contributions, where N is the number of visits to the tip. The expectation of N , and hence that of τ , is O(log(δ −2 t)), and therefore it will introduce a negligible error into the above computations.
Proof of Lemma 6.2. As explained in in the proof of Lemma 6.1, it suffices to consider the case of a cone. Moreover, since, in the notation of that proof, γ δ δ −2 t ∈ ∪ k {σ k (y)} if and only ifγ δ δ −2 t ∈ ∪ k {σ k (y)}, it is in fact sufficient to consider the case of a plane, where it is immediate from the local Central limit theorem.

])
D ,ϕ t f ; g = E (1,t) ( f, g) + E (2,t) x,y,t f (y), g(x) − R ,ϕ x,y,t g(y) P (x, y, t) dxdy, and Q t (x) = P (x, y, t)dy is the probability that γ t started at x did not stop by time t.
For a compact set K ⊂ and K ε its small neighborhood, and h ∈ D(E ,ϕ ) smooth, we have x,y,t h(x), h(y) − R ,ϕ x,y,t h(x) P (x, y, t)dxdy as is readily seen by identifying h with a functionh using a local trivialization near each x and using Taylor approximation ofh and (10.1-10.2). Taking a supremum over all such K , we conclude that h ∈ H 1,ϕ ( ). Finally, we refine D(E ,ϕ ) ⊂ H 1,ϕ ( ) to D(E ,ϕ ) ⊂ H 1,ϕ ∂ ( ). We do the trivial line bundle case; the general one only differs by heavier notation. Let 0 ≤ φ ≤ 1 be a smooth function compactly supported in a Euclidean neighborhood of x ∈ ∂ D \ S and identically equal to 1 in a smaller neighborhood. If f ∈ D(E ), then we have E (2,t) and both terms remain bounded as t → 0. Hence, also φ f ∈ D(E ). Identifying supp φ f with a subdomain of the upper half-plane H equipped with Dirichlet boundary conditions, due to (10.2), we have 1