# Harmonic functions on Tutte's embeddings and linearized Monge-Ampère equation

Mikhail Basok, Dmitry Chelkak, Benoît Laslier, Marianna Russkikh

<figure><img src="image_1.png" /></figure>

## To cite this version:

Mikhail Basok, Dmitry Chelkak, Benoît Laslier, Marianna Russkikh. Harmonic functions on Tutte's embeddings and linearized Monge-Ampère equation. 2025. hal-05416195

HAL Id: hal-05416195

https://hal.science/hal-05416195v1

Preprint submitted on 15 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
HARMONIC FUNCTIONS ON TUTTE'S EMBEDDINGS
AND LINEARIZED MONGE-AMPÈRE EQUATION

MIKHAIL BASOK<sup>A</sup>, DMITRY CHELKAK<sup>B</sup>, BENOÎT LASLIER<sup>C</sup>, AND MARIANNA RUSSKIKH<sup>D</sup>

ABSTRACT. We prove convergence of solutions of Dirichlet problems and Green's functions on Tutte's harmonic embeddings to those of the linearized Monge-Ampère equation $\mathcal{L}_{\varphi}h = 0$. The potential $\varphi$ appears as the limit of piecewise linear potentials associated with the embeddings and the only assumption that we use is the uniform convexity of $\varphi$. Even if $\varphi$ is quadratic, this setup significantly generalizes known results for discrete harmonic functions on orthodiagonal tilings. Motivated by potential applications to the analysis of 2d lattice models on irregular graphs, we also study the situation in which the limits are harmonic in a different complex structure.

# 1. Introduction and main results

**1.1. Motivation.** Convergence of discrete harmonic functions in the small mesh size limit is a classical topic dating back to the first works on numerical methods [Lyu26; CFL28]. In two dimensions, this question is intrinsically related to the notion of discrete holomorphic functions and their convergence, with important contributions made in the middle of the last century by Ferrand [Lel55, Chapter V] who greatly developed this theory on the square grid and Duffin [Duf68] who, in particular, suggested a generalization to the so-called rhombic tilings of the complex plane. From the combinatorial perspective, two related classical setups for which the notions of discrete holomorphic and harmonic functions make perfect sense are Brooks-Stone-Smith-Tutte's square (or, more generally, rectangular) tilings [Bro+40] and Tutte's barycentric, or harmonic, embeddings [Tut63]. One of the sources of the modern interest to discrete complex analysis is the 2d statistical physics; e.g., see [Smi10]. In particular, the framework of rhombic lattices was important for the study of scaling limits of free fermionic models with Baxter's Z-invariant weights; e.g., see [Mer01; Ken02; CS12; BTR17].

Recent advances of the probabilistic approach to the Liouville Quantum Gravity (LQG) in 2d (e.g., see [GHS23] and references therein) motivate a development of discrete complex analysis techniques for possibly very irregular planar grids that appear as embeddings of given (random) planar graphs into the plane. In particular, applications of Tutte's embeddings and square tilings in the LQG context recently appeared in [GMS21; BGS25]. Going in this direction, the last three authors of this paper suggested a general framework of the so-called t-embeddings of graphs equipped with a bipartite dimer model [CLR23], also known as Coulomb gauges [Ken+21], which simultaneously generalizes Tutte's embeddings and the so-called s-embeddings of planar graphs carrying the Ising model. Shortly afterwards, it was understood that t-embeddings should be viewed as piecewise linear surfaces in the Minkowski space $\mathbb{R}^{2,2}$, called t-surfaces below, which provide a correct description of the complex structure arising in the small mesh size limit [CR24; Che24; CLR21; BNR25].

<footnote><sup>A</sup> UNIVERSITY OF HELSINKI, DEPARTMENT OF MATHEMATICS AND STATISTICS, 00500 HELSINKI, FINLAND.</footnote>

<footnote><sup>B</sup> UNIVERSITY OF MICHIGAN, DEPARTMENT OF MATHEMATICS, ANN ARBOR, MI 48109-1043, USA.</footnote>

<footnote><sup>C</sup> UNIVERSITÉ PARIS CITÉ AND SORBONNE UNIVERSITÉ, CNRS, LABORATOIRE DE PROBABILITÉS, STATISTIQUE ET MODÉLISATION, F-75013 PARIS, FRANCE AND DÉPARTEMENT DE MATHÉMATIQUES ET APPLICATIONS, ÉCOLE NORMALE SUPÉRIEURE, UNIVERSITÉ PSL, CNRS, 75005 PARIS, FRANCE.</footnote>

<footnote><sup>D</sup> UNIVERSITY OF NOTRE DAME, DEPARTMENT OF MATHEMATICS, 255 HURLEY BLDG, NOTRE DAME, IN 46556, UNITED STATES OF AMERICA.</footnote>

<footnote>E-mail: mikhail.basok@helsinki.fi, dchelkak@umich.edu, laslier@lpsm.paris, mrusskik@nd.edu.</footnote>
In this paper we apply the formalism of [CLR23] to Tutte's embeddings and prove the convergence of discrete harmonic functions to solutions of the linearized Monge-Ampère equation under no regularity assumptions on the graphs in question. To see that linearized Monge-Ampère equations provide the right language for this problem, the reader may observe that in the degenerate case when the potential $\varphi$ is piecewise linear, the elliptic operator $L_\varphi$ given by (1.9) generates the Brownian motion on the cable system (metric graph) of the corresponding Tutte's embedding.

To the best of our knowledge, recent literature on the subject was mostly focused on the very special case – so-called orthodiagonal tilings – which appears when Tutte's embeddings of a graph and of its dual form non-self-intersecting quads. In this setup, the convergence to usual harmonic functions (that is, solutions of the equation $-\Delta h = 0$, which corresponds to the quadratic potential $\varphi(w) = \frac{1}{2}|w|^2$ in (1.9)) under no regularity assumptions on the graphs was obtained in [GJN20]; see also references therein and [BP25; BG25]. However, this passage from Tutte's embeddings to orthodiagonal ones is quite restrictive. In particular, it generally does not work if one starts with a given harmonic embedding $\Gamma_\delta$ and then superimposes $\Gamma_\delta$ with the harmonic embedding $\Gamma_\delta^*$ of the dual graph. (Even if edges of $\Gamma_\delta^*$ 'stay close' to the corresponding edges of $\Gamma_\delta$, the quads formed by them typically self-intersect.) In this paper we aim to provide a much more general framework that can be applied to the study of 2d lattice models on irregular weighted graphs given in advance and harmonically embedded into the plane, which avoids this constraint and also does not require that the potential $\varphi$ is quadratic.

## 1.2. Tutte's embeddings $\Gamma_\delta$, potentials $\Phi_\delta$, and random walks on $\Gamma_\delta$.
Throughout this paper we work with a sequence $(\Gamma_\delta)_{\delta \to 0}$ of Tutte's embeddings of finite planar weighted graphs, each of which covers a fixed simply connected domain $U \subset \mathbb{C}$. Let us emphasize that the parameter $\delta$ is not the 'mesh size' of $\Gamma_\delta$ and can be much bigger than a 'typical size' of its faces; see equivalent properties (CONV), (LIP), (RW), and Remark 1.3 below. Each $\Gamma_\delta$ is a planar graph drawn in the complex plane such that its edges are non-intersecting straight segments. We denote by $v$ the vertices of this graph and by $\mathcal{H}_\delta(v) \in \mathbb{C}$ the coordinates of these vertices in the complex plane. In what follows we sometimes identify $v$ and $\mathcal{H}_\delta(v)$ if no confusion arises. Each edge $(vv')$ has a pre-assigned weight, or conductance, $c_{vv'} > 0$. The embedding $\mathcal{H}$ is called Tutte's, barycentric, or harmonic if the identity $\sum_{v' \sim v} c_{vv'}(\mathcal{H}(v') - \mathcal{H}(v)) = 0$ holds for each inner vertex of the graph.

Given a harmonic embedding $\mathcal{H}_\delta$, one can define a function $\mathcal{H}_\delta^*$ on the dual graph by setting

$$ (1.1) \qquad \mathcal{H}_\delta^*(v_2^*) - \mathcal{H}_\delta^*(v_1^*) = ic_{v_1v_2}(\mathcal{H}_\delta(v_2) - \mathcal{H}_\delta(v_1)) $$

if $(v_1^*v_2^*)$ is the dual edge to $(v_1v_2)$ and $v_1^*$ is to the right from $(v_1v_2)$. The harmonicity of $\mathcal{H}_\delta$ implies that the sum of these increments around each inner vertex of $\Gamma$ vanishes and hence $\mathcal{H}_\delta^*$ – which is called the dual Tutte's embedding – is well defined up to a global additive constant. Let us emphasize that this dual embedding can have overlaps if $U$ is not convex; see Section 2 for details.

Let us now define piecewise constant maps $\Psi_\delta: U \to \mathbb{C}$ by setting

$$ (1.2) \qquad \Psi_\delta(w) := \mathcal{H}_\delta^*(v^*) \quad \text{if } w \in \mathbb{C} \text{ belongs to the face of } \Gamma_\delta \text{ corresponding to the dual vertex } v^*. $$

(If $w \in \mathbb{C}$ is a vertex of $\Gamma_\delta$ or lies on an edge, we choose one of the incident faces arbitrarily.) A very important property of $\Psi_\delta$ is the existence of a piecewise linear function $\Phi_\delta: U \to \mathbb{R}$ such that

$$ \Psi_\delta = 2\partial_{\bar{w}}\Phi_\delta = \partial_x\Phi_\delta + i\partial_y\Phi_\delta \quad \text{for all } w = x + iy \text{ lying on faces of } \Gamma_\delta. $$

One can construct $\Phi_\delta$ explicitly by defining its values on vertices of $\Gamma_\delta$ via the identity

$$ (1.3) \qquad \Phi_\delta(v_2) - \Phi_\delta(v_1) = \operatorname{Re} \left[ \overline{\mathcal{H}_\delta^*(v_1^*)(\mathcal{H}_\delta(v_2) - \mathcal{H}_\delta(v_1))} \right] = \operatorname{Re} \left[ \overline{\mathcal{H}_\delta^*(v_2^*)(\mathcal{H}_\delta(v_2) - \mathcal{H}_\delta(v_1))} \right] $$

and noting that $\Phi_\delta$ restricted to vertices of a face $v^*$ of $\Gamma_\delta$ coincides with a linear function with gradient $\mathcal{H}_\delta^*(v^*)$. We call $\Phi_\delta$ the potential corresponding to $\Gamma_\delta$. Since $\mathcal{H}_\delta^*$ and hence $\Psi_\delta$ are defined up to a global additive constant, the potential $\Phi_\delta$ is defined by $\Gamma_\delta$ up to a linear function.
It follows from (1.1) and $c_{v_1v_2}>0$ that the piecewise linear function $\Phi_\delta$ is convex near each edge of $\Gamma_\delta$. Therefore, $\Phi_\delta$ is a convex function on each convex subset of $U$. The only property of potentials $\Phi_\delta$ that we use in our paper is their uniform convexity above scales $\delta \to 0$; see also Remark 1.3 below.

**Property (CONV).** We say that the potentials $(\Phi_\delta)_{\delta \to 0}$ have property (CONV) on the domain $U \subset \mathbb{C}$ if there exist constants $\lambda > 0$ and $C > 0$ such that the inequalities

$$ (1.4) \qquad \lambda |w_2 - w_1|^2 \le \Phi_\delta(w_2) - 2\Phi_\delta(\frac{1}{2}(w_1+w_2)) + \Phi_\delta(w_1) \le \lambda^{-1} |w_2 - w_1|^2 $$

hold for all straight segments $[w_1; w_2] \subset U$ such that $|w_2 - w_1| \ge C\delta$.

**Remark 1.1.** The well-studied ‘nicest’ setup for discrete complex analysis beyond the square grid is that of rhombic lattices or isoradial graphs introduced by Duffin in [Duf68]; see also [CS11]. In this case one can easily deduce from (1.1) that the restriction of $\Phi_\delta$ onto the vertices of $\Gamma_\delta$ exactly equals the quadratic potential $\Phi_\delta(v) = \frac{1}{2}|\mathcal{H}_\delta(v)|^2$. In a more general setup of orthodiagonal tilings [GJN20; BP25; BG25], one has $\Psi_\delta(w) = w + O(\delta)$ and hence $\Phi_\delta(w) = \frac{1}{2}|w|^2 + O(\delta|w|)$, where $\delta \to 0$ denotes the size of the largest tile. As for the converse, let us emphasize that the assumption $\Psi_\delta(w) = w + O(\delta)$ does not put the analysis back to the setup of [GJN20; BP25; BG25] even if $\delta$ is comparable with the lengths of all the edges of $\Gamma_\delta$ since the quads $\mathcal{H}_\delta(v_1)\mathcal{H}_\delta^*(v_1)\mathcal{H}_\delta(v_2)\mathcal{H}_\delta^*(v_2)$ may self-intersect. See also Remark 1.3 and Remark 1.4 below.

In terms of mappings $\Psi_\delta: U \to \mathbb{C}$, property (CONV) can be formulated as follows.

**Property (LIP).** We say that the maps $(\Psi_\delta)_{\delta \to 0}$ defined by (1.2) have property (LIP) on the domain $U \subset \mathbb{C}$ if there exist constants $\varkappa > 0$ and $C > 0$ such that

$$ (1.5) \qquad \mathrm{Re} \frac{\Psi_\delta(w_2) - \Psi_\delta(w_1)}{w_2 - w_1} \ge \varkappa \quad \text{and} \quad \left| \frac{\Psi_\delta(w_2) - \Psi_\delta(w_1)}{w_2 - w_1} \right| \le \varkappa^{-1} $$

for all straight segments $[w_1; w_2] \subset U$ such that $|w_2 - w_1| \ge C\delta$.

**Remark 1.2.** It is not hard to see that property (CONV) for potentials $\Phi_\delta$ on $U$ implies property (LIP) for their gradients $\Psi_\delta$ on each subdomain $U' \subset U$ and vice versa; we refer the reader to Lemma 5.1 for the proof of this elementary fact. (It is sufficient to assume that $U'$ lies in an $O(\delta)$-interior of $U$.) More importantly, property (LIP) is also equivalent to the $\kappa$-Lipschitzness, $\kappa = \kappa(\varkappa) < 1$, of the so-called origami map in the terminology of [CLR23], which makes the techniques developed in [CLR23] applicable to the setup of this paper; see Section 2.2 below.

Given a harmonic embedding $\Gamma_\delta$, we define the continuous time random walk $X_t$ on its vertices so that the jump rates from a vertex $v$ to its neighbors $v'$ are proportional to $c_{vv'}$, with the normalization chosen so that $|X_t|^2 - t$ is a martingale. (Note also that $X_t$ is a $\mathbb{C}$-valued martingale since the embedding is harmonic.) The normalization factors

$$ (1.6) \qquad \mu_\delta(v) := \sum_{v' \sim v} c_{vv'} |\mathcal{H}_\delta(v') - \mathcal{H}_\delta(v)|^2 $$

give us a standard invariant measure for $X_t$ if one excludes the effect of the boundary of $\Gamma_\delta$ from the consideration. Let us formulate two important properties of these random walks.

**Property (RW).** We say that a sequence $(\Gamma_\delta)_{\delta \to 0}$ of Tutte's embeddings has property (RW) on a domain $U' \subset U$ if there exists constants $c, C > 0$ such that the following holds:

(a) The random walk on $\Gamma_\delta$ is uniformly elliptic starting from the scale $\delta$: for each $\alpha \in \mathbb{C}$ such that $|\alpha| = 1$ and each vertex $v \in \Gamma_\delta \cap U'$ one has $\mathrm{Var}(\mathrm{Re}[\bar{\alpha}X_\tau^v(B(v, C\delta))]) \ge c\delta^2$, where $X_\tau^v$ denotes the random walk started at $v$ and $\tau(B(v, C\delta))$ is the first exist time of $X_\tau^v$ from the ball $B(v, C\delta)$.

(b) The invariant measure $\mu_\delta$ is comparable to the Lebesgue measure starting from the scale $\delta$: for each $x \in U'$ one has $c\delta^2 \le \sum_{v \in B(x, C\delta)} \mu_\delta(v) \le c^{-1}\delta^2$.
Our first result shows the equivalence of the property (RW) of random walks on Tutte's embeddings $\Gamma_{\delta}$ to (equivalent) 'geometric' properties (CONV)/(LIP) of $\Gamma_{\delta}$ discussed above.

**Theorem 1.** (i) Properties (CONV) and/or (LIP) of Tutte's embeddings $\Gamma_{\delta}$ on a domain $U \subset \mathbb{C}$ imply property (RW) of the random walks on $\Gamma_{\delta}$ on each subdomain $U' \in U$.

(ii) Vice versa, property (RW) on $U'$ implies properties (CONV), (LIP) on each subdomain $U'' \in U'$. Moreover, it is enough to assume that $U'$, resp. $U''$, lies in an $O(\delta)$-interior of $U$, resp. $U'$.

The fact that (LIP) implies (RW) follows from the framework developed in [CLR23]; in particular, [CLR23, Proposition 6.4] gives the uniform ellipticity of the random walk $X_t$ above scale $\delta$. We recap this material in Section 2.3 and also sketch a self-contained proof of the ellipticity property in Remark 4.3. We prove the converse statement in Section 5; see Theorem 5.6.

**1.3. Main convergence results.** Let $\Omega \in U$ be a smooth simply connected domain. Both smoothness and simply connectedness assumptions on $\Omega$ are not essential for our results and made only for simplicity; see also Remark 1.5 below. Put $\Omega_{\delta}$ to be the subgraph of $\Gamma_{\delta}$ spanned by those vertices of $\Gamma_{\delta}$ that lie inside $\Omega$. (If this intersection has more than one connected components, we keep only the component that covers the bulk of $\Omega$.) Let $\bar{\Omega}_{\delta}$ be the union of $\Omega_{\delta}$ with neighbors of these vertices and $\partial\Omega_{\delta} := \bar{\Omega}_{\delta} \setminus \Omega_{\delta}$. Given a function $H_{\delta}: \bar{\Omega}_{\delta} \to \mathbb{R}$ on vertices of $\bar{\Omega}_{\delta}$ and $v \in \Omega_{\delta}$ we define

$$ (1.7) \qquad [\mathcal{L}_{\delta}H_{\delta}](v) = (\mu_{\delta}(v))^{-1} \sum_{v' \sim v} c_{vv'}(H_{\delta}(v) - H_{\delta}(v')). $$

A function $H_{\delta}$ is called discrete harmonic in $\Omega_{\delta}$ if $[\mathcal{L}_{\delta}H_{\delta}](v) = 0$ for all $v \in \Omega_{\delta}$.

Assume now that potentials $(\Phi_{\delta})_{\delta \to 0}$ converge, uniformly on compacts, to a function $\varphi: U \to \mathbb{R}$. It follows from (1.4) and (1.5) that the limit $\varphi$ is uniformly convex and $C^{1,1}$-smooth; the latter means that $\varphi$ is continuously differentiable and has uniformly Lipschitz gradient $\psi = 2\varphi_{\bar{w}} = \varphi_x + i\varphi_y$. Let

$$ (1.8) \qquad A_{\varphi} = \begin{pmatrix} \varphi_{yy} & -\varphi_{xy} \\ -\varphi_{xy} & \varphi_{xx} \end{pmatrix} = \begin{pmatrix} \mathrm{Re}(\psi_w - \psi_{\bar{w}}) & \mathrm{Im}(\psi_w - \psi_{\bar{w}}) \\ -\mathrm{Im}(\psi_w + \psi_{\bar{w}}) & \mathrm{Re}(\psi_w + \psi_{\bar{w}}) \end{pmatrix}. $$

Note that $A_{\varphi}$ is defined almost everywhere due to Rademacher's theorem and is uniformly elliptic, i.e., $\lambda \le A_{\varphi} \le \lambda^{-1}$ due to (1.4). Consider the differential operator

$$ (1.9) \qquad \mathcal{L}_{\varphi}h = -\varphi_{yy}h_{xx} + 2\varphi_{xy}h_{xy} - \varphi_{xx}h_{yy} = -\mathrm{div}(A_{\varphi}\nabla h). $$

The equation $\mathcal{L}_{\varphi}h = 0$ is known under the name linearized Monge-Ampère equation [CG97]. Even though the matrix $A_{\varphi}$ is only measurable, the fact that the operator $\mathcal{L}_{\varphi}$ can also be written in the divergence form allows one to speak about weak solutions, that is, functions $h \in W_{\mathrm{loc}}^{1,2}(\Omega)$ such that $\int_{\Omega}(A_{\varphi}\nabla h \cdot \nabla \phi)(w)d^2w = 0$ for all smooth test functions $\phi \in C_{0}^{\infty}(\Omega)$.

Another important feature coming from the divergence expression for $\mathcal{L}_{\varphi}$ is the notion of a harmonic conjugate. Define $\ast = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}$ and note that for every weak solution $h$ there is a function $h^* \in W_{\mathrm{loc}}^{1,2}(\Omega)$ satisfying $\nabla h^* = \ast A_{\varphi}\nabla h$. Such $h^*$ is called an $A_{\varphi}$-harmonic conjugate of $h$; see [AIM09, Chapter 16] and Section 3 for more details. Note that $h^*$ is a weak solution of the 'dual' equation

$$ (1.10) \qquad \mathcal{L}_{\varphi}^{\times}h^* = \mathrm{div}(\ast A_{\varphi}^{-1} \ast \nabla h^*) = 0. $$

The standard theory of elliptic equations (see [Ken+00] and references therein) guarantees that weak solutions $h, h^*$ are continuous in $\Omega$. Moreover, the Dirichlet problem $\mathcal{L}_{\varphi}h = 0$ in $\Omega$, $h|_{\partial\Omega} = g$ has a unique solution $h \in W_{\mathrm{loc}}^{1,2}(\Omega) \cap C(\bar{\Omega})$ for each continuous function $g: \partial\Omega \to \mathbb{R}$ [LSW63; Ken+00].

Note also that, as in (1.1), with each discrete harmonic function $H_{\delta}$ on a Tutte embedding $\Gamma_{\delta}$ one can associate its harmonic conjugate $H_{\delta}^*$, which is a discrete harmonic function on the dual embedding. Our first convergence theorem deals with solutions of discrete Dirichlet problems.

**Theorem 2.** In the setup described above, let $g$ be a continuous function defined in a neighborhood of $\partial\Omega$ and $h$ be the solution of the Dirichlet problem $\mathcal{L}_{\varphi}h = 0$, $h|_{\partial\Omega} = g$. Let $H_{\delta}: \bar{\Omega}_{\delta} \to \mathbb{R}$ be the solution of the same Dirichlet problem for $\mathcal{L}_{\delta}$, i.e., the unique discrete harmonic function in $\Omega_{\delta}$ such
that $h_\delta(v) = g(v)$ for all $v \in \partial\Omega_\delta$. Then, the functions $H_\delta$ converge to $h$ uniformly in $\bar{\Omega}$ as $\delta \to 0$. Moreover, discrete harmonic conjugates $H_\delta^*$ can be chosen such that they converge to an $A_\varphi$-harmonic conjugate $h^*$ of $h$ uniformly on compacts of $\Omega$.

**Remark 1.3.** Instead of assuming that $(\Gamma_\delta)_{\delta \to 0}$ have property (CONV), one can start with an arbitrary sequence of Tutte's embeddings $\Gamma^{(n)}$ such that the potentials $\Phi^{(n)}$ converge to a uniformly convex function $\varphi: U \to \mathbb{R}$. This implies the existence of parameters $\delta^{(n)} \to 0$ such that $\Phi^{(n)}$ have property (CONV) with twice smaller $\lambda$ starting from the scale $\delta^{(n)}$, and one can rename $\Phi^{(n)}$ by $\Phi_{\delta^{(n)}}$.

We also have a similar convergence result for Green's functions. Namely, let $G_\Omega^\delta(\cdot, v_0)$ be the discrete Green function in $\Omega_\delta$, i.e., the unique function that is discrete harmonic in $\Omega_\delta \setminus \{v_0\}$, has zero boundary conditions at $\partial\Omega_\delta$, and such that $\mu_\delta(v_0)[\mathcal{L}_\delta G_\Omega^\delta(\cdot, v_0)](v_0) = 1$. Equivalently, $\mu_\delta(v_0)G_\Omega^\delta(v, v_0)$ is the expected time spent at $v_0 \in \Omega_\delta$ by the random walk $X_t^v$ on $\Gamma_\delta$ started at $v$ before it exists $\Omega_\delta$. Let $G_\Omega(\cdot, w_0)$ be the Green's function of the elliptic operator $\mathcal{L}_\varphi$, that is, the unique positive function that belongs to the space $W_{\text{loc}}^{1,2}(\Omega \setminus \{w_0\}) \cap C(\bar{\Omega} \setminus \{w_0\})$, vanishes at $\partial\Omega$, and satisfies the identity

$$ \int_{\Omega} A_{\varphi} \nabla_{w} G(w, w_0) \cdot \nabla \phi(w) \, d^2 w = \phi(w_0) $$

for all smooth compactly supported test functions $\phi \in C_0^\infty(\Omega)$.

**Theorem 3.** In the setup described above, discrete Green's functions $G_\Omega^\delta$ converge, as $\delta \to 0$, to the Green's function $G_\Omega$, uniformly on each compact $K \subset (\bar{\Omega} \times \bar{\Omega}) \setminus \text{diag}$.

We prove Theorem 2 and Theorem 3 in Section 4 first under an additional smoothness assumption $\varphi \in C^3$ in Section 4.3 and then in the general case $\varphi \in C^{1,1}$. A simplification comes from the fact that for $\varphi \in C^3$ each 'very weak' solution of the equation $\mathcal{L}_\varphi h = 0$, i.e., a function $h$ such that $\int_{\Omega} h(w)[\mathcal{L}_\varphi \phi](w) d^2 w = 0$ for all $\phi \in C_0^\infty(\Omega)$, is automatically a weak solution. In the general case $\varphi \in C^{1,1}$, we rely upon the a priori regularity theory developed in [CLR23] and an idea somewhat similar to the one used in [Mah25] in the Ising model context in order to show that each discrete harmonic function on $\Gamma_\delta$ can be locally approximated by discrete harmonic functions on the same graph with bounded gradients. However, we do not know whether or not all discrete harmonic functions on $\Gamma_\delta$ have uniformly bounded gradients. In the positive direction, we prove that (a) the convergence in Theorem 2 and 3 also holds for $\nabla H_\delta$ and $\nabla G_\Omega^\delta(\cdot, v_0)$ if the graphs $(\Gamma_\delta)_{\delta \to 0}$ satisfy a very mild additional regularity assumption (EXP-FAT), and (b) the same convergence always holds for the 'gradients defined on a mesoscopic scale'. See Section 4.4 and Section 4.5 for details.

**Remark 1.4.** In the special setup of orthodiagonal tilings, Theorem 2 has been obtained in [GJN20]; see also references therein and [BG25, Theorem B] for an alternative proof. Moreover, in [BG25, Theorem A] the authors show that in this case [BG25, Theorem B] implies the convergence of trajectories of the random walks on $\Gamma_\delta$ to trajectories of the Brownian motion, up to a possible time-reparametrization. We believe that Theorem 2 supplemented by Theorem 3 is sufficient to prove the convergence of random walks $X_t$ on $\Gamma_\delta$ to the driftless diffusion defined by the elliptic operator $\mathcal{L}_\varphi$ together with the time-parametrizations. However, we do not explore this question in our paper.

**Remark 1.5.** In applications of convergence results to 2d lattice models, it is often important to drop the smoothness assumption on $\partial\Omega$ and to use the Carathéodory topology for the convergence $\Omega_\delta \to \Omega$ as $\delta \to 0$. The reason for such a generalization stems from the fact that it allows to use a simple compactness argument in order to deduce that the convergence, e.g., of functions $H_\delta$ to $h$ in Theorem 2 is uniform with respect to the domain $\Omega$; e.g., see a discussion in [CS11, Section 3.2]. In some setups it is also important to consider harmonic measures of boundary arcs or fractal interfaces instead of solutions of Dirichlet problems with continuous boundary data $g$. Such generalizations of Theorem 2 and Theorem 3 can be deduced from the uniform ellipticity property (RW) of random walks on $\Gamma_\delta$, which implies the so-called weak-Beurling estimate similarly to [CS11].
**1.4. Convergence to harmonic functions in a non-trivial metric and t-surfaces in $\mathbb{R}^{2,2}$**

In general, there exists no coordinate change $w = w(\zeta)$ under which the operator $\mathcal{L}_{\varphi}$ becomes proportional to the standard Laplacian in the new complex coordinate $\zeta$. Still, the situations in which such a coordinate change exists are of special interest from the 2d statistical physics perspective: in this case, one recovers the conformal invariance (in the complex structure given by $\zeta$) of the limits of discrete harmonic functions on $\Gamma_{\delta}$. Our next result, Theorem 4 given below, provides a few equivalent conditions that characterize the existence of such $\zeta$.

To formulate this theorem we need to introduce a geometric object originating from the main tool that we use in our paper: *t-surfaces* $\Theta^{\delta} \subset \mathbb{R}^{2,2}$ corresponding to Tutte's embeddings $\Gamma_{\delta}$; see Section 2 for details. As $\delta \to 0$, these discrete surfaces $\Theta^{\delta}$ converge to the (locally) space-like surface

$$ (1.11) \qquad \Theta = \left\{ \left( \frac{1}{2}(w + \psi(w)); \frac{1}{2}(\bar{\psi}(\bar{w}) - \bar{w}) \right) \in \mathbb{C}^{1,1} \cong \mathbb{R}^{2,2} \mid w \in U \right\}, $$

where $\psi = 2\varphi_{\bar{w}} = \varphi_x + i\varphi_y$ and the quadratic form in the Minkowski space $\mathbb{C}^{1,1} \cong \mathbb{R}^{2,2}$ is given by

$$ (1.12) \qquad \|(x_1 + ix_2; x_3 + ix_4)\|_{1,1}^2 = \|(x_1, x_2; x_3, x_4)\|_{2,2}^2 = x_1^2 + x_2^2 - x_3^2 - x_4^2. $$

**Theorem 4.** (i) The intrinsic metric on the surface $\Theta$ induced by the metric (1.12) in $\mathbb{C}^{1,1} \cong \mathbb{R}^{2,2}$ coincides with the Riemannian metric on the domain $U \subset \mathbb{C} \cong \mathbb{R}^2$ defined by the Hessian $D^2\varphi$.

(ii) The following conditions are equivalent:

(a) The convex potential $\varphi$ solves the Monge-Ampère equation $\det D^2\varphi = \text{const}$ or, equivalently, its gradient $\psi = 2\varphi_{\bar{w}} = \varphi_x + i\varphi_y$ preserves the area up to a global multiplicative constant.

(b) $\Theta$ is a maximal space-like surface in $\mathbb{R}^{2,2}$ (i.e., its mean curvature vanishes at every point).

(c) There exists a coordinate change $w = w(\zeta)$ such that we have $\mathcal{L}h = 0$ if and only if $h$ is harmonic in $\zeta$. If this holds, then $\zeta$ is a conformal coordinate on $\Theta$.

(d) The two differential equations $\mathcal{L}_{\varphi}h = 0$ and $\mathcal{L}_{\varphi}^{\times}h = 0$ are equivalent.

(e) The equation $\mathcal{L}_{\varphi}\psi = 0$ holds (in the weak sense).

**Remark 1.6.** Assertion (d) means that random walks on the weighted planar graphs in question and on their duals ‘behave in a similar way’ on large scales. (In general, the operators $\mathcal{L}_{\varphi}$ and $\mathcal{L}_{\varphi}^{\times}$ have the same diffusion parts but different drifts; see Section 6 for details.) One can try to verify this condition in various setups of interest as a precursor of the conformal invariance of 2d lattice models such as LERW/UST before choosing particular embeddings of these graphs into the complex plane.

**1.5. Organization of the paper.** We recall the combinatorial construction of t-embeddings and t-surfaces associated with Tutte's embedding and the modern discrete complex analysis framework developed in [CLR23] in Section 2. Section 3 is devoted to basic properties of elliptic equations in 2d. We prove main convergence theorems in Section 4 basing upon classical Dirichlet energy and Caccioppoli estimates. Section 5 contains the proof of the equivalence of properties (CONV)/(LIP) and property (RW). Theorem 4 is proved in Section 6.

**Acknowledgments.** The work of Mikhail Basok was supported by Research Council of Finland grants 333932 and 363549. The work of Benoît Laslier was partially supported by the DIMERS project ANR-18-CE40-0033. Marianna Russkikh was partially supported by a Simons Foundation Travel Support for Mathematicians MPS-TSM00007877. Part of this research was performed while the authors were visiting the Institute for Pure and Applied Mathematics (IPAM), which is supported by the National Science Foundation (Grant No. DMS-1925919).

## 2. T-surface of a harmonic embedding

Given a weighted planar graph $(\Gamma_{\delta}, c)$ embedded into the complex plane harmonically, there are two discrete surfaces that can be associated with it. The first is the graph of the convex potential $\Phi_{\delta}$ associated with the embedding and discussed in the introduction; this surface belongs to the Euclidean
<figure><img src="image_3.png" /><figcaption>FIGURE 1. From left to right: an example of a harmonic embedding $\mathcal{H}$ of a weighted graph $\Gamma$, dual harmonic embedding $\mathcal{H}^*$ of $\Gamma^*$ (with modified outer vertex), and the superposition graph $\Gamma \cup \Gamma^*$. Note that we do not fix the embedding of $\Gamma \cup \Gamma^*$ into $\mathbb{C}$.</figcaption></figure>

space $\mathbb{R}^3$. The second is called a *t-surface*: this is a polygonal surface $\Theta_\delta$ embedded into $\mathbb{R}^{2,2} \cong \mathbb{C}^{1,1}$, which is obtained as the lift to this Minkowski space of the *t-embedding* associated with the harmonic embedding as described in [CLR23, Section 8.1]. The structure of $\Theta_\delta$ plays a crucial role in the discrete complex analysis technique developed in [CLR23] upon which our approach is based. In this section we focus on constructing this t-surface and reviewing basic facts about it. For shortness, we drop the index $\delta$ until the end of this section.

**2.1. Definition of the t-surface.** Let $(\Gamma, c)$ be a weighted finite planar graph with a chosen outer face. We denote by $\Gamma^*$ the dual graph of $\Gamma$ modified as follows: replace the vertex $v_{\text{out}}^*$ of $\Gamma^*$ corresponding to the outer face of $\Gamma$ with $n = \deg v_{\text{out}}^*$ vertices $v_{\text{out},1}, \dots, v_{\text{out},n}$ of degree 1; see Figure 1 for an example. Note that we still have a bijection between the edges of $\Gamma$ and $\Gamma^*$.

We define the superposition graph $\Gamma \cup \Gamma^*$ to be the bipartite graph whose vertices correspond to vertices of $\Gamma$, vertices of $\Gamma^*$ and midpoints of edges of $\Gamma$, and whose edges correspond to half-edges of $\Gamma$ and $\Gamma^*$ (see Figure 1). We call the two classes of vertices black and white and say that a vertex of $\Gamma \cup \Gamma^*$ is black if it is a vertex of $\Gamma$ or $\Gamma^*$ and white otherwise.

Throughout the rest of the section we will assume that $\Gamma$ is embedded harmonically. Denote by $\mathcal{H}: \Gamma \to \mathbb{C}$ the map that sends a vertex of $\Gamma$ to its position in the complex plane. Let $\mathcal{H}^*: \Gamma^* \to \mathbb{C}$ be the dual embedding, that is, for each edge $v_1v_2$ and the corresponding dual edge $v_1^*v_2^*$ we have

$$ \mathcal{H}^*(v_2^*) - \mathcal{H}^*(v_1^*) = ic_{v_1v_2}(\mathcal{H}(v_2) - \mathcal{H}(v_1)) $$

if $v_1^*$ is on the right of $v_1v_2$.

Let us emphasize that, unlike in Figure 1, the map $\mathcal{H}^*$ is not necessary a proper embedding (unless the boundary polygon of $\mathcal{H}$ is required to be convex). Even more importantly, $\mathcal{H}$ and $\mathcal{H}^*$ are not necessarily 'aligned', i.e., it may be not possible to shift $\mathcal{H}^*$ so that each vertex $\mathcal{H}^*(v^*)$ belongs to the corresponding face of $\Gamma$. In other words, embeddings $\mathcal{H}$ and $\mathcal{H}^*$ do not induce an embedding of the superposition graph $\Gamma \cup \Gamma^*$ in general. Below we will show that they nevertheless define an embedding of the corner graph of $\Gamma$, which coincides with the dual graph of the superposition graph $\Gamma \cup \Gamma^*$ away from the boundary. Let us define this corner graph formally:

**Definition 2.1.** The corner graph $\mathcal{V}$ of $\Gamma$ is the graph whose vertices correspond to pairs $(v, v^*)$ of incident vertices of $\Gamma$ and $\Gamma^*$. Two vertices $(v_1, v_1^*)$ and $(v_2, v_2^*)$ are linked by an edge of $\mathcal{V}$ if either $v_1 = v_2$ and $v_1^* \sim v_2^*$ or $v_1 \sim v_2$ and $v_1^* = v_2^*$; see Figure 2 for an example.
<figure><img src="image_4.png" /><figcaption>FIGURE 2. From left to right: an example of a harmonic embedding $\mathcal{H}$ of a weighted graph $\Gamma$, the corner graph $\mathcal{V}$ of $\Gamma$, and the t-embedding of $\mathcal{V}$ constructed out of $\mathcal{H}$.</figcaption></figure>

Note that there is a bijection between the edges of $\mathcal{V}$ and the edges of $\Gamma \cup \Gamma^*$. Also, there is a bijection between faces of $\mathcal{V}$ and all vertices of $\Gamma \cup \Gamma^*$ except the black boundary ones. The latter correspondence allows us to color the faces of $\mathcal{V}$ in black and white accordingly to the coloring of vertices of $\Gamma \cup \Gamma^*$; see Figure 2 for an example. We denote by $B(\mathcal{V})$ and $W(\mathcal{V})$ the sets of black and white faces of $\mathcal{V}$, respectively. The next definition follows [CLR23, Section 8.1]:

**Definition 2.2.** T-embedding $\mathcal{T}: \mathcal{V} \to \mathbb{C}$ of the corner graph is the map $\mathcal{T}(v, v^*) = \frac{1}{2}(\mathcal{H}(v) + \mathcal{H}^*(v^*))$.

An example of a t-embedding $\mathcal{T}$ is depicted in Figure 2. It is easy to see that for each black face $b \in B(\mathcal{V})$ the polygon $2\mathcal{T}(b)$ is a translation of the face of $\Gamma$ or $\Gamma^*$ corresponding to $b$. If $w \in W(\mathcal{V})$ is a white face of $\mathcal{V}$, then $2\mathcal{T}(w)$ is a rectangle whose sides are translations of the edge of $\Gamma$ and the edge of $\Gamma^*$ corresponding to $w$, and whose orientation agrees with the orientation of the neighboring black faces. From these definitions, it is easy to see that $\mathcal{T}$ is always locally proper; moreover, it is globally proper if we assume that the boundary of $\mathcal{H}(\Gamma)$ is convex.

Following the terminology of [CLR23], each t-embedding $\mathcal{T}$ has an *origami map* $\mathcal{O}$ associated with it. Informally speaking, the mapping $\mathcal{O}$, defined up to rotations and translations, is obtained from the map $\mathcal{T}$ by folding the plane along all the edges of $\mathcal{T}$ consequently (one can prove that this procedure is consistent). Note that the images of a given face of $\mathcal{V}$ in $\mathcal{T}$ and in $\mathcal{O}$ are isometric to each other, and that these images have the same orientation for white faces and the opposite ones for black faces. In the setup of this paper one can use an explicit definition

$$ (2.1) \qquad \mathcal{O}: \mathcal{V} \to \mathbb{C}, \qquad \mathcal{O}(v, v^*) = \frac{1}{2}(\overline{\mathcal{H}^*(v^*)} - \overline{\mathcal{H}(v)}); $$

see also equation (2.9) below and [CLR23, Section 8.1] for more details.

The pair of maps $\mathcal{T}$ and $\mathcal{O}$ gives rise to an embedding of the graph $\mathcal{V}$ into the Minkowski space $\mathbb{C}^{1,1} \cong \mathbb{R}^{2,2}$. We call this embedding a *t-surface* and denote it by $\Theta$, that is,

$$ (2.2) \qquad \Theta = (\mathcal{T}, \mathcal{O}): \mathcal{V} \to \mathbb{C}^{1,1}. $$

It is convenient to view $\mathcal{T}$ and $\mathcal{O}$ as piecewise linear functions on $\Theta$, which in its turn can be viewed as a polygonal surface in $\mathbb{C}^{1,1}$. It is easy to see that on each face of $\Theta$ we have

$$ (2.3) \qquad |d\mathcal{T}| = |d\mathcal{O}|; $$

see also equation (2.9) below. Recall that the function $\Phi$, defined by (1.3), is the convex piecewise linear potential associated with $\Gamma$ and that $\Psi = 2\partial_{\bar{w}}\Phi = \partial_x\Phi + i\partial_y\Phi$.
**Lemma 2.3.** (i) If $b \in B(\mathcal{V})$ corresponds to a vertex $v$ of $\Gamma$ and $p \in \Theta(b)$, then $(T - \overline{\mathcal{O}})(p) = \mathcal{H}(v)$. (ii) If $u \in W(\mathcal{V})$ corresponds to an edge of $\Gamma$ and $p \in \Theta(u)$, then $(T - \overline{\mathcal{O}})(p)$ lies on the image of this edge of $\Gamma$ in the harmonic embedding $\mathcal{H}$.

(iii) Let $b \in B(\mathcal{V})$ corresponds to a face of $\Gamma$. Assume that $p \in \Theta(b)$ and let $w = (T - \overline{\mathcal{O}})(p)$. Then, $w$ belongs to the image of this face of $\Gamma$ in $\mathcal{H}$ and one has

$$ (2.4) \qquad \mathcal{T}(p) = \frac{1}{2}(w + \Psi(w)), \quad \mathcal{O}(p) = \frac{1}{2}(\overline{\Psi(w)} - w). $$

*Proof.* This follows directly from the definition of $\mathcal{T}$ and $\mathcal{O}$; see also [CLR23, Section 8.1]. □

Let $\mathcal{W} = \mathcal{T} - \overline{\mathcal{O}}$. According to Lemma 2.3, the mapping $\mathcal{W}$ projects $\Theta$ onto the union of faces of the harmonic embedding $\mathcal{H}$ of $\Gamma$. As already mentioned above, the t-embedding $\mathcal{T}$ is not always globally proper (it may have overlaps if the boundary polygon $\partial[\mathcal{H}(\Gamma)]$ of the harmonic embedding $\mathcal{H}$ is not convex), which means that $\mathcal{T}$ is not necessarily an injective map from the t-surface $\Theta$ to $\mathbb{C}$. However, $\mathcal{T}$ is always locally one-to-one. Using (2.3), it is easy to make this statement quantitative:

**Lemma 2.4.** Let $p \in \Theta$ be an arbitrary point on the t-surface and $2r \le \text{dist}(\mathcal{W}(p), \partial[\mathcal{H}(\Gamma)])$. There is an open (in the topology induced on $\Theta$ from $\mathbb{C}^2$) connected subset $p \in B_{\Theta}(p, r) \subset \Theta$ such that $\mathcal{T}: B_{\Theta}(p, r) \to \mathcal{B}(\mathcal{T}(p), r)$ is a bijection. Moreover, $\mathcal{W}(B_{\Theta}(p, r)) \subset \mathcal{B}(\mathcal{W}(p), 2r)$.

*Proof.* Let $B_{\Theta}(p, r)$ be the connected component of $\mathcal{T}^{-1}(\mathcal{B}(\mathcal{T}(p), r))$ that contains $p$. For a point $w \in \mathcal{B}(\mathcal{T}(p), r)$, let $\ell_w$ be the straight segment connecting $\mathcal{T}(p)$ and $w$ and let $\ell'_w \subset \ell$ be the longest subsegment that contains $\mathcal{T}(p)$ and does not intersect $\mathcal{T}(\partial\Theta)$. Denote by $\gamma_w \subset \mathcal{T}$ the lift of $\ell'_w$ onto the t-surface. From the fact that $\mathcal{T}$ is locally proper it is easy to see that $\mathcal{T}$ is one-to-one on $\gamma_w$. It follows from 2.3 that $\text{diam}(\mathcal{W}(\gamma_w)) \le 2r$. If $\ell'_w \subsetneq \ell_w$, we would also have $\mathcal{W}(\gamma_w) \cap \mathcal{H}(\partial\Gamma) \ne \emptyset$, a contradiction. Let $B_{\Theta}(p, r)$ denote the union of paths $\gamma_w$ as above. By construction, $\mathcal{T}$ is a bijection between $B_{\Theta}(p, r)$ and $\mathcal{B}(w, r)$. As $\text{diam}(\mathcal{W}(\gamma_w)) \le 2r$, we have $\mathcal{W}(B_{\Theta}(p, r)) \subset \mathcal{B}(\mathcal{W}(p), 2r)$. □

**2.2. Equivalence of property (LIP) and property Lip($\kappa$, $\delta$) of the t-surface $\Theta$.** According to Lemma 2.3 we can view $\mathcal{T}$ and $\mathcal{O}$ as pullbacks of the functions $w \mapsto \frac{1}{2}(w + \Psi)$ and $w \mapsto \frac{1}{2}(\overline{\Psi(w)} - w)$ from the plane into which $\Gamma$ is embedded by $\mathcal{H}$ to the corresponding t-surface $\Theta$. (Note that, contrary to $\Psi$, the maps $\mathcal{T}$ and $\mathcal{O}$ are well-defined and continuous on $\Theta$.) If the potential $\Phi$ was smooth and uniformly convex, then one could easily prove that $\mathcal{O}$ is, locally, a $\kappa$-Lipschitz function of $\mathcal{T}$ for some $\kappa < 1$ that depends only on the constant $\lambda$ in (1.4) or, equivalently, on the constant $\varkappa$ in (1.5); see also Lemma 5.1. Moreover, $\Phi$ is uniformly convex if and only if this is true. In our setup $\mathcal{O}$ is locally a linear function of $\mathcal{T}$, which means that this $\kappa$-Lipschitzness property cannot hold at very small distances. However, it is reasonable to expect that property (CONV)/(LIP) implies that $\mathcal{O}$ is a $\kappa$-Lipschitz function of $\mathcal{T}$ starting from the scale $\delta$. The following definition appeared in [CLR23]:

**Definition 2.5.** Let $\delta > 0$ and $0 < \kappa < 1$ be given. We say that a t-surface $\Theta = (\mathcal{T}, \mathcal{O})$ satisfies property Lip($\kappa$, $\delta$) if for each set $U_{\Theta} \subset \Theta$ such that $\mathcal{T}$ restricted to $U_{\Theta}$ is one-to-one and $\mathcal{T}(U_{\Theta})$ is convex and for each $p_1, p_2 \in U$ one has

$$ (2.5) \qquad |\mathcal{O}(p_2) - \mathcal{O}(p_1)| \le \kappa |\mathcal{T}(p_2) - \mathcal{T}(p_1)| \quad \text{if } |\mathcal{T}(p_2) - \mathcal{T}(p_1)| \ge \delta. $$

The next lemma shows that the property Lip($\kappa$, $\delta$) is equivalent to property (LIP) in the same sense as (LIP) is equivalent (see Lemma 5.1) to property (CONV), i.e., $\mathcal{O}(\delta)$-away from the boundary.

**Lemma 2.6.** Let $\Theta$ and $\Psi$ be the t-surface and the map $\Psi = 2\partial_w\Phi$, respectively, corresponding to $\mathcal{H}$. (i) Assume that $\Theta$ satisfy property Lip($\kappa$, $\delta$). Then, $\Psi$ has property (LIP) with $\varkappa = \frac{1-\kappa}{1+\kappa}$ and $C = 2$ provided that the points $w_1, w_2$ in (1.5) are $8(1-\kappa)^{-1}\delta$-away from the boundary of $\mathcal{H}$.

(ii) Vice versa, if $\Psi$ has property (LIP), then $\Theta$ satisfy Lip($\kappa$, $C'\delta$) for each $\kappa > \frac{1-\kappa}{1+\kappa}$ provided that the points $p_1, p_2$ in (2.5) are $C'\delta$-away from $\partial\Theta$, where $C'$ depends on $\varkappa$, $\kappa$, and $C$ in (LIP) only.
Proof. (i) Let $r = 4(1 - \kappa)^{-1}\delta$, recall the notation $\mathcal{W} = \mathcal{T} - \overline{\mathcal{O}}$ and note that Lemma 2.3 yields
$$ (2.6) \qquad \Psi \circ \mathcal{W} = \mathcal{T} + \overline{\mathcal{O}} $$

if we restrict these maps to the black faces of $\Theta$ corresponding to faces of $\Gamma$. Let $p_1 \in \Theta$, $z_1 = \mathcal{T}(p_1)$, and $w_1 = \mathcal{W}(p_1)$ be such that $\mathrm{dist}(w_1, \partial[\mathcal{H}(\Gamma)]) \ge 2r$ and consider the neighborhood $B_\Theta = B_\Theta(p_1, r)$ from Lemma 2.4; recall that $\mathcal{T}: B_\Theta \to B(z_1, r)$ is a bijection and that $\mathcal{W}(B_\Theta) \subset B(w_1, 2r)$. It is also easy to see that the Lipschitz property (2.5) implies that $\mathcal{W}(B_\Theta) \supset B(w_1, 4\delta)$ since the image of the circle $|z - z_1| = r = 4(1 - \kappa)^{-1}\delta$ under the map $\mathcal{W} \circ \mathcal{T}^{-1}$ stays at least $4\delta$-away from $w_1$ and encircles $w_1$ due to topological reasons.

Let $w_2 \in B(w_1, 4\delta) \setminus B(w_1, 2\delta)$; without true loss of generality we can also assume that $w_1$ and $w_2$ lie in the interiors of faces of $\mathcal{H}$. Denote $p_2 = \mathcal{W}^{-1}(w_2) \in B_\Theta$ and $z_2 = \mathcal{T}(p_2) \in B(z_1, r)$. Note that $|z_2 - z_1| \ge \delta$ as the mapping $\mathcal{W} \circ \mathcal{T}^{-1}$ is 2-Lipschitz. It is now easy to see from (2.6) and (2.5) that property (1.5) with $\varkappa = (1 - \kappa)/(1 + \kappa)$ holds for these $w_1$ and $w_2$. The general case follows by linearity: split the segment $[w_1; w_2]$ into subsegments of length between $2\delta$ and $4\delta$.

(ii) Let $r = C'\delta$, where $C' > C(1 - \kappa)^{-1}$ is a large enough constant. First, assume that $p_1 \in \Theta$ belongs to a black face corresponding to a face of $\Gamma$ and $\mathrm{dist}(\mathcal{T}(p_1), \partial[\mathcal{H}(\Gamma)]) \ge 2r$. Let $p_2 \in B_\Theta = B_\Theta(p_1, r)$ be another such point with $|\mathcal{T}(p_2) - \mathcal{T}(p_1)| \ge \frac{1}{2}r$. Denote $w_1 = \mathcal{W}(p_1)$, $w_2 = \mathcal{W}(p_2)$ as above, and recall that $w_2 \in \mathcal{W}(B_\Theta) \subset B(w_1, 2r)$. It follows from (2.6) and the upper bound in (1.5) that the image of the circle $|w - w_1| = C\delta$ under the map $\mathcal{T} \circ \mathcal{W}^{-1}$ lies inside $B(p_1, C(1 - \varkappa)^{-1}\delta)$. Therefore, $|w_2 - w_1| \ge C\delta$. It is now easy to deduce from (2.4) and the lower bound in (1.5) that

$$ \left| \frac{\mathcal{O}(p_2) - \mathcal{O}(p_1)}{\mathcal{T}(p_2) - \mathcal{T}(p_1)} \right| = \left| \frac{(\Psi(w_2) - \Psi(w_1)) - (w_2 - w_1)}{(w_2 - w_1) + (\Psi(w_2) - \Psi(w_1))} \right| \le \frac{1 - \varkappa}{1 + \varkappa}. $$

It remains to consider points $p_1, p_2 \in \Theta$ that does not necessarily belong to black faces of $\Theta$ that correspond to faces of $\Gamma$ but still satisfy $p_2 \in B_\Theta(p_1, r)$ and $|\mathcal{T}(p_2) - \mathcal{T}(p_1)| \ge \frac{1}{2}r$. Replacing these points by nearby points on such faces and using (1.5) again, it is easy to see that

$$ |\mathcal{O}(p_2) - \mathcal{O}(p_1)| \le \frac{1 - \varkappa}{1 + \varkappa} |\mathcal{T}(p_2) - \mathcal{T}(p_1)| + 2C(1 + \varkappa^{-1})\delta. $$

Therefore, for each $\kappa > \frac{1-\varkappa}{1+\varkappa}$ one can find a large enough constant $C'$ such that the estimate (2.5) holds without additional assumptions on $p_1, p_2$ provided that $|\mathcal{T}(p_2) - \mathcal{T}(p_1)| \ge C'\delta$. $\square$

**2.3. Corollaries of property Lip($\kappa, \delta$) for random walks and harmonic functions on $\Gamma$.** In order to simplify the notation, from now onwards we usually identify vertices $v$ of an abstract weighted graph $(\Gamma, c)$ with their positions $\mathcal{H}(v)$ in the complex plane if no confusion arises.

**Proposition 2.7.** Let $\Gamma$ be a harmonic embedding and the $t$-surface $\Theta$ be defined as above. Let $X_t$ be the continuous time simple random walk on $\Gamma$ parameterized so that $|X_t|^2 - t$ is a local martingale. Denote by $X_t^v$ the random walk started at $X_0^v = v$ and stopped at $\partial\Gamma$. Assume that $\Theta$ has property Lip($\kappa, \delta$) for some $0 < \kappa < 1$ and $\delta > 0$ (or, equivalently, $\Gamma$ has properties (CONV)/(LIP)). Then there exists a constant $C = C(\kappa) > 0$ such that for each $t \ge C\delta^2$, each vertex $v$ of $\Gamma$ such that $\mathrm{dist}(v, \partial\Gamma) \ge \sqrt{t}$, and each $\theta \in \mathbb{R}$ we have

$$ (2.7) \qquad \mathrm{Var}(\mathrm{Re}(e^{i\theta} X_t^v)) \ge C^{-1}t. $$

The same assertion holds for the dual graph $\Gamma^*$.

Proof. As explained in [CLR23, Section 8.1], $\Gamma$ and $\Gamma^*$ are T-graphs obtained from $\Theta = (\mathcal{T}, \mathcal{O})$ as projections $\mathcal{T} - \overline{\mathcal{O}}$ and $\mathcal{T} + \overline{\mathcal{O}}$ (cf. Lemma 2.3). The ellipticity estimate (2.7) for the random walk on such T-graphs then follows from [CLR23, Proposition 6.4]. To make this entirely clear, note that since (2.7) is local and $\Theta$ satisfies Lip($\kappa, \delta$), it is enough to assume that the projection $\Theta \mapsto \mathcal{T}$ is one-to-one (cf. the proof of Lemma 2.6), that is, $\mathcal{T}$ is a proper t-embedding as defined in [CLR23]; then $\mathcal{O}$ is a concrete instance of its origami map. If we switch colors of faces of $\mathcal{T}$, then $\pm \overline{\mathcal{O}}$ will represent
the origami of the obtained t-embedding (see [CLR23, Section 2.2]), to which we can apply [CLR23, Proposition 6.4] to get the desired estimates for $\mathcal{T} - \bar{\mathcal{O}}$ and $\mathcal{T} + \bar{\mathcal{O}}$. $\square$

Remark 2.1. Instead of referring to [CLR23, Proposition 6.4] one can give a self-contained proof of Proposition 2.7 using the notion of discrete extremal length; see Remark 4.3 for further comments.

The following lemmas are standard corollaries of Proposition 2.7; see [CLR23, Sections 6.3, 6.4]):

**Lemma 2.8 (crossing estimates).** Under the assumptions of Proposition 2.7 there exists a constant $C = C(\kappa) > 0$ such that for each $r \ge C\delta$, each vertex $v$ of $\Gamma$ such that $\mathrm{dist}(v, \partial\Gamma) \ge Cr$, and for each $\theta_0 \in \mathbb{R}$ we have

$$ \mathbb{P} \left[ \begin{array}{l} \text{random walk } X_t^v \text{ exists the disc } B(v,r) \\ \text{through the arc } \{v+e^{i\theta}r, |\theta-\theta_0| \le \pi/4\} \end{array} \right] \ge C^{-1}. $$

The same assertion holds for the dual graph $\Gamma^*$.

**Lemma 2.9 (Harnack inequality).** Under the assumptions of Proposition 2.7 there exists a constant $C = C(\kappa) > 0$ such that for each $r \ge C\delta$, each vertex $v$ of $\Gamma$ such that $\mathrm{dist}(v, \partial\Gamma) \ge 2r$ is covered by faces of $\Gamma$, and for each non-negative discrete harmonic function $H$ on $B(v, 2r)$ we have

$$ \max_{v \in B(v,r)} H(v) \le C \min_{v \in B(v,r)} H(v). $$

The same assertion holds for the dual graph $\Gamma^*$.

**Lemma 2.10 (Hölder-type decay of oscillations).** Under the assumptions of Proposition 2.7 there exist constants $\beta = \beta(\kappa) > 0$ and $C = C(\kappa) > 0$ such that for each $R \ge r \ge C\delta$, each vertex $v$ of $\Gamma$ such that $\mathrm{dist}(v, \partial\Gamma) \ge CR$, and for each discrete harmonic function $H$ on $B(v, R)$ we have

$$ \mathrm{osc}_{B(v,r)} H \le C \left(\frac{r}{R}\right)^{\beta} \mathrm{osc}_{B(v,R)} H. $$

The same assertion holds for the dual graph $\Gamma^*$.

Note that property Lip($\kappa$, $\delta$) implies item (a) in property (RW) due to Lemma 2.8. Let us show that Lip($\kappa$, $\delta$) also implies item (b) in the latter property.

**Lemma 2.11.** Let $\Gamma$ be a harmonic embedding and the t-surface $\Theta$ be defined as above. Assume that $\Theta$ has property Lip($\kappa$, $\delta$) for some $0 < \kappa < 1$ and $\delta > 0$ (or, equivalently, $\Gamma$ has properties (CONV)/(LIP)). There exists a constant $C = C(\kappa) > 0$ such that for each $r \ge C\delta$ and each vertex $v$ of $\Gamma$ such that $\mathrm{dist}(v, \partial\Gamma) \ge Cr$ we have

$$ C^{-1}r^2 \le \sum_{v_1, v_2 \in \Gamma \cap B(v,r): \ v_1 \sim v_2} c_{v_1 v_2} |v_1 - v_2|^2 \le Cr^2. $$

*Proof.* Since the required estimate is local, it is enough to assume that the projection $\Theta \mapsto \mathcal{T}$ is one-to-one; see Lemma 2.4). Given an edge $v_1v_2$ of $\Gamma$ denote by $u(v_1v_2) \in W(\mathcal{V})$ the corresponding white face of $\mathcal{V}$. It follows from Definition 2.2 that

$$ c_{v_1v_2}|v_1 - v_2|^2 = 4\mathrm{Area}(\mathcal{T}(u(v_1v_2))). $$

The claim now follows from [CLR23, Eq. (6.1) and Lemma 6.3]. $\square$

**2.4. T-holomorphic functions.** Let $\Gamma$ be a weighted planar graph and $\Gamma^*$ be its dual. There is a standard way to define a discrete holomorphic function on the superposition of $\Gamma$ and $\Gamma^*$ by declaring it to be a pair of $H: \Gamma \to \mathbb{R}$ and $H^*: \Gamma^* \to \mathbb{R}$ that are harmonic conjugate to each other. Assume now that $\Gamma$ is embedded harmonically and $\mathcal{T}$ is the corresponding t-embedding; see Section 2.1. In this case the aforementioned discrete holomorphic functions correspond to t-white-holomorphic functions on $\mathcal{T}$ as defined in [CLR23]. It is convenient to introduce an alternative definition of the latter using the terminology developed in Section 2.1.
We begin by fixing an instance of the *origami square root function* $\eta$ associated with the t-embedding $\mathcal{T}$ and the *origami map* $\mathcal{O}$; see [CLR23, Section 2.2]. The function $\eta$ is defined on the set $B(\mathcal{V}) \cup W(\mathcal{V})$ of faces of $\mathcal{T}$. Recall that $\mathcal{V}$ is the corner graph and there is a correspondence between $B(\mathcal{V})$ on the one side and vertices of $\Gamma$ and those of $\Gamma^*$ on the other, as well as between $W(\mathcal{V})$ and edges of $\Gamma$, which are also in a one-to-one correspondence with edges of $\Gamma^*$. (See Definition 2.1 and the discussion below it.) Define the function $\eta$ as follows:

$$ (2.8) \qquad \begin{array}{ll} \eta_b = 1 & \text{if } b \in \Gamma \subset B(\mathcal{V}); \\ \eta_b = i & \text{if } b \in \Gamma^* \subset B(\mathcal{V}); \end{array} \qquad \eta_u = \pm i \frac{\overline{v_1 - v_2}}{|v_1 - v_2|} \quad \text{if } u \in W(\mathcal{V}) \text{ corresponds to an edge } v_1 v_2 \text{ of } \Gamma, $$

where the $\pm$ signs are chosen arbitrary. This function helps to relate the gradients of $\mathcal{T}$ and $\mathcal{O}$ viewed as piecewise linear mappings on the t-surface $\Theta = (\mathcal{T}, \mathcal{O})$ corresponding to $\mathcal{T}$:

$$ (2.9) \qquad d\mathcal{O}(p) = \begin{cases} \eta_u^2 d\mathcal{T}(p) & \text{if } p \text{ belongs to a white face } \Theta(u), \ u \in W(\mathcal{V}), \\ \bar{\eta}_b^2 d\overline{\mathcal{T}}(p) & \text{if } p \text{ belongs to a black face } \Theta(b), \ b \in B(\mathcal{V}). \end{cases} $$

Given a function $F : B(\mathcal{V}) \to \mathbb{C}$, let $Fd\mathcal{T}$ be the 1-form defined on edges of $\Theta$, which equals $F(b) d\mathcal{T}$ on each edge incident to a face $\Theta(b)$, $b \in B(\mathcal{V})$, of the t-surface $\Theta$.

**Definition 2.12.** A function $F : B(\mathcal{V}) \to \mathbb{C}$ is called t-white-holomorphic if for each $b \in B(\mathcal{V})$ we have $F(b) \in \eta_b \mathbb{R}$ and the 1-form $Fd\mathcal{T}$ defined on edges of $\Theta$ is closed.

The assertion that the 1-form $Fd\mathcal{T}$ is closed implies a non-trivial relation on $F$ for each $u \in W(\mathcal{V})$. Namely, let $u = u(v_1v_2)$ corresponds to an edge $v_1v_2$ of $\Gamma$. As usual, we denote by $v_1^*v_2^*st$ the dual edge of $\Gamma^*$ oriented such that $v_1^*$ is on the right of $v_1v_2$. By definition of $\mathcal{T}$, the fact that the form $Fd\mathcal{T}$ is closed around $u$ can be written as

$$ (-F(v_1) + F(v_2)(\mathcal{H}^*(v_2^*) - \mathcal{H}^*(v_1^*)) + (F(v_1^*) - F(v_2^*)))(\mathcal{H}(v_2) - \mathcal{H}(v_1)) = 0, $$

which is equivalent to saying that

$$ F(v_2^*) - F(v_1^*) = ic_{v_1v_2}(F(v_2) - F(v_1)) $$

since $\mathcal{H}^*(v_2^*) - \mathcal{H}^*(v_1^*) = ic_{v_1v_2}(\mathcal{H}(v_2) - \mathcal{H}(v_1))$ by definition of the dual harmonic embedding. Taking into account that $F(b)$ is purely real when $b \in \Gamma$ and purely imaginary when $b \in \Gamma^*$ we conclude that the functions $H = F|_{\Gamma}$ and $H^* = -iF|_{\Gamma^*}$ are harmonic conjugate to each other. Thus, we obtained the definition of discrete holomorphicity that we have started with. We address the reader to [CLR23, Section 8.1] for more comments.

By definition given above, a t-white-holomorphic function $F$ is a pair $(H, iH^*)$ of harmonic functions, whereas, by the analogy with the continuous holomorphicity, the actual holomorphic function should correspond to the sum $H + iH^*$. This expression does not make sense 'as is' because $H$ and $H^*$ are defined on different sets, but we can give it a meaning by extending $F$ to white faces of $\mathcal{V}$ as follows. Recall that white face of $\mathcal{V}$ are rectangles.

**Definition 2.13.** For each white face of $\mathcal{V}$ choose a diagonal splitting it into two triangles. Denote by $W_{\text{spl}}(\mathcal{V})$ the set of these triangles. Given a t-white-holomorphic function $F$, define the function $F^\circ : W_{\text{spl}}(\mathcal{V}) \to \mathbb{C}$ so that for each $u \in W_{\text{spl}}(\mathcal{V})$ incident to a black face $b \in B(\mathcal{V})$ we have

$$ F(b) = \text{Pr}(F^\circ(u), \eta_b \mathbb{R}). $$

Note that each triangle $u \in W_{\text{spl}}(\mathcal{V})$ has exactly two black faces $v \in \Gamma$ and $v^* \in \Gamma^*$ incident to it. By definition, $F(v) \in \eta_v \mathbb{R} = \mathbb{R}$ and $F(v^*) \in \eta_{v^*} \mathbb{R} = i\mathbb{R}$. Thus, we simply have

$$ F^\circ(u) = F(v) + F(v^*) \quad \text{if} \quad v \sim u \sim v^*. $$

We call the values of $F^\circ$ 'true complex values' of a t-white-holomorphic function $F$.
**Lemma 2.14.** Let $F: B \to \mathbb{C}$ be a $t$-white-holomorphic function and $F^\circ$ be its extension to $W_{\text{spl}}(\mathcal{V})$ defined above. Define a piecewise constant 1-form $F^\circ d\mathcal{T} + \overline{F^\circ d\mathcal{O}}$ on the $t$-surface $\Theta$ as follows:

$$ (F^\circ d\mathcal{T} + \overline{F^\circ d\mathcal{O}})(p) = \begin{cases} F^\circ(u) d\mathcal{T}(p) + \overline{F^\circ(u) d\mathcal{O}(p)} & \text{if } p \in \Theta(u), u \in W_{\text{spl}}(\mathcal{V}), \\ 2F(b) d\mathcal{T}(p) = F^\circ(u) d\mathcal{T}(p) + \overline{F^\circ(u) d\mathcal{O}(p)} & \text{if } p \in \Theta(b), b \in B(\mathcal{V}), \end{cases} $$

where in the second line one can take any white triangle $u$ incident to the face $b$. This piecewise constant form is well-defined, does not depend on the choices of $u \sim b$ due to (2.9), and is closed.

*Proof.* See [CLR23, Proposition 3.7 and Section 5]. □

In Section 4.4 and Section 4.5 we also use the notion of t-black-holomorphic functions on $\Theta$. Similarly to Definition 2.12, a t-black-holomorphic function $G: W(\mathcal{V}) \to \mathbb{C}$ is originally defined on white (rectangular) faces of $\Theta$ so that $G(u) \in \eta_u \mathbb{R}$ for each $u$ and the differential 1-form $Gd\mathcal{T}$ defined on edges of $\Theta$ as $G(u)d\mathcal{T}$ is closed. It follows from (2.8) that the imaginary part of the primitive of this form is constant on vertices $(vv^*) \in \mathcal{V}$ that share the same vertex $v \in \Gamma$ and similarly for the real part and $v^* \in \Gamma^*$. This allows one to consider the primitives

$$ H(v) := \int^v \mathrm{Im}[Gd\mathcal{T}], \quad H^*(v^*) := -\int^{v^*} \mathrm{Re}[Gd\mathcal{T}]. $$

By construction $H$ and $H^*$ are harmonically conjugate to each other. In particular, $H$ is harmonic on $\Gamma$ and $H^*$ is harmonic on $\Gamma^*$. In other words, t-black-holomorphic functions on $W(\mathcal{V})$ are, up to the multiple $i$, gradients of harmonic functions on $\Gamma$.

In order to define 'true complex values' $G^\bullet$ of $G$, we triangulate each black face of $\Theta$, extend $H$ and $H^*$ linearly to thus obtained triangles $b \in B_{\text{spl}}(\mathcal{V})$, and define $G^\bullet(b)$ to be the gradient of $H$ (if $b$ is a part of a face of $\Gamma^*$) or gradient of $H^*$ (if $b$ is a part of a face of $\Gamma$), viewed as a complex number and multiplied by $i$. Similarly to Lemma 2.14, this allows one to introduce a closed differential 1-form $G^\bullet d\mathcal{T} + \overline{G^\bullet d\mathcal{O}}$ on $\Theta$. We address the reader to [CLR23, Section 8.1] for more details.

### 3. Elliptic equations and closed 1-forms in 2D

In this section we collect some facts about linearized Monge-Ampère equation that we will need later. This is a classical topic at the crossroad of the theory of elliptic PDE and the optimal transport problem. We do not attempt to summarize the literature and only mention the foundational paper [CG97] and the monograph [Fig17] where the link with the transport problem is addressed. For our purposes it is enough to restrict the discussion to linearized Monge-Ampère equations with uniformly convex potentials. Such equations are uniformly elliptic which makes the standard theory of elliptic PDE applicable. To this end, we address the reader to [Ken+00] and references therein for the standard properties of elliptic operators in divergence form with measurable coefficients, and the book [AIM09] as a standard reference to the interplay between Beltrami and elliptic equations in 2D.

We begin with the notion of $A$-harmonic conjugate that can be related with any elliptic operator in 2D in divergence form. Let $\Omega \subset \mathbb{C}$ be a simply connected domain and $A$ be a $2 \times 2$ matrix whose coefficients are bounded real-valued functions in $\Omega$. In what follows we will use the notation $w = x + iy$ for the coordinate in $\mathbb{C}$. Consider the differential operator

$$ \mathcal{L}h = -\mathrm{div}(A\nabla h). $$

We assume that the operator $\mathcal{L}$ has uniformly bounded coefficients and that $\mathcal{L}$ is strongly elliptic: that is, there exists a constant $\lambda > 0$ such that for each $w \in \Omega$ and $\nu \in \mathbb{R}^2$ we have

$$ (3.1) \qquad \lambda|\nu|^2 \le \nu \cdot A(w)\nu \le \lambda^{-1}|\nu|^2. $$

In what follows we do not make any smoothness assumption on the coefficients of $A$. In particular, the equation $\mathcal{L}h = 0$ is understood in a weak form, that is, we say that $\mathcal{L}h = 0$ if the gradient of $h$ is
locally square integrable and for each test function $\xi \in C_0^\infty(\Omega)$ we have

$$ \int_{\Omega} (\nabla \xi \cdot A \nabla h) dxdy = 0. $$

Denote by $*$ the Hodge star operator in $\mathbb{R}^2$, that is,

$$ * = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}. $$

Classically [AIM09, Chapter 16.1.3], with every solution $h$ of the equation $\mathcal{L}h = 0$ we can associate its $A$-harmonic conjugate $h^*$ defined by

$$ (3.2) \qquad \nabla h^* = *A\nabla h; $$

note that the equation $\mathcal{L}h = 0$ is equivalent to $*A\nabla h$ being curl-free, that is to $*A\nabla h$ being the gradient of a function. Note also that

$$ \operatorname{div}(*A^{-1} * \nabla h^*) = - \operatorname{div}(*\nabla h) = 0 $$

hence $h^*$ belongs to the kernel of the adjoint operator

$$ \mathcal{L}^* h = \operatorname{div}(*A^{-1} * \nabla h). $$

For example, if $A = \mathrm{Id}_{2 \times 2}$, then $\mathcal{L}$ is the (positively defined) Euclidean Laplacian and $h^*$ is the usual harmonic conjugate of a harmonic function $h$. In this case $f = h + ih^*$ is holomorphic and the relation (3.2) is nothing but the Cauchy–Riemann equation.

Let us now focus on elliptic operators $\mathcal{L}_\varphi$ corresponding to the linearized Monge–Ampère equation. Recall that in this case the matrix $A$ takes the form

$$ (3.3) \qquad A = \begin{pmatrix} \varphi_{yy} & -\varphi_{xy} \\ -\varphi_{xy} & \varphi_{xx} \end{pmatrix} = \begin{pmatrix} \operatorname{Re}(\psi_w - \psi_{\bar{w}}) & \operatorname{Im}(\psi_w - \psi_{\bar{w}}) \\ -\operatorname{Im}(\psi_w + \psi_{\bar{w}}) & \operatorname{Re}(\psi_w + \psi_{\bar{w}}) \end{pmatrix} = \begin{pmatrix} v_y & -u_y \\ -v_x & u_x \end{pmatrix}, $$

where $\psi = u + iv = 2\partial_{\bar{w}}\varphi$ and $\varphi$ is a convex function. Note that $A$ is symmetric. The boundedness of coefficients of $A$ corresponds to Lipschitzness of $\psi$ and the ellipticity condition (3.1) corresponds to the existence of $a > 0$ and $\mu < 1$, such that

$$ (3.4) \qquad \psi_w \ge a, \quad |\psi_{\bar{w}}| \le \mu \psi_w $$

holds almost everywhere (note that $\psi$ is differentiable almost everywhere by Rademacher's theorem). The aforementioned properties of $\psi$ are equivalent to $\varphi$ being uniformly convex, that is, to the existence of $\lambda > 0$ such that for each segment $[w_1, w_2] \subset \Omega$ we have

$$ (3.5) \qquad \lambda|w_1 - w_2|^2 \le \varphi(w_2) - 2\varphi(\frac{1}{2}(w_1 + w_2)) + \varphi(w_1) \le \lambda^{-1}|w_1 - w_2|^2. $$

**Lemma 3.1.** Assume that $\varphi$ is a uniformly convex function in a domain $\Omega$ and $\psi = 2\partial_{\bar{w}}\varphi$. Assume that $h, h^*$ are continuous real-valued functions in $\Omega$ such that the differential form $h d\psi + ih^* dw$ is closed ‘in the weak sense’, i.e., that its integral over each piece-wise smooth contractible loop vanishes. Then the following holds:

(i) $h$ is a very weak solution of the equation $\mathcal{L}_\varphi h = 0$, i.e. for each test function $\xi \in C_0^2(\Omega)$ we have $\int_{\Omega} (h \cdot \mathcal{L}_\varphi \xi) dxdy = 0$.

(ii) Assume additionally that $h$ is differentiable almost everywhere and the gradient of $h$ is locally square integrable. Then $h^*$ is an $A$-harmonic conjugate of $h$.

*Proof.* Let $\xi \in C_0^1(\Omega)$ be a continuously differentiable real-valued test function. Integrating over level lines of $\xi$ we obtain the identity

$$ (3.6) \qquad \int_{\Omega} (h \, d\psi + ih^* \, dw) \wedge d\xi = 0. $$
Taking the real and the imaginary parts separately and plugging $\psi = \varphi_x + i\varphi_y$, $w = x + iy$ in, this can be rewritten as

$$ \int_{\Omega} (h(\varphi_{xx}\xi_y - \varphi_{xy}\xi_x) + h^*\xi_x) dxdy = 0 = \int_{\Omega} (h(\varphi_{xy}\xi_y - \varphi_{yy}\xi_x) + h^*\xi_y) dxdy. $$

Assume now that $\xi$ is twice differentiable. Using the first equation with $\xi$ replaced by $\xi_y$, the second with $\xi$ replaced by $\xi_x$, and subtracting the two, we see that

$$ (3.7) \qquad \int (h \cdot \mathcal{L}_{\varphi}\xi) dxdy = \int_{\Omega} (h \cdot (\varphi_{xx}\xi_{yy} - 2\varphi_{xy}\xi_{xy} + \varphi_{yy}\xi_{xx})) dxdy = 0 $$

for each real-valued test function $\xi \in C_0^2(\Omega)$, as required in (i).

To prove (ii), assume first that $h$ and $\psi$ are smooth. In this case, for each $\xi \in C_0^1(\Omega)$ we can use integration by parts and write

$$ (3.8) \qquad \int_{\Omega} hd\psi \wedge d\xi = - \int_{\Omega} \xi dh \wedge d\psi. $$

Both sides of the equation (3.8) are continuous with respect to the norm $\|h\|_{L^2(U)} + \|\nabla h\|_{L^2(U)} + \|\nabla\psi\|_{L^\infty(U)}$, where $U \subset \Omega$ is any open set containing $\text{supp} \, \xi$. Thus, (3.8) is actually valid even when $\psi$ is only Lipschitz and $\nabla h$ is locally square integrable, as in this case we can approximate $\psi$ and $h$ by smooth functions and use the equality (3.8) for them. We conclude that under the assumption made in (ii) the relation (3.6) reads as

$$ (3.9) \qquad \int_{\Omega} \xi dh \wedge h\psi = i \int_{\Omega} h^* dw \wedge d\xi. $$

A direct calculation using the explicit form (3.3) of $A$ shows that a $A$-harmonic conjugate $\tilde{h}^*$ of $h$ satisfies the same equation which implies that $h^* - \tilde{h}^*$ is locally constant, hence $h^*$ is also an $A$-harmonic conjugate. $\square$

One can ask if the property of $h$ being a very weak solution of $\mathcal{L}_{\varphi}h = 0$ in the sense of Lemma 3.1(i) implies that $h$ is differentiable. In the next lemma we give a simple proof of this fact for $\varphi \in C^3(\Omega)$.

**Lemma 3.2.** Assume that $\varphi \in C^3(\Omega)$ and $h \in C(\Omega)$ is such that for each $\xi \in C_0^2(\Omega)$ we have $\int_{\Omega}(h \cdot \mathcal{L}_{\varphi}\xi) dxdy = 0$. Then $h \in C^2(\Omega)$ and $\mathcal{L}_{\varphi}h = 0$.

*Proof.* Consider a ball $B = B(p, r)$ such that $\bar{B} \subset \Omega$ and let $h_0$ be the solution of the equation $\mathcal{L}_{\varphi}h_0 = 0$ in $B$ with continuous boundary values $h_0 = h$ on $\partial B$. Since we assume that $\varphi \in C^3$, the matrix $A$ has $C^1$-smooth coefficients. Therefore, this strong solution exists and is $C^2$-smooth inside $B$. Thus, it remains to prove that $h = h_0$ in $B$ or that $\int_B((h - h_0) \cdot \phi) dxdy = 0$ for all $\phi \in C_0^\infty(B)$.

To this end, note that the difference $h - h_0$ is continuous in $\bar{B}$, has zero boundary values on $\partial B$ and satisfies (3.7) for each twice differentiable compactly supported test function $\xi$ in $B$. Now let $\xi$ be the (strong) solution of the nonhomogeneous equation $\mathcal{L}_{\varphi}\xi = \phi$ in $B$ with zero boundary values. As $\phi \in C_0^\infty(B)$, classical Schauder estimates [GT01, Chapter 6] imply that $\xi$ is $C^2$-smooth inside $B$ and has uniformly bounded gradient in $B$ including near $\partial B$. Finally, for sufficiently small $\varepsilon > 0$ denote $\xi_\varepsilon = \eta_\varepsilon\xi$, where $\eta_\varepsilon : B \to [0, 1]$ is a smooth function such that $\eta_\varepsilon = 1$ in $B(p, (1-2\varepsilon)r)$, $\eta_\varepsilon = 0$ outside $B(p, (1-\varepsilon)r)$, and $|D^2\eta_\varepsilon| = O(\varepsilon^{-2})$. Since $\xi_\varepsilon \in C_0^2(B)$, by applying (3.7) to $h - h_0$ and $\xi_\varepsilon$ we see that

$$ (3.10) \qquad \int_B ((h - h_0) \cdot \phi) dxdy = \int_B ((h - h_0) \cdot \mathcal{L}_{\varphi}\xi) dxdy = \int_B ((h - h_0) \cdot \mathcal{L}_{\varphi}[(1 - \eta_\varepsilon)\xi]) dxdy. $$

Note that $[\mathcal{L}_{\varphi}(1 - \eta_\varepsilon)\xi](w) \neq 0$ only if $\varepsilon \le \text{dist}(w, \partial B) \le 2\varepsilon$. For such $w$ we have $[\mathcal{L}_{\varphi}\xi](w) = 0$ and $\xi(w) = O(\varepsilon)$ as the gradient of $\xi$ is uniformly bounded in $\bar{B}$. Hence, $\mathcal{L}_{\varphi}[(1 - \eta_\varepsilon)\xi](w) = O(\varepsilon^{-1})$, which allows us to estimate the integral (3.10) by the maximum of $|h - h_0|$ in the $\varepsilon$-neighborhood of $\partial B$. Since $h - h_0$ is continuous in $\bar{B}$, taking the limit as $\varepsilon \to 0$ completes the proof of $h = h_0$. $\square$
Classically (see [LSW63, Section 6] or [TKB13] and references therein), the operator $\mathcal{L}_{\varphi}$ admits a Green's function in a domain $\Omega$, that is, a function $G(w_1, w_2)$ that solves the equation

$$ \mathcal{L}_{\varphi} G(\cdot, w_2) = \delta_{w_2} $$

(where $\delta_{w_2}$ is the delta measure at $w_2$) for each $w_2 \in \Omega$ and vanishes along the boundary of $\Omega$. The standard theory of elliptic equations imply that $G$ is symmetric, positive on $\Omega \times \Omega \setminus \text{diag}$, and continuous on $\bar{\Omega} \times \bar{\Omega} \setminus \text{diag}$. Moreover, any non-negative continuous function $\tilde{G}(\cdot, w_2)$ in $\bar{\Omega} \setminus \{w_2\}$ that vanishes at the boundary of $\Omega$ and satisfies the equation $\mathcal{L}_{\varphi}\tilde{G}(\cdot, w_2) = 0$ coincides with $G(\cdot, w_2)$ up to a multiplicative constant. Indeed, standard application of Harnack inequality (see e.g. [LSW63, Section 7]) implies that such $\tilde{G}(w_1, w_2)$ has the same growth as $G(w_1, w_2)$ when $w_1 \to w_2$, which implies the assertion due to a classical result from [SW66]. The multiplicative constant can be fixed by computing the monodromy of the A-harmonic conjugate:

**Lemma 3.3.** Let $G$ be the Green's function of the operator $\mathcal{L}_{\varphi}$. Fix an arbitrary $w_2 \in \Omega$ and let $G^*(w_1, w_2)$ be an A-harmonic conjugate of $G(w_1, w_2)$ in the first variable. Then $G^*(w_1, w_2)$ is additively multivalued and its monodromy around $w_2$ in the counterclockwise direction is equal to $-1$.

*Proof.* By the standard continuity arguments (e.g., see [LSW63, Section 5]) it is enough to assume that $\varphi$ is smooth, thus, by Schauder estimates, $G$ is smooth as well. Fix a smooth function $\xi: \mathbb{R} \to \mathbb{R}$ with a small support that is constant in a vicinity of 0. Let $m$ be the monodromy of $G^*(\cdot, w_2)$. Then, using the definition of A-harmonic conjugate and the definition of the Hodge star $*$, and integrating along concentric circles around $w_2$ we can compute

$$ -m\xi(0) = \int_{\Omega} \nabla G^*(w, w_2) \cdot *\nabla \xi(|w - w_0|) d^2w = \int_{\Omega} A\nabla G(w, w_2) \cdot \nabla \xi(|w - w_0|) d^2w = \xi(0) $$

which implies $m = -1$. $\square$

<figure><img src="image_5.png" /></figure>

## 4. Convergence of discrete harmonic functions on $\Gamma_{\delta}$ as $\delta \to 0$

The goal of this section is to prove Theorem 2 and Theorem 3. Our strategy can be described as follows. First we observe that property (LIP), Lemma 2.6 and Proposition 2.7 imply that the random walk on $\Gamma_{\delta}$ is uniformly elliptic starting from the scale $\delta$. This immediately implies (see Lemma 2.10) that the family $(H_{\delta})$ is precompact in $C^0$ (and in fact $C^{\alpha}$ for a small $\alpha > 0$) topology and that subsequential limits have correct boundary conditions. Further, discrete Dirichlet energy of $H_{\delta}$ can be bounded on compacts due to a discrete Caccioppoli estimate which, in its turn, allows to estimate harmonic conjugate functions $H_{\delta}^*$ and to conclude that the family $(H_{\delta}^*)$ is precompact in $C^0$ topology as well. Passing to a subsequence and applying Lemma 2.14 we obtain a pair of continuous functions $h, h^*$ such that the 1-form $h \, d\psi + ih^* \, dw$ is closed. This allows us to recover the coefficients of the elliptic equation that $h$ satisfies by using results from Section 3.

### 4.1. Oscillations estimate via the Dirichlet energy.
In this section we recall a well-known equicontinuity estimate for discrete harmonic functions via their Dirichlet energies; see Lemma 4.3 and Remark 4.2 below. In its simplest form, this estimate dates back to Lusternik [Lyu26]; see also [Sko13; GJN20] and references therein.

In what follows we assume that $\Gamma_{\delta}$ is a family of graphs embedded harmonically and satisfying property (LIP). Given a set $U \subset \mathbb{C}$ we denote by $U_{\delta}$ the set of vertices of $\Gamma_{\delta}$ that belong to $U$ and by $\partial U_{\delta}$ the set of vertices of $\Gamma_{\delta}$ that do not belong to $U$ but have neighbors in $U_{\delta}$. We put $\bar{U}_{\delta} = U_{\delta} \cup \partial U_{\delta}$ and denote by $E(U_{\delta})$ the set of edges of $\Gamma_{\delta}$ with both ends in $U_{\delta}$. Also, we denote by $E(U_{\delta}^*)$ the set of dual to $E(U_{\delta})$ edges of $\Gamma_{\delta}^*$ and by $U_{\delta}^*$ the set of all vertices of $\Gamma_{\delta}^*$ incident to these edges.
**Definition 4.1.** Let $\delta > 0$ and $U \subset \mathbb{C}$ be fixed. Given a function $H$ on $U_\delta$ we define its Dirichlet energy over $U$ as

$$ \mathcal{E}_U^\delta(H) = \sum_{v_1 v_2 \in E(U_\delta)} c_{v_1 v_2} (H(v_1) - H(v_2))^2. $$

Similarly, given a function $H^*$ on $U_\delta^*$ we set $\mathcal{E}_U^{\delta,*}(H^*) = \sum_{v_1^* v_2^* \in E(U_\delta^*)} c_{v_1 v_2} (H^*(v_1^*) - H^*(v_2^*))^2$.

**Lemma 4.2.** Assume that $U \subset \mathbb{C}$ is bounded and that there exists $r_0 > 0$ such that the $r_0$-neighborhood of $U$ is covered by faces of $\Gamma_\delta$ for all $\delta > 0$ small enough. Let $g$ be a continuously differentiable function defined in a $\delta$-neighborhood of $U$ (and hence on all $U_\delta$). Then,

$$ \mathcal{E}_U^\delta(g) \le \|\nabla g\|_\infty^2 \cdot \sum_{v_1 v_2 \in E(U_\delta)} c_{v_1 v_2} |v_1 - v_2|^2. $$

*Proof.* This follows from the trivial estimate $|g(v_2) - g(v_1)| \le \|\nabla g\|_\infty \cdot |v_2 - v_1|$. $\square$

Given two subsets $A, B \subset U_\delta$, we denote by $\mathrm{EL}_U^\delta(A \leftrightarrow B)$ the (edge) extremal length of paths connecting $A$ and $B$ inside $U_\delta$; this quantity is also known as the effective resistance between $A$ and $B$ in the electrical network $U_\delta$. Similarly, given $A^*, B^* \subset U_\delta^*$, we denote by $\mathrm{EL}_U^{\delta,*}(A^* \leftrightarrow B^*)$ the extremal length of paths connecting $A^*$ and $B^*$ inside $U_\delta^*$. We address the reader to [LP16, Chapter 2] or to [Che16, Section 6], [BP25, Section 2.2] and references therein for the definition and basic properties of the discrete extremal length.

For $K \subset \mathbb{C}$ and a function $H^* : K_\delta^* \to \mathbb{R}$ define

$$ \mathrm{osc}_K H^* = \max_{K_\delta^*} H^* - \min_{K_\delta^*} H^*. $$

**Lemma 4.3.** Let $U \subset \mathbb{C}$ be an open simply connected set and assume that $U$ is covered by faces of $\Gamma_\delta$ for all $\delta > 0$ small enough. For each compact subset $K \subset U$ there exist constants $C, \delta_0 > 0$ such that for each $\delta \le \delta_0$ and every discrete harmonic in $U_\delta$ function $H : \bar{U}_\delta \to \mathbb{R}$ we have

$$ (4.1) \qquad \mathrm{osc}_K H^* \le C \sqrt{\mathcal{E}_U^\delta(H)}, $$

where $H^*$ is the harmonic conjugate to $H$. (Note that $H^*$ is uniquely defined, up to an additive constant, in the bulk of $U$ if $\delta$ is sufficiently small.)

*Proof.* It is enough to prove that for each $p \in U$ there exist $C, r > 0$ such that (4.1) holds for $K = \bar{B}(p, r)$ and small enough $\delta$. Fix $p$ and let $r = \frac{1}{2}\mathrm{dist}(p, \partial U)$. Put $A(p, r, 2r) = B(p, 2r) \setminus \bar{B}(p, r)$. Using the maximum principle for the function $H^*$ it is easy to show the existence of two paths $\gamma_M^*$ and $\gamma_m^*$ inside the discrete annulus $A(p, r, 2r)_\delta^*$ connecting its boundary components and such that

$$ H^*|_{\gamma_M^*} \ge \max_{B(p,r)_\delta^*} H^*, \qquad H^*|_{\gamma_m^*} \le \min_{B(p,r)_\delta^*} H^*. $$

The basic properties of the extremal length/effective resistance immediately imply that

$$ (\mathrm{osc}_{B(p,r)} H^*)^2 \le \mathcal{E}_{B(p,2r)}^{\delta,*}(H^*) \cdot \mathrm{EL}_{A(p,r,2r)}^{\delta,*}(\gamma_M^* \leftrightarrow \gamma_m^*) \le \mathcal{E}_{B(p,2r)}^{\delta}(H) \cdot \mathrm{EL}_{A(p,r,2r)}^{\delta,*}(\gamma_M^* \leftrightarrow \gamma_m^*), $$

where we also used the fact that the Dirichlet energies of the harmonic function $H$ and of its harmonic conjugate $H^*$ are equal. We can now bound $\mathrm{EL}_{A(p,r,2r)}^{\delta,*}(\gamma_M^* \leftrightarrow \gamma_m^*)$ from above by the extremal length of circuits separating two boundary components of the annulus. Classically (see [FF56, Theorem 1]), the latter equals the inverse of the extremal length $\mathrm{EL}_{A(p,r,2r)}^{\delta}(\partial B(p, r) \leftrightarrow \partial B(p, 2r))$, where $\partial B(p, r)$ denotes the set of vertices in $A(p, r, 2r)_\delta$ incident to vertices from $\bar{B}(p, r)_\delta$, and the same for $\partial B(p, 2r)$. Therefore,

$$ (\mathrm{osc}_{B(p,r)} H^*)^2 \le \mathcal{E}_{B(p,2r)}^{\delta}(H) \cdot (\mathrm{EL}_{A(p,r,2r)}^{\delta}(\partial B(p, r) \leftrightarrow \partial B(p, 2r)))^{-1}. $$
To estimate $\mathrm{EL}^\delta_{A(p,r,2r)}(\partial B(p,r) \leftrightarrow \partial B(p,2r))$ from below we can consider the metric given by the discrete gradient of the function $g_p(v) = \log|v-p|$:

$$ \mathrm{EL}^\delta_{A(p,r,2r)}(\partial B(p,r) \leftrightarrow \partial B(p,2r)) \geq \left( \min_{\partial B(p,2r)} g_p - \max_{\partial B(p,r)} g_p \right)^2 \cdot (\varepsilon^\delta_{A(p,r,2r)}(g_p))^{-1}. $$

This allows us to conclude that

$$ (4.2) \qquad (\mathrm{osc}_{B(p,r)} H^*)^2 \leq (\log 2)^{-2} \cdot \varepsilon^\delta_{A(p,r,2r)}(g_p) \cdot \varepsilon^\delta_{B(p,2r)}(H). $$

Due to Lemma 4.2 and Lemma 2.11, $\varepsilon^\delta_{A(p,r,2r)}(g_p)$ is bounded from above by a constant that depends on the constant in property (LIP) only (provided that $\delta$ is small enough). The claim easily follows. $\square$

**Remark 4.1.** Same arguments as in the proof of Lemma 4.3 can be used to estimate the oscillations of a harmonic function $H$ on the primal graph $\Gamma_\delta$. The subtle difference is that on the last step we need to estimate the dual extremal length $\mathrm{EL}^{\delta,*}_{A(p,r,2r)}(\partial B(p,r) \leftrightarrow \partial B(p,2r))$ of the annulus which is roughly the conformal modulus of $\Psi_\delta(A(p,r,2r))$. This conformal modulus is however comparable with the conformal modulus of $A(p,r,2r)$ as follows from Lip($\kappa, \delta$) that the t-surface $\Theta_\delta$ of $\Gamma_\delta$ satisfies (see Lemma 2.6).

**Remark 4.2.** A simple modification of the arguments from Lemma 4.3 imply that for every compact $K \subset U$ and a discrete harmonic function $H$ on $U_\delta$ one has

$$ (4.3) \qquad \omega_K(t, H) \leq C |\log t|^{-1/2} \cdot \sqrt{\varepsilon_U^\delta(H)} \quad \text{for all } t \geq C\delta, $$

where $\omega_K^H(t)$ is the modulus of continuity of $H$ on $K$ and a constant $C$ depends on $K, U$ and the constants in property (LIP) only. In fact, Proposition 2.7 gives a stronger Hölder-type estimate for $\omega_K^H(t)$ once we know that the oscillations of $H$ are bounded. However, as explained in the next remark, the estimate (4.3) can be bootstrapped to obtain a simple self-contained proof of Proposition 2.7.

**Remark 4.3.** Note that the arguments used above to derive the inequality (4.3) can be applied to any (not necessary harmonic) function $H$ satisfying the following property: for each $v \in U_\delta$ there exist two paths $\gamma_+, \gamma_-$ connecting $v$ to $\partial U_\delta$ such that $H(v_+) \geq H(v) - \mathrm{cst} \cdot \delta$ for all $v_+ \in \gamma_+$ and $H(v_-) \leq H(v) + \mathrm{cst} \cdot \delta$ for all $v_- \in \gamma_-$. This observation can be used to give an independent proof of the ellipticity estimate from Proposition 2.7; let us sketch it here. Pick $v_0 \in U_\delta$ and consider the square $Q = v_0 + [-L\delta, L\delta]^2$, where $L$ is a sufficiently large constant. Assume that a non-negative smooth function $f$ equals 1 at the middle of the right side of $Q$ (that is, at $v_0 + L\delta$) and vanishes in an $O(\delta)$-neighborhood of other three sides. Let a function $H$ be harmonic in $Q_\delta$ and equal to $f$ outside $Q$. Note that the Dirichlet energy $\varepsilon^\delta_{B(v_0,2L\delta)}(H) \leq \varepsilon^\delta_{B(v_0,2L\delta)}(f)$ is bounded by an absolute constant and one can choose $f$ so that $H$ satisfies the aforementioned property on the existence of paths $\gamma_\pm$ emanating from each $v$ and running to $\partial B(v_0, 2L\delta)$. Then, the estimate (4.3) implies that $H$ is bounded away from zero in a small vicinity of the point $v_0 + L\delta$ which implies a uniform lower bound on the probability that a random walk started in this vicinity exits $Q$ through its right side. Finally, one can show that such a weaker version form of Lemma 2.8 is ultimately equivalent to Proposition 2.7.

**4.2. Discrete Caccioppoli estimate.** Caccioppoli's theorem is a classical result in the theory of elliptic equations asserting that the $L^2$ norm over a ball $B(p,r)$ of the gradient of a solution to an elliptic equation can be controlled via the $L^2$ norm of the solution itself over a larger ball $B(p,2r)$. This admits a simple analogue for discrete harmonic functions.

**Proposition 4.4.** Let $\Gamma_\delta$ be a weighted planar graph embedded harmonically into plane. There exists an absolute constant $C > 0$ such that the following holds. Let $r > 0$ and $p \in \mathbb{C}$ be such that $B(p,2r)$ is covered by faces of $\Gamma_\delta$. Then for each harmonic in $B(p,2r)_\delta$ function $H: \bar{B}(p,2r)_\delta \to \mathbb{R}$ we have
the following estimate of its Dirichlet energy on $B(p, r)$ (see Definition 4.1):

$$ \mathcal{E}_{B(p,r)}^{\delta}(H) \leq \frac{C}{r^2} \sum_{v_1 v_2 \in E(B(p, 2r)_{\delta})} c_{v_1 v_2} |v_1 - v_2|^2 (H(v_1) + H(v_2))^2. $$

*Proof.* Let $\phi: \mathbb{C} \to [0, 1]$ be a smooth function equal to 1 over $B(p, r)$ and vanishing near the boundary and outside of $B(p, 2r)$. Since $H$ is harmonic in $B(p, 2r)$, the discrete integration by parts gives

$$ \sum_{v_1 v_2 \in E(\Gamma_{\delta})} c_{v_1 v_2} (H(v_1) - H(v_2)) (\phi(v_1)^2 H(v_1) - \phi(v_2)^2 H(v_2)) = 0. $$

Simple algebraic manipulations allow to rewrite this identity as

$$ \sum_{v_1 v_2 \in E(\Gamma_{\delta})} c_{v_1 v_2} (H(v_1) - H(v_2))^2 (\phi(v_1)^2 + \phi(v_2)^2) = \sum_{v_1 v_2 \in E(\Gamma_{\delta})} c_{v_1 v_2} (H(v_1)^2 - H(v_2)^2) (\phi(v_2)^2 - \phi(v_1)^2). $$

Taking the square and applying Cauchy-Schwartz's inequality to the right-hand side we get the upper bound

$$ \sum_{v_1 v_2 \in E(\Gamma_{\delta})} c_{v_1 v_2} (H(v_1) - H(v_2))^2 (\phi(v_1) + \phi(v_2))^2 \times \sum_{v_1 v_2 \in E(\Gamma_{\delta})} c_{v_1 v_2} (H(v_1) + H(v_2))^2 (\phi(v_2) - \phi(v_1))^2. $$

Since $(\phi(v_1) + \phi(v_2))^2 \leq 2(\phi(v_1)^2 + \phi(v_2)^2)$ and $\phi = 1$ on $B(p, r)$, this gives the estimate

$$ \begin{aligned} 2\mathcal{E}_{B(p,r)}^{\delta}(H) &\leq \sum_{v_1 v_2 \in E(\Gamma_{\delta})} c_{v_1 v_2} (H(v_1) - H(v_2))^2 (\phi(v_1)^2 + \phi(v_2)^2) \\ &\leq 2 \sum_{v_1 v_2 \in E(\Gamma_{\delta})} c_{v_1 v_2} (H(v_1) + H(v_2))^2 (\phi(v_2) - \phi(v_1))^2, \end{aligned} $$

which implies that

$$ \mathcal{E}_{B(p,r)}^{\delta}(H) \leq \|\nabla\phi\|_{\infty}^{2} \sum_{v_1 v_2 \in E(B(p, 2r)_{\delta})} c_{v_1 v_2} |v_1 - v_2|^2 (H(v_1) + H(v_2))^2. $$

The proposition follows.

□

**Corollary 4.5.** Assume that $(\Gamma_{\delta})_{\delta>0}$ is a family of graphs embedded harmonically into $\mathbb{C}$ and satisfying property (LIP). Assume that $r > 0$ and $p \in \mathbb{C}$ are such that the disc $B(p, 3r)$ is covered by the faces of $\Gamma_{\delta}$ for each sufficiently small $\delta > 0$. There exist $C, \delta_0 > 0$ such that for each $\delta \leq \delta_0$ and each harmonic on $B(p, 2r)_{\delta}$ function $H: \bar{B}(p, 2r)_{\delta} \to \mathbb{R}$ we have

$$ \mathcal{E}_{B(p,r)}^{\delta} \leq C \max_{v \in B(p, 2r)_{\delta}} H(v)^2. $$

*Proof.* This follows from Lemma 2.11 and Proposition 4.4.

□

**4.3. Proofs of Theorem 2 and Theorem 3 on the assumption $\varphi \in C^3$.** We begin with the proof of Theorem 2 and then discuss the only additional ingredient – a priori uniform bound on the values of the Green function – that is needed to prove Theorem 3 along the same lines.

*Proof of Theorem 2 assuming that $\varphi \in C^3$.* Note that the functions $H_{\delta}$ are uniformly bounded by maximum principle. Corollary 4.5 implies that for each compact $K \subset \Omega$ there is a constant $C = C(K)$ such that

$$ (4.4) \qquad \mathcal{E}_{K}^{\delta}(H_{\delta}) \leq C(K). $$

Lemma 4.3 implies that harmonic conjugate functions $H_{\delta}^*$ of $H_{\delta}$ can be chosen to be uniformly bounded on compacts of $\Omega$. Combining property (LIP), Lemma 2.6 and Proposition 2.7 we conclude that the family $(H_{\delta}, H_{\delta}^*)_{\delta}$ is precompact in the topology of uniform convergence on every compact in $\Omega$. Let
$(h, h^*)$ be a subsequential limit of this family. Denote $f = h+ih^*$ and let $z = \frac{1}{2}(w+\psi)$, $\vartheta = \frac{1}{2}(\bar{\psi}-\bar{w})$. By Lemma 2.14 we know that the 1-form

$$f dz + \bar{f} d\bar{\vartheta} = h d\psi + ih^* dw$$

is closed ‘in the weak sense’: we have $\int_{\gamma} (f dz + \bar{f} d\bar{\vartheta}) = 0$ for every piece-wise smooth contractible loop $\gamma \subset \Omega$. Applying Lemma 3.1 and Lemma 3.2 (here we use that $\varphi \in C^3$) we conclude that $h \in C^2(\Omega)$, we have $\mathcal{L}_{\varphi}h = 0$ and $h^*$ is an A-harmonic conjugate of $h$.

Using Lemma 2.8 it is easy to see that $h$ is continuous up to $\bar{\Omega}$, we have $h|_{\partial\Omega} = g$ and $H_\delta$ converge to $h$ uniformly in $\bar{\Omega}$ along the corresponding subsequence. It follows that $h$ solves the Dirichlet problem for $\mathcal{L}_{\varphi}$ with boundary conditions $g$, thus $h$ satisfies the assumptions of the theorem. By the uniqueness of the solution of the Dirichlet problem we conclude that $H_\delta$ converge to $h$. $\square$

Among the arguments used in the proof if Theorem 2 given above, the only missing ingredient in the setup of Theorem 3 is an a priori uniform bound on the values of discrete harmonic functions under consideration, from which the proof starts. The next lemma provides this missing ingredient.

**Lemma 4.6.** In the setup of Theorem 3 for every compact $K \subset \bar{\Omega} \times \bar{\Omega} \setminus \text{diag}$ there exists a constant $C > 0$ that depends on $K$ and the constants in property (LIP) only such that for each $\delta$ small enough we have

$$\max_{(v_1,v_2) \in K} G_{\Omega}^{\delta}(v_1, v_2) \leq C.$$

*Proof.* We can assume without loss of generality that $\text{dist}(v_2, \partial\Omega) \geq C\delta$ for some $C > 0$ big enough: the desired estimate in the other case follows from the monotonicity with respect to the domain. Below we use the symbol $O(\cdot)$ to denote an estimate with a constant that depends on the constants from property (LIP) only.

Let $X_t$ be the continuous time random walk on $\Gamma_\delta$ with jump rates from $v$ to $v_1$ proportional to $c_{vv_1}$ and parameterized such that $|X_t|^2 - t$ is a local martingale. Recall that (see (1.6) and the discussion above it)

$$\mu_{\delta}(v) = \sum_{v_1 \sim v} c_{vv_1} |v - v_1|^2$$

is an invariant (in the bulk) measure for $X_t$. By Dynkin's formula we have

$$(4.5) \qquad G_{\Omega}^{\delta}(v_1, v_2) = \frac{1}{\mu_{\delta}(v_2)} \mathbb{E} \int_{0}^{\tau_{\Omega}^{\delta}(v_1)} \mathbb{1}[X_{t}^{v_1} = v_2] dt,$$

where $X_t^v$ stands for the random walk started at $v$ and $\tau_{\Omega}^{\delta}(v)$ is the first time when $X_t^v$ exits $\Omega$.

Put $r = \frac{1}{4}\text{dist}(v_2, \partial\Omega)$ and assume that $|v_1 - v_2| \geq 2r$. Using that $|X_t|^2 - t$ is a local martingale and crossing estimates from Lemma 2.8 it is easy to see that

$$\mathbb{E} \int_{0}^{\tau_{\Omega}(v_1)} \mathbb{1}[X_{t}^{v_1} \in B(v_2, r)] = O(r^2).$$

Due to (4.5), this gives

$$\sum_{v \in B(v_2, r)_{\delta}} \mu_{\delta}(v) G_{\Omega}(v_1, v) = O(r^2).$$

By Harnack inequality (Lemma 2.9) there exist $c > 0$ such that for each $v \in B(v_2, r)$ we have

$$G_{\Omega}^{\delta}(v_1, v) \geq c \cdot G_{\Omega}^{\delta}(v_1, v_2).$$

From Lemma 2.11 we also know that $\sum_{v \in B(p, r)_{\delta}} \mu_{\delta}(v) \geq cr^2$ for some $c > 0$. Thus, $G_{\Omega}^{\delta}(v_1, v_2) = O(1)$. We conclude that $G_{\Omega}^{\delta}(v_1, v_2)$ is uniformly bounded when $|v_1 - v_2| \geq \max(C\delta, \frac{1}{2}\text{dist}(v_2, \partial\Omega))$. This implies the estimate over an arbitrary compact $K \subset \bar{\Omega} \times \bar{\Omega} \setminus \text{diag}$ by Harnack inequality. $\square$
*Proof of Theorem 3 assuming* $\varphi \in C^3$. Using Lemma 4.6 and Lemma 2.10, as well as Lemma 2.8 to estimate $G_{\Omega}^{\delta}$ near the boundary, we conclude that the family $(G_{\Omega}^{\delta})_{\delta>0}$ is precompact in the topology of uniform convergence on compacts in $\overline{\Omega} \times \overline{\Omega} \setminus \text{diag}$. It remains to prove that every subsequential limit coincides with the Green's function $G$ of $\mathcal{L}$. To this end we employ the characterization of $G$ discussed before Lemma 3.3.

Let $\hat{G}$ be such a subsequential limit. Clearly, $\hat{G}$ is non-negative, continuous in $\overline{\Omega} \times \overline{\Omega} \setminus \text{diag}$ and vanishes along the boundary of $\Omega$. Now fix $v_2 \in \Omega$ and put $G_{\delta}(v) = G_{\Omega}^{\delta}(v, v_2)$ for simplicity. Corollary 4.5 and Lemma 4.3 allow us to define harmonic conjugates $G_{\delta}^{*}$ such that they are uniformly bounded on compact subsets of the universal cover of $\Omega \setminus \{v_2\}$. Arguing as in the proof of Theorem 2 above we conclude that $\hat{G}(\cdot, w_0) \in C^2(\Omega \setminus \{w_0\})$, we have $\mathcal{L}_{\varphi}\hat{G}(\cdot, w_0) = 0$ in $\Omega \setminus \{w_0\}$ and $G_{\delta}^{*}(\cdot, w_0)$ converge to an A-harmonic conjugate of $\hat{G}(\cdot, w_0)$ uniformly on compacts of the universal cover of $\Omega \setminus \{w_0\}$. In particular, we conclude that A-harmonic conjugates of $\hat{G}(\cdot, w_0)$ have monodromy $-1$ around $w_0$. Due to Lemma 3.3, this characterizes $\hat{G}$ as the Green's function $G$. $\square$

**4.4.** $C^1$ convergence under additional regularity assumption (EXP-FAT). The main convergence results of our paper, Theorem 2 and Theorem 3 claim the convergence of discrete harmonic functions $H_{\delta}$ to solutions of the equation $\mathcal{L}_{\varphi}h = 0$ under no local 'non-degeneracy' assumptions on the embeddings $\Gamma_{\delta}$. At the same time, in applications of such results to the dimer model one often needs to prove convergence of the *gradients* of $H_{\delta}$. (For instance, see [BLQ23, Section 1.7] for a short discussion of the link between dimer model observables in this setup and discrete Green's functions.) In this section we strengthen Theorem 2 and Theorem 3 to the convergence of gradients of $H_{\delta}$ under a mild additional regularity assumption (EXP-FAT) on $\Gamma_{\delta}$ that originated in [CLR23]. We say that a face of $\Gamma_{\delta}$ is $\rho$-fat if it can be triangulated so that each triangle contains a disc of radius $\rho$.

**Assumption (EXP-FAT) (cf. [CLR23, Assumption 5.9]).** Assume that a family of harmonic embeddings $(\Gamma_{\delta})_{\delta \to 0}$ and an open set $U \subset \mathbb{C}$ covered by each of $\Gamma_{\delta}$ be given. We say that $\Gamma_{\delta}$ have property ExpFat in $U$ if there exists a function $\delta'(\delta)$ such that $\delta'(\delta) \to 0$ as $\delta \to 0$ and the following holds:

if one removes from $\Gamma_{\delta}$ all $\delta \exp(-\delta'(\delta)\delta^{-1})$-fat faces, then all vertex-connected components of $\Gamma_{\delta}$ that are fully contained in $U$ have diameters less than $\delta'(\delta)$.

Let us begin with a short motivating discussion. Assume that the potential $\varphi$ is smooth enough and that a real-valued function $h$ is a solution of the equation $\mathcal{L}_{\varphi}h = 0$. A direct computation shows that this is equivalent to the fact that the 1-form $\text{Im}[h_w d\varphi]$ is closed. Note also that the 1-form $\text{Re}[h_w dw] = \frac{1}{2}dh$ is always closed. These real-valued 1-forms are the imaginary and the real parts of the complex-valued 1-form

$$ (4.6) \qquad h_w dz - h_{\bar{w}} d\theta, \quad z = \frac{1}{2}(w + \psi), \quad \theta = \frac{1}{2}(\bar{\psi} - \bar{w}), $$

which is therefore also closed. Let $\alpha \in \mathbb{T}$ and note that

$$ h_w dz - h_{\bar{w}} d\theta = \bar{\alpha} (\text{Re}[\alpha h_w] d\psi_{\alpha} + i \text{Im}[\alpha h_w] dw_{\alpha}), \quad w_{\alpha} = z + \alpha^2 \theta, \quad \psi_{\alpha} = z - \alpha^2 \theta. $$

Repeating the proof of Lemma 3.1 one can see from here that each projection $\text{Re}[\alpha h_w]$ of the gradient $h_w$ of $h$ satisfies an elliptic PDE (which depends on $\alpha$) in the coordinate $w_{\alpha}$.

The framework developed in [CLR23] provides a 'discrete version' of the computations made above when a smooth potential $\varphi$ is replaced with a piecewise linear potential $\Phi_{\delta}$ and $h$ with a discrete harmonic function $H_{\delta}$ on $\Gamma_{\delta}$. Splitting each face of $\Gamma_{\delta}$ into triangles and extending $H_{\delta}$ inside these triangles linearly we can view $H_{\delta}$ as a piecewise linear function and its gradient $F_{\delta} := \partial_{w}H_{\delta}$ as a piecewise constant function. (Note the advantage of the discrete setup: $F_{\delta}$ is well-defined pointwise for each discrete potential $\Phi_{\delta}$.) The discrete analog of the 1-form (4.6) is then the 1-form $F_{\delta}d\mathcal{T}_{\delta} - \bar{F}_{\delta}d\mathcal{O}_{\delta}$ on the t-surface $\Theta_{\delta}$ associated with $\Gamma_{\delta}$. (See [CLR23, Section 4.2 and Section 5] for an extension of this 1-form to the whole t-surface from black faces.) As above, this gives a certain discrete harmonicity property for each of the functions $\text{Re}[\alpha F_{\delta}], \alpha \in \mathbb{T}$, in the coordinate $\mathcal{T}_{\delta}+\alpha^2\mathcal{O}_{\delta}$; see [CLR23, Section 4.3]
for details. Importantly, if $\Gamma_\delta$ has property (CONV)/(LIP), with no additional assumptions on $\Gamma_\delta$, then the corresponding random walks on T-graphs $\mathcal{T}_\delta + \alpha^2\mathcal{O}_\delta$ are uniformly elliptic starting from the scale $\delta$. This leads to the following:

**Proposition 4.7.** Let a harmonic embedding $\Gamma_\delta$ cover an open set $U \subset \mathbb{C}$, and $H_\delta: \bar{U}_\delta \to \mathbb{R}$ be a discrete harmonic function. Triangulate faces of $\Gamma_\delta$ arbitrarily and extend $H_\delta$ to $U$ by linearity.

(i) The gradient $\nabla H_\delta$ of $H_\delta$ satisfies the maximum principle: $\max_U |\nabla H_\delta| = \max_{\partial U} |\nabla H_\delta|$, where we take the maximum of $|\nabla H_\delta|$ over all incident faces if $\partial U$ passes through a vertex of $\Gamma_\delta$.

(ii) Assume in addition that $\Gamma_\delta$ has property (LIP). There exists constants $C, \beta > 0$ (depending on constants in (LIP) only) such that for each $R \ge r \ge C\delta$ and $w \in U$ such that $B(w, R) \subset U$ we have

$$ \mathrm{osc}_{B(w,r)} \nabla H_{\delta} \le C(r/R)^{\beta} \mathrm{osc}_{B(w,r)} \nabla H_{\delta}, $$

where $\mathrm{osc}_B F = \sup_{w_1,w_2 \in B} |F(w_1) - F(w_2)|$. The exponent $\beta$ depends only on the constant $\lambda$ in (1.4).

(iii) In the same setup, there exist constants $\beta_0, C_0 > 0$ depending on the constants in property (LIP) only such that the following holds. Let $r \ge C_0\delta$ and $w \in U$ be such that $B(w, r) \subset U$. Then,

$$ \text{either } \max_{B(v, \frac{1}{2}r)} |\nabla H_{\delta}| \le C_0 r^{-1} \mathrm{osc}_{B(v,r)} H_{\delta}, $$

$$ \text{or } \max_{B(v, \frac{3}{4}r)} |\nabla H_{\delta}| \ge \exp(\beta_0 r \delta^{-1}) C_0 r^{-1} \mathrm{osc}_{B(w,r)} H_{\delta}. $$

*Proof.* Note that the gradient of $H_\delta$ can be seen as a t-black-holomorphic function on the t-embedding $\mathcal{T}_\delta = \frac{1}{2}(\mathcal{T}_\delta + \mathcal{O}_\delta)$; see the discussion at the end of Section 2.4 and [CLR23, Proposition 8.2(ii)]. Recall also that the mapping from $\Gamma_\delta$ to $\mathcal{T}_\delta$ is bi-Lipschitz on convex sets starting from the scale $\delta$. Item (i) follows from [CLR23, Proposition 4.17 and Proposition 5.7]. Items (ii) and (iii) follow from [CLR23, Proposition 6.13 and Theorem 6.17]. □

Given the alternative provided by item (iii) in Proposition 4.7, we can now benefit from the fact that the gradient of a discrete harmonic function $H_\delta$ on a $\rho$-fat face of $\Gamma_\delta$ cannot be bigger than $\rho^{-1} \max |H_\delta|$. Together with the maximal principle for $\nabla H_\delta$ from item (i), this allows to control $\nabla H_\delta$ if $\Gamma_\delta$ does not have clusters of not $\exp(-o_\delta \to 0(1)\delta^{-1})$-fat faces that percolate to the boundary.

**Theorem 4.8.** In the setup of Theorem 2, assume additionally that harmonic embeddings $\Gamma_\delta$ satisfy Assumption (EXP-FAT), where $\delta$ is taken from property (CONV). There exists a constant $\beta > 0$ that depends only on the constant $\lambda$ in (1.4) such that the following holds: if we extend discrete harmonic functions $H_\delta$ linearly to $\Omega \in U$ using the triangulations of faces from Assumption (EXP-FAT), then $h \in C^{1,\beta}(\Omega)$ and the gradients $\nabla H_\delta$ converge to $\nabla h$ uniformly on compact subsets of $\Omega$. Moreover, for each compact $K \subset \Omega$ there exists a constant $C = C(K, \lambda) > 0$ such that $\|h\|_{C^{1,\beta}(K)} \le C\|h\|_{C(\bar{\Omega})}$.

Similarly, if $\Gamma_\delta$ satisfy Assumption (EXP-FAT) in the setup of Theorem 3, then the gradients of $G_\delta^\Omega$ converge to the gradients of $G_\Omega$ uniformly on compact subsets of $(\Omega \times \Omega) \setminus \mathrm{diag}$.

**Remark 4.4.** Let us emphasize that in Theorem 4.8 we do not impose any additional regularity assumption on the potential $\varphi \in C^{1,1}$ besides the uniform convexity (1.4). In particular, one can use Theorem 4.8 to prove that $h \in C_{\mathrm{loc}}^{1,\beta}$ for each (weak) solution $h \in W_{\mathrm{loc}}^{1,2}$ of the equation $\mathcal{L}_\varphi h = 0$ with uniformly convex potential $\varphi$ by discretizing the operator $\mathcal{L}_\varphi$; see [GN11] for related results.

*Proof.* Consider the setup of Theorem 2 first. The family $(H_\delta)_{\delta>0}$ of solutions of discrete Dirichlet problems in $\Omega_\delta$ is uniformly bounded in $\bar{\Omega}$ by the maximum principle. Combining Assumption (EXP-FAT) with the alternative provided by item (iii) and the maximum principle from item (i) of Proposition 4.7 we conclude that the gradients $\nabla H_\delta$ are uniformly bounded on compact subsets of $\Omega$. Moreover, the Hölder-type estimate from the item (ii) of the same proposition and a standard pre-compactness argument implies that any subsequential limit $h$ of $H_\delta$ belongs to $C^{1,\beta}(\Omega)$ and that $\nabla H_\delta$ converge to $\nabla h$ uniformly on compacts along the corresponding subsequence. Since all estimates in Proposition 4.7 are uniform in $\delta$, we also get the desired estimate $\|h\|_{C^{1,\beta}(K)} \le C\|h\|_{C(\Omega)}$.
It remains to prove that $h$ is the (weak) solution of the Dirichlet problem for $\mathcal{L}_\varphi$ even if $\varphi \notin C^3$. To this end we repeat the proof of Theorem 2 given in Section 4.3 and note that we do not need to use Lemma 3.2 anymore since under Assumption (EXP-FAT) we already know that each subsequential limit $h$ of functions $H_\delta$ is continuously differentiable inside $\Omega$.

Similar arguments apply in the setup of Theorem 3: discrete Green's functions $G_\delta^\Omega$ are uniformly bounded on compact subsets of $(\Omega \times \Omega) \setminus \text{diag}$ due to Lemma 4.6 which relies upon the property (CONV) only and does not require any additional regularity of $\varphi$. As above, this implies that the gradients of discrete Green's functions are uniformly bounded on compacts. (Note that in Assumption (EXP-FAT) we insist that each point of $\Omega$ is surrounded by a small contour consisting of $\delta \exp(-\delta'\delta^{-1})$-fat faces.) The proof finishes as before, with no need in Lemma 3.2 as we now know that subsequential limits of $G_\delta^\Omega$ are continuously differentiable. $\square$

*Remark 4.5.* Note that in the setup of Theorem 2 we only used the fact that each compact set $K \subset \Omega$ is separated from $\partial\Omega$ by a cycle of $\delta \exp(-\delta'(\delta)\delta^{-1})$-fat faces of $\Gamma_\delta$ for all sufficiently small $\delta \le \delta_0(K)$. This weaker form of Assumption (EXP-FAT) is not enough in the setup of Theorem 3 as we also need the existence of a small circuit of $\delta \exp(-\delta'(\delta)\delta^{-1})$-fat faces separating $K$ from the pole of $G_\delta^\Omega$.

**4.5. Proof of Theorem 2 and Theorem 3 for general uniformly convex potentials** $\varphi \in C^{1,1}$. Recall that the only place in the proofs of Theorem 2 and Theorem 3 in Section 4.3 where we used the additional regularity assumption $\varphi \in C^3$ is Lemma 3.2, which is needed to prove that subsequential limits of discrete harmonic functions on $\Gamma_\delta$ are differentiable. Generalizing Lemma 3.2 to general uniformly convex potentials $\varphi \in C^{1,1}$ appears to be out of our reach. In this section we provide an alternative proof of the differentiability of all such limits – see Proposition 4.9 below – which can be used instead of Lemma 3.2 in the general case.

Since the statement is local, we can assume that we are given a sequence of uniformly bounded discrete harmonic functions in a disc $B(w_0, R) \subset U$ such that $H_\delta \to h \in C(B(w_0, R))$ as $\delta \to 0$, uniformly on compacts. Our goal is to prove that $h$ is continuously differentiable and, moreover, that

$$ (4.7) \qquad \|h\|_{C^{1,\beta}(B(w_0,r))} \le C(\lambda, r/R) \cdot \|h\|_{C(B(w_0,R))} \quad \text{for all } r < R, $$

where the exponent $\beta > 0$ depends only on the constant $\lambda$ in (1.4). Put $\tilde{r} = \frac{1}{2}(r+R)$. The idea of the proof is to modify $\Gamma_\delta$ outside the disc $B(w_0, \tilde{r})$ and to construct approximations $\tilde{H}_\delta^\varepsilon$ of given functions $H_\delta$ such that the gradients $\nabla \tilde{H}_\delta^\varepsilon$ are uniformly bounded in $B(w_0, \tilde{r})$ due to the alternative from item (iii) in Proposition 4.7 applied on the *modified graph* $\tilde{\Gamma}_\delta$; see also Remark 4.5 above.

Let

$$ \tilde{V}_\delta = \left\{ v_\delta \in \mathbb{C} : \begin{array}{l} \text{either } v_\delta \in \Gamma_\delta \text{ and } |v_\delta - w_0| < \tilde{r} \\ \text{or } v_\delta \in \delta\mathbb{Z}^2 \text{ and } \tilde{r} \le |v_\delta - w_0| \le R \end{array} \right\} $$

and denote by $\tilde{P}_\delta$ the upper convex hull of the set $\{(v_\delta; \Phi_\delta(v_\delta)), v_\delta \in V_\delta\} \subset \mathbb{C} \times \mathbb{R}$. Since $\Phi_\delta$ is convex, $\tilde{P}_\delta$ is the supergraph of a convex function $\tilde{\Phi}_\delta$ that is defined on the convex hull of $V_\delta$ and is linear on faces of a certain graph whose vertices belong to $V_\delta$. We denote this graph by $\tilde{\Gamma}_\delta$ and define the weight of an edge of $\Gamma_\delta$ as the jump of the gradient of $\tilde{\Phi}_\delta$ between two adjacent faces divided by the length of this edge. (In other words, the weighted graph $\tilde{\Gamma}_\delta$ and the modified potential $\tilde{\Phi}_\delta$ are related to each other in the usual way; e.g., see (1.1).) It follows from (1.4) that for each small enough $\varepsilon > 0$ there exists $\delta_0(\varepsilon) > 0$ such that for all $\delta \le \delta_0(\varepsilon)$

* all edges of $\tilde{\Gamma}_\delta$ in the disc $B(w_0, R - \varepsilon)$ have length at most $C\delta$,

* all faces of $\tilde{\Gamma}_\delta$ in the annulus $B(w_0, R - \varepsilon) \setminus B(w_0, \tilde{r} + \varepsilon)$ are $c\delta$-fat,

* $\tilde{\Phi}_\delta$ equals $\Phi_\delta$ (and hence $\tilde{\Gamma}_\delta$ coincides with $\Gamma_\delta$ as weighted graphs) in the disc $B(w_0, \tilde{r} - \varepsilon)$,

where the constants $C, c > 0$ depend only on the constants in property (CONV). We are now in the position to establish the smoothness of all subsequential limits of discrete harmonic functions $H_\delta$.
**Proposition 4.9.** There exist an exponent $\beta > 0$ depending only on the constant $\lambda$ in (1.4) such that the estimate (4.7) holds for each (subsequential, uniform on compacts) limit $h \in C(B(w_0, r))$ of discrete harmonic functions $H_\delta$ on $\Gamma_\delta$.

*Proof.* Without loss of generality, assume that $\|h\|_{C(B(w_0, R))} = 1$ and recall that the functions $H_\delta$ satisfy a Hölder-type regularity estimate from Lemma 2.10, which is uniform in $\delta$. Passing to the limit $\delta \to 0$, this gives the estimate $\omega_h(\rho) \le C\rho^\beta/(R-\tilde{R})^\beta$ for the modulus of continuity of the function $h$ in each smaller disc $B(w_0, \tilde{R})$, where $C, \beta > 0$ depend on constants in (1.4) only.

Let $H_\delta^\varepsilon$ be the solution of the discrete Dirichlet problem on $\Gamma_\delta$ in the disc $B(w_0, \tilde{r} + 2\varepsilon)$ with boundary values $h$ and $\tilde{H}_\delta^\varepsilon$ be the solution of the same boundary value problem on the modified graph $\tilde{\Gamma}_\delta$. Recall that these two graphs coincide in $B(w_0, \tilde{r} - \varepsilon)$ if $\delta \le \delta_0(\varepsilon)$. We claim that

$$ (4.8) \qquad \max_{B(w_0, \tilde{r}-\varepsilon)} |\tilde{H}_\delta^\varepsilon - H_\delta^\varepsilon| \le C(\varepsilon/\sqrt{\varepsilon})^\beta + \omega_h(\sqrt{\varepsilon}) \quad \text{for all } \delta \le \delta_0(\varepsilon). $$

Due to the maximum principle, it is enough to prove this estimate for vertices $v^\delta$ near $\partial B(w_0, \tilde{r} - \varepsilon)$. Each of the two values $\tilde{H}_\delta^\varepsilon(v_\delta)$ and $H_\delta^\varepsilon(v_\delta)$ can be written as the expectation of the boundary value at the exit point of the corresponding random walk from $B(w_0, \tilde{r} + 2\varepsilon)$. Since both random walks are uniformly elliptic (see Lemma 2.8), the probability that at least one of these hitting points does not belong to a $\sqrt{\varepsilon}$-vicinity of $v_\delta$ is bounded by $O((\varepsilon/\sqrt{\varepsilon})^\beta)$, which implies the desired upper bound.

Since $H_\delta^\varepsilon \to h$ as $\delta \to 0$ uniformly on compacts, it follows from the maximum principle that $H_\delta^\varepsilon \to h$ in $B(w_0, \tilde{r})$, uniformly in $\varepsilon$. Choose $\varepsilon = \varepsilon(\delta) \to 0$ so that $\delta \le \delta_0(\varepsilon)$. It follows from (4.8) that

$$ \tilde{H}_\delta^{\varepsilon(\delta)} \to h \quad \text{as} \quad \delta \to 0 \text{ uniformly on compact subsets of } B(w_0, \tilde{r}). $$

By construction, the modified graph $\tilde{\Gamma}_\delta$ contains a circuit of $c\delta$-fat faces near the boundary of the disc $B(w_0, \tilde{r} + 2\varepsilon)$. Using Proposition 4.7 similarly to the proof of Theorem 4.8, we conclude that the gradients $\nabla \tilde{H}_\delta^\varepsilon$ are uniformly bounded in the disc $B(w_0, \tilde{r})$. Passing to the limit $\delta \to 0$ in the Hölder-type estimate for the gradients $\nabla \tilde{H}_\delta^\varepsilon$ provided by item (ii) in Proposition 4.7 gives (4.7). $\square$

*Proof of Theorem 2 and Theorem 3 for $\varphi \in C^{1,1}$. As already mentioned in the beginning of this section, to prove Theorem 2 and Theorem 3 in the general case one simply repeats the proof given in Section 4.3 with Lemma 3.2 replaced by Proposition 4.9.* $\square$

**Remark 4.6.** The approximation trick used above does not allow us to conclude that the gradients $\nabla H_\delta$ are always uniformly bounded. However, it allows one to prove that $H_\delta$ are ‘differentiable at a mesoscopic scale’; the result recently obtained in [Pec25] for harmonic functions on orthodiagonal tilings. To this end, note that in the arguments given above one can set $\tilde{r} = R - \delta^\alpha$, $\alpha < 1$, and $\varepsilon = C_0\delta$ with a big enough constant $C_0$ depending on constants in property (CONV) only. If one replaces $\sqrt{\varepsilon}$ in the proof of (4.8) with $\delta^{\frac{1}{2}(1+\alpha)}$, the right-hand side of this estimate reads as $O(\delta^{\frac{1}{2}\beta(1-\alpha)})$. Thus,

$$ \begin{aligned} H_\delta(w_2) - H_\delta(w_1) &= H_\delta^\varepsilon(w_2) - H_\delta^\varepsilon(w_1) + O(\delta^{\frac{1}{2}\beta(1-\alpha)}) \\ &= \text{Re} \left[ \partial_w H_\delta^\varepsilon(w_1)(w_2 - w_1) \right] + O(|w_2 - w_1|^{1+\beta}) \end{aligned} $$

if $|w_2 - w_1| \ge \delta^{\frac{1}{2}\beta(1-\alpha)/(1+\beta)}$ and $w_1, w_2 \in B(w_0, r)$, where the constant in $O(\cdot)$ depends on $r/R$.

We do not know whether or not the gradients of uniformly bounded discrete harmonic functions are always uniformly bounded if no additional regularity assumption like (EXP-FAT) is imposed on $\Gamma_\delta$ and leave this as an open question.

## 5. Equivalence of properties (CONV)/(LIP) and (RW)

We begin this Section with an (elementary) proof of the equivalence of properties (CONV) and (LIP). Recall that the property (LIP) is also equivalent to a similar Lipschitz property Lip($\kappa, \delta$) of the t-surfaces $\Theta_\delta$; see Definition 2.5 and Lemma 2.6.
The proof of the property (RW) on t-surfaces satisfying Lip(κ, δ) is contained in Proposition 2.7 and Lemma 2.11. The main content of this section is the converse statement: the fact that (RW) implies property Lip(κ, δ) of the underlying t-surfaces.

**Lemma 5.1.** (i) Assume that the potential $\Phi_{\delta}$ satisfies property (CONV) with constant $\lambda > 0$ on an open set $U \subset \mathbb{C}$. Then its gradient $\Psi_{\delta} = 2\partial_{w}\Phi_{\delta}$ satisfies property (LIP) with constant $\varkappa = \frac{1}{8}\lambda$ on each open set $U' \subset U$ such that $\text{dist}(U', \partial U) \ge 4C\delta$, where $C$ is the constant from (1.4).

(ii) Vice versa, if the map $\Psi_{\delta} = 2\partial_{w}\Phi_{\delta}$ satisfies property (LIP) with constant $\varkappa > 0$ on a set $U \subset \mathbb{C}$, then the potential $\Phi_{\delta}$ satisfies property (CONV) with constant $\lambda = \frac{1}{4}\varkappa$ on the same set $U$.

*Proof.* We start with item (ii). By linearity, it is enough to check (1.4) for segments $[w_1; w_2] \subset U'$ of length $4r = |w_2 - w_1| \in [2C\delta, 4C\delta]$, where $C$ is the constant from (1.5). Rotating and translating the graph $\Gamma_{\delta}$, we can assume that $w_1 = -2r$ and $w_2 = 2r$. We have

$$ (5.1) \qquad \Phi_{\delta}(2r) - 2\Phi_{\delta}(0) + \Phi_{\delta}(-2r) = \int_{-r}^{r} \mathrm{Re}[\Psi_{\delta}(t+r) - \Psi_{\delta}(t-r)]dt. $$

Due to (1.5), we have $2\varkappa r \le \mathrm{Re}[\Psi_{\delta}(t+r) - \Psi_{\delta}(t-r)] \le 2\varkappa^{-1}r$ for all $t$. This gives (1.4) with $\lambda = \frac{1}{4}\varkappa$.

Let us now pass to item (i). Again, it is sufficient to check (1.5) for segments $[w_1; w_2] \subset U'$ of length $2r = |w_2 - w_1| \in [C\delta, 2C\delta]$, where $C$ is the constant from (1.4). Rotating and translating the graph $\Gamma_{\delta}$, we can assume that $w_1 = -r$, $w_2 = r$ and, moreover, that the set $U$ contains the twice bigger disc $\bar{B}(0, 2r)$. Since $\Psi_{\delta}$ is convex, the function $\mathrm{Re}\,\Phi_{\delta}(x)$ is an increasing function of $x \in \mathbb{R}$ and the identity (5.1) gives

$$ \lambda r^2 \le \Phi_{\delta}(r) - 2\Phi_{\delta}(0) + \Phi_{\delta}(-r) \le r \cdot \mathrm{Re}[\Psi_{\delta}(r) - \Psi_{\delta}(-r)] \le \Phi_{\delta}(2r) - 2\Phi_{\delta}(0) + \Phi_{\delta}(-2r) \le 16\lambda^{-1}r^2. $$

To conclude the proof of (1.4), it remains to prove a similar upper bound for $|\mathrm{Im}[\Psi_{\delta}(r) - \Psi_{\delta}(-r)]|$. To this end, note that the convexity of $\Phi_{\delta}$ implies that

$$ r^{-1}(\Phi_{\delta}(2ir) - \Phi_{\delta}(r)) \ge -\mathrm{Re}\,\Psi_{\delta}(r) + 2\mathrm{Im}\,\Psi_{\delta}(r) \quad \text{and} \quad r^{-1}(\Phi_{\delta}(-2ir) - \Phi_{\delta}(-r)) \ge \mathrm{Re}\,\Psi_{\delta}(r) - 2\mathrm{Im}\,\Psi_{\delta}(r). $$

Therefore,

$$ 2 \mathrm{Im}[\Psi_{\delta}(r) - \Psi_{\delta}(-r)] \le \mathrm{Re}[\Psi_{\delta}(r) - \Psi_{\delta}(-r)] + r^{-1}(\Phi_{\delta}(2ir) - \Phi_{\delta}(0) + \Phi_{\delta}(-2ir)) \le 32\lambda^{-1}r. $$

The upper bound for $2 \mathrm{Im}[\Psi_{\delta}(-r) - \Psi_{\delta}(r)]$ is similar and we obtain (1.5) with $\varkappa = \frac{1}{8}\lambda$. $\square$

We now pass to proving that property (RW) of random walks on an open set $U$ implies property Lip(κ, δ) of the corresponding t-surfaces on each subset $U' \in U$. In fact, it is sufficient to assume that $\text{dist}(U', \partial U) \ge C\delta$ with a sufficiently large constant $C > 0$ that depends on constants in property (RW) only. We start with the following observation:

**Lemma 5.2.** The assertion of Lemma 4.3 is still valid if instead of property (LIP) we assume that the graphs $\Gamma_{\delta}$ satisfy the upper bound on the invariant measure $\mu_{\delta}$ from property (RW).

*Proof.* By examining the proof of Lemma 4.3 one can see that the only way in which the property (LIP) is used there is the estimate of the Dirichlet energy of a smooth function using the upper bound from Lemma 2.11. The latter upper bound is nothing but the estimate of the invariant measure $\mu_{\delta}$ from property (RW). Thus, we can replace Lemma 2.11 by property (RW) keeping the proof unchanged. $\square$

An immediate corollary of this observation is

**Lemma 5.3.** Assume that $\Gamma_{\delta}$ have property (RW) on an open set $U \subset \mathbb{C}$. There exists a constant $C > 0$ depending on the constants in (RW) only such that one has

$$ |\Psi_{\delta}(w_1) - \Psi_{\delta}(w_2)| \le C|w_1 - w_2| \quad \text{if} \quad |w_1 - w_2| \ge C\delta $$

for each segment $[w_1; w_2]$ lying in the $C\delta$-interior of $U$.
*Proof.* Due to linearity, it is sufficient to consider the case $C_0\delta \le |w_1 - w_2| \le 2C_0\delta$ for a large constant $C_0 > 0$. Let $B = B(w_1, 2|w_1 - w_2|)$ and $\mathcal{H}_\delta : \Gamma_\delta \to \mathbb{C}$ be the map that sends each vertex to its position in the complex plane. Provided that $C_0$ is big enough, it follows from property (RW) that

$$ \varepsilon_U^\delta(\mathcal{H}_\delta) \le \frac{1}{2} \sum_{v \in B_\delta} \mu_\delta(v) \le 2c^{-1}|w_1 - w_2|^2. $$

Applying Lemma 4.3 and Lemma 5.2 to the real and the imaginary parts of $\mathcal{H}_\delta$ we conclude that $|\Psi_\delta(w_1) - \Psi_\delta(w_2)| \le C|w_1 - w_2|$, where $C$ depends only on constants in the property (RW). $\square$

Recall the notation $W_\delta = \mathcal{T}_\delta - \bar{\mathcal{O}}_\delta$ introduce in Section 2 and denote $W_\delta^* = \mathcal{T}_\delta + \bar{\mathcal{O}}_\delta$. It follows from Lemma 2.3) that $W_\delta$ projects $\Theta_\delta$ onto the harmonic embedding $\Gamma_\delta$ while $W_\delta^*$ projects $\Theta_\delta$ onto $\Gamma_\delta^*$ embedded via the conjugate harmonic embedding.

**Corollary 5.4.** Assume that $\Gamma_\delta$ have property (RW) on an open set $U \subset \mathbb{C}$. There exists a constant $C > 0$ depending on the constants in (RW) only such that one has

$$ |\mathcal{T}_\delta(p_1) - \mathcal{T}_\delta(p_2)| \le C|W_\delta(p_1) - W_\delta(p_2)| \quad \text{if} \quad |W_\delta(p_1) - W_\delta(p_2)| \ge C\delta. $$

for all pairs of points $p_1, p_2 \in \Theta_\delta$ such that the segment $[W_\delta(p_1); W_\delta(p_2)]$ lies in the $C\delta$-interior of $U$.

*Proof.* This directly follows from Lemma 5.3 as $\mathcal{T}_\delta(p) = W_\delta(p) + W_\delta^*(p) = W_\delta(p) + \Psi_\delta(W_\delta(p))$. $\square$

Given an edge $v_1v_2$ of $\Gamma_\delta$ denote by $u(v_1v_2)$ the corresponding white face of the corner graph $\mathcal{V}_\delta$ (see Definition 2.1 and the discussion below it). Recall that $\mu_\delta(v_1) = 4 \sum_{v_2 \sim v_1} \text{Area}[\mathcal{T}_\delta(u(v_1v_2))]$. For $\alpha \in \mathbb{T}$, define

$$ \begin{align} \mu_\delta^\alpha(v_1) &= 2 \sum_{v_2 \sim v_1} \text{Area}[(\mathcal{T}_\delta - \alpha^2 \mathcal{O}_\delta)(u(v_1v_2))] \tag{5.2} \\ &= \sum_{v_2 \sim v_1} c_{v_1v_2} (\text{Im}[\bar{\alpha} \eta_{u(v_1v_2)}])^2 |v_1 - v_2|^2 = \sum_{v_2 \sim v_1} c_{v_1v_2} (\text{Re}[\bar{\alpha}(v_1 - v_2)])^2. \nonumber \end{align} $$

Above, $\eta$ is the origami square root function defined by (2.8), the second equality follows from (2.9) and the last equality follows from the definition of $\eta$. Note that we always have $\mu_\delta^\alpha \le \mu_\delta$. As we will see in the next proposition, the ellipticity estimate from property (RW) implies that these measures are in fact uniformly comparable to each other starting from the scale $\delta$.

**Proposition 5.5.** Assume that $\Gamma_\delta$ have property (RW) on an open set $U \subset \mathbb{C}$. There exist constants $c, C > 0$ depending on the constants in (RW) only such that for each $\alpha \in \mathbb{T}$ one has

$$ \sum_{v \in B(w,r)} \mu_\delta^\alpha(v) \ge cr^2. $$

for each $w \in U$ and $r \ge C\delta$ such that the disc $B(w,r)$ lies in the $C\delta$-interior of $U$.

*Proof.* Denote $Z = \sum_{v \in B(w,r)} \mu_\delta(v)$ and $\nu = Z^{-1}\mu_\delta \mathbb{1}_{B(w,r)}$. Recall that $Z \ge cr^2$ as we assume property (RW). Let $X_t^\nu$ be the random walk on $\Gamma_\delta$ stopped upon leaving $B(w,r)$ with $X_0$ distributed according to $\nu$. Provided that $C$ is large enough, Proposition 2.7 and a standard large deviation estimate (e.g., see [CLR23, Proposition 6.1]) imply that for each $\alpha \in \mathbb{T}$ one has

$$ (5.3) \qquad \mathrm{Var}(\mathrm{Re}[\bar{\alpha}(X_{C\delta^2}^\nu - X_0^\nu)]) \ge \delta^2 \sum_{v \in B(w,r/2)} \nu(v) \ge \frac{1}{4} c^2 \delta^2, $$

where the last estimate is a part of the property (RW).

Now let $X_t^{v_1}$ be the random walk on $\Gamma_\delta$ started from a fixed vertex $v_1 \in \Gamma_\delta$. Using the last equality in (5.2) it is easy to see that

$$ \frac{d}{dt} \mathrm{Var}(\mathrm{Re}[\bar{\alpha} X_t^{v_1}])|_{t=0} = \frac{\mu_\delta^\alpha(v_1)}{\mu_\delta(v_1)}. $$
Applying this identity to $X_t^\nu$ we obtain the following estimate:

$$ (5.4) \quad \mathrm{Var}(\mathrm{Re}[\bar{\alpha}(X_{C\delta^2}^\nu - X_0^\nu)]) \le \int_0^{C\delta^2} \sum_{v \in B(w,r)} \frac{\mu_\delta^\alpha(v)}{\mu_\delta(v)} \nu(v) dt = C\delta^2 Z^{-1} \sum_{v \in B(w,r)} \mu_\delta^\alpha(v), $$

where the first inequality follows from the fact that the measure $\nu$ is sub-stationary (as $\mu$ is stationary). Combining (5.3), (5.4), and the estimate $Z \ge cr^2$ one sees that $\sum_{v \in B(w,r)} \mu_\delta^\alpha(v) \ge \frac{1}{4}c^2C^{-1}r^2$. $\square$

We are now in the position to prove the remaining part of Theorem 1.

**Theorem 5.6.** Assume that $\Gamma_\delta$ have property (RW) on an open set $U$ and $U' \subset U$. Then, for all sufficiently small $\delta$ the $t$-surfaces $\Theta_\delta$ associated with $\Gamma_\delta$ have property Lip($\kappa$, $\delta$) on $U'$, where $k < 1$ depends only on constants in (RW). Moreover, it is sufficient to assume that $U'$ lies in the $C\delta$-interior of $U$ with sufficiently large constant $C > 0$.

*Proof.* Let $R_0 > 0$ be a sufficiently large constant that will be chosen at the end of the proof and consider a point $p \in \Theta_\delta$ such that $B(W_\delta(p), 2R_0\delta) \subset U$. It follows from Lemma 2.4 that

$$ D := B(\mathcal{T}_\delta(p), R_0\delta) \subset \mathcal{T}_\delta(B(W_\delta(p), 2R_0\delta)) $$

and that $\mathcal{T}_\delta$ is one-to-one if restricted to the connected component of $\mathcal{T}_\delta^{-1}(D)$ that contains $p$; note also that it contains the disc $B(W_\delta(p), 2C^{-1}R_0\delta)$ due to Corollary 5.4. In what follows, we replace $\Gamma_\delta$ by this connected component and assume that $\mathcal{T}_\delta$ is a global bijection from $\Theta_\delta$ onto the disc $D$.

Since we now assume that $\mathcal{T}_\delta$ is one-to-one, we can identify points $p \in \Theta_\delta$ with their images $z = \mathcal{T}_\delta(p) \in D$ and view all maps defined on $\Theta_\delta$ as maps on $D$ in order to simplify the notation. For example, in this new notation we have $\mathcal{T}_\delta(z) = z$ and $W_\delta(z) = z - \overline{O_\delta(z)}$. Corollary 5.4 implies that

$$ (5.5) \quad B(W_\delta(z), C^{-1}r\delta) \subset W_\delta(B(z, r\delta)) $$

whenever $r \ge C$ and $B(z, r\delta) \subset D$ (where the constant $C > 0$ is taken from Corollary 5.4).

Put $z_0 = \mathcal{T}_\delta(p)$ and fix $R \in [\frac{1}{2}R_0, R_0 - 2C]$ and $z_1 \in D$ such that $|z_1 - z_0| = R\delta$. Let $\varepsilon \in [0, 1]$ and $\alpha \in \mathbb{T}$ be defined by the equality

$$ \mathcal{O}_\delta(z_1) - \mathcal{O}_\delta(z_0) = (1 - \varepsilon)\bar{\alpha}^2(z_1 - z_0). $$

Our goal is to prove that $\varepsilon$ cannot be too close to 0. For $t \in [0, 1]$, denote $z_t = (1 - t)z_0 + tz_1$. Since $|d\mathcal{O}_\delta(z_t)| = |z_1 - z_0|dt$ almost everywhere and

$$ \int_0^1 \mathrm{Re} \left[ 1 - \alpha^2 \frac{d}{dt} \frac{\mathcal{O}_\delta(z_t)}{z_1 - z_0} \right] dt = \mathrm{Re} \left[ 1 - \alpha^2 \frac{\mathcal{O}_\delta(z_1) - \mathcal{O}_\delta(z_0)}{z_1 - z_0} \right] = \varepsilon, $$

Jensen's inequality for the concave function $[0, 2] \ni x \mapsto \sqrt{x^2 + (1 - (1 - x)^2)} = \sqrt{2x}$ implies that

$$ \left| \frac{(z_t - \alpha^2 \mathcal{O}_\delta(z_t)) - (z_0 - \alpha^2 \mathcal{O}(z_0))}{z_1 - z_0} \right| \le \int_0^1 \left| 1 - \alpha^2 \frac{d}{dt} \frac{\mathcal{O}_\delta(z_t)}{z_1 - z_0} \right| dt \le \sqrt{2\varepsilon} $$

for all $t \in [0, 1]$.

As $z \mapsto \mathcal{O}_\delta(z)$ is a 1-Lipschitz function, the mapping $D \ni z \mapsto z - \alpha^2 \mathcal{O}_\delta(z)$ has no overlaps; see also [CLR23, proof of Proposition 4.3]. Let $Q$ be the $2C\delta$-neighborhood of the segment $[z_0; z_1]$. It follows from the last estimate that

$$ (5.6) \quad \mathrm{Area}[(\mathcal{T}_\delta - \alpha^2 \mathcal{O}_\delta)(Q)] \le \pi(2C\delta + \sqrt{2\varepsilon}|z_1 - z_0|)^2 = \pi(2C + \sqrt{2\varepsilon}R)^2\delta^2. $$

On the other hand, Corollary 5.4 implies that $|W_\delta(z_1) - W_\delta(z_0)| \ge C^{-1}|z_1 - z_0|$. Splitting this segment into $\lfloor |W_\delta(z_1) - W_\delta(z_0)|/(2C^2\delta) \rfloor$ pieces of length $2C\delta$, using inclusion (5.5) for discs centered at the midpoints of these pieces, definition (5.2) of $\mu_\delta^\alpha$, and Proposition 5.5, we obtain the lower bound

$$ (5.7) \quad \mathrm{Area}[(\mathcal{T}_\delta - \alpha^2 \mathcal{O}_\delta)(Q)] \ge \frac{1}{4} \sum_{v \in W_\delta(Q')} \mu_\delta^\alpha(v) \ge \frac{1}{4} \left\lfloor \frac{|z_1 - z_0|}{2C^2\delta} \right\rfloor \cdot c\pi\delta^2, $$
where $Q'$ stands for the $C\delta$-neighborhood of $[z_0; z_1]$. Recall that $|z_1 - z_0| = R\delta$ and $\frac{1}{2}R_0 \le R \le R_0$. Provided that the constant $R_0$ is chosen large enough, combining (5.6) and (5.7) one obtains a lower bound on $\varepsilon$ that depends only on the constants in property (RW). $\square$

**6. Harmonicity in the intrinsic metric of the surface $\Theta$: proof of Theorem 4**

We begin with a few definitions and preliminary lemmas. Recall that $\mathcal{L}_{\varphi}h = -\operatorname{div}(A_{\varphi}\nabla h)$ where the matrix $A_{\varphi}$ is given by (1.8) and $\mathcal{L}_{\varphi}^{\times}h = \operatorname{div}(*A_{\varphi}^{-1}*\nabla h)$ is the 'dual' operator that describes the limits of A-harmonic conjugate functions to solutions of the equation $\mathcal{L}_{\varphi}h = 0$; see (1.10). It is easy to see that $\mathcal{L}_{\varphi}^{\times}\psi = 0$; in fact, the function $-i\psi$ is an A-harmonic conjugate of $w$.

**Lemma 6.1.** Let $Jac_w(\psi) = |\psi_w|^2 - |\psi_{\bar{w}}|^2$ denote the Jacobian of the mapping $\psi: U \to \mathbb{C} \cong \mathbb{R}^2$. The following equation holds in the weak sense:

$$ \mathcal{L}_{\varphi}\psi = -2\partial_{\bar{w}} Jac_w(\psi). $$

*Proof.* This is a straightforward computation based upon the identity

$$ \int_U h \mathcal{L}_{\varphi} g \, dxdy = 2 \int_U (h_w \ h_{\bar{w}}) \begin{pmatrix} -\psi_{\bar{w}} & \bar{\psi}_{\bar{w}} \\ \psi_w & -\bar{\psi}_w \end{pmatrix} \begin{pmatrix} g_w \\ g_{\bar{w}} \end{pmatrix} \, dxdy $$

applied to $h = \psi$ and a smooth compactly supported test function $g \in C_0^{\infty}(U)$. $\square$

Consider now the surface $\Theta = \{(\frac{1}{2}(w+\psi(w)); \frac{1}{2}(\overline{\psi(w)}-\overline{\bar{w}})) | w \in U\}$ from Theorem 4. In general, $\Theta$ is not smooth because $\psi$ is not. However $\Theta$ has tangent space almost everywhere by Rademacher's theorem since $\psi$ is Lipschitz. The inner product on $\mathbb{C}^{1,1}$ induces inner products in these tangent spaces. By a slight abuse of the terminology we regard this family of inner products as a Riemannian metric on $\Theta$. It is easy to see that this metric is positively definite as the uniform convexity of the potential $\varphi$ implies that $|\psi_{\bar{w}}| = 2|\varphi_{w\bar{w}}| < 2\varphi_{w\bar{w}} = \psi_w$; see also (3.4).

As we endowed $\Theta$ with a Riemannian metric, we can consider its conformal parametrization. More precisely, given a domain $D \subset \mathbb{C}$ we say that a map $D \ni \zeta \mapsto (z(\zeta), \theta(\zeta)) \in \Theta$ is a conformal parametrization of $\Theta$ if the mapping $\zeta \mapsto w(\zeta) = z(\zeta) + \theta(\zeta)$ is quasiconformal and we have

$$ (6.1) \qquad z_{\zeta}\bar{z}_{\zeta} = \theta_{\zeta}\bar{\theta}_{\zeta} $$

almost everywhere. (This definition makes sense since quasiconformal mappings are almost everywhere differentiable and the mappings $w \mapsto z$ and $w \mapsto \theta$ are Lipschitz.) The condition (6.1) is equivalent to requiring that the first fundamental form on $\Theta$ is (almost everywhere) diagonal in the coordinate $\zeta$.

**Lemma 6.2.** In the coordinate $w = x+iy$ the first fundamental form on $\Theta$ coincides with the Hessian of the potential $\varphi$, that is, the first fundamental form is given by

$$ \varphi_{xx} \, dx^2 + \varphi_{yy} \, dy^2 + 2\varphi_{xy} \, dxdy = \det A_{\varphi} \cdot (dx \ dy) A_{\varphi}^{-1} \begin{pmatrix} dx \\ dy \end{pmatrix}. $$

*Proof.* A straightforward computation expresses the first fundamental form on $\Theta$ as

$$ \frac{1}{4}(|\psi_w + 1|^2|dw|^2 + 2 \operatorname{Re}[(\psi_w + 1)\bar{\psi}_w(dw)^2]) - \frac{1}{4}(|\psi_w - 1|^2|dw|^2 + 2 \operatorname{Re}[(\psi_w - 1)\bar{\psi}_w](dw)^2), $$ which equals $\operatorname{Re}[\psi_w |dw|^2 + \bar{\psi}_w(dw)^2]$. Note that $\psi_w = \frac{1}{2}(\varphi_{xx} + \varphi_{yy})$ and $\bar{\psi}_w = \frac{1}{2}(\varphi_{xx} - \varphi_{yy}) - i\varphi_{xy}$. $\square$

**Corollary 6.3.** Denote $\rho = (\det A_{\varphi})^{1/2} = (\det D^2\varphi)^{1/2}$ and assume that $D \ni \zeta \mapsto (z(\zeta), \theta(\zeta))$ is a conformal parametrization of $\Theta$. Then,

$$ \operatorname{Jac}_{\zeta}(w) \cdot \mathcal{L}_{\varphi}h = -\operatorname{div}_{\zeta}(\rho \nabla_{\zeta}h), $$

$$ \operatorname{Jac}_{\zeta}(w) \cdot \mathcal{L}_{\varphi}^{\times}h = -\operatorname{div}_{\zeta}(\rho^{-1}\nabla_{\zeta}h) $$

for each $h$ with locally square summable gradient, where equalities are understood in the weak sense. (Note that $\rho$ is bounded away from zero and infinity since we assume that $\varphi$ is uniformly convex.)
*Proof.* By definition of a conformal coordinate, the matrix $A_{\varphi}^{-1} \det A_{\varphi}$ of the first fundamental form on $\Theta$ becomes a scalar function in the coordinate $\zeta$. Therefore, $A_{\varphi}^{-1} = \rho^{-1}(\det J)^{-1}J^{\top}J$, where $J \in \mathbb{R}^{2 \times 2}$ is the Jacobian matrix of the coordinate change $w(\zeta)$; in particular, $Jac_{\zeta}(w) = \det J$. Taking the inverse, we see that $A_{\varphi} = (\det J) \cdot \rho J^{-1}(J^{-1})^{\top}$ and we also have $\ast A_{\varphi}^{-1} \ast = -(\det J) \cdot \rho^{-1}J^{-1}(J^{-1})^{\top}$ since $\ast J^{\top} = (\det J) \cdot J^{-1} \ast$. The claim now follows by examining how quadratic forms of $\mathcal{L}_{\varphi}$ and $\mathcal{L}_{\varphi}^{\times}$ change under this coordinate change. $\square$

One says that a surface $\Theta \subset \mathbb{C}^{1,1} \cong \mathbb{R}^{2,2}$ is maximal if it is smooth and its mean curvature vector vanishes; we emphasize that the mean curvature is calculated with respect to the Minkovski metric (1.12). As in Euclidean spaces, vanishing of the mean curvature means that $\Theta$ is locally a stationary point of the area functional; however, for space-like surfaces in $\mathbb{R}^{2,2}$ such stationary points correspond to maxima rather than minima, hence the name. A direct computation similar to the classical Euclidean setup shows that $\Theta$ is maximal if and only if its conformal parametrization $(z(\zeta), \theta(\zeta))$ is harmonic, i.e.,

$$ (6.2) \qquad z_{\zeta\bar{\zeta}} = \theta_{\zeta\bar{\zeta}} = 0. $$

(Note that this is equivalent to saying that $w$ and $\psi$ are harmonic in $\zeta$.) In particular, maximal surfaces are real analytic. We are now in the position to prove Theorem 4.

*Proof of Theorem 4.* We have already proved (i) in Lemma 6.2, let us prove (ii).

First, assume that assertion (a) holds, i.e., that $\det D^2\varphi = \rho^2$ is constant. Since $\mathcal{L}_{\varphi}w = 0$ and $\mathcal{L}_{\varphi}^{\times}\psi = 0$, Corollary 6.3 implies that $w$ and $\psi$ are harmonic functions of $\zeta$, that is, (6.2) holds and the surface $\Theta$ is maximal as required in (b).

Second, let us prove that (b) implies (c). Applying Corollary 6.3 again and using the fact that $w$ is harmonic in $\zeta$ we can write

$$ 0 = \mathrm{Jac}_{\zeta}(w)\mathcal{L}_{\varphi}w = \mathrm{div}_{\zeta}(\rho\nabla_{\zeta}w) = \nabla_{\zeta}\rho \cdot \nabla_{\zeta}w, $$

which is only possible if $\nabla_{\zeta}\rho = 0$ since $\rho$ is a real-valued function and the map $\zeta \mapsto w$ preserves orientation. Therefore, both $\mathcal{L}_{\varphi}$ and $\mathcal{L}_{\varphi}^{\times}$ are proportional to the Laplacian in the coordinate $\zeta$.

Third, let us prove that (c) implies (d). Let $w(\zeta)$ be the coordinate change from (c); note that we do not assume in advance that $\zeta$ is a conformal parametrization of $\Theta$. However, as we have $\mathcal{L}_{\varphi}[\zeta] = 0$ and $\mathcal{L}_{\varphi}[\zeta^2] = 0$, the operator $\mathcal{L}_{\varphi}$ has to be proportional to the Laplacian in the coordinate $\zeta$. Therefore, the Jacobi matrix of $w(\zeta)$ diagonalizes $A_{\varphi}$, which implies that $\zeta$ is a conformal parametrization of $\Theta$. Using Corollary 6.3 and $\mathcal{L}_{\varphi}[\zeta] = 0$ again, we see that $\nabla_{\zeta}\rho = 0$. Hence, $\rho$ is constant and the operators $\mathcal{L}_{\varphi}$ and $\mathcal{L}_{\varphi}^{\times}$ are multiples of each other.

Finally, (d) trivially implies (e) since $\mathcal{L}_{\varphi}\psi = 0$ and it is easy to see that (e) implies (a). Indeed, it follows from Lemma 6.1 that $\partial_{\bar{w}}Jac_w(\psi) = -\frac{1}{2}\mathcal{L}_{\varphi}\psi = 0$, which means that the real-valued function $Jac_w(\psi) = \det A_{\varphi} = \det D^2\varphi$ is constant. $\square$

References

[AIM09] Kari Astala, Tadeusz Iwaniec, and Gaven Martin. *Elliptic partial differential equations and quasiconformal mappings in the plane*. Vol. 48. Princeton Mathematical Series. Princeton University Press, Princeton, NJ, 2009, pp. xviii+677. ISBN: 978-0-691-13777-3.

[BG25] Ahmed Bou-Rabee and Ewain Gwynne. “Random walk on sphere packings and Delaunay triangulations in arbitrary dimension”. In: *Proc. Lond. Math. Soc. (3)* 131.3 (2025), Paper No. e70084. ISSN: 0024-6115. DOI: 10.1112/plms.70084. URL: https://doi.org/10.1112/plms.70084.

[BGS25] Federico Bertacco, Ewain Gwynne, and Scott Sheffield. “Scaling limits of planar maps under the Smith embedding”. In: *Ann. Probab.* 53.3 (2025), pp. 1138–1196. ISSN: 0091-1798. DOI: 10.1214/24-aop1731. URL: https://doi.org/10.1214/24-aop1731.
[BLQ23] Nathanaël Berestycki, Marcin Lis, and Wei Qian. "Free boundary dimers: random walk representation and scaling limit". English. In: *Probab. Theory Relat. Fields* 186.3-4 (2023), pp. 735-812. ISSN: 0178-8051. DOI: 10.1007/s00440-023-01203-x.

[BNR25] Tomas Berggren, Matthew Nicoletti, and Marianna Russkikh. "Perfect t-embeddings of doubly periodic Aztec diamonds". In: *ArXiv e-prints* (2025). arXiv: 2508.04938 [math.PR]. URL: https://arxiv.org/abs/2508.04938.

[BP25] Ilia Binder and David Pechersky. "Orthodiagonal Maps, Tilings of Rectangles, and their Convergence to Conformal Maps". In: *ArXiv e-prints* (2025). arXiv: 2407.20851 [math.CV]. URL: https://arxiv.org/abs/2407.20851.

[Bro+40] R. L. Brooks et al. "The dissection of rectangles into squares". In: *Duke Math. J.* 7 (1940), pp. 312-340. ISSN: 0012-7094. URL: http://projecteuclid.org/euclid.dmj/1077492259.

[BTR17] Cédric Boutillier, Béatrice de Tilière, and Kilian Raschel. "The Z-invariant massive Laplacian on isoradial graphs". In: *Invent. Math.* 208.1 (2017), pp. 109-189. ISSN: 0020-9910. DOI: 10.1007/s00222-016-0687-z. URL: https://doi.org/10.1007/s00222-016-0687-z.

[CFL28] R. Courant, K. Friedrichs, and H. Lewy. "Über die partiellen Differenzengleichungen der mathematischen Physik". In: *Math. Ann.* 100.1 (1928), pp. 32-74. ISSN: 0025-5831. DOI: 10.1007/BF01448839. URL: https://doi.org/10.1007/BF01448839.

[CG97] Luis A Caffarelli and Cristian E Gutiérrez. "Properties of the solutions of the linearized Monge-Ampère equation". In: *American Journal of Mathematics* 119.2 (1997), pp. 423-465.

[Che16] Dmitry Chelkak. "Robust discrete complex analysis: a toolbox". In: *Ann. Probab.* 44.1 (2016), pp. 628-683. ISSN: 0091-1798. DOI: 10.1214/14-A0P985. URL: https://doi.org/10.1214/14-A0P985.

[Che24] Dmitry Chelkak. "Ising model and s-embeddings of planar graphs". In: *Ann. Sci. Éc. Norm. Supér.* (4) 57.5 (2024), pp. 1271-1346. ISSN: 0012-9593.

[CLR21] Dmitry Chelkak, Benoît Laslier, and Marianna Russkikh. "Bipartite dimer model: perfect t-embeddings and Lorentz-minimal surfaces". In: *ArXiv e-prints* (2021). arXiv: 2109.06272 [math.PR]. URL: https://arxiv.org/abs/2109.06272.

[CLR23] Dmitry Chelkak, Benoît Laslier, and Marianna Russkikh. "Dimer model and holomorphic functions on t-embeddings of planar graphs". In: *Proc. Lond. Math. Soc.* (3) 126.5 (2023), pp. 1656-1739. ISSN: 0024-6115. DOI: 10.1112/plms.12516. URL: https://doi.org/10.1112/plms.12516.

[CR24] Dmitry Chelkak and Sanjay Ramassamy. "Fluctuations in the Aztec diamonds via a space-like maximal surface in Minkowski 3-space". In: *Confluentes Math.* 16 (2024), pp. 1-17.

[CS11] Dmitry Chelkak and Stanislav Smirnov. "Discrete complex analysis on isoradial graphs". In: *Adv. Math.* 228.3 (2011), pp. 1590-1630. ISSN: 0001-8708. DOI: 10.1016/j.aim.2011.06.025. URL: https://doi.org/10.1016/j.aim.2011.06.025.

[CS12] Dmitry Chelkak and Stanislav Smirnov. "Universality in the 2D Ising model and conformal invariance of fermionic observables". In: *Invent. Math.* 189.3 (2012), pp. 515-580. ISSN: 0020-9910. DOI: 10.1007/s00222-011-0371-2. URL: https://doi.org/10.1007/s00222-011-0371-2.

[Duf68] R. J. Duffin. "Potential theory on a rhombic lattice". English. In: *J. Comb. Theory* 5 (1968), pp. 258-272. ISSN: 0097-3165. DOI: 10.1016/S0021-9800(68)80072-9.

[FF56] L. R. Ford Jr. and D. R. Fulkerson. "Maximal flow through a network". In: *Canadian J. Math.* 8 (1956), pp. 399-404. ISSN: 0008-414X. DOI: 10.4153/CJM-1956-045-5. URL: https://doi.org/10.4153/CJM-1956-045-5.

[Fig17] Alessio Figalli. *The Monge-Ampere equation and its applications*. European Mathematical Society Zürich, 2017.
[GHS23] Ewain Gwynne, Nina Holden, and Xin Sun. "Mating of trees for random planar maps and Liouville quantum gravity: a survey". In: Topics in statistical mechanics. Vol. 59. Panor. Synthèses. Soc. Math. France, Paris, 2023, pp. 41–120.

[GJN20] Ori Gurel-Gurevich, Daniel C. Jerison, and Asaf Nachmias. "The Dirichlet problem for orthodiagonal maps". In: Adv. Math. 374 (2020), pp. 107–379. ISSN: 0001-8708. DOI: 10.1016/j.aim.2020.107379. URL: https://doi.org/10.1016/j.aim.2020.107379.

[GMS21] Ewain Gwynne, Jason Miller, and Scott Sheffield. "The Tutte embedding of the mated-CRT map converges to Liouville quantum gravity". In: Ann. Probab. 49.4 (2021), pp. 1677–1717. ISSN: 0091-1798. DOI: 10.1214/20-aop1487. URL: https://doi.org/10.1214/20-aop1487.

[GN11] Cristian E. Gutiérrez and Truyen Nguyen. "Interior gradient estimates for solutions to the linearized Monge-Ampère equation". In: Adv. Math. 228.4 (2011), pp. 2034–2070. ISSN: 0001-8708. DOI: 10.1016/j.aim.2011.06.035. URL: https://doi.org/10.1016/j.aim.2011.06.035.

[GT01] David Gilbarg and Neil S. Trudinger. Elliptic Partial Differential Equations of Second Order. 2nd ed. Classics in Mathematics. Springer Berlin, Heidelberg, 2001, pp. xviii+677. ISBN: 978-3-540-41160-4.

[Ken+00] Carlos Kenig et al. "A new approach to absolute continuity of elliptic measure, with applications to non-symmetric equations". In: Advances in mathematics 153.2 (2000), pp. 231–298.

[Ken+21] Richard Kenyon et al. "Dimers and circle patterns". In: Annales Scientifiques de l'École Normale Supérieure (2021).

[Ken02] R. Kenyon. "The Laplacian and Dirac operators on critical planar graphs". In: Invent. Math. 150.2 (2002), pp. 409–439. ISSN: 0020-9910. DOI: 10.1007/s00222-002-0249-4. URL: https://doi.org/10.1007/s00222-002-0249-4.

[Lel55] Jacqueline Lelong-Ferrand. Représentation conforme et transformations à intégrale de Dirichlet bornée. Gauthier-Villars, Paris, 1955, pp. viii+259.

[LP16] Russell Lyons and Yuval Peres. Probability on trees and networks. Vol. 42. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, New York, 2016, pp. xv+699. ISBN: 978-1-107-16015-6. DOI: 10.1017/9781316672815. URL: https://doi.org/10.1017/9781316672815.

[LSW63] Walter Littman, Guido Stampacchia, and Hans F Weinberger. "Regular points for elliptic equations with discontinuous coefficients". In: Annali della Scuola Normale Superiore di Pisa-Scienze Fisiche e Matematiche 17.1-2 (1963), pp. 43–77.

[Lyu26] L. A. Lyusternik. "Über einige Anwendungen der direkten Methoden in der Variationsrechnung." German. In: Rec. Math. Moscou 33 (1926), pp. 173–201.

[Mah25] Rémy Mahfouf. "Crossing estimates for the Ising model on general s-embeddings". In: Proc. Lond. Math. Soc. (3) 131.4 (2025), Paper No. e70092. ISSN: 0024-6115. DOI: 10.1112/plms.70092. URL: https://doi.org/10.1112/plms.70092.

[Mer01] Christian Mercat. "Discrete Riemann surfaces and the Ising model". In: Comm. Math. Phys. 218.1 (2001), pp. 177–216. ISSN: 0010-3616. DOI: 10.1007/s002200000348. URL: https://doi.org/10.1007/s002200000348.

[Pec25] David Pechersky. "A Polynomial Rate of Convergence for the Dirichlet Problem on Orthodiagonal Maps". In: ArXiv e-prints (2025). arXiv: 2503.20284 [math.PR]. URL: https://arxiv.org/abs/2503.20284.

[Sko13] M. Skopenkov. "The boundary value problem for discrete analytic functions". In: Adv. Math. 240 (2013), pp. 61–87. ISSN: 0001-8708. DOI: 10.1016/j.aim.2013.03.002. URL: https://doi.org/10.1016/j.aim.2013.03.002.
[Smi10] Stanislav Smirnov. "Discrete complex analysis and probability". In: *Proceedings of the International Congress of Mathematicians. Volume I*. Hindustan Book Agency, New Delhi, 2010, pp. 595-621.

[SW66] James Serrin and HF Weinberger. "Isolated singularities of solutions of linear elliptic equations". In: *American Journal of Mathematics* 88.1 (1966), pp. 258-272.

[TKB13] Justin L Taylor, Seick Kim, and Russell Murray Brown. "The Green function for elliptic systems in two dimensions". In: *Communications in Partial Differential Equations* 38.9 (2013), pp. 1574-1600.

[Tut63] William Thomas Tutte. "How to draw a graph". In: *Proceedings of the London Mathematical Society* 3.1 (1963), pp. 743-767.