# The Z-Dirac and massive Laplacian operators in the Z-invariant Ising model

Béatrice de Tilière

► To cite this version:

Béatrice de Tilière. The Z-Dirac and massive Laplacian operators in the Z-invariant Ising model. Electronic Journal of Probability, 2021, 26, 26 (53), pp.1-86. 10.1214/21-EJP601. hal-03876780

HAL Id: hal-03876780

https://hal.science/hal-03876780

Submitted on 28 May 2024

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The Z-Dirac and massive Laplacian operators in the Z-invariant Ising model

Béatrice de Tilière

► To cite this version:

Béatrice de Tilière. The Z-Dirac and massive Laplacian operators in the Z-invariant Ising model. Electronic Journal of Probability, 2021, 26, 26 (53), pp.1-86. 10.1214/21-EJP601. hal-03876780

HAL Id: hal-03876780

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

Submitted on 28 May 2024

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Electron. J. Probab. 26 (2021), article no. 53, 1–86.
ISSN: 1083-6489 https://doi.org/10.1214/21-EJP601

# The $Z$-Dirac and massive Laplacian operators in the $Z$-invariant Ising model

Béatrice de Tilière*

**Abstract**

Consider an elliptic parameter $k$; we introduce a family of $Z^u$-Dirac operators $(K(u))_{u \in C}$, relate them to the $Z$-massive Laplacian of [11], and extend to the full $Z$-invariant case the results of Kenyon [45] on discrete holomorphic and harmonic functions, which correspond to the case $k = 0$. We prove through combinatorial identities, how and why the $Z^u$-Dirac and $Z$-massive Laplacian operators appear in the $Z$-invariant Ising model, considering the case of infinite and finite isoradial graphs. More precisely, consider the dimer model on the Fisher graph $G^F$ arising from a $Z$-invariant Ising model. We express coefficients of the inverse Fisher Kasteleyn operator as a function of the inverse $Z^u$-Dirac operator and also as a function of the $Z$-massive Green function; in particular this proves a (massive) random walk representation of important observables of the Ising model. We prove that the squared partition function of the Ising model is equal, up to a constant, to the determinant of the $Z$-massive Laplacian operator with specific boundary conditions, the latter being the partition function of rooted spanning forests. To show these results, we relate the inverse Fisher Kasteleyn operator and that of the dimer model on the bipartite graph $G^Q$ arising from the XOR-Ising model, and we prove matrix identities between the Kasteleyn matrix of $G^Q$ and the $Z^u$-Dirac operator, that allow to reach inverse matrices as well as determinants.

**Keywords:** Ising model; massive Laplacian and Dirac operators; discrete massive harmonic and holomorphic functions; dimer model; spanning forests and spanning trees; $Z$-invariance.

**MSC2020 subject classifications:** 82B20; 82B23; 33E05; 05A19.

Submitted to EJP on February 15, 2018, final version accepted on March 3, 2021.

Supersedes arXiv:1801.00207v2.

# 1 Introduction

This paper is inspired by three sets of results suggesting connections between the Ising model on a planar graph $G$ and (massive) random walks on $G$ and its dual $G^*$.

<footnote>*PSL University Paris-Dauphine, CNRS, UMR 7534, CEREMADE, 75016 Paris, France. Institut universitaire de France. E-mail: detiliere@ceremade.dauphine.fr</footnote>
<figure><img src="image_2.png" /></figure>

Electron. J. Probab. 26 (2021), article no. 53, 1–86.
ISSN: 1083-6489 https://doi.org/10.1214/21-EJP601

# The Z-Dirac and massive Laplacian operators in the Z-invariant Ising model

Béatrice de Tilière*

## Abstract

Consider an elliptic parameter $k$; we introduce a family of $Z^u$-Dirac operators $(K(u))_{u \in C}$, relate them to the $Z$-massive Laplacian of [11], and extend to the full $Z$-invariant case the results of Kenyon [45] on discrete holomorphic and harmonic functions, which correspond to the case $k = 0$. We prove through combinatorial identities, how and why the $Z^u$-Dirac and $Z$-massive Laplacian operators appear in the $Z$-invariant Ising model, considering the case of infinite and finite isoradial graphs. More precisely, consider the dimer model on the Fisher graph $G^F$ arising from a $Z$-invariant Ising model. We express coefficients of the inverse Fisher Kasteleyn operator as a function of the inverse $Z^u$-Dirac operator and also as a function of the $Z$-massive Green function; in particular this proves a (massive) random walk representation of important observables of the Ising model. We prove that the squared partition function of the Ising model is equal, up to a constant, to the determinant of the $Z$-massive Laplacian operator with specific boundary conditions, the latter being the partition function of rooted spanning forests. To show these results, we relate the inverse Fisher Kasteleyn operator and that of the dimer model on the bipartite graph $G^Q$ arising from the XOR-Ising model, and we prove matrix identities between the Kasteleyn matrix of $G^Q$ and the $Z^u$-Dirac operator, that allow to reach inverse matrices as well as determinants.

**Keywords:** Ising model; massive Laplacian and Dirac operators; discrete massive harmonic and holomorphic functions; dimer model; spanning forests and spanning trees; Z-invariance.

**MSC2020 subject classifications:** 82B20; 82B23; 33E05; 05A19.

Submitted to EJP on February 15, 2018, final version accepted on March 3, 2021.

Supersedes arXiv:1801.00207v2.

## 1 Introduction

This paper is inspired by three sets of results suggesting connections between the Ising model on a planar graph $G$ and (massive) random walks on $G$ and its dual $G^*$.

<footnote>*PSL University Paris-Dauphine, CNRS, UMR 7534, CEREMADE, 75016 Paris, France. Institut universitaire de France. E-mail: detiliere@ceremade.dauphine.fr</footnote>
The Z-Dirac and massive Laplacian operators in the Z-invariant Ising model

• Messikh [60] observes that large deviation estimates of a massive random walk occur when computing the correlation length of the super-critical Ising model on $Z^2$; a result later proved by Beffara and Duminil-Copin [7] using the FK-Ising observable of [68] away from the critical point.

• In Smirnov and Chelkak-Smirnov's proof of conformal invariance of the critical Z-invariant Ising model [67, 68, 18], the key discrete tools are observables – spin or FK (see also [39]) – that are holomorphic. Discrete holomorphic functions in turn are naturally related to harmonic functions [33, 59, 45, 17]. The paper [58] addresses the question of the massive version of these observables and that of proving convergence to massive SLE's. In particular, they give ideas of proofs in the case of the massive harmonic explorer and of loop erased random walks, see [19] for a complete proof. As mentioned in the previous point, a massive version of the FK-Ising observable is fruitfully used in the paper [7].

• We prove, through combinatorial constructions, that the squared partition function of the critical Z-invariant Ising model is equal, up to a multiplicative constant, to the partition function of spanning trees [29, 30]. An abstract proof of this identity is given in the toroidal Z-invariant case in [10, 12].

The main contribution of this paper is to provide a unified framework for all of the above, which holds in the full Z-invariant case, in the infinite and finite cases. Our main results are obtained as a combination of intermediate steps that are interesting in their own respect. We nevertheless feel that, before listing statements leading to the principal Ising results, we should convey the main ideas.

Let us first be more precise about operators underlying our "inspiration" papers. Large deviation estimates of massive random walks are related to the massive Green function, the latter being the inverse of the massive Laplacian operator. By definition discrete holomorphic functions are in the kernel of the Dirac operator, which is a Kasteleyn matrix/operator of the double graph $G^D$ [45]; harmonic functions are in the kernel of the Laplacian operator. The spanning tree partition function is equal to the determinant of the Laplacian operator [50]. Summarizing, a central role is played by the Dirac operator (at criticality) and the (massive) Laplacian in the (super) critical Ising model.

Our first contribution is to introduce one of the missing pieces of the puzzle, namely the full Z-invariant version of the (critical) Dirac operator of [45], referred to as the $Z^u$-Dirac operator, $u$ being a natural free parameter disappearing at criticality. This is the subject of Section 3, as well as its connections to the Z-invariant massive Laplacian of [11] and the study of the corresponding dimer model on the double graph $G^D$.

To study the Ising model, we use Fisher's correspondence [35] relating the low (or high) temperature expansion of the model [51, 52] to the dimer model on the Fisher graph $G^F$ with associated Kasteleyn matrix/operator $K^F$. The partition function of the dimer model is the Pfaffian of $K^F$, and the Boltzmann/Gibbs measures are explicitly expressed using coefficients of $K^F$ and its inverse $(K^F)^{-1}$ [72, 44, 23, 47, 9]. This means that knowing the determinant of $K^F$ and its inverse amounts to fully understanding the partition function of the Ising model and probabilities of its low (or high) temperature expansion. Notably, coefficients of the inverse Kasteleyn operator $(K^F)^{-1}$ are also related to other important observables of the Ising model as the spin-Ising observable of [18], see [69, 16], and the FK-Ising observable of [68, 18] which is also the fermionic spinor observable of [39] (up to normalization), see [32, 63].

Consider the dimer model on $G^F$ arising from a Z-invariant Ising model. Our main contribution is to prove matrix identities relating the Kasteleyn operator $K^F$ and the $Z^u$-Dirac operator and also the Z-massive Laplacian of [11]. The strength of these identities is that they allow to reach inverse operators and also, after some extra
• Messikh [60] observes that large deviation estimates of a massive random walk occur when computing the correlation length of the super-critical Ising model on $Z^2$; a result later proved by Beffara and Duminil-Copin [7] using the FK-Ising observable of [68] away from the critical point.

• In Smirnov and Chelkak-Smirnov's proof of conformal invariance of the critical Z-invariant Ising model [67, 68, 18], the key discrete tools are observables – spin or FK (see also [39]) – that are holomorphic. Discrete holomorphic functions in turn are naturally related to harmonic functions [33, 59, 45, 17]. The paper [58] addresses the question of the massive version of these observables and that of proving convergence to massive SLE's. In particular, they give ideas of proofs in the case of the massive harmonic explorer and of loop erased random walks, see [19] for a complete proof. As mentioned in the previous point, a massive version of the FK-Ising observable is fruitfully used in the paper [7].

• We prove, through combinatorial constructions, that the squared partition function of the critical Z-invariant Ising model is equal, up to a multiplicative constant, to the partition function of spanning trees [29, 30]. An abstract proof of this identity is given in the toroidal Z-invariant case in [10, 12].

The main contribution of this paper is to provide a unified framework for all of the above, which holds in the full Z-invariant case, in the infinite and finite cases. Our main results are obtained as a combination of intermediate steps that are interesting in their own respect. We nevertheless feel that, before listing statements leading to the principal Ising results, we should convey the main ideas.

Let us first be more precise about operators underlying our "inspiration" papers. Large deviation estimates of massive random walks are related to the massive Green function, the latter being the inverse of the massive Laplacian operator. By definition discrete holomorphic functions are in the kernel of the Dirac operator, which is a Kasteleyn matrix/operator of the double graph $G^D$ [45]; harmonic functions are in the kernel of the Laplacian operator. The spanning tree partition function is equal to the determinant of the Laplacian operator [50]. Summarizing, a central role is played by the Dirac operator (at criticality) and the (massive) Laplacian in the (super) critical Ising model.

Our first contribution is to introduce one of the missing pieces of the puzzle, namely the full Z-invariant version of the (critical) Dirac operator of [45], referred to as the $Z^u$-Dirac operator, $u$ being a natural free parameter disappearing at criticality. This is the subject of Section 3, as well as its connections to the Z-invariant massive Laplacian of [11] and the study of the corresponding dimer model on the double graph $G^D$.

To study the Ising model, we use Fisher's correspondence [35] relating the low (or high) temperature expansion of the model [51, 52] to the dimer model on the Fisher graph $G^F$ with associated Kasteleyn matrix/operator $K^F$. The partition function of the dimer model is the Pfaffian of $K^F$, and the Boltzmann/Gibbs measures are explicitly expressed using coefficients of $K^F$ and its inverse $(K^F)^{-1}$ [72, 44, 23, 47, 9]. This means that knowing the determinant of $K^F$ and its inverse amounts to fully understanding the partition function of the Ising model and probabilities of its low (or high) temperature expansion. Notably, coefficients of the inverse Kasteleyn operator $(K^F)^{-1}$ are also related to other important observables of the Ising model as the spin-Ising observable of [18], see [69, 16], and the FK-Ising observable of [68, 18] which is also the fermionic spinor observable of [39] (up to normalization), see [32, 63].

Consider the dimer model on $G^F$ arising from a Z-invariant Ising model. Our main contribution is to prove matrix identities relating the Kasteleyn operator $K^F$ and the $Z^u$-Dirac operator and also the Z-massive Laplacian of [11]. The strength of these identities is that they allow to reach inverse operators and also, after some extra
work, determinants. As a consequence, in the finite and infinite cases, we express coefficients of $(K^F)^{-1}$ using the inverse $Z^u$-Dirac operator and also using the Z-massive Green function; this is the subject of Section 5, see also the corresponding part of the introduction. In essence, this proves that the contour Ising Boltzmann/Gibbs measures can be computed from (massive) random walks (with specific boundary conditions in the finite case). In the finite case we also prove that the squared Ising partition function is equal, up to an explicit constant, to the determinant of the massive Laplacian, that is to the partition function of rooted spanning forests, see Corollary 4.5 and also Theorem 1.7 of the introduction. Comments on how these results connect to our "inspiration" and other papers are given at the end of this section.

Section 2 contains preliminaries. Section 4 contains the main intermediate step: we consider the dimer model on the bipartite graph $G^Q$ arising from the XOR-Ising model [74] constructed from two independent Z-invariant Ising models [32, 13]. We prove matrix identities relating its Kasteleyn matrix/operator $K^Q$ and the $Z^u$-Dirac operator. In Section 5, building on the work of Dubédat [32], we express coefficients of the inverse operator $(K^F)^{-1}$ using coefficients of the inverse operator $(K^Q)^{-1}$; this result holds for the dimer model on the Fisher graph $G^F$ arising from any 2d-Ising model, not necessarily Z-invariant. Note that this result can also be derived from the paper [16], see Remark 1.9. Combining this with the results of Section 4 then allows us to deduce the Ising results. In Section 6, we specify some of our results in two important cases: the Z-invariant critical case, and the full Z-invariant case when the underlying graph is $Z^2$.

To give detailed statements, let us be more precise about Z-invariant models [64, 43], fully developed by Baxter [4, 5, 6], see also [65, 2, 3]. A Z-invariant model is naturally defined on an isoradial graph $G = (V, E)$; parameters are chosen so that the partition function only changes by a constant when performing a star-triangle transformation of the underlying graph, i.e., they are required to satisfy the Yang-Baxter equations. The solution to this set of equations for the Ising model has, given the embedding of the graph, a free elliptic parameter $k$, such that $(k')^2 := 1 - k^2 \in (0, \infty)$, and the coupling constants $J$ are [6]:

$$ \forall e \in E, \quad J_e = \frac{1}{2} \ln \left( \frac{1 + \operatorname{sn}(\theta_e | k)}{\operatorname{cn}(\theta_e | k)} \right), $$

where $\operatorname{sn}, \operatorname{cn}$ are two of the Jacobi elliptic trigonometric functions, and $\bar{\theta}_e = \theta_e \frac{\pi}{2K}$ is an angle associated to the edge $e$ in the isoradial embedding. When $k = 0$, i.e. $k' = 1$, the elliptic functions $\operatorname{sn}, \operatorname{cn}$, are the trigonometric functions $\operatorname{sin}, \operatorname{cos}$, and the Ising model is critical [55, 22, 57]. As $(k')^2$ varies from 0 to $\infty$, the coupling constants range from $\infty$ to 0 [12] thus covering the whole range of inverse temperatures. In the paper [11], we introduce Z-invariant rooted spanning forests with associated operator the massive Laplacian $\Delta^m$; when $k = 0$, we recover the critical Laplacian of [45]. We also prove an explicit local expression for its inverse, the Z-massive Green function $G^m$, using the discrete massive exponential function [11]. We are now ready to give a detailed overview of this paper.

**Section 3: $Z^u$-Dirac and Z-massive Laplacian operators** Fix an elliptic parameter $k$, and let $T(k)$ be the torus $C/(4KZ + 4iK'Z)$, where $K, K'$ are defined in Section 2.5.4. Denote by $C(k)$ the component $\mathbb{R}/4KZ$ of $T(k)$; whenever no confusion occurs, we will omit the argument $k$ in $C(k)$. We introduce a family of $Z^u$-Dirac operators $(K(u))_{u \in C}$ on the double graph $G^D = (W \cup B, E^D)$ associated to pairs of dual directed spanning trees, extending to the full Z-invariant case the Dirac operator $\bar{\partial}$ of [45], corresponding to $k = 0$. In the finite case, we introduce a family of operators $(K^\partial(u))_{u \in C'}$, with boundary conditions tuned for the Ising model, where $C'$ is a subset of $C$ defined in Equation (3.4). Although these operators play a key role in the Z-invariant Ising model,
work, determinants. As a consequence, in the finite and infinite cases, we express coefficients of $(K^F)^{-1}$ using the inverse $Z^u$-Dirac operator and also using the Z-massive Green function; this is the subject of Section 5, see also the corresponding part of the introduction. In essence, this proves that the contour Ising Boltzmann/Gibbs measures can be computed from (massive) random walks (with specific boundary conditions in the finite case). In the finite case we also prove that the squared Ising partition function is equal, up to an explicit constant, to the determinant of the massive Laplacian, that is to the partition function of rooted spanning forests, see Corollary 4.5 and also Theorem 1.7 of the introduction. Comments on how these results connect to our "inspiration" and other papers are given at the end of this section.

Section 2 contains preliminaries. Section 4 contains the main intermediate step: we consider the dimer model on the bipartite graph $G^Q$ arising from the XOR-Ising model [74] constructed from two independent Z-invariant Ising models [32, 13]. We prove matrix identities relating its Kasteleyn matrix/operator $K^Q$ and the $Z^u$-Dirac operator. In Section 5, building on the work of Dubédat [32], we express coefficients of the inverse operator $(K^F)^{-1}$ using coefficients of the inverse operator $(K^Q)^{-1}$; this result holds for the dimer model on the Fisher graph $G^F$ arising from any 2d-Ising model, not necessarily Z-invariant. Note that this result can also be derived from the paper [16], see Remark 1.9. Combining this with the results of Section 4 then allows us to deduce the Ising results. In Section 6, we specify some of our results in two important cases: the Z-invariant critical case, and the full Z-invariant case when the underlying graph is $Z^2$.

To give detailed statements, let us be more precise about Z-invariant models [64, 43], fully developed by Baxter [4, 5, 6], see also [65, 2, 3]. A Z-invariant model is naturally defined on an isoradial graph $G = (V, E)$; parameters are chosen so that the partition function only changes by a constant when performing a star-triangle transformation of the underlying graph, i.e., they are required to satisfy the Yang-Baxter equations. The solution to this set of equations for the Ising model has, given the embedding of the graph, a free elliptic parameter $k$, such that $(k')^2 := 1 - k^2 \in (0, \infty)$, and the coupling constants $J$ are [6]:

$$ \forall e \in E, \quad J_e = \frac{1}{2} \ln \left( \frac{1 + \operatorname{sn}(\theta_e | k)}{\operatorname{cn}(\theta_e | k)} \right), $$

where $\operatorname{sn}, \operatorname{cn}$ are two of the Jacobi elliptic trigonometric functions, and $\bar{\theta}_e = \theta_e \frac{\pi}{2K}$ is an angle associated to the edge $e$ in the isoradial embedding. When $k = 0$, i.e. $k' = 1$, the elliptic functions $\operatorname{sn}, \operatorname{cn}$, are the trigonometric functions $\operatorname{sin}, \operatorname{cos}$, and the Ising model is critical [55, 22, 57]. As $(k')^2$ varies from 0 to $\infty$, the coupling constants range from $\infty$ to 0 [12] thus covering the whole range of inverse temperatures. In the paper [11], we introduce Z-invariant rooted spanning forests with associated operator the massive Laplacian $\Delta^m$; when $k = 0$, we recover the critical Laplacian of [45]. We also prove an explicit local expression for its inverse, the Z-massive Green function $G^m$, using the discrete massive exponential function [11]. We are now ready to give a detailed overview of this paper.

**Section 3: $Z^u$-Dirac and Z-massive Laplacian operators** Fix an elliptic parameter $k$, and let $T(k)$ be the torus $C/(4KZ + 4iK'Z)$, where $K, K'$ are defined in Section 2.5.4. Denote by $C(k)$ the component $\mathbb{R}/4KZ$ of $T(k)$; whenever no confusion occurs, we will omit the argument $k$ in $C(k)$. We introduce a family of $Z^u$-Dirac operators $(K(u))_{u \in C}$ on the double graph $G^D = (W \cup B, E^D)$ associated to pairs of dual directed spanning trees, extending to the full Z-invariant case the Dirac operator $\bar{\partial}$ of [45], corresponding to $k = 0$. In the finite case, we introduce a family of operators $(K^\partial(u))_{u \in C'}$, with boundary conditions tuned for the Ising model, where $C'$ is a subset of $C$ defined in Equation (3.4). Although these operators play a key role in the Z-invariant Ising model,
they are interesting in their own respect. In the specific case $k = 0$, results we obtain
can be found in [45, 71, 43]. We prove, see also Theorem 3.6:

**Theorem 1.1.**

* Infinite case. Let $u \in C$, then the $Z^u$-Dirac operator $K(u)$, the Z-massive Laplacian $\Delta^m$ and the dual $\Delta^{m,*}$ of [11] satisfy the following identity:

$$
\overline{K(u)}^t K(u) = k' \begin{pmatrix} \Delta^m & 0 \\ 0 & \Delta^{m,*} \end{pmatrix}.
$$

* Finite case. Let $u \in C'$, then the $Z^u$-Dirac operators $K(u)$, $K^\partial(u)$, the Z-massive Laplacian $\Delta^{m,\partial}(u)$ and the dual $\Delta^{m,*}$ satisfy the following identity:

$$
\overline{K^{\partial}(u)}^t K(u) = k' \begin{pmatrix} \Delta^{m,\partial}(u) & Q(u) \\ 0 & \Delta^{m,*} \end{pmatrix}.
$$

A function $F \in C^B$ is said to be $Z^u$-holomorphic if $K(u)F = 0$. As a consequence of Theorem 1.1, if $F$ is $Z^u$-holomorphic, then $F|_{V}$ is $Z$-massive harmonic on $G$ and $F|_{V^*}$ is $Z$-massive harmonic on $G^*$, thus explaining the part “Dirac” in “$Z^u$-Dirac operator”.

Note that in the finite case, the Z-massive Laplacian $\Delta^{m,\partial}(u)$ arising in the context of the Ising model might have negative masses and conductances along the boundary, see Equation (3.6) where it is explicitly defined. As a consequence, this operator is not necessarily positive definite.

In the infinite case, Theorem 1.1 yields the following relations for inverse operators,
see also Corollary 3.15; the statement in the finite case is given in Corollary 3.17.

**Corollary 1.2 (Infinite case).** For every $u \in C$, consider the operator $K(u)^{-1}$ mapping $C^W$ to $C^B$ whose coefficients are defined by, for every $v, f, w$ as in Figure 14,

$$
\begin{align*}
K(u)_{v,w}^{-1} &= e^{-i \frac{\alpha_f + \beta_f}{2}} (k')^{-1} \mathrm{sc}(\theta_f)^{\frac{1}{2}} \left( [\mathrm{dn}(u_{\alpha_f}) \mathrm{dn}(u_{\beta_f})]^{\frac{1}{2}} G_{v,v_2}^m - [\mathrm{dn}(u_{\alpha_f+2K}) \mathrm{dn}(u_{\beta_f+2K})]^{\frac{1}{2}} G_{v,v_1}^m \right) \\
K(u)_{f,w}^{-1} &= -ie^{-i \frac{\alpha_f + \beta_f}{2}} (k')^{-1} \mathrm{sc}(\theta_f^*)^{\frac{1}{2}} \left( [\mathrm{dn}((u_{\beta_f})^*) \mathrm{dn}((u_{\alpha_f+2K})^*)]^{\frac{1}{2}} G_{f,f_2}^{m,*} + \right. \\
&\qquad \left. - [\mathrm{dn}((u_{\beta_f-2K})^*) \mathrm{dn}((u_{\alpha_f})^*)]^{\frac{1}{2}} G_{f,f_1}^{m,*} \right),
\end{align*}
$$

where $G^m$ and $G^{m,*}$ are the Z-massive and dual Z-massive Green functions of [11]. Then $K(u)^{-1}$ is the unique inverse of the $Z^u$-Dirac operator $K(u)$ decreasing to zero at infinity.

This gives, in Theorem 3.19, an explicit local expression for a Gibbs measure of the
dimer model on the double graph $G^D$, where the locality property is inherited from that
of the $Z$-massive Green functions of [11]. Using the KPW-Temperley bijection [71, 49],
probabilities of pairs of dual directed spanning trees are computed using the Green
function of a massive, non-directed random walk. Apart from the locality property which
is specific, a similar result is obtained by Chhita [20] in the case of $\mathbb{Z}^2$ with a specific
choice of weights.

In Theorem 3.8 and Corollary 3.10, we restrict to the finite case and prove relations
on determinants; we show,

**Theorem 1.3.** Let $M_0$ be a dimer configuration of $G^{D,'}$. Then, for every $u \in C'$, we have

$$
\begin{align*}
| \det K(u) | &= (k')^{\frac{|V^*|}{2}} \left( \prod_{w \in W} \mathrm{sc}(\theta_w) \right)^{\frac{1}{2}} \left( \prod_{e=wx \in M_0} [\mathrm{dn}(u_{\alpha_e}) \mathrm{dn}(u_{\beta_e})]^{\frac{1}{2}} \right) \det \Delta^{m,*}, \\
| \det K^\partial(u) | &= (k')^{\frac{|V|-1}{2}} \left( \prod_{w \in W} \mathrm{cs}(\theta_w) \right)^{\frac{1}{2}} \left( \prod_{e=wx \in M_0} [k' \mathrm{nd}(u_{\alpha_e}) \mathrm{nd}(u_{\beta_e})]^{\frac{1}{2}} \right) \det \Delta^{m,\partial}(u).
\end{align*}
$$
they are interesting in their own respect. In the specific case $k = 0$, results we obtain
can be found in [45, 71, 43]. We prove, see also Theorem 3.6:

**Theorem 1.1.**

* Infinite case. Let $u \in C$, then the $Z^u$-Dirac operator $K(u)$, the Z-massive Laplacian $\Delta^m$ and the dual $\Delta^{m,*}$ of [11] satisfy the following identity:

$$
\overline{K(u)}^t K(u) = k' \begin{pmatrix} \Delta^m & 0 \\ 0 & \Delta^{m,*} \end{pmatrix}.
$$

* Finite case. Let $u \in C'$, then the $Z^u$-Dirac operators $K(u)$, $K^\partial(u)$, the Z-massive Laplacian $\Delta^{m,\partial}(u)$ and the dual $\Delta^{m,*}$ satisfy the following identity:

$$
\overline{K^{\partial}(u)}^t K(u) = k' \begin{pmatrix} \Delta^{m,\partial}(u) & Q(u) \\ 0 & \Delta^{m,*} \end{pmatrix}.
$$

A function $F \in C^B$ is said to be $Z^u$-holomorphic if $K(u)F = 0$. As a consequence of Theorem 1.1, if $F$ is $Z^u$-holomorphic, then $F|_{V}$ is $Z$-massive harmonic on $G$ and $F|_{V^*}$ is $Z$-massive harmonic on $G^*$, thus explaining the part “Dirac” in “$Z^u$-Dirac operator”.

Note that in the finite case, the Z-massive Laplacian $\Delta^{m,\partial}(u)$ arising in the context of the Ising model might have negative masses and conductances along the boundary, see Equation (3.6) where it is explicitly defined. As a consequence, this operator is not necessarily positive definite.

In the infinite case, Theorem 1.1 yields the following relations for inverse operators,
see also Corollary 3.15; the statement in the finite case is given in Corollary 3.17.

**Corollary 1.2 (Infinite case).** For every $u \in C$, consider the operator $K(u)^{-1}$ mapping $C^W$ to $C^B$ whose coefficients are defined by, for every $v, f, w$ as in Figure 14,

$$
\begin{align*}
K(u)_{v,w}^{-1} &= e^{-i \frac{\alpha_f + \beta_f}{2}} (k')^{-1} \mathrm{sc}(\theta_f)^{\frac{1}{2}} \left( \left[ \mathrm{dn}(u_{\alpha_f}) \mathrm{dn}(u_{\beta_f}) \right]^{\frac{1}{2}} G_{v,v_2}^m - \left[ \mathrm{dn}(u_{\alpha_f+2K}) \mathrm{dn}(u_{\beta_f+2K}) \right]^{\frac{1}{2}} G_{v,v_1}^m \right) \\
K(u)_{f,w}^{-1} &= -i e^{-i \frac{\alpha_f + \beta_f}{2}} (k')^{-1} \mathrm{sc}(\theta_f^*)^{\frac{1}{2}} \left( \left[ \mathrm{dn}((u_{\beta_f})^*) \mathrm{dn}((u_{\alpha_f+2K})^*) \right]^{\frac{1}{2}} G_{f,f_2}^{m,*} + \right. \\
&\qquad \left. - \left[ \mathrm{dn}((u_{\beta_f-2K})^*) \mathrm{dn}((u_{\alpha_f})^*) \right]^{\frac{1}{2}} G_{f,f_1}^{m,*} \right),
\end{align*}
$$

where $G^m$ and $G^{m,*}$ are the Z-massive and dual Z-massive Green functions of [11]. Then $K(u)^{-1}$ is the unique inverse of the $Z^u$-Dirac operator $K(u)$ decreasing to zero at infinity.

This gives, in Theorem 3.19, an explicit local expression for a Gibbs measure of the
dimer model on the double graph $G^D$, where the locality property is inherited from that
of the $Z$-massive Green functions of [11]. Using the KPW-Temperley bijection [71, 49],
probabilities of pairs of dual directed spanning trees are computed using the Green
function of a massive, non-directed random walk. Apart from the locality property which
is specific, a similar result is obtained by Chhita [20] in the case of $\mathbb{Z}^2$ with a specific
choice of weights.

In Theorem 3.8 and Corollary 3.10, we restrict to the finite case and prove relations
on determinants; we show,

**Theorem 1.3.** Let $M_0$ be a dimer configuration of $G^{D,'}$. Then, for every $u \in C'$, we have

$$
\begin{align*}
| \det K(u) | &= (k')^{\frac{|V^*|}{2}} \left( \prod_{w \in W} \mathrm{sc}(\theta_w) \right)^{\frac{1}{2}} \left( \prod_{e=wx \in M_0} [\mathrm{dn}(u_{\alpha_e}) \mathrm{dn}(u_{\beta_e})]^{\frac{1}{2}} \right) \det \Delta^{m,*}, \\
| \det K^\partial(u) | &= (k')^{\frac{|V|-1}{2}} \left( \prod_{w \in W} \mathrm{cs}(\theta_w) \right)^{\frac{1}{2}} \left( \prod_{e=wx \in M_0} [k' \mathrm{nd}(u_{\alpha_e}) \mathrm{nd}(u_{\beta_e})]^{\frac{1}{2}} \right) \det \Delta^{m,\partial}(u).
\end{align*}
$$
As a consequence, the partition function of pairs of dual directed spanning trees is equal, up to a constant, to the partition function of rooted spanning forests. In the critical case, $k = 0$, this is an easy consequence of Temperley's bijection [71], but the correspondence does not extend when $k \neq 0$. The main tools of the proof are gauge equivalences on bipartite adjacency matrices and on adjacency matrices of digraphs, see also Appendix A.

The $Z^u$-Dirac operator is equivalent to a model of directed spanning trees. In Proposition 3.13, we prove that the latter is Z-invariant, thus explaining the part "$Z^u$" of the terminology "$Z^u$-Dirac operator".

**Section 4: Kasteleyn operator of the graph $G^Q$ and $Z^u$-Dirac operator** We consider the dimer model on the graph $G^Q$ arising from the Z-invariant XOR-Ising model, with Kasteleyn matrix $K^Q$. The main result of this section, and one of the main result of this paper, is Theorem 4.2 proving the following relations between the matrix $K^Q$ and the $Z^u$-Dirac operator. The matrices $S(u)$ and $T(u)$ are defined in Section 4.1 and the statement is as follows.

**Theorem 1.4.**

* Infinite case. Let $u \in C$, then the Kasteleyn matrix $K^Q$, the $Z^u$-Dirac operator $K(u)$ and the matrices $S(u)$, $T(u)$ are related by the following identity:

$$K^Q T(u) = S(u) K(u).$$

* Finite case. Let $u \in C'$, then the Kasteleyn matrix $\underline{K}^Q$, the $Z^u$-Dirac operator $\underline{K}^\theta(u)$ and the matrices $S(u)$, $T(u)$ are related by the following identity:

$$\underline{K}^Q T(u) = S(u) \underline{K}^\theta(u).$$

The matrix relations established in Theorem 1.4 have quite remarkable consequences. They allow to establish a one parameter family of relations between inverse operators which turn out to be very useful for some specific choices of $u$. Also, after some extra work, they permit to relate determinants in the finite case, thus partition functions. More precisely, in the infinite case, Theorem 1.4 yields the following relations on inverse matrices, see also Corollary 4.12; the statement in the finite case is the subject of Corollary 4.14.

**Corollary 1.5 (Infinite case).** For every $u \in C$, for every $w, v, f$, and every $w, b, b'$ as in Figure 20,

$$(K^Q)_{w,b}^{-1} \operatorname{cn}(u_{\beta_f}) - i (K^Q)_{w,b'}^{-1} \operatorname{sn}(u_{\beta_f}) \operatorname{dn}(u_{\alpha_f}) = \frac{e^{i \frac{\bar{\beta}_f - \bar{\beta}_i}{2}}}{\Lambda(u_{\alpha_f}, u_{\beta_f})} [\operatorname{cn}(u_{\beta_i}) K(u)_{v,w}^{-1} - i \operatorname{sn}(u_{\beta_i}) K(u)_{f,w}^{-1}],$$

$$\text{where } \Lambda(u_{\alpha}, u_{\beta}) = [\operatorname{sn} \theta \operatorname{cn} \theta \operatorname{nd}(u_{\alpha}) \operatorname{nd}(u_{\beta})]^{\frac{1}{2}}, \text{ and } \theta = u_{\alpha} - u_{\beta}.$$

Specifying the value of the parameter $u$ allows to express coefficients of the inverse Kasteleyn operator $(K^Q)^{-1}$ using the inverse $Z^u$-Dirac operator; combining this with Theorem 1.2 yields an expression using the Z-massive Green function of [11]. Note that since the dimer model on $G^Q$ is also related to a free-fermion, zero-field 6-vertex model [62, 76, 32], see also Section 2.3.3, the following corollary also relates this 6-vertex model Boltzmann measure to the inverse $Z^u$-Dirac operator and massive Green function. We obtain, see also Corollary 4.16 and Corollary 4.18 for the finite case,
As a consequence, the partition function of pairs of dual directed spanning trees is equal, up to a constant, to the partition function of rooted spanning forests. In the critical case, $k = 0$, this is an easy consequence of Temperley's bijection [71], but the correspondence does not extend when $k \neq 0$. The main tools of the proof are gauge equivalences on bipartite adjacency matrices and on adjacency matrices of digraphs, see also Appendix A.

The $Z^u$-Dirac operator is equivalent to a model of directed spanning trees. In Proposition 3.13, we prove that the latter is Z-invariant, thus explaining the part "$Z^u$" of the terminology "$Z^u$-Dirac operator".

**Section 4: Kasteleyn operator of the graph $G^Q$ and $Z^u$-Dirac operator** We consider the dimer model on the graph $G^Q$ arising from the Z-invariant XOR-Ising model, with Kasteleyn matrix $K^Q$. The main result of this section, and one of the main result of this paper, is Theorem 4.2 proving the following relations between the matrix $K^Q$ and the $Z^u$-Dirac operator. The matrices $S(u)$ and $T(u)$ are defined in Section 4.1 and the statement is as follows.

**Theorem 1.4.**

* Infinite case. Let $u \in C$, then the Kasteleyn matrix $K^Q$, the $Z^u$-Dirac operator $K(u)$ and the matrices $S(u)$, $T(u)$ are related by the following identity:

$$K^Q T(u) = S(u) K(u).$$

* Finite case. Let $u \in C'$, then the Kasteleyn matrix $\underline{K}^Q$, the $Z^u$-Dirac operator $\underline{K}^\theta(u)$ and the matrices $S(u)$, $T(u)$ are related by the following identity:

$$\underline{K}^Q T(u) = S(u) \underline{K}^\theta(u).$$

The matrix relations established in Theorem 1.4 have quite remarkable consequences. They allow to establish a one parameter family of relations between inverse operators which turn out to be very useful for some specific choices of $u$. Also, after some extra work, they permit to relate determinants in the finite case, thus partition functions. More precisely, in the infinite case, Theorem 1.4 yields the following relations on inverse matrices, see also Corollary 4.12; the statement in the finite case is the subject of Corollary 4.14.

**Corollary 1.5 (Infinite case).** For every $u \in C$, for every $w, v, f$, and every $w, b, b'$ as in Figure 20,

$$(K^Q)_{w,b}^{-1} \operatorname{cn}(u_{\beta_f}) - i (K^Q)_{w,b'}^{-1} \operatorname{sn}(u_{\beta_f}) \operatorname{dn}(u_{\alpha_f}) = \frac{e^{i \frac{\bar{\beta}_f - \bar{\beta}_i}{2}}}{\Lambda(u_{\alpha_f}, u_{\beta_f})} \left[ \operatorname{cn}(u_{\beta_i}) K(u)_{v,w}^{-1} - i \operatorname{sn}(u_{\beta_i}) K(u)_{f,w}^{-1} \right],$$

$$\text{where } \Lambda(u_\alpha, u_\beta) = [\operatorname{sn}\theta \operatorname{cn}\theta \operatorname{nd}(u_\alpha) \operatorname{nd}(u_\beta)]^{\frac{1}{2}}, \text{ and } \theta = u_\alpha - u_\beta.$$

Specifying the value of the parameter $u$ allows to express coefficients of the inverse Kasteleyn operator $(K^Q)^{-1}$ using the inverse $Z^u$-Dirac operator; combining this with Theorem 1.2 yields an expression using the Z-massive Green function of [11]. Note that since the dimer model on $G^Q$ is also related to a free-fermion, zero-field 6-vertex model [62, 76, 32], see also Section 2.3.3, the following corollary also relates this 6-vertex model Boltzmann measure to the inverse $Z^u$-Dirac operator and massive Green function. We obtain, see also Corollary 4.16 and Corollary 4.18 for the finite case,
**Corollary 1.6 (Infinite case).** For every w, b of $G^Q$ as in Figure 22,

$$ (K^Q)_{w,b}^{-1} = \frac{e^{\frac{\beta_f - \beta_i}{2}}}{[\operatorname{cn}(\theta_f) \operatorname{sn}(\theta_f) \operatorname{nd}(\theta_f)]^{\frac{1}{2}}} \left( \operatorname{cn}\left(\frac{\beta_f - \beta_i}{2}\right) K(\beta_f)_{v,w}^{-1} - i \operatorname{sn}\left(\frac{\beta_f - \beta_i}{2}\right) K(\beta_f)_{f,w}^{-1} \right) $$

$$ (K^Q)_{w,b}^{-1} = e^{-i \frac{\beta_i + \alpha_f}{2}} (k')^{-1} \left( \frac{\operatorname{cn}\left(\frac{\beta_f - \beta_i}{2}\right)}{\operatorname{cn}(\theta_f)} [\operatorname{dn}(\theta_f) G_{v,v_2}^m - k' G_{v,v_1}^m] + \right. \\ \left. - \frac{\operatorname{sn}\left(\frac{\beta_f - \beta_i}{2}\right)}{\operatorname{sn}(\theta_f)} [\operatorname{dn}(\theta_f) G_{f,f_2}^{m,*} - G_{f,f_1}^{m,*}] \right). $$

This proves, in an alternative way, an explicit local expression for a Gibbs measure of the dimer model on the graph $G^Q$ [12], where the locality property is seen as directly inherited from that of the Z-massive Green function.

Again, using Theorem 1.4 and additional combinatorial arguments, we prove in Theorem 4.3 that the determinants of $K^Q$ and of the $Z^u$-Dirac operator are equal, up to an explicit constant. By [32], the determinant of $K^Q$ is equal up to a constant, to the squared partition function of the Ising model. Combining this with Theorem 1.3 gives, see also Corollary 4.5 for the explicit value of $C(u)$,

**Theorem 1.7.** For every $u \in C''$,

$$ [Z_{\text{Ising}}^{+}(G, J)]^2 = C(u) |\det \Delta^{m,\theta}(u)|, $$

where $C''$ is a subset of $C$ defined in Equation (3.19).

When $k=0$, we essentially recover the result of [30]; see Section 6.1. The above proves that the squared partition function of the Z-invariant Ising model is equal, up to a multiplicative constant, to that of rooted directed spanning forests. This relates on the level of partition functions two classical, apparently very different, models of statistical mechanics.

**Section 5: Dimer model on the Fisher graph $G^F$ and the Kasteleyn matrix $K^Q$**

Consider the dimer model on $G^F$ with Kasteleyn matrix $K^F$ arising from an Ising model with coupling constants $J$, not necessarily Z-invariant, and the corresponding dimer model on the bipartite graph $G^Q$, with (real) Kasteleyn matrix $\tilde{K}^Q$. Note that the entries of $\tilde{K}^Q$ are fully determined by the coupling constants $J$. Building on the work of Dubédat [32] and proving additional matrix relations, we express coefficients of the inverse operator $(K^F)^{-1}$ using coefficients of the inverse operator $(\tilde{K}^Q)^{-1}$. Partitioning vertices of $G^F$ as $A \cup B$ as in [32], we obtain, see also Theorem 5.4,

**Theorem 1.8 (Finite and infinite cases).** Using the notation of Figure 26, there are four cases to consider:

1. For every $a \in A$ and every $b \in B$ such that, when the graph $G^F$ is moreover finite, $b$ is not a boundary vertex:

$$ (K^F)_{a,b}^{-1} = \frac{1}{1 + e^{-4J_e^f}} [(\tilde{K}^Q)_{w,b}^{-1} + (\tilde{K}^Q)_{w,b'}^{-1} \varepsilon_{b',b} e^{-2J_e^f}]. $$

2. When the graph $G^F$ is finite, for every $a \in A$ and every boundary vertex $b$ of $B$, we have

$$ (K^F)_{a,b}^{-1} = (\tilde{K}^Q)_{w,b}^{-1}. $$

3. For every $a, a \in A$,

$$ (K^F)_{a,a}^{-1} = -\frac{1}{2} (\tilde{K}^Q)_{w,b}^{-1} \varepsilon_{b,a} + \kappa_{a,a}, $$

where $\kappa_{a,a} = 0$ if $a$ and $a$ do not belong to the same decoration, and to $\pm \frac{1}{4}$ if they do.
**Corollary 1.6 (Infinite case).** For every w, b of $G^Q$ as in Figure 22,

$$ (K^Q)_{w,b}^{-1} = \frac{e^{\frac{\beta_f - \beta_i}{2}}}{[\operatorname{cn}(\theta_f) \operatorname{sn}(\theta_f) \operatorname{nd}(\theta_f)]^{\frac{1}{2}}} \left( \operatorname{cn}\left(\frac{\beta_f - \beta_i}{2}\right) K(\beta_f)_{v,w}^{-1} - i \operatorname{sn}\left(\frac{\beta_f - \beta_i}{2}\right) K(\beta_f)_{f,w}^{-1} \right) \\ (K^Q)_{w,b}^{-1} = e^{-i \frac{\beta_i + \alpha_f}{2}} (k')^{-1} \left( \frac{\operatorname{cn}\left(\frac{\beta_f - \beta_i}{2}\right)}{\operatorname{cn}(\theta_f)} [\operatorname{dn}(\theta_f) G_{v,v_2}^m - k' G_{v,v_1}^m] + \right. \\ \left. - \frac{\operatorname{sn}\left(\frac{\beta_f - \beta_i}{2}\right)}{\operatorname{sn}(\theta_f)} [\operatorname{dn}(\theta_f) G_{f,f_2}^{m,*} - G_{f,f_1}^{m,*}] \right). $$

This proves, in an alternative way, an explicit local expression for a Gibbs measure of the dimer model on the graph $G^Q$ [12], where the locality property is seen as directly inherited from that of the Z-massive Green function.

Again, using Theorem 1.4 and additional combinatorial arguments, we prove in Theorem 4.3 that the determinants of $K^Q$ and of the $Z^u$-Dirac operator are equal, up to an explicit constant. By [32], the determinant of $K^Q$ is equal up to a constant, to the squared partition function of the Ising model. Combining this with Theorem 1.3 gives, see also Corollary 4.5 for the explicit value of $C(u)$,

**Theorem 1.7.** For every $u \in C''$,

$$ [Z_{\text{Ising}}^{+}(G, J)]^2 = C(u) |\det \Delta^{m,\theta}(u)|, $$

where $C''$ is a subset of $C$ defined in Equation (3.19).

When $k=0$, we essentially recover the result of [30]; see Section 6.1. The above proves that the squared partition function of the Z-invariant Ising model is equal, up to a multiplicative constant, to that of rooted directed spanning forests. This relates on the level of partition functions two classical, apparently very different, models of statistical mechanics.

**Section 5: Dimer model on the Fisher graph $G^F$ and the Kasteleyn matrix $K^Q$**

Consider the dimer model on $G^F$ with Kasteleyn matrix $K^F$ arising from an Ising model with coupling constants J, not necessarily Z-invariant, and the corresponding dimer model on the bipartite graph $G^Q$, with (real) Kasteleyn matrix $\tilde{K}^Q$. Note that the entries of $\tilde{K}^Q$ are fully determined by the coupling constants J. Building on the work of Dubédat [32] and proving additional matrix relations, we express coefficients of the inverse operator $(K^F)^{-1}$ using coefficients of the inverse operator $(\tilde{K}^Q)^{-1}$. Partitioning vertices of $G^F$ as $A \cup B$ as in [32], we obtain, see also Theorem 5.4,

**Theorem 1.8 (Finite and infinite cases).** Using the notation of Figure 26, there are four cases to consider:

1. For every $a \in A$ and every $b \in B$ such that, when the graph $G^F$ is moreover finite, $b$ is not a boundary vertex:

$$ (K^F)_{a,b}^{-1} = \frac{1}{1 + e^{-4J_e^f}} [(\tilde{K}^Q)_{w,b}^{-1} + (\tilde{K}^Q)_{w,b'}^{-1} \varepsilon_{b',b} e^{-2J_e^f}]. $$

2. When the graph $G^F$ is finite, for every $a \in A$ and every boundary vertex $b$ of $B$, we have

$$ (K^F)_{a,b}^{-1} = (\tilde{K}^Q)_{w,b}^{-1}. $$

3. For every $a, a \in A$,

$$ (K^F)_{a,a}^{-1} = -\frac{1}{2} (\tilde{K}^Q)_{w,b}^{-1} \varepsilon_{b,a} + \kappa_{a,a}, $$

where $\kappa_{a,a} = 0$ if $a$ and $a$ do not belong to the same decoration, and to $\pm \frac{1}{4}$ if they do.
4. For every $\underline{b}, b \in B$,

$$ (K^F)_{\underline{b},b}^{-1} = -\epsilon_{b,a_1} (K^F)_{a_1,b}^{-1} + \epsilon_{b,a_2} (K^F)_{a_2,b}^{-1}, $$

where $(K^F)_{a_1,b}^{-1}, (K^F)_{a_2,b}^{-1}$ are given by Case 1.

As a consequence, the Boltzmann/Gibbs measures of the dimer model on the non-bipartite graph $G^F$ can be computed using the inverse Kasteleyn operator of the bipartite graph $G^Q$. Note that in the finite case, we do not need positivity of the coupling constants $J$. When $J < 0$, the dimer model on the Fisher graph $G^F$ has positive weights 1 and $e^{-2J_e}$ on edges, and is related to a bipartite dimer model with some negative weights, see also Remark 5.5.

**Remark 1.9.**

* As mentioned in [32], bosonization identities somehow prove the existence of such linear relations, but working them out requires more work, which is the subject of the above theorem.

* As pointed out by one of the referees whom I thank, these relations can also be derived from the paper [16]. Following his/her suggestion, let us explain how. Using the notation of [16] we have that $K^F$ is $\hat{F} = -iU_F^*F U_F$. The definition of the matrix $F$ implies, see [16, Section 3.1],

$$ F^{-1} = \begin{pmatrix} I & 0 \\ -JB^* & J \end{pmatrix} \begin{pmatrix} C^{-1} & 0 \\ 0 & -J \end{pmatrix} \begin{pmatrix} I & -BJ \\ 0 & J \end{pmatrix}. $$

Then, [16, Lemma 3.4] expresses coefficients of $C^{-1}$ using those of $D^{-1}$, where $\hat{D} = -iU_C^*D U_C$ is essentially the operator $\hat{K}^Q$. Putting this together gives Theorem 1.8.

Next we restrict to the Z-invariant case. The coefficient $(K^F)_{a,a}^{-1}$ is equal, up to an additive constant, to the coefficient $(K^Q)_{w,b}^{-1}$, and is thus expressed using the inverse $Z^u$-Dirac operator using Corollary 1.6 in the infinite case, and Corollary 4.18 in the finite case. The same holds for $(K^F)_{a,b}^{-1}$ when $b$ is a boundary vertex. The coefficient $(K^F)_{b,b}^{-1}$ is a simple linear combination of two coefficients $(K^F)_{a_1,b'}$ and $(K^F)_{a_2,b'}$, so we are left with expressing the coefficient $(K^F)_{a,b}^{-1}$. Choosing a specific value of $u$ in Corollary 1.5, and using Corollary 1.2 gives, see also Corollary 5.8 for the finite case,

**Corollary 1.10 (Infinite case).** Let $u^f = \frac{\alpha_f + \beta_f}{2} + K$. Then,

$$ \begin{aligned} (K^F)_{a,b}^{-1} &= q_{b,w} e^{i \frac{\beta_f - \beta_i}{2} (k') \frac{1}{2} \frac{\operatorname{cn}\left(\frac{K-\theta_f}{2}\right) (1+(k')^{-1} \operatorname{dn}(\theta_f))}{2[\operatorname{cn}(\theta_f) \operatorname{sn}(\theta_f)]^{\frac{1}{2}}}} \left( \operatorname{cn}(\mathbf{u}_{\beta_i}^f) K(\mathbf{u}^f)_{\underline{v},w}^{-1} - i \operatorname{sn}(\mathbf{u}_{\beta_i}^f) K(\mathbf{u}^f)_{\underline{f},w}^{-1} \right) \\ &= q_{b,w} e^{-i \frac{\alpha_f + \beta_i}{2} \frac{1}{2} \frac{\operatorname{cn}\left(\frac{K-\theta_f}{2}\right) (1+(k')^{-1} \operatorname{dn}(\theta_f))}{2}} \times \\ & \qquad \left( \frac{\operatorname{cn}(\mathbf{u}_{\beta_i}^f)}{\operatorname{cn}(\theta_f)} (G_{\underline{v},v_2}^m - G_{\underline{v},v_1}^m) - \frac{\operatorname{sn}(\mathbf{u}_{\beta_i}^f)}{\operatorname{sn}(\theta_f)} \left( \operatorname{nd}\left(\frac{K-\theta_f}{2}\right) G_{\underline{f},f_2}^{m,*} - \operatorname{nd}\left(\frac{K+\theta_f}{2}\right) G_{\underline{f},f_1}^{m,*} \right) \right). \end{aligned} $$

*Connection to previously known results.* Apart from allowing to compute the contour Ising Boltzmann/Gibbs measures, coefficients of the inverse Kasteleyn matrix $(K^F)^{-1}$ are important observables of the Ising model: $(K^F)_{\underline{b},b}^{-1}$ is related to the spin-Ising observable [69, 18], see for example [16]. Dubédat [32] proves that $(K^F)_{a,a}^{-1}$ is the fermionic spinor observable of [39] and, referring to Nienhuis-Knops [63], mentions that it is also the FK-Ising observable of [67, 18] (up to normalization).

As a consequence, in the specific case $k=0$ (the critical case), Theorem 1.8, Corollary 1.10 and Corollary 1.6 are deeply related to the discrete part of [18] proving that these observables are holomorphic, and integrating the square to obtain close to harmonic functions, see also Section 6.1. Our results have two important features: they
4. For every $\underline{b}, b \in B$,

$$ (K^F)_{\underline{b},b}^{-1} = -\varepsilon_{b,a_1} (K^F)_{a_1,b}^{-1} + \varepsilon_{b,a_2} (K^F)_{a_2,b}^{-1}, $$

where $(K^F)_{a_1,b}^{-1}, (K^F)_{a_2,b}^{-1}$ are given by Case 1.

As a consequence, the Boltzmann/Gibbs measures of the dimer model on the non-bipartite graph $G^F$ can be computed using the inverse Kasteleyn operator of the bipartite graph $G^Q$. Note that in the finite case, we do not need positivity of the coupling constants $J$. When $J < 0$, the dimer model on the Fisher graph $G^F$ has positive weights 1 and $e^{-2J_e}$ on edges, and is related to a bipartite dimer model with some negative weights, see also Remark 5.5.

**Remark 1.9.**

* As mentioned in [32], bosonization identities somehow prove the existence of such linear relations, but working them out requires more work, which is the subject of the above theorem.

* As pointed out by one of the referees whom I thank, these relations can also be derived from the paper [16]. Following his/her suggestion, let us explain how. Using the notation of [16] we have that $K^F$ is $\hat{F} = -iU_F^*F U_F$. The definition of the matrix $F$ implies, see [16, Section 3.1],

$$ F^{-1} = \begin{pmatrix} I & 0 \\ -JB^* & J \end{pmatrix} \begin{pmatrix} C^{-1} & 0 \\ 0 & -J \end{pmatrix} \begin{pmatrix} I & -BJ \\ 0 & J \end{pmatrix}. $$

Then, [16, Lemma 3.4] expresses coefficients of $C^{-1}$ using those of $D^{-1}$, where $\hat{D} = -iU_C^*D U_C$ is essentially the operator $\hat{K}^Q$. Putting this together gives Theorem 1.8.

Next we restrict to the Z-invariant case. The coefficient $(K^F)_{\underline{a},a}^{-1}$ is equal, up to an additive constant, to the coefficient $(K^Q)_{\underline{w},b}^{-1}$, and is thus expressed using the inverse $Z^u$-Dirac operator using Corollary 1.6 in the infinite case, and Corollary 4.18 in the finite case. The same holds for $(K^F)_{\underline{a},b}^{-1}$ when $b$ is a boundary vertex. The coefficient $(K^F)_{\underline{b},b}^{-1}$ is a simple linear combination of two coefficients $(K^F)_{\underline{a}_1,b'}$ and $(K^F)_{\underline{a}_2,b'}$, so we are left with expressing the coefficient $(K^F)_{\underline{a},b}^{-1}$. Choosing a specific value of $u$ in Corollary 1.5, and using Corollary 1.2 gives, see also Corollary 5.8 for the finite case,

**Corollary 1.10 (Infinite case).** Let $u^f = \frac{\alpha_f + \beta_f}{2} + K$. Then,

$$ (K^F)_{\underline{a},b}^{-1} = q_{b,w} e^{i \frac{\beta_f - \beta_i}{2} (k') \frac{1}{2} \frac{\operatorname{cn}\left(\frac{K-\theta_f}{2}\right) (1+(k')^{-1} \operatorname{dn}(\theta_f))}{2[\operatorname{cn}(\theta_f) \operatorname{sn}(\theta_f)]^{\frac{1}{2}}}} \left( \operatorname{cn}(\mathbf{u}_{\beta_i}^f) K(\mathbf{u}^f)_{\underline{v},w}^{-1} - i \operatorname{sn}(\mathbf{u}_{\beta_i}^f) K(\mathbf{u}^f)_{\underline{f},w}^{-1} \right) \\ = q_{b,w} e^{-i \frac{\alpha_f + \beta_i}{2} \frac{1}{2} \frac{\operatorname{cn}\left(\frac{K-\theta_f}{2}\right) (1+(k')^{-1} \operatorname{dn}(\theta_f))}{2}} \times \\ \times \left( \frac{\operatorname{cn}(\mathbf{u}_{\beta_i}^f)}{\operatorname{cn}(\theta_f)} (G_{\underline{v},v_2}^m - G_{\underline{v},v_1}^m) - \frac{\operatorname{sn}(\mathbf{u}_{\beta_i}^f)}{\operatorname{sn}(\theta_f)} \left( \operatorname{nd}\left(\frac{K-\theta_f}{2}\right) G_{\underline{f},f_2}^{m,*} - \operatorname{nd}\left(\frac{K+\theta_f}{2}\right) G_{\underline{f},f_1}^{m,*} \right) \right). $$

*Connection to previously known results.* Apart from allowing to compute the contour Ising Boltzmann/Gibbs measures, coefficients of the inverse Kasteleyn matrix $(K^F)^{-1}$ are important observables of the Ising model: $(K^F)_{\underline{b},b}^{-1}$ is related to the spin-Ising observable [69, 18], see for example [16]. Dubédat [32] proves that $(K^F)_{\underline{a},a}^{-1}$ is the fermionic spinor observable of [39] and, referring to Nienhuis-Knops [63], mentions that it is also the FK-Ising observable of [67, 18] (up to normalization).

As a consequence, in the specific case $k=0$ (the critical case), Theorem 1.8, Corollary 1.10 and Corollary 1.6 are deeply related to the discrete part of [18] proving that these observables are holomorphic, and integrating the square to obtain close to harmonic functions, see also Section 6.1. Our results have two important features: they
prove that in the infinite and finite cases, these observables have an exact explicit expression involving Green functions, and also that these expressions not only hold at criticality but in the full Z-invariant regime. Note that since our results hold for all $k$, by taking $k$ tending to 0 at an appropriate speed, one may obtain interesting features of the near critical Ising model.

The paper [56] gives a *non-backtracking random walk* representation of the inverse Kac-Ward operator, the latter being connected to the inverse Kasteleyn operator. In this paper, we give a (*massive*) random walk representation of the inverse Kasteleyn operator where, in the finite case, this random walk has some vortices along the boundary. In the critical case, and for one choice of $u$ (namely $u = i\infty$), part of the relation of Theorem 1.4 was obtained in [21]. Let us end this introduction with a comment on the paper [7] based on an observation by Messikh [60] about the occurrence of large deviation estimates of a massive random walk in the correlation length of the super-critical Ising model on $Z^2$. The proof consists in showing that, in the super-critical regime, spin correlations are approximated by the FK-Ising observable, and then analyzing the latter. By Theorem 1.8 and Corollary 1.6, the latter a directly related to the massive Green function, thus explaining the occurrence of the massive random walk, see also Section 6.2 specifying our results to the case where $G = Z^2$.

A word on the parameter $u$. As noted above, the parameter $u$ of the $Z^u$-Dirac operator disappears at criticality, and so some motivation on where it comes from may be useful. When writing this paper, the parameter $u$ naturally arose from elliptic trigonometric identities; it provided an additional degree of freedom which turned out to be very useful: for example in Corollary 1.6, choosing specific values of $u$ enables us to express the operator $(K^Q)^{-1}$ as a function of the operator $K(u)^{-1}$. While in the revision process of this paper, the paper [8] came out and allows to give more insight in the periodic case. By [12], we know that spectral curves of $Z$-invariant Ising models with elliptic weights are in correspondence with genus 1 Harnack curves with central symmetry. From [8] we know that genus 1 Harnack curves arise from bipartite dimer models with Fock's elliptic weights [36]. Indeed, it turns out that the $Z^u$-Dirac operator is gauge equivalent to Fock's elliptic Kasteleyn operator on the double graph $G^D$ as shown in [8, Section 8.2]. Now Fock's elliptic operator has a natural additional parameter (denoted $t$ in [8] instead of $u$ here) which is in bijection with the hole of the amoeba of the spectral curve. In the critical case, corresponding to genus 0 Harnack curves, there is no hole in the amoeba and thus no extra parameter. Note that, although the construction is not explicit, Goncharov and Kenyon [37], already proved a correspondence between bipartite dimer models and Harnack curves with marked points on the ovals. One last point on this topic: as observed after the statement of Theorem 1.1, every $Z^u$-holomorphic function is $Z$-massive harmonic; this in particular raises the question of thoroughly studying a massive version of discrete complex analysis in the spirit of the papers [59, 17].

## 2 Preliminaries

This section contains all the preliminaries required for this paper. We give the definitions of the Ising model, the dimer model per se, the dimer model on decorated graphs arising from the Ising model, from the XOR-Ising model and from pairs of dual directed spanning trees; we also define the rooted directed spanning forests model. We end with isoradial graphs, $Z$-invariance and the $Z$-invariant versions of the above models. In the whole of the paper, planar embedded graphs are supposed to be simply connected.
prove that in the infinite and finite cases, these observables have an exact explicit expression involving Green functions, and also that these expressions not only hold at criticality but in the full Z-invariant regime. Note that since our results hold for all $k$, by taking $k$ tending to 0 at an appropriate speed, one may obtain interesting features of the near critical Ising model.

The paper [56] gives a *non-backtracking random walk* representation of the inverse Kac-Ward operator, the latter being connected to the inverse Kasteleyn operator. In this paper, we give a (*massive*) random walk representation of the inverse Kasteleyn operator where, in the finite case, this random walk has some vortices along the boundary. In the critical case, and for one choice of $u$ (namely $u = i\infty$), part of the relation of Theorem 1.4 was obtained in [21]. Let us end this introduction with a comment on the paper [7] based on an observation by Messikh [60] about the occurrence of large deviation estimates of a massive random walk in the correlation length of the super-critical Ising model on $Z^2$. The proof consists in showing that, in the super-critical regime, spin correlations are approximated by the FK-Ising observable, and then analyzing the latter. By Theorem 1.8 and Corollary 1.6, the latter a directly related to the massive Green function, thus explaining the occurrence of the massive random walk, see also Section 6.2 specifying our results to the case where $G = Z^2$.

A word on the parameter $u$. As noted above, the parameter $u$ of the $Z^u$-Dirac operator disappears at criticality, and so some motivation on where it comes from may be useful. When writing this paper, the parameter $u$ naturally arose from elliptic trigonometric identities; it provided an additional degree of freedom which turned out to be very useful: for example in Corollary 1.6, choosing specific values of $u$ enables us to express the operator $(K^Q)^{-1}$ as a function of the operator $K(u)^{-1}$. While in the revision process of this paper, the paper [8] came out and allows to give more insight in the periodic case. By [12], we know that spectral curves of $Z$-invariant Ising models with elliptic weights are in correspondence with genus 1 Harnack curves with central symmetry. From [8] we know that genus 1 Harnack curves arise from bipartite dimer models with Fock's elliptic weights [36]. Indeed, it turns out that the $Z^u$-Dirac operator is gauge equivalent to Fock's elliptic Kasteleyn operator on the double graph $G^D$ as shown in [8, Section 8.2]. Now Fock's elliptic operator has a natural additional parameter (denoted $t$ in [8] instead of $u$ here) which is in bijection with the hole of the amoeba of the spectral curve. In the critical case, corresponding to genus 0 Harnack curves, there is no hole in the amoeba and thus no extra parameter. Note that, although the construction is not explicit, Goncharov and Kenyon [37], already proved a correspondence between bipartite dimer models and Harnack curves with marked points on the ovals. One last point on this topic: as observed after the statement of Theorem 1.1, every $Z^u$-holomorphic function is $Z$-massive harmonic; this in particular raises the question of thoroughly studying a massive version of discrete complex analysis in the spirit of the papers [59, 17].

## 2 Preliminaries

This section contains all the preliminaries required for this paper. We give the definitions of the Ising model, the dimer model per se, the dimer model on decorated graphs arising from the Ising model, from the XOR-Ising model and from pairs of dual directed spanning trees; we also define the rooted directed spanning forests model. We end with isoradial graphs, $Z$-invariance and the $Z$-invariant versions of the above models. In the whole of the paper, planar embedded graphs are supposed to be simply connected.
## 2.1 The Ising model

Consider a finite, planar, simple graph $G = (V, E)$. Suppose that edges of $G$ are assigned positive *coupling constants* $J = (J_e)_{e \in E}$. The *Ising model on G with free boundary conditions* is defined as follows. A *spin configuration* is a function on vertices of $G$ taking values in $\{-1, 1\}$. The probability on the set of spin configurations $\{-1, 1\}^V$ is given by the *Ising Boltzmann measure* $P_{\text{Ising}}$, defined by:

$$ \forall \sigma \in \{-1, 1\}^V, \quad P_{\text{Ising}}(\sigma) = \frac{1}{Z_{\text{Ising}}(G, J)} \exp\left(\sum_{e=vv' \in E} J_e \sigma_v \sigma_{v'}\right), $$

where $Z_{\text{Ising}}(G, J) = \sum_{\sigma \in \{-1, 1\}^V} \exp\left(\sum_{e=vv' \in E} J_e \sigma_v \sigma_{v'}\right)$ is the normalizing constant known as the *Ising partition function*.

From now on, we suppose that the planar graph $G$ is embedded. *Boundary vertices* of $G$ are vertices on the boundary of the unbounded face of $G$. The Ising model with + *boundary conditions* has the additional restriction that boundary vertices have +1 spin. Denote by $P_{\text{Ising}}^{+}$ and $Z_{\text{Ising}}^{+}(G, J)$ the corresponding Boltzmann measure and partition function<fnref n="1" />.

Denote by $\bar{G}^* = (\bar{V}^*, \bar{E}^*)$ the dual graph of $G$, and by o the vertex of $\bar{G}^*$ corresponding to the unbounded face of $G$. Consider also the *restricted dual graph* $G^* = (V^*, E^*)$ obtained from $\bar{G}^*$ by removing the vertex o and all of its incident edges. A *polygon configuration* of $G^*$ is a subset of edges such that every vertex has even degree; let $\mathcal{P}(G^*)$ denote the set of polygon configurations of $G^*$. Then, the *low temperature expansion* (LTE) of the Ising partition function with + boundary conditions is [51, 52]:

$$ Z_{\text{Ising}}^{+}(G, J) = \left(\prod_{e \in E} e^{J_e}\right) \sum_{P \in \mathcal{P}(G^*)} \prod_{e^* \in P} e^{-2J_e^*}. \qquad (2.1) $$

Polygon configurations of this expansion separate clusters of ±1 spins of the Ising model.

In this paper we consider the case where the graph is finite or infinite. The definition of the Boltzmann measure does not hold in the infinite case but extends naturally, and this will be clarified as we go along.

In the finite case, we consider the Ising model with + boundary conditions. It will be crucial to use the *boundary trick* of Chelkak and Smirnov [18] consisting in adding one extra vertex with +1 spin on every boundary edge of the graph. This has no effect on the Ising model, but the graph gains geometric freedom along the boundary, which will be key to handling boundary terms in Theorem 4.2. In order not to introduce too many graphs and confuse the reader, from now on we let $G$ be the graph we started from with the extra vertex on every boundary edge, then $\bar{G}^*$ is its dual graph and $G^*$ its restricted dual. Figure 1 provides an example of: a graph $G$, its restricted dual $G^*$, a spin configuration with + boundary conditions and the corresponding low temperature polygon configuration of $G^*$.

In the infinite case, we suppose that the embedded graph together with its faces cover the whole plane. So as not to have too many notation, and since it will be clear from the setting, we also denote by $G$ the infinite graph; the dual graph is denoted $G^*$.

## 2.2 The dimer model

Throughout the paper, we use the dimer model defined on three decorated versions of the graph $G$. Prior to defining these decorated graphs, we recall the definition of the

<footnote n="1">¹Note that the Ising model with + boundary conditions on the graph $G$ can be seen as the Ising model with free boundary conditions on the graph $G'$ obtained from $G$ by merging all boundary edges and vertices into a single vertex.</footnote>
## 2.1 The Ising model

Consider a finite, planar, simple graph $G = (V, E)$. Suppose that edges of $G$ are assigned positive coupling constants $J = (J_e)_{e \in E}$. The Ising model on $G$ with free boundary conditions is defined as follows. A spin configuration is a function on vertices of $G$ taking values in $\{-1, 1\}$. The probability on the set of spin configurations $\{-1, 1\}^V$ is given by the Ising Boltzmann measure $P_{\text{Ising}}$, defined by:

$$ \forall \sigma \in \{-1, 1\}^V, \quad P_{\text{Ising}}(\sigma) = \frac{1}{Z_{\text{Ising}}(G, J)} \exp\left(\sum_{e=vv' \in E} J_e \sigma_v \sigma_{v'}\right), $$

where $Z_{\text{Ising}}(G, J) = \sum_{\sigma \in \{-1, 1\}^V} \exp\left(\sum_{e=vv' \in E} J_e \sigma_v \sigma_{v'}\right)$ is the normalizing constant known as the *Ising partition function*.

From now on, we suppose that the planar graph $G$ is embedded. *Boundary vertices* of $G$ are vertices on the boundary of the unbounded face of $G$. The Ising model with + *boundary conditions* has the additional restriction that boundary vertices have +1 spin. Denote by $P_{\text{Ising}}^+$ and $Z_{\text{Ising}}^+(G, J)$ the corresponding Boltzmann measure and partition function<fnref n="1" />.

Denote by $\bar{G}^* = (\bar{V}^*, \bar{E}^*)$ the dual graph of $G$, and by o the vertex of $\bar{G}^*$ corresponding to the unbounded face of $G$. Consider also the *restricted dual graph* $G^* = (V^*, E^*)$ obtained from $\bar{G}^*$ by removing the vertex o and all of its incident edges. A *polygon configuration* of $G^*$ is a subset of edges such that every vertex has even degree; let $\mathcal{P}(G^*)$ denote the set of polygon configurations of $G^*$. Then, the *low temperature expansion* (LTE) of the Ising partition function with + boundary conditions is [51, 52]:

$$ Z_{\text{Ising}}^{+}(G, J) = \left(\prod_{e \in E} e^{J_e}\right) \sum_{P \in \mathcal{P}(G^*)} \prod_{e^* \in P} e^{-2J_e^*}. \qquad (2.1) $$

Polygon configurations of this expansion separate clusters of ±1 spins of the Ising model.

In this paper we consider the case where the graph is finite or infinite. The definition of the Boltzmann measure does not hold in the infinite case but extends naturally, and this will be clarified as we go along.

In the finite case, we consider the Ising model with + boundary conditions. It will be crucial to use the *boundary trick* of Chelkak and Smirnov [18] consisting in adding one extra vertex with +1 spin on every boundary edge of the graph. This has no effect on the Ising model, but the graph gains geometric freedom along the boundary, which will be key to handling boundary terms in Theorem 4.2. In order not to introduce too many graphs and confuse the reader, from now on we let $G$ be the graph we started from with the extra vertex on every boundary edge, then $\bar{G}^*$ is its dual graph and $G^*$ its restricted dual. Figure 1 provides an example of: a graph $G$, its restricted dual $G^*$, a spin configuration with + boundary conditions and the corresponding low temperature polygon configuration of $G^*$.

In the infinite case, we suppose that the embedded graph together with its faces cover the whole plane. So as not to have too many notation, and since it will be clear from the setting, we also denote by $G$ the infinite graph; the dual graph is denoted $G^*$.

## 2.2 The dimer model

Throughout the paper, we use the dimer model defined on three decorated versions of the graph $G$. Prior to defining these decorated graphs, we recall the definition of the

<footnote n="1">¹Note that the Ising model with + boundary conditions on the graph $G$ can be seen as the Ising model with free boundary conditions on the graph $G'$ obtained from $G$ by merging all boundary edges and vertices into a single vertex.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1: An example of a graph G (black), of its restricted dual G* (grey), of a spin configuration with + boundary conditions on G and its corresponding polygon configuration on G* (red). Vertices of G are pictured as bullets – black ones represent vertices of the original graph and grey ones are the additional vertices on boundary edges – vertices of G* are pictured as diamonds.</figcaption></figure>

dimer model per se, as well as that of the Kasteleyn matrix. We also recall the founding results that we will use.

Consider a planar, simple, embedded graph $G = (V, E)$. A dimer configuration of $G$, also known as a perfect matching, is a subset of edges such that every vertex is incident to exactly one edge of this subset. Denote by $M(G)$ the set of dimer configurations of the graph $G$. Suppose that a positive weight function $\nu$ is assigned to edges of $G$.

### 2.2.1 Finite case

Suppose that the graph $G$ is finite, and that $|V|$ is even. Then, the probability of occurrence of a dimer configuration, chosen with respect to the dimer Boltzmann measure $P_{\text{dimer}}$, is given by

$$ \forall M \in \mathcal{M}(G), \quad P_{\text{dimer}}(M) = \frac{\prod_{e \in M} \nu_e}{Z_{\text{dimer}}(G, \nu)}, $$

where $Z_{\text{dimer}}(G, \nu) = \sum_{M \in \mathcal{M}(G)} \prod_{e \in M} \nu_e$ is the normalizing constant known as the dimer partition function.

The main tool used to study the dimer model is the Kasteleyn matrix [40, 42, 72], it is defined as follows. A face-cycle is a cycle of $G$ bounding a bounded face of the graph. A Kasteleyn orientation is an orientation of the edges such that every face-cycle is clockwise odd, meaning that when traveling clockwise around a face-cycle, the number of co-oriented edges is odd. By the results of [42], a Kasteleyn orientation always exists for planar graphs. A Kasteleyn matrix, denoted by $K$, is a weighted, directed, adjacency matrix of the graph $G$ associated to the weight function $\nu$ and to a Kasteleyn orientation. More precisely, rows and columns of the matrix $K$ are indexed by vertices of $G$, and non-zero coefficients of $K$ are defined by,

$$ \forall \text{ edge } (x, y) \text{ of } G, \quad K_{x,y} = \varepsilon_{x,y} \nu_{xy}, $$

where

$$ \varepsilon_{x,y} = \begin{cases} 1 & \text{if } xy \in E \text{ and } x \to y \\ -1 & \text{if } xy \in E \text{ and } x \leftarrow y. \end{cases} $$
<figure><img src="image_3.png" /><figcaption>Figure 1: An example of a graph G (black), of its restricted dual G* (grey), of a spin configuration with + boundary conditions on G and its corresponding polygon configuration on G* (red). Vertices of G are pictured as bullets – black ones represent vertices of the original graph and grey ones are the additional vertices on boundary edges – vertices of G* are pictured as diamonds.</figcaption></figure>

dimer model per se, as well as that of the Kasteleyn matrix. We also recall the founding results that we will use.

Consider a planar, simple, embedded graph $G = (V, E)$. A dimer configuration of $G$, also known as a perfect matching, is a subset of edges such that every vertex is incident to exactly one edge of this subset. Denote by $\mathcal{M}(G)$ the set of dimer configurations of the graph $G$. Suppose that a positive weight function $\nu$ is assigned to edges of $G$.

### 2.2.1 Finite case

Suppose that the graph $G$ is finite, and that $|V|$ is even. Then, the probability of occurrence of a dimer configuration, chosen with respect to the dimer Boltzmann measure $\mathbb{P}_{\text{dimer}}$, is given by

$$ \forall M \in \mathcal{M}(G), \quad \mathbb{P}_{\text{dimer}}(M) = \frac{\prod_{e \in M} \nu_e}{Z_{\text{dimer}}(G, \nu)}, $$

where $Z_{\text{dimer}}(G, \nu) = \sum_{M \in \mathcal{M}(G)} \prod_{e \in M} \nu_e$ is the normalizing constant known as the dimer partition function.

The main tool used to study the dimer model is the Kasteleyn matrix [41, 42, 72], it is defined as follows. A face-cycle is a cycle of $G$ bounding a bounded face of the graph. A Kasteleyn orientation is an orientation of the edges such that every face-cycle is clockwise odd, meaning that when traveling clockwise around a face-cycle, the number of co-oriented edges is odd. By the results of [42], a Kasteleyn orientation always exists for planar graphs. A Kasteleyn matrix, denoted by $K$, is a weighted, directed, adjacency matrix of the graph $G$ associated to the weight function $\nu$ and to a Kasteleyn orientation. More precisely, rows and columns of the matrix $K$ are indexed by vertices of $G$, and non-zero coefficients of $K$ are defined by,

$$ \forall \text{ edge } (x, y) \text{ of } G, \quad K_{x,y} = \varepsilon_{x,y} \nu_{xy}, $$

where

$$ \varepsilon_{x,y} = \begin{cases} 1 & \text{if } xy \in E \text{ and } x \to y \\ -1 & \text{if } xy \in E \text{ and } x \leftarrow y. \end{cases} $$
Note that the matrix $K$ is skew symmetric.

When the graph $G$ is bipartite, the set of vertices can be split into $V = W \cup B$, where $W$ represents the set of white vertices, $B$ the set of black ones, and vertices in $W$ are only adjacent to vertices in $B$. Suppose that $|W| = |B|$ for otherwise $G$ has no dimer configurations. The Kasteleyn matrix $K$ is naturally block diagonal with two 0 blocks corresponding to rows/columns indexed by $W/W$ or $B/B$. It thus suffices to consider the bipartite, weighted, directed, adjacency matrix of the graph $G$, denoted by $\tilde{K}$. It has rows indexed by white vertices of $G$ and column by black ones. Non-zero coefficients are defined by:

$$ \forall \text{ edge } w b \text{ of } G, \quad \tilde{K}_{w,b} = \varepsilon_{w,b} \nu_{wb}. $$

Note that the bipartite Kasteleyn matrix can also be defined as minus the transpose of the above matrix $\tilde{K}$; rows are then indexed by black vertices and columns by white ones. Actually both bipartite Kasteleyn matrices are considered in this paper.

The two founding results of the dimer model are: an explicit expression for the partition function [41, 42, 72] and for the dimer Boltzmann measure [44]. Here are their statements.

**Theorem 2.1** ([41, 42, 72]). The dimer partition function of the graph $G$ with weight function $\nu$ is equal to:

$$ Z_{\text{dimer}}(G, \nu) = |\text{Pf} K|. $$

When the graph $G$ is moreover bipartite, we have:

$$ Z_{\text{dimer}}(G, \nu) = |\det \tilde{K}|. $$

**Theorem 2.2** ([44]). The probability of occurrence of a subset $\mathcal{E} = \{e_1 = x_1 y_1, \dots, e_l = x_l y_l\}$ of edges of $G$, chosen with respect to the dimer Boltzmann measure $P_{\text{dimer}}$ is equal to

$$ P_{\text{dimer}}(e_1, \dots, e_l) = \left( \prod_{j=1}^{l} K_{x_i, y_i} \right) \text{Pf}(K^{-1})_{\mathcal{E}}^t, $$

where $(K^{-1})_{\mathcal{E}}$ is the sub-matrix of the inverse Kasteleyn matrix $K^{-1}$ whose rows and columns are indexed by vertices $x_1, y_1, \dots, x_l, y_l$.

When the graph $G$ is moreover bipartite, the subset of edges $\mathcal{E}$ is written as $\mathcal{E} = \{e_1 = w_1 b_1, \dots, e_l = w_l b_l\}$, and we also have,

$$ P_{\text{dimer}}(e_1, \dots, e_l) = \left( \prod_{j=1}^{l} \tilde{K}_{w_i, b_i} \right) \det(\tilde{K}^{-1})_{\mathcal{E}}, $$

where $(\tilde{K}^{-1})_{\mathcal{E}}$ is the sub-matrix of the inverse bipartite Kasteleyn matrix $\tilde{K}^{-1}$ whose rows are indexed by black vertices $b_1, \dots, b_l$ and columns by white vertices $w_1, \dots, w_l$.

### 2.2.2 Infinite case

Suppose that the graph $G$ is infinite. The dimer Boltzmann measure is not well defined and is replaced by the notion of Gibbs measure. A Gibbs measure is a probability measure on $M(G)$ satisfying the DLR-conditions: when one fixes a dimer configuration in an annular region, then perfect matchings inside and outside of the annulus are independent; moreover, the probability of a dimer configuration in the finite region separated by the annulus is proportional to the product of the edge-weights.
Note that the matrix $K$ is skew symmetric.

When the graph $G$ is bipartite, the set of vertices can be split into $V = W \cup B$, where $W$ represents the set of white vertices, $B$ the set of black ones, and vertices in $W$ are only adjacent to vertices in $B$. Suppose that $|W| = |B|$ for otherwise $G$ has no dimer configurations. The Kasteleyn matrix $K$ is naturally block diagonal with two 0 blocks corresponding to rows/columns indexed by $W/W$ or $B/B$. It thus suffices to consider the bipartite, weighted, directed, adjacency matrix of the graph $G$, denoted by $\tilde{K}$. It has rows indexed by white vertices of $G$ and column by black ones. Non-zero coefficients are defined by:

$$ \forall \text{ edge } w b \text{ of } G, \quad \tilde{K}_{w,b} = \varepsilon_{w,b} \nu_{wb}. $$

Note that the bipartite Kasteleyn matrix can also be defined as minus the transpose of the above matrix $\tilde{K}$; rows are then indexed by black vertices and columns by white ones. Actually both bipartite Kasteleyn matrices are considered in this paper.

The two founding results of the dimer model are: an explicit expression for the partition function [41, 42, 72] and for the dimer Boltzmann measure [44]. Here are their statements.

**Theorem 2.1** ([41, 42, 72]). The dimer partition function of the graph $G$ with weight function $\nu$ is equal to:

$$ Z_{\text{dimer}}(G, \nu) = |\text{Pf} K|. $$

When the graph $G$ is moreover bipartite, we have:

$$ Z_{\text{dimer}}(G, \nu) = |\det \tilde{K}|. $$

**Theorem 2.2** ([44]). The probability of occurrence of a subset $\mathcal{E} = \{e_1 = x_1 y_1, \dots, e_l = x_l y_l\}$ of edges of $G$, chosen with respect to the dimer Boltzmann measure $P_{\text{dimer}}$ is equal to

$$ P_{\text{dimer}}(e_1, \dots, e_l) = \left( \prod_{j=1}^{l} K_{x_i, y_i} \right) \text{Pf}(K^{-1})_{\mathcal{E}}^t, $$

where $(K^{-1})_{\mathcal{E}}$ is the sub-matrix of the inverse Kasteleyn matrix $K^{-1}$ whose rows and columns are indexed by vertices $x_1, y_1, \dots, x_l, y_l$.

When the graph $G$ is moreover bipartite, the subset of edges $\mathcal{E}$ is written as $\mathcal{E} = \{e_1 = w_1 b_1, \dots, e_l = w_l b_l\}$, and we also have,

$$ P_{\text{dimer}}(e_1, \dots, e_l) = \left( \prod_{j=1}^{l} \tilde{K}_{w_i, b_i} \right) \det(\tilde{K}^{-1})_{\mathcal{E}}, $$

where $(\tilde{K}^{-1})_{\mathcal{E}}$ is the sub-matrix of the inverse bipartite Kasteleyn matrix $\tilde{K}^{-1}$ whose rows are indexed by black vertices $b_1, \dots, b_l$ and columns by white vertices $w_1, \dots, w_l$.

### 2.2.2 Infinite case

Suppose that the graph $G$ is infinite. The dimer Boltzmann measure is not well defined and is replaced by the notion of Gibbs measure. A Gibbs measure is a probability measure on $M(G)$ satisfying the DLR-conditions: when one fixes a dimer configuration in an annular region, then perfect matchings inside and outside of the annulus are independent; moreover, the probability of a dimer configuration in the finite region separated by the annulus is proportional to the product of the edge-weights.
Consider a Kasteleyn orientation of the graph $G$ and the corresponding Kasteleyn matrix $K$, then $K$ can also be seen as an operator acting on $C^V$:

$$ \forall F \in C^V, \quad (KF)_x = \sum_{y \sim x} K_{x,y} F_y. $$

When $G$ is bipartite, the bipartite Kasteleyn matrix $\tilde{K}$ is an operator mapping $C^W$ to $C^B$:

$$ \forall F \in C^W, \quad (\tilde{K}F)_b = \sum_{w \sim b} \tilde{K}_{b,w} F_w. $$

Explicit expressions for Gibbs measures typically involve inverse operators of a family of gauge equivalent Kasteleyn operators, see [47] for definitions and details. An inverse operator $L$ of a Kasteleyn operator $K$ is asked to satisfy the following conditions:

* $KL = \mathrm{Id}$ or $LK = \mathrm{Id}$,

* For every $x$, $L_{x,y} \to 0$ as $y$ tends to infinity.

Existence of an inverse Kasteleyn operator and explicit expressions for coefficients are proved for: $Z^2$-periodic bipartite graphs using Fourier techniques [23, 47]; $Z^2$-periodic (non-bipartite) Fisher graphs [9, 32]; non-periodic, bipartite isoradial graphs, bipartite quadri-tiling graphs, and Fisher graphs, all with specific weights arising from $Z$-invariance [45, 10, 12], see Sections 2.3.2, 2.3.3, 2.5 for definitions; coefficients of the inverse then have the remarkable property of being local. We refer to the original papers for the explicit expressions.

Uniqueness is established when the graph $G$ is $Z^2$-periodic [66, 9], and when the graph is non-periodic in the setting of the papers [45, 12, 8]. When the inverse Kasteleyn operator exists and is unique, it is denoted by $K^{-1}$. Note that uniqueness and the fact that the product $(KK^{-1})K = K(K^{-1}K)$ is associative implies that if $K^{-1}$ is a right, resp. left, inverse it is also a left, resp. right, inverse [24].

Consider the $\sigma$-field generated by cylinder sets of $M(G)$. In all of the above cases, there is an explicit expression for a Gibbs measure $P_{\mathrm{dimer}}$ on $(M(G), \mathcal{F})$ whose probabilities on cylinder sets is given by the formulas of Theorem 2.2 with $K^{-1}$ being the inverse Kasteleyn operator above. When the graph $G$ is moreover $Z^2$-periodic, this Gibbs measure is obtained as weak limit of the Boltzmann measures on the toroidal exhaustion $(G_n)_{n \ge 1}$, where $G_n = G/nZ^2$. We refer to the original papers for an exact statement, see also Theorem 3.19 which has the same form.

## 2.3 Dimer models on decorated graphs

In this paper, an important role is played by the dimer model on the double graph $G^D$, a model in correspondence with random pairs of dual directed spanning trees [71, 14, 49]. Furthermore, we consider two dimer representations of the Ising model. The first is related to the LTE of the Ising model [51, 52], while the second arises from the XOR-Ising model, built from two independent copies of the Ising model [32, 13]. The two corresponding dimer models live on the Fisher graph $G^F$ and the bipartite graph $G^Q$, respectively. The three graphs $G^D$, $G^F$ and $G^Q$ are decorated versions of the graph $G$.

In the next three sections, we define these decorated graphs and the mappings considered. We treat the case where the graph $G$ is infinite or finite. In the finite case, the graph $G$, the dual graph $\bar{G}^*$ and the restricted dual $G^*$ are those defined in Section 2.1, where recall that $G$ has an additional vertex on every boundary edge, and that $o$ denotes the vertex of $\bar{G}^*$ corresponding to the unbounded face of $G$. In the infinite case, the dual graph is $G^*$. Figures illustrate the finite case; a local picture of the infinite case is obtained by looking at the interior of the finite case.
Consider a Kasteleyn orientation of the graph $G$ and the corresponding Kasteleyn matrix $K$, then $K$ can also be seen as an operator acting on $C^V$:

$$ \forall F \in C^V, \quad (KF)_x = \sum_{y \sim x} K_{x,y} F_y. $$

When $G$ is bipartite, the bipartite Kasteleyn matrix $\tilde{K}$ is an operator mapping $C^W$ to $C^B$:

$$ \forall F \in C^W, \quad (\tilde{K}F)_b = \sum_{w \sim b} \tilde{K}_{b,w} F_w. $$

Explicit expressions for Gibbs measures typically involve inverse operators of a family of gauge equivalent Kasteleyn operators, see [47] for definitions and details. An inverse operator $L$ of a Kasteleyn operator $K$ is asked to satisfy the following conditions:

* $KL = \mathrm{Id}$ or $LK = \mathrm{Id}$,

* For every $x$, $L_{x,y} \to 0$ as $y$ tends to infinity.

Existence of an inverse Kasteleyn operator and explicit expressions for coefficients are proved for: $Z^2$-periodic bipartite graphs using Fourier techniques [23, 47]; $Z^2$-periodic (non-bipartite) Fisher graphs [9, 32]; non-periodic, bipartite isoradial graphs, bipartite quadri-tiling graphs, and Fisher graphs, all with specific weights arising from $Z$-invariance [45, 10, 12], see Sections 2.3.2, 2.3.3, 2.5 for definitions; coefficients of the inverse then have the remarkable property of being *local*. We refer to the original papers for the explicit expressions.

Uniqueness is established when the graph $G$ is $Z^2$-periodic [66, 9], and when the graph is non-periodic in the setting of the papers [45, 12, 8]. When the inverse Kasteleyn operator exists and is unique, it is denoted by $K^{-1}$. Note that uniqueness and the fact that the product $(KK^{-1})K = K(K^{-1}K)$ is associative implies that if $K^{-1}$ is a right, resp. left, inverse it is also a left, resp. right, inverse [24].

Consider the $\sigma$-field generated by cylinder sets of $M(G)$. In all of the above cases, there is an explicit expression for a Gibbs measure $P_{\mathrm{dimer}}$ on $(M(G), \mathcal{F})$ whose probabilities on cylinder sets is given by the formulas of Theorem 2.2 with $K^{-1}$ being the inverse Kasteleyn operator above. When the graph $G$ is moreover $Z^2$-periodic, this Gibbs measure is obtained as weak limit of the Boltzmann measures on the toroidal exhaustion $(G_n)_{n \ge 1}$, where $G_n = G/nZ^2$. We refer to the original papers for an exact statement, see also Theorem 3.19 which has the same form.

## 2.3 Dimer models on decorated graphs

In this paper, an important role is played by the dimer model on the *double graph* $G^D$, a model in correspondence with random pairs of dual directed spanning trees [71, 14, 49]. Furthermore, we consider two dimer representations of the Ising model. The first is related to the LTE of the Ising model [51, 52], while the second arises from the XOR-Ising model, built from two independent copies of the Ising model [32, 13]. The two corresponding dimer models live on the *Fisher graph* $G^F$ and the bipartite graph $G^Q$, respectively. The three graphs $G^D$, $G^F$ and $G^Q$ are decorated versions of the graph $G$.

In the next three sections, we define these decorated graphs and the mappings considered. We treat the case where the graph $G$ is infinite or finite. In the finite case, the graph $G$, the dual graph $\bar{G}^*$ and the restricted dual $G^*$ are those defined in Section 2.1, where recall that $G$ has an additional vertex on every boundary edge, and that $\bar{0}$ denotes the vertex of $\bar{G}^*$ corresponding to the unbounded face of $G$. In the infinite case, the dual graph is $G^*$. Figures illustrate the finite case; a local picture of the infinite case is obtained by looking at the interior of the finite case.
### 2.3.1 Dimers on the double graph $G^D$ and Temperley's bijection

The double graph is denoted by $G^D = (V^D, E^D)$. It is defined as follows, see also Figure 2 (left).

**Infinite case.** Embed the dual graph $G^*$ so that edges of the primal and the dual intersect at a single point. The double graph is obtained by superimposing $G$ and $G^*$ and adding an extra vertex at the crossing of each primal and dual edge.

**Finite case.** It is constructed similarly to the infinite case from the superimposition of $G$ and the dual graph $\bar{G}^*$. Edges incident to the vertex o are then removed.

In the infinite and fine cases, the double graph $G^D$ is bipartite and face-cycles are quadrangles. The set of black vertices of $G^D$, denoted by $B$, consists of vertices of $G$ and $G^*$; the set of white vertices of $G^D$, denoted by $W$, consists of vertices at the crossing of edges of $G$ and $G^*$ in the infinite case, and of $G$ and $\bar{G}^*$ in the finite case. White vertices are in bijection with edges of the graph $G$, or equivalently with edges of the dual graph. We thus have, $V^D = B \cup W$, where $B = V \cup V^*$ and $W \leftrightarrow E$.

Suppose again that $G$ is finite, fix a vertex $r$ of $G$ amongst the additional vertices on boundary edges, and let $V^r = V \setminus \{r\}$. Denote by $G^{D,r}$ the graph obtained from $G^D$ by removing the vertex $r$ and all edges incident to it. The graph $G^{D,r}$ is also bipartite; its set of black vertices is $B^r$, where $B^r = V^r \cup V^*$ and its set of white vertices is $W^r = W \leftrightarrow E$, see Figure 2 (right) for an example. Note that $G^{D,r}$ has the same number of black and white vertices: $|B^r| = |W^r|$.

**Bijection between pairs of dual directed spanning trees and dimers** Suppose that $G$ is finite. Prior to stating the bijection, we need a few definitions. A tree of $G$ is an acyclic connected subset of edges. A spanning tree is a tree spanning all vertices of the graph. Let $v$ be a vertex of $G$, then a $v$-directed spanning tree ($v$-dST) is obtained from a spanning tree by directing all edges towards the vertex $v$, referred to as the *root*; with such an orientation, every vertex has exactly one outgoing edge except the root which has none. Given a spanning tree, the set of dual edges of the edges absent in the spanning tree form a spanning tree of the dual graph $\bar{G}^*$, known as the *dual spanning tree*.

Consider the fixed boundary vertex $r$ of $G$ as above. Denote by $\mathcal{T}(G)$ the set of $r$-dST of $G$, by $\mathcal{J}^o(\bar{G}^*)$ the set of $o$-dST of $\bar{G}^*$, and by $\mathcal{T}^{r,o}(G, \bar{G}^*)$ the set of pairs of dual directed spanning trees (dST-pairs) of $G$ and $\bar{G}^*$ such that the primal tree is rooted at $r$ and the dual tree is rooted at $o$, see Figure 2 (left) for an example.

The result of Temperley [71], extended by [14] to general non-directed graphs and by [49] to the directed case, proves a weight preserving bijection between dimer configurations of the double graph $G^{D,r}$ and dST-pairs of $\mathcal{T}^{r,o}(G, \bar{G}^*)$. It relies on the following bijection between edges of $G^{D,r}$ and directed edges of $G$ and $\bar{G}^*$. Let $w \in W^r = W$, $x \in V^r \cup V^*$ such that $wx$ is an edge of $G^{D,r}$, then

$$wx \leftrightarrow \begin{cases} (v, v') \text{ of } G & \text{if } x = v \in V^r \text{ and } v' \text{ is s.t. } w \text{ belongs to the edge } (v, v') \\ (f, f') \text{ of } \bar{G}^* & \text{if } x = f \in V^* \text{ and } f' \text{ is s.t. } w \text{ belongs to the edge } (f, f'). \end{cases} \quad (2.2)$$

Note that there are no directed edges of $G$ exiting the vertex $r$, and no directed edges of $\bar{G}^*$ exiting the vertex $o$. Using this bijection, a subset of edges of $G^{D,r}$ corresponds to a subset of directed edges of $G$ and $\bar{G}^*$; Temperley's bijection states that subsets defining dimer configurations are in correspondence with subsets defining dST-pairs of $\mathcal{T}^{r,o}(G, \bar{G}^*)$. An example is provided in Figure 2, the vertex $o$ is represented in a spread-out way, i.e., the dotted line should be thought of as being the single vertex $o$.

Let $c$ be a weight function on edges of $G^{D,r}$ and $\tilde{c}$ a weight function on directed edges of $G, \bar{G}^*$. The relation between $c$ and $\tilde{c}$ which makes Temperley's bijection weight
### 2.3.1 Dimers on the double graph $G^D$ and Temperley's bijection

The double graph is denoted by $G^D = (V^D, E^D)$. It is defined as follows, see also Figure 2 (left).

**Infinite case.** Embed the dual graph $G^*$ so that edges of the primal and the dual intersect at a single point. The double graph is obtained by superimposing $G$ and $G^*$ and adding an extra vertex at the crossing of each primal and dual edge.

**Finite case.** It is constructed similarly to the infinite case from the superimposition of $G$ and the dual graph $\bar{G}^*$. Edges incident to the vertex o are then removed.

In the infinite and fine cases, the double graph $G^D$ is bipartite and face-cycles are quadrangles. The set of black vertices of $G^D$, denoted by $B$, consists of vertices of $G$ and $G^*$; the set of white vertices of $G^D$, denoted by $W$, consists of vertices at the crossing of edges of $G$ and $G^*$ in the infinite case, and of $G$ and $\bar{G}^*$ in the finite case. White vertices are in bijection with edges of the graph $G$, or equivalently with edges of the dual graph. We thus have, $V^D = B \cup W$, where $B = V \cup V^*$ and $W \leftrightarrow E$.

Suppose again that $G$ is finite, fix a vertex $r$ of $G$ amongst the additional vertices on boundary edges, and let $V^r = V \setminus \{r\}$. Denote by $G^{D,r}$ the graph obtained from $G^D$ by removing the vertex $r$ and all edges incident to it. The graph $G^{D,r}$ is also bipartite; its set of black vertices is $B^r$, where $B^r = V^r \cup V^*$ and its set of white vertices is $W^r = W \leftrightarrow E$, see Figure 2 (right) for an example. Note that $G^{D,r}$ has the same number of black and white vertices: $|B^r| = |W^r|$.

**Bijection between pairs of dual directed spanning trees and dimers** Suppose that $G$ is finite. Prior to stating the bijection, we need a few definitions. A tree of $G$ is an acyclic connected subset of edges. A spanning tree is a tree spanning all vertices of the graph. Let $v$ be a vertex of $G$, then a $v$-directed spanning tree ($v$-dST) is obtained from a spanning tree by directing all edges towards the vertex $v$, referred to as the *root*; with such an orientation, every vertex has exactly one outgoing edge except the root which has none. Given a spanning tree, the set of dual edges of the edges absent in the spanning tree form a spanning tree of the dual graph $\bar{G}^*$, known as the *dual spanning tree*.

Consider the fixed boundary vertex $r$ of $G$ as above. Denote by $\mathcal{T}(G)$ the set of $r$-dST of $G$, by $\mathcal{J}^o(\bar{G}^*)$ the set of $o$-dST of $\bar{G}^*$, and by $\mathcal{T}^{r,o}(G, \bar{G}^*)$ the set of pairs of dual directed spanning trees (dST-pairs) of $G$ and $\bar{G}^*$ such that the primal tree is rooted at $r$ and the dual tree is rooted at $o$, see Figure 2 (left) for an example.

The result of Temperley [71], extended by [14] to general non-directed graphs and by [49] to the directed case, proves a weight preserving bijection between dimer configurations of the double graph $G^{D,r}$ and dST-pairs of $\mathcal{T}^{r,o}(G, \bar{G}^*)$. It relies on the following bijection between edges of $G^{D,r}$ and directed edges of $G$ and $\bar{G}^*$. Let $w \in W^r = W$, $x \in V^r \cup V^*$ such that $wx$ is an edge of $G^{D,r}$, then

$$wx \leftrightarrow \begin{cases} (v, v') \text{ of } G & \text{if } x = v \in V^r \text{ and } v' \text{ is s.t. } w \text{ belongs to the edge } (v, v') \\ (f, f') \text{ of } \bar{G}^* & \text{if } x = f \in V^* \text{ and } f' \text{ is s.t. } w \text{ belongs to the edge } (f, f'). \end{cases} \quad (2.2)$$

Note that there are no directed edges of $G$ exiting the vertex $r$, and no directed edges of $\bar{G}^*$ exiting the vertex $o$. Using this bijection, a subset of edges of $G^{D,r}$ corresponds to a subset of directed edges of $G$ and $\bar{G}^*$; Temperley's bijection states that subsets defining dimer configurations are in correspondence with subsets defining dST-pairs of $\mathcal{T}^{r,o}(G, \bar{G}^*)$. An example is provided in Figure 2, the vertex $o$ is represented in a spread-out way, i.e., the dotted line should be thought of as being the single vertex $o$.

Let $c$ be a weight function on edges of $G^{D,r}$ and $\tilde{c}$ a weight function on directed edges of $G, \bar{G}^*$. The relation between $c$ and $\tilde{c}$ which makes Temperley's bijection weight
<figure><img src="image_4.png" /><figcaption>Figure 2: Left: double graph $G^D$ of the graph $G$ of Figure 1 (grey lines), r-directed spanning tree of $G$ (turquoise) and dual o-directed spanning tree of $\bar{G}^*$ (purple). Right: graph $G^{D,r}$ (grey) and the dimer configuration (red) in bijection with the pair of dual directed spanning trees of the left figure.</figcaption></figure>

preserving naturally arises from the bijection between edges of $G^{D,r}$ and directed edges of $G, \bar{G}^*$. Let $w \in W^r, x \in V^r \cup V^*$, such that $wx$ is an edge of $G^{D,r}$. Using the notation of (2.2), we have

$$c_{wx} = \begin{cases} \tilde{c}_{v,v'} & \text{if } x = v \in V^r \\ \tilde{c}_{f,f'} & \text{if } x = f \in V^*, \end{cases} \qquad (2.3)$$

and $\tilde{c}_{r,v} = 0$ for every vertex $v \in V^r$ adjacent to $r$, $\tilde{c}_{o,f'} = 0$ for every vertex $f' \in V^*$ adjacent to $o$.

**Model on pairs of dual directed spanning trees** Suppose that directed edges of $G, \bar{G}^*$ are assigned the weight function $\tilde{c}$. Consider the Boltzmann measure on dST-pairs, denoted $P_{\text{dST-pairs}}^{r,o}$, defined by

$$\forall (T, T^*) \in \mathcal{T}^{r,o}(G, \bar{G}^*), \quad P_{\text{dST-pairs}}^{r,o}(T, T^*) = \frac{(\prod_{(v,v') \in T} \tilde{c}_{v,v'})(\prod_{(f,f') \in T^*} \tilde{c}_{f,f'})}{Z_{\text{dST-pairs}}^{r,o}(G, \bar{G}^*)},$$

where $Z_{\text{dST-pairs}}^{r,o}((G, \bar{G}^*), \tilde{c}) = \sum_{(T, T^*) \in \mathcal{T}^{r,o}(G, \bar{G}^*)} (\prod_{(v,v') \in T} \tilde{c}_{v,v'})(\prod_{(f,f') \in T^*} \tilde{c}_{f,f'})$, is the dST-pairs partition function. As a consequence of the KPW-Temperley bijection [71, 49], we have

$$Z_{\text{dST-pairs}}^{r,o}((G, \bar{G}^*), \tilde{c}) = Z_{\text{dimer}}(G^{D,r}, c).$$

There is also a natural correspondence between the dST-pairs Boltzmann measure $P_{\text{dST-pairs}}^{r,o}$ and the dimer Boltzmann measure $P_{\text{dimer}}^{D}$ on $G^{D,r}$ with weight function $c$.

Note that if $\tilde{c} \equiv 1$ on edges of $\bar{G}^*$, resp. on edges of $G$, then $Z_{\text{dST-pairs}}^{r,o}((G, \bar{G}^*), \tilde{c})$ is equal to the partition function $Z_{\text{dST}}^r(G, \tilde{c})$ of r-directed spanning trees of $G$, resp. $Z_{\text{dST}}^o(\bar{G}^*, \tilde{c})$ of o-directed spanning trees of $\bar{G}^*$.

### 2.3.2 Dimers on the Fisher graph $G^F$ and the LTE of the Ising model

The Fisher graph is denoted by $G^F = (V^F, E^F)$. It is constructed as follows [35, 32], see Figure 3 for an example.
<figure><img src="image_4.png" /><figcaption>Figure 2: Left: double graph $G^D$ of the graph $G$ of Figure 1 (grey lines), r-directed spanning tree of $G$ (turquoise) and dual o-directed spanning tree of $\bar{G}^*$ (purple). Right: graph $G^{D,r}$ (grey) and the dimer configuration (red) in bijection with the pair of dual directed spanning trees of the left figure.</figcaption></figure>

preserving naturally arises from the bijection between edges of $G^{D,r}$ and directed edges of $G, \bar{G}^*$. Let $w \in W^r, x \in V^r \cup V^*$, such that $wx$ is an edge of $G^{D,r}$. Using the notation of (2.2), we have

$$c_{wx} = \begin{cases} \tilde{c}_{v,v'} & \text{if } x = v \in V^r \\ \tilde{c}_{f,f'} & \text{if } x = f \in V^*, \end{cases} \qquad (2.3)$$

and $\tilde{c}_{r,v} = 0$ for every vertex $v \in V^r$ adjacent to $r$, $\tilde{c}_{o,f'} = 0$ for every vertex $f' \in V^*$ adjacent to $o$.

**Model on pairs of dual directed spanning trees** Suppose that directed edges of $G, \bar{G}^*$ are assigned the weight function $\tilde{c}$. Consider the *Boltzmann measure on dST-pairs*, denoted $P_{\text{dST-pairs}}^{r,o}$, defined by

$$\forall (T, T^*) \in \mathcal{T}^{r,o}(G, \bar{G}^*), \quad P_{\text{dST-pairs}}^{r,o}(T, T^*) = \frac{(\prod_{(v,v') \in T} \tilde{c}_{v,v'})(\prod_{(f,f') \in T^*} \tilde{c}_{f,f'})}{Z_{\text{dST-pairs}}^{r,o}(G, \bar{G}^*)},$$

where $Z_{\text{dST-pairs}}^{r,o}((G, \bar{G}^*), \tilde{c}) = \sum_{(T, T^*) \in \mathcal{T}^{r,o}(G, \bar{G}^*)} (\prod_{(v,v') \in T} \tilde{c}_{v,v'})(\prod_{(f,f') \in T^*} \tilde{c}_{f,f'})$, is the *dST-pairs partition function*. As a consequence of the KPW-Temperley bijection [71, 49], we have

$$Z_{\text{dST-pairs}}^{r,o}((G, \bar{G}^*), \tilde{c}) = Z_{\text{dimer}}(G^{D,r}, c).$$

There is also a natural correspondence between the dST-pairs Boltzmann measure $P_{\text{dST-pairs}}^{r,o}$ and the dimer Boltzmann measure $P_{\text{dimer}}^{D}$ on $G^{D,r}$ with weight function $c$.

Note that if $\tilde{c} \equiv 1$ on edges of $\bar{G}^*$, resp. on edges of $G$, then $Z_{\text{dST-pairs}}^{r,o}((G, \bar{G}^*), \tilde{c})$ is equal to the partition function $Z_{\text{dST}}^r(G, \tilde{c})$ of r-directed spanning trees of $G$, resp. $Z_{\text{dST}}^o(\bar{G}^*, \tilde{c})$ of o-directed spanning trees of $\bar{G}^*$.

### 2.3.2 Dimers on the Fisher graph $G^F$ and the LTE of the Ising model

The *Fisher graph* is denoted by $G^F = (V^F, E^F)$. It is constructed as follows [35, 32], see Figure 3 for an example.
**Infinite case.** Start from the dual graph $G^*$ and replace every vertex of $G^*$ by a *decoration* made of triangles, where each of the triangles corresponds to an edge incident to this vertex, then join the triangles in a circular way.

**Finite case.** Start from the dual graph $\bar{G}^*$ and do the same procedure as in the infinite case. Then, remove the decoration of the vertex o as well as all edges of $\bar{G}^*$ incident to this decoration.

In both the infinite and finite case, the Fisher graph consists of *internal edges*, which are edges of the decorations, and *external edges* which are in bijection with edges of $G^*$ and will often be identified with them. Each decoration has a dual vertex in its center, giving a way of identifying decorations and vertices of $G^*$.

<figure><img src="image_5.png" /><figcaption>Figure 3: The Fisher graph $G^F$ for the LTE expansion of the Ising model on $G$ with + boundary conditions (black); one of the $2^{13}$ dimer configurations corresponding to the polygon configuration of Figure 1.</figcaption></figure>

**Mapping between LTE polygon configurations and dimers** Suppose that $G$ is finite. Fisher [35] introduces a mapping between polygon configurations of $G^*$ and dimer configurations of the corresponding Fisher graph $G^F$. To a given polygon configuration of $G^*$, there corresponds $2^{|V^*|}$ dimer configurations of $G^F$: edges of the polygon configuration are exactly the external edges of the dimer configurations and given these external edges, there is exactly two ways of filling each decoration so as to have a dimer configuration [35], see Figure 3 for an example. This mapping naturally extends when the graph $G$ is infinite.

We consider polygon configurations arising from the LTE expansion of the Ising model on $G$ with + boundary conditions and coupling constants $J$. In order for this correspondence to be weight preserving up to a multiplicative constant, the dimer weight function $\mu^j$ on edges of $G^F$ is defined to be, see Equation (2.1), for every edge $e$ of $G^F$,

$$ \mu_e^j = \begin{cases} 1 & \text{if } e \text{ is an internal edge} \\ e^{-2J_e} & \text{if } e \text{ is an external edge arising from a dual edge } e^* \text{ of } G^*. \end{cases} $$

Let $P_F^\text{dimer}$ and $Z_\text{dimer}(G^F, \mu^j)$ be the corresponding dimer Boltzmann measure and partition function. Then, as a consequence of Fisher's correspondence we have,

$$ Z_{\text{Ising}}^{+}(G, J) = 2^{-|V^*|} \left( \prod_{e \in E} e^{J_e} \right) Z_{\text{dimer}}(G^F, \mu^j). \quad (2.4) $$
**Infinite case.** Start from the dual graph $G^*$ and replace every vertex of $G^*$ by a *decoration* made of triangles, where each of the triangles corresponds to an edge incident to this vertex, then join the triangles in a circular way.

**Finite case.** Start from the dual graph $\bar{G}^*$ and do the same procedure as in the infinite case. Then, remove the decoration of the vertex o as well as all edges of $\bar{G}^*$ incident to this decoration.

In both the infinite and finite case, the Fisher graph consists of *internal edges*, which are edges of the decorations, and *external edges* which are in bijection with edges of $G^*$ and will often be identified with them. Each decoration has a dual vertex in its center, giving a way of identifying decorations and vertices of $G^*$.

<figure><img src="image_5.png" /><figcaption>Figure 3: The Fisher graph $G^F$ for the LTE expansion of the Ising model on $G$ with + boundary conditions (black); one of the $2^{13}$ dimer configurations corresponding to the polygon configuration of Figure 1.</figcaption></figure>

**Mapping between LTE polygon configurations and dimers** Suppose that $G$ is finite. Fisher [35] introduces a mapping between polygon configurations of $G^*$ and dimer configurations of the corresponding Fisher graph $G^F$. To a given polygon configuration of $G^*$, there corresponds $2^{|V^*|}$ dimer configurations of $G^F$: edges of the polygon configuration are exactly the external edges of the dimer configurations and given these external edges, there is exactly two ways of filling each decoration so as to have a dimer configuration [35], see Figure 3 for an example. This mapping naturally extends when the graph $G$ is infinite.

We consider polygon configurations arising from the LTE expansion of the Ising model on $G$ with + boundary conditions and coupling constants $J$. In order for this correspondence to be weight preserving up to a multiplicative constant, the dimer weight function $\mu^j$ on edges of $G^F$ is defined to be, see Equation (2.1), for every edge $e$ of $G^F$,

$$ \mu_e^j = \begin{cases} 1 & \text{if } e \text{ is an internal edge} \\ e^{-2J_e} & \text{if } e \text{ is an external edge arising from a dual edge } e^* \text{ of } G^*. \end{cases} $$

Let $P_F^\text{dimer}$ and $Z_\text{dimer}(G^F, \mu^j)$ be the corresponding dimer Boltzmann measure and partition function. Then, as a consequence of Fisher's correspondence we have,

$$ Z_{\text{Ising}}^{+}(G, J) = 2^{-|V^*|} \left( \prod_{e \in E} e^{J_e} \right) Z_{\text{dimer}}(G^F, \mu^j). \quad (2.4) $$
### 2.3.3 Dimers on the bipartite graph $G^Q$ and the XOR-Ising model

The quadri-tiling graph is denoted by $G^Q = (V^Q, E^Q)$, where the name comes from the paper [28]. In both the finite and infinite cases, we start from the preceding definition of the double graph $G^D$. Recall that face-cycles of $G^D$ are quadrangles consisting of two black and two white vertices, then add the edges joining opposite black vertices in quadrangles.

**Infinite case.** The graph $G^Q$ is the dual of the modified graph $G^D$.

**Finite case.** The graph $G^Q$ is the restricted dual of the modified graph $G^D$, see Figure 4.

Vertices of $G^Q$ are partitioned as $V^Q = B \cup W$, and the bipartite coloring is fixed as in Figure 4. Black, resp. white, vertices of $G^Q$ are denoted by b, resp. w, with or with sub/super-scripts.

In the infinite case, the graph $G^Q$ consists of *quadrangles* that are joined by *external edges*. Quadrangles are in bijection with edges of G, or equivalently edges of $G^*$, or equivalently white vertices of $G^D$: each quadrangle has a white vertex of $G^D$ in its interior, two of its edges are "parallel" to an edge of G and the two other edges are "parallel" to the dual edge of $G^*$. Face-cycles of $G^Q$ other than quadrangles either have a vertex of G or a vertex of $G^*$ in their interior.

In the finite case, the description is similar away from the boundary. Along the boundary "quadrangles" in bijection with boundary edges of G, or equivalently with boundary white vertices of $G^D$, are actually reduced to single edges "parallel" to boundary edges of G. We refer to those degenerate quadrangles as *boundary quadrangles* of $G^Q$, keeping in mind that they actually are *edges*. Note that some quadrangle edges of $G^Q$ are boundary edges of $G^Q$ (in the sense that they belong to the boundary of the unbounded face) but still belong to "full" quadrangles; as such they are *not* boundary quadrangle edges.

<figure><img src="image_6.png" /><figcaption>Figure 4: The quadri-tiling graph $G^Q$: take the double graph $G^D$ of Figure 2 (left) and add edges joining opposite black vertices in quadrangle-faces; the restricted dual of this graph is $G^Q$.</figcaption></figure>

We consider the dimer model on the bipartite graph $G^Q$ arising from the XOR-Ising model [74], also known as the *polarization* of the Ising model [38], obtained by taking the product of the spins of two independent Ising models. There are two mappings leading to the dimer model on $G^Q$, both of them are rather long to describe so that we refer to the original papers: [32] based on results of [40, 75, 34, 76] for the first approach, and [13] based on results of [62, 76] for the second one. Note that the last part of the above constructions establishes a correspondence between the dimer model on $G^Q$ and a free-fermion, zero-field 6-vertex model on the medial graph [62, 76, 32]. In
### 2.3.3 Dimers on the bipartite graph $G^Q$ and the XOR-Ising model

The quadri-tiling graph is denoted by $G^Q = (V^Q, E^Q)$, where the name comes from the paper [28]. In both the finite and infinite cases, we start from the preceding definition of the double graph $G^D$. Recall that face-cycles of $G^D$ are quadrangles consisting of two black and two white vertices, then add the edges joining opposite black vertices in quadrangles.

**Infinite case.** The graph $G^Q$ is the dual of the modified graph $G^D$.

**Finite case.** The graph $G^Q$ is the restricted dual of the modified graph $G^D$, see Figure 4.

Vertices of $G^Q$ are partitioned as $V^Q = B \cup W$, and the bipartite coloring is fixed as in Figure 4. Black, resp. white, vertices of $G^Q$ are denoted by b, resp. w, with or with sub/super-scripts.

In the infinite case, the graph $G^Q$ consists of *quadrangles* that are joined by *external edges*. Quadrangles are in bijection with edges of G, or equivalently edges of $G^*$, or equivalently white vertices of $G^D$: each quadrangle has a white vertex of $G^D$ in its interior, two of its edges are "parallel" to an edge of G and the two other edges are "parallel" to the dual edge of $G^*$. Face-cycles of $G^Q$ other than quadrangles either have a vertex of G or a vertex of $G^*$ in their interior.

In the finite case, the description is similar away from the boundary. Along the boundary "quadrangles" in bijection with boundary edges of G, or equivalently with boundary white vertices of $G^D$, are actually reduced to single edges "parallel" to boundary edges of G. We refer to those degenerate quadrangles as *boundary quadrangles* of $G^Q$, keeping in mind that they actually are *edges*. Note that some quadrangle edges of $G^Q$ are boundary edges of $G^Q$ (in the sense that they belong to the boundary of the unbounded face) but still belong to "full" quadrangles; as such they are *not* boundary quadrangle edges.

<figure><img src="image_6.png" /><figcaption>Figure 4: The quadri-tiling graph $G^Q$: take the double graph $G^D$ of Figure 2 (left) and add edges joining opposite black vertices in quadrangle-faces; the restricted dual of this graph is $G^Q$.</figcaption></figure>

We consider the dimer model on the bipartite graph $G^Q$ arising from the XOR-Ising model [74], also known as the *polarization* of the Ising model [38], obtained by taking the product of the spins of two independent Ising models. There are two mappings leading to the dimer model on $G^Q$, both of them are rather long to describe so that we refer to the original papers: [32] based on results of [40, 75, 34, 76] for the first approach, and [13] based on results of [62, 76] for the second one. Note that the last part of the above constructions establishes a correspondence between the dimer model on $G^Q$ and a free-fermion, zero-field 6-vertex model on the medial graph [62, 76, 32]. In
order not to lengthen this exposition we refer the reader to, for example, Sections 5.1 and 5.2 of [13] for a summary of this correspondence.

The dimer weight function $\nu^j$ on $G^Q$ is defined by, for every edge $e$ of $G^Q$,

$$ \nu_e^j = \begin{cases} 1 & \text{if } e \text{ is an external edge} \\ 1 & \text{if } e \text{ is a boundary quadrangle edge in the finite case} \\ \tanh(2J_e) & \text{if } e \text{ is a quadrangle edge/non-boundary quadrangle edge in the infinite/finite case, "parallel" to an edge } e \text{ of } G \\ \cosh^{-1}(2J_e) & \text{if } e \text{ is a quadrangle edge, "parallel" to a dual edge } e^* \text{ of an edge } e \text{ of } G. \end{cases} $$

When the graph $G^Q$ is finite, we let $P_{\text{dimer}}^Q$ and $Z_{\text{dimer}}(G^Q, \nu^j)$ be the corresponding dimer Boltzmann measure and partition function. As a consequence of [32], see also Corollary 5.3, we have:

$$ Z_{\text{dimer}}(G^Q, \nu^j)^2 = 2^{|V^*|} \prod_{e^* \in E^*} (1 + e^{-4J_e}) Z_{\text{dimer}}(G^Q, \nu^j). $$

Combining this with Equation (2.4) for the Ising partition function, and denoting by $E^\theta$ the set of boundary edges of the graph $G$, we obtain

$$ \begin{aligned} [Z_{\text{Ising}}^{+}(G, J)]^2 &= 2^{-|V^*|} \left(\prod_{e \in E} e^{2J_e}\right) \left(\prod_{e^* \in E^*} (1 + e^{-4J_e})\right) Z_{\text{dimer}}(G^Q, \nu^j) \\ &= 2^{-|V^*|+|E^*|} \left(\prod_{e \in E^\theta} e^{2J_e}\right) \left(\prod_{e^* \in E^*} \cosh(2J_e)\right) Z_{\text{dimer}}(G^Q, \nu^j) \\ &= 2^{|V|-1} \left(\prod_{e \in E^\theta} \frac{e^{2J_e}}{2}\right) \left(\prod_{e^* \in E^*} \cosh(2J_e)\right) Z_{\text{dimer}}(G^Q, \nu^j), \end{aligned} \quad (2.5) $$

where in the last line we used that $|E| = |E^*| + |E^\theta|$ and Euler's formula: $|E| = |V| + |V^*| - 1$.

## 2.4 Rooted directed spanning forests and directed spanning trees

We also need the model of *rooted directed spanning forests* on the graphs $G$ and $G^*$, resp. $G^\Gamma$ and $G^*$, in the infinite, resp. finite, case. So as to include both the primal and the dual graphs, we now define this model on a simple graph $G = (V, E)$.

Suppose that vertices are assigned non-negative masses, denoted $m = (m_x)_{x \in V}$, and that directed edges have positive conductances, denoted $\rho$, meaning that every directed edge $(x, x')$ has conductance $\rho_{x, x'}$.

A *rooted directed spanning forest* (rdSF) of $G$ is a subset of edges spanning all vertices of the graph, such that each connected component is a directed tree rooted at a vertex of $G$, denoted $x_T$. Let $\mathcal{F}(G)$ denote the set of rdSF of the graph $G$.

Suppose that $G$ is finite and consider the *Boltzmann measure on rdSF*, denoted $P_{\text{rdSF}}$, defined by:

$$ \forall F \in \mathcal{F}(G), \quad P_{\text{rdSF}}(F) = \frac{\prod_{T \in F} (m_{x_T} \prod_{(x,x') \in T} \rho_{x,x'})}{Z_{\text{rdSF}}(G, \rho, m)}, $$

where $Z_{\text{rdSF}}(G, \rho, m) = \sum_{F \in \mathcal{F}(G)} \prod_{T \in F} m_{x_T} \prod_{(x,x') \in T} \rho_{x,x'}$ is the *rdSF partition function*. Whenever conductances are symmetric, i.e., $\rho_{x,x'} = \rho_{x',x}$, we will remove the “d” in rdSF.

As a consequence of the directed version of Kirchhoff's matrix-tree theorem [50, 73], the rdSF partition function is computed using the massive Laplacian operator/matrix as follows. The *massive Laplacian operator* $\Delta^m: \mathbb{C}^V \to \mathbb{C}^V$ is defined by:

$$ \forall F \in \mathbb{C}^V, \quad (\Delta^m F)_x = \sum_{x' \sim x} \rho_{x,x'} (F_x - F_{x'}) + m_x F_x. $$
order not to lengthen this exposition we refer the reader to, for example, Sections 5.1 and 5.2 of [13] for a summary of this correspondence.

The dimer weight function $\nu^j$ on $G^Q$ is defined by, for every edge $e$ of $G^Q$,

$$ \nu_e^j = \begin{cases} 1 & \text{if } e \text{ is an external edge} \\ 1 & \text{if } e \text{ is a boundary quadrangle edge in the finite case} \\ \tanh(2J_e) & \text{if } e \text{ is a quadrangle edge/non-boundary quadrangle edge} \\ \cosh^{-1}(2J_e) & \text{in the infinite/finite case, "parallel" to an edge } e \text{ of } G \\ \cosh^{-1}(2J_e) & \text{if } e \text{ is a quadrangle edge, "parallel" to a dual edge } e^* \\ \cosh^{-1}(2J_e) & \text{of an edge } e \text{ of } G. \end{cases} $$

When the graph $G^Q$ is finite, we let $P_{\text{dimer}}^Q$ and $Z_{\text{dimer}}(G^Q, \nu^j)$ be the corresponding dimer Boltzmann measure and partition function. As a consequence of [32], see also Corollary 5.3, we have:

$$ Z_{\text{dimer}}(G^Q, \nu^j)^2 = 2^{|V^*|} \prod_{e^* \in E^*} (1 + e^{-4J_e}) Z_{\text{dimer}}(G^Q, \nu^j). $$

Combining this with Equation (2.4) for the Ising partition function, and denoting by $E^\theta$ the set of boundary edges of the graph $G$, we obtain

$$ \begin{aligned} [Z_{\text{Ising}}^{+}(G, J)]^2 &= 2^{-|V^*|} \left(\prod_{e \in E} e^{2J_e}\right) \left(\prod_{e^* \in E^*} (1 + e^{-4J_e})\right) Z_{\text{dimer}}(G^Q, \nu^j) \\ &= 2^{-|V^*|+|E^*|} \left(\prod_{e \in E^\theta} e^{2J_e}\right) \left(\prod_{e^* \in E^*} \cosh(2J_e)\right) Z_{\text{dimer}}(G^Q, \nu^j) \\ &= 2^{|V|-1} \left(\prod_{e \in E^\theta} \frac{e^{2J_e}}{2}\right) \left(\prod_{e^* \in E^*} \cosh(2J_e)\right) Z_{\text{dimer}}(G^Q, \nu^j), \end{aligned} \quad (2.5) $$

where in the last line we used that $|E| = |E^*| + |E^\theta|$ and Euler's formula: $|E| = |V| + |V^*| - 1$.

## 2.4 Rooted directed spanning forests and directed spanning trees

We also need the model of *rooted directed spanning forests* on the graphs $G$ and $G^*$, resp. $G^\Gamma$ and $G^*$, in the infinite, resp. finite, case. So as to include both the primal and the dual graphs, we now define this model on a simple graph $G = (V, E)$.

Suppose that vertices are assigned non-negative masses, denoted $m = (m_x)_{x \in V}$, and that directed edges have positive conductances, denoted $\rho$, meaning that every directed edge $(x, x')$ has conductance $\rho_{x, x'}$.

A *rooted directed spanning forest* (rdSF) of $G$ is a subset of edges spanning all vertices of the graph, such that each connected component is a directed tree rooted at a vertex of $G$, denoted $x_T$. Let $\mathcal{F}(G)$ denote the set of rdSF of the graph $G$.

Suppose that $G$ is finite and consider the *Boltzmann measure on rdSF*, denoted $P_{\text{rdSF}}$, defined by:

$$ \forall F \in \mathcal{F}(G), \quad P_{\text{rdSF}}(F) = \frac{\prod_{T \in F} (m_{x_T} \prod_{(x,x') \in T} \rho_{x,x'})}{Z_{\text{rdSF}}(G, \rho, m)}, $$

where $Z_{\text{rdSF}}(G, \rho, m) = \sum_{F \in \mathcal{F}(G)} \prod_{T \in F} m_{x_T} \prod_{(x,x') \in T} \rho_{x,x'}$ is the *rdSF partition function*. Whenever conductances are symmetric, i.e., $\rho_{x,x'} = \rho_{x',x}$, we will remove the “d” in rdSF.

As a consequence of the directed version of Kirchhoff's matrix-tree theorem [50, 73], the rdSF partition function is computed using the massive Laplacian operator/matrix as follows. The *massive Laplacian operator* $\Delta^m: \mathbb{C}^V \to \mathbb{C}^V$ is defined by:

$$ \forall F \in \mathbb{C}^V, \quad (\Delta^m F)_x = \sum_{x' \sim x} \rho_{x,x'} (F_x - F_{x'}) + m_x F_x. $$
The operator $\Delta^m$ is represented by a matrix, also denoted $\Delta^m$, whose non-zero coefficients are given by:

$$ \Delta_{x,x'}^m = \begin{cases} -\rho_{x,x'} & \text{if } (x, x') \text{ is an edge of } G \\ m_x + \sum_{x' \sim x} \rho_{x,x'} & \text{if } x' = x \text{ is a vertex of } G. \end{cases} $$

A function $F \in C^V$ is said to be massive-harmonic, if $\Delta^m F = 0$.

Consider the graph $G_\dagger$ constructed from $G$ by adding a cemetery vertex $\dagger$ and an edge $(x, \dagger)$ for every vertex $x$ such that $m_x \neq 0$. Define the modified weight function $\rho^m$ on (directed) edges of $G_\dagger$ by,

$$ \forall \text{ edge } (x, x') \text{ of } G_\dagger, \quad \rho_{x,x'}^m = \begin{cases} \rho_{x,x'} & \text{if } x' \neq \dagger \\ m_x & \text{if } x' = \dagger. \end{cases} $$

There is a natural weight-preserving bijection between $\mathcal{T}^\dagger(G_\dagger)$ and $\mathcal{F}(G)$: a $\dagger$-directed spanning tree of $G_\dagger$ corresponds to the rooted directed spanning forest of $G$ obtained by replacing every edge $(x, \dagger)$ of the dST by a root of the rdSF.

Denote by $\Delta_\dagger$ the (non-massive) Laplacian matrix of $G_\dagger$ with weight function $\rho^m$ on the edges. Then, $\Delta^m$ is the Laplacian matrix $\Delta_\dagger$ from which one has removed the row and column corresponding to the cemetery $\dagger$ and thus, by Kirchhoff's matrix-tree theorem [50, 73], the determinant of $\Delta^m$ counts $\rho^m$ weighted $\dagger$-dST of $G_\dagger$. Using the bijection between $\dagger$-dST of $G_\dagger$ and rdSF of $G$, we thus have,

**Theorem 2.3 ([50, 73]).**

$$ Z_{\text{rdSF}}(G, \rho, m) = Z_{\text{dST}}^{\dagger}(G_{\dagger}, \rho^m) = \det(\Delta^m). $$

When there is at least one vertex with positive mass, the massive Green function, denoted $G^m$, is the inverse of the massive Laplacian $\Delta^m$. Since in the remainder of the paper, graphs are written with the letter $G$ with or without superscripts, we believe that the notation $G^m$ will not create confusion. The massive Green function is naturally related to the expected number of visits of the network random walk associated to the conductances $\rho$ and masses $m$, see for example Appendix D of [11], where a number of facts are recalled.

## 2.5 Isoradial graphs and Z-invariance

Sections 3, 4 and 5.3 use Z-invariant models defined on isoradial graphs. We recall these notions, related concepts and more specifically give the definitions of the Z-invariant versions of the Ising model on $G$, of the dimer model on the decorated graphs $G^F$ and $G^Q$ and of rooted spanning forests on $G$ or $G^*$.

### 2.5.1 Isoradial graphs, diamond graphs and angles

Isoradial graphs naturally appear when considering a discrete version of the Cauchy-Riemann equations, see [33] and [59, 45, 17]; they also arise in Z-invariant models when solving the corresponding Yang-Baxter equations [6, 25]; the name isoradial comes from the paper [45].

Suppose that $G$ is an infinite, planar graph. Then $G$ is said to be isoradial if it can be embedded in the plane in such a way that every face is inscribable in a circle of the same radius, and such that the circumcircles are in the interior of the faces. We consider $G$ as an embedded graph and take the common radius to be 2. Note that the dual $G^*$ of an isoradial graph is also isoradial, an embedding of $G^*$ is obtained by taking as vertices the circumcenters of the circles.
The operator $\Delta^m$ is represented by a matrix, also denoted $\Delta^m$, whose non-zero coefficients are given by:

$$ \Delta_{x,x'}^m = \begin{cases} -\rho_{x,x'} & \text{if } (x, x') \text{ is an edge of } G \\ m_x + \sum_{x' \sim x} \rho_{x,x'} & \text{if } x' = x \text{ is a vertex of } G. \end{cases} $$

A function $F \in C^V$ is said to be massive-harmonic, if $\Delta^m F = 0$.

Consider the graph $G_\dagger$ constructed from $G$ by adding a cemetery vertex $\dagger$ and an edge $(x, \dagger)$ for every vertex $x$ such that $m_x \neq 0$. Define the modified weight function $\rho^m$ on (directed) edges of $G_\dagger$ by,

$$ \forall \text{ edge } (x, x') \text{ of } G_\dagger, \quad \rho_{x,x'}^m = \begin{cases} \rho_{x,x'} & \text{if } x' \neq \dagger \\ m_x & \text{if } x' = \dagger. \end{cases} $$

There is a natural weight-preserving bijection between $\mathcal{T}^\dagger(G_\dagger)$ and $\mathcal{F}(G)$: a $\dagger$-directed spanning tree of $G_\dagger$ corresponds to the rooted directed spanning forest of $G$ obtained by replacing every edge $(x, \dagger)$ of the dST by a root of the rdSF.

Denote by $\Delta_\dagger$ the (non-massive) Laplacian matrix of $G_\dagger$ with weight function $\rho^m$ on the edges. Then, $\Delta^m$ is the Laplacian matrix $\Delta_\dagger$ from which one has removed the row and column corresponding to the cemetery $\dagger$ and thus, by Kirchhoff's matrix-tree theorem [50, 73], the determinant of $\Delta^m$ counts $\rho^m$ weighted $\dagger$-dST of $G_\dagger$. Using the bijection between $\dagger$-dST of $G_\dagger$ and rdSF of $G$, we thus have,

**Theorem 2.3 ([50, 73]).**

$$ Z_{\text{rdSF}}(G, \rho, m) = Z_{\text{dST}}^{\dagger}(G_{\dagger}, \rho^m) = \det(\Delta^m). $$

When there is at least one vertex with positive mass, the massive Green function, denoted $G^m$, is the inverse of the massive Laplacian $\Delta^m$. Since in the remainder of the paper, graphs are written with the letter $G$ with or without superscripts, we believe that the notation $G^m$ will not create confusion. The massive Green function is naturally related to the expected number of visits of the network random walk associated to the conductances $\rho$ and masses $m$, see for example Appendix D of [11], where a number of facts are recalled.

## 2.5 Isoradial graphs and Z-invariance

Sections 3, 4 and 5.3 use Z-invariant models defined on isoradial graphs. We recall these notions, related concepts and more specifically give the definitions of the Z-invariant versions of the Ising model on $G$, of the dimer model on the decorated graphs $G^F$ and $G^Q$ and of rooted spanning forests on $G$ or $G^*$.

### 2.5.1 Isoradial graphs, diamond graphs and angles

Isoradial graphs naturally appear when considering a discrete version of the Cauchy-Riemann equations, see [33] and [59, 45, 17]; they also arise in Z-invariant models when solving the corresponding Yang-Baxter equations [6, 25]; the name isoradial comes from the paper [45].

Suppose that $G$ is an infinite, planar graph. Then $G$ is said to be isoradial if it can be embedded in the plane in such a way that every face is inscribable in a circle of the same radius, and such that the circumcircles are in the interior of the faces. We consider $G$ as an embedded graph and take the common radius to be 2. Note that the dual $G^*$ of an isoradial graph is also isoradial, an embedding of $G^*$ is obtained by taking as vertices the circumcenters of the circles.
This definition also holds when the graph is finite. Recall that in this case, the notation $G$ is used for the graph having an additional vertex on each boundary edge. We now fix the isoradial embedding of $G$ when the original graph (the one without the additional vertices) is isoradial. This is done in the same way as in [18]: each additional boundary vertex corresponds to a boundary edge $xy$ of the original graph and we embed this additional vertex in the middle of the arc joining $x$ and $y$, see Figure 5 (left and right).

<figure><img src="image_7.png" /><figcaption>Figure 5: Left: original isoradial graph (black lines) with circumcircles (grey) and dual vertices embedded as circumcenters (diamonds). Right: isoradial graph $G$ (black lines) with the additional boundary vertices (grey bullets); diamond graph $G^{\diamond}$ (grey lines); rhombus (yellow) assigned to an edge $e = (v, v')$ with the corresponding half-angle $\tilde{\theta}_e$ and rhombus vectors $2e^{i\tilde{\alpha}_e}, 2e^{i\tilde{\beta}_e}$; a boundary rhombus pair of $\mathcal{R}^{\theta}$ (light grey).</figcaption></figure>

In the infinite case, the *diamond graph*, denoted $G^{\diamond}$, is constructed from an isoradial graph $G$ and its dual $G^*$ as follows: its vertex set is $V \cup V^*$, the vertices of $G$ and $G^*$; and each dual vertex is joined to all vertices bounding the face it corresponds to. Since the graph $G$ is isoradial, faces of the diamond graph $G^{\diamond}$ are side-length-2 rhombi.

There is a bijection between rhombi of $G^{\diamond}$ and pairs $e, e^*$ of primal and dual edges, the latter being the two diagonals of the rhombi. To every edge $e$, one assigns an *angle* $\tilde{\theta}_e \in (0, \frac{\pi}{2})$ defined to be the half-angle of the corresponding rhombus at the edge $e$. We furthermore ask that $\tilde{\theta}_e \in (\varepsilon, \frac{\pi}{2} - \varepsilon)$, for some $\varepsilon > 0$. The rhombus angle of the dual edge $e^*$ is $\tilde{\theta}_{e^*} = \frac{\pi}{2} - \tilde{\theta}_e := \tilde{\theta}_e^*$. To a directed edge $e = (v, v')$ of $G$ we further assign two rhombus vectors $2e^{i\tilde{\alpha}_e}, 2e^{i\tilde{\beta}_e}$ of $G^{\diamond}$, such that $2e^{i\tilde{\alpha}_e}$ is on the right of the edge $(v, v')$, see Figure 5 (right). The angles $\tilde{\alpha}_e$ and $\tilde{\beta}_e$ are defined so that $\frac{\tilde{\beta}_e - \tilde{\alpha}_e}{2} = \tilde{\theta}_e$. Whenever no confusion occurs, we remove the subscript $e$ from the notation. In the finite case, the *diamond graph*, also denoted $G^{\diamond}$, is constructed in a similar way from $G$ and its restricted dual $G^*$. One then adds the missing half-rhombi along the boundary; they may overlap but this causes no problem, see Figure 5 (right). Angles and rhombus vectors assigned to edges are defined in the same way.

Because of the additional vertex on each boundary edge and because of our choice of embedding, rhombi along the boundary of $G^{\diamond}$ come in pairs, both having the same rhombus half-angle; let us denote by $\mathcal{R}^{\theta}$ the set of boundary rhombus pairs, an instance is highlighted in Figure 5 (right, light grey).
This definition also holds when the graph is finite. Recall that in this case, the notation $G$ is used for the graph having an additional vertex on each boundary edge. We now fix the isoradial embedding of $G$ when the original graph (the one without the additional vertices) is isoradial. This is done in the same way as in [18]: each additional boundary vertex corresponds to a boundary edge $xy$ of the original graph and we embed this additional vertex in the middle of the arc joining $x$ and $y$, see Figure 5 (left and right).

<figure><img src="image_7.png" /><figcaption>Figure 5: Left: original isoradial graph (black lines) with circumcircles (grey) and dual vertices embedded as circumcenters (diamonds). Right: isoradial graph $G$ (black lines) with the additional boundary vertices (grey bullets); diamond graph $G^{\diamond}$ (grey lines); rhombus (yellow) assigned to an edge $e = (v, v')$ with the corresponding half-angle $\tilde{\theta}_e$ and rhombus vectors $2e^{i\tilde{\alpha}_e}, 2e^{i\tilde{\beta}_e}$; a boundary rhombus pair of $\mathcal{R}^{\theta}$ (light grey).</figcaption></figure>

In the infinite case, the *diamond graph*, denoted $G^{\diamond}$, is constructed from an isoradial graph $G$ and its dual $G^*$ as follows: its vertex set is $V \cup V^*$, the vertices of $G$ and $G^*$; and each dual vertex is joined to all vertices bounding the face it corresponds to. Since the graph $G$ is isoradial, faces of the diamond graph $G^{\diamond}$ are side-length-2 rhombi.

There is a bijection between rhombi of $G^{\diamond}$ and pairs $e, e^*$ of primal and dual edges, the latter being the two diagonals of the rhombi. To every edge $e$, one assigns an *angle* $\tilde{\theta}_e \in (0, \frac{\pi}{2})$ defined to be the half-angle of the corresponding rhombus at the edge $e$. We furthermore ask that $\tilde{\theta}_e \in (\epsilon, \frac{\pi}{2} - \epsilon)$, for some $\epsilon > 0$. The rhombus angle of the dual edge $e^*$ is $\tilde{\theta}_{e^*} = \frac{\pi}{2} - \tilde{\theta}_e := \tilde{\theta}_e^*$. To a directed edge $e = (v, v')$ of $G$ we further assign two rhombus vectors $2e^{i\tilde{\alpha}_e}, 2e^{i\tilde{\beta}_e}$ of $G^{\diamond}$, such that $2e^{i\tilde{\alpha}_e}$ is on the right of the edge $(v, v')$, see Figure 5 (right). The angles $\tilde{\alpha}_e$ and $\tilde{\beta}_e$ are defined so that $\frac{\tilde{\beta}_e - \tilde{\alpha}_e}{2} = \tilde{\theta}_e$. Whenever no confusion occurs, we remove the subscript $e$ from the notation. In the finite case, the *diamond graph*, also denoted $G^{\diamond}$, is constructed in a similar way from $G$ and its restricted dual $G^*$. One then adds the missing half-rhombi along the boundary; they may overlap but this causes no problem, see Figure 5 (right). Angles and rhombus vectors assigned to edges are defined in the same way.

Because of the additional vertex on each boundary edge and because of our choice of embedding, rhombi along the boundary of $G^{\diamond}$ come in pairs, both having the same rhombus half-angle; let us denote by $\mathcal{R}^{\theta}$ the set of boundary rhombus pairs, an instance is highlighted in Figure 5 (right, light grey).
2.5.2 Isoradial embeddings of the decorated graphs $G^D$ and $G^Q$

Consider an isoradial graph $G$, its dual $G^*$ in the infinite case and its restricted dual $G^*$ in the finite case. The double graph $G^D$ is embedded so that the black vertices are those of $G$ and $G^*$ and the white vertices are at the crossing of the diagonals of the rhombi of $G^\circ$, see Figure 6 (left).

<figure><img src="image_8.png" /><figcaption>Figure 6: Left: isoradial embedding of the graph $G^D$ (black lines). Right: isoradial embedding of the graph $G^Q$ (black lines). In both cases is also pictured the diamond graph $G^{\diamond/2}$ (grey lines), a boundary rhombus pair of $\mathcal{R}^\theta$ and the root pair (light grey).</figcaption></figure>

In the infinite case and in the finite non-boundary case, consider the embedding of the bipartite graph $G^Q$ where external edges have length-0 and their endpoints become a single vertex in the middle of the rhombus edges of $G^\circ$, see Figure 6 (right, inner vertices); then, inner quadrangles of $G^Q$ are rectangles. Note that although external edges are embedded as single vertices, they still consist of two vertices joined by (a length-0) edge, i.e., the combinatorics of the graph does not change.

When the graph $G^Q$ is *finite*, the procedure along the boundary is different. Consider a boundary rhombus pairs of $\mathcal{R}^\theta$, the following notation will be used throughout the paper and is illustrated in Figure 7 below. Let $v^\ell, v^c, v^r$ be the vertices of $G$ in cw order and let $f^c$ be the vertex of $G^*$; note that $v^c$ is the additional vertex on the edge $v^\ell v^r$ of the original graph. Denote by $w^\ell, w^r$ the white vertices of the double graph $G^D$ and by $w^\ell, b^\ell, w^c, b^r, w^r$ the vertices of $G^Q$. Then, taking the same convention as in the infinite case for the embedding gives Figure 7 (left); but it turns out that the appropriate embedding to obtain Theorem 4.2 is that of Figure 7 (center), see also Figure 6 (right), where the boundary quadrangle edge $b^\ell w^c$ has length-0 and is "replaced" by the external edge $w^\ell b^\ell$. This change of embedding preserves the combinatorics of the graph; it has the effect of exchanging the colors of the bipartite coloring of $G^Q$ in the left rhombus of the rhombus pair.

We will often be using the fact that vertices/edges of the boundary rhombus pairs of $\mathcal{R}^\theta$ encode: boundary vertices/edges of $G$, boundary vertices of the restricted dual $G^*$, where a *boundary vertex* of $G^*$ is defined to be a vertex adjacent to the vertex $o$ in $\bar{G}^*$; boundary quadrangle vertices/edges of $G^Q$, where recall that boundary quadrangles are degenerate and reduced to edges in bijection with boundary edges of $G$. We will use the notation $\{v^\ell \in \mathcal{R}^\theta\}$ for the set of vertices of type $v^\ell$ belonging to boundary rhombus pairs, and similarly for other vertices or edges of $\mathcal{R}^\theta$.
### 2.5.2 Isoradial embeddings of the decorated graphs $G^D$ and $G^Q$

Consider an isoradial graph $G$, its dual $G^*$ in the infinite case and its restricted dual $G^*$ in the finite case. The double graph $G^D$ is embedded so that the black vertices are those of $G$ and $G^*$ and the white vertices are at the crossing of the diagonals of the rhombi of $G^\circ$, see Figure 6 (left).

<figure><img src="image_8.png" /><figcaption>Figure 6: Left: isoradial embedding of the graph $G^D$ (black lines). Right: isoradial embedding of the graph $G^Q$ (black lines). In both cases is also pictured the diamond graph $G^{\diamond/2}$ (grey lines), a boundary rhombus pair of $\mathcal{R}^\theta$ and the root pair (light grey).</figcaption></figure>

In the infinite case and in the finite non-boundary case, consider the embedding of the bipartite graph $G^Q$ where external edges have length-0 and their endpoints become a single vertex in the middle of the rhombus edges of $G^\circ$, see Figure 6 (right, inner vertices); then, inner quadrangles of $G^Q$ are rectangles. Note that although external edges are embedded as single vertices, they still consist of two vertices joined by (a length-0) edge, i.e., the combinatorics of the graph does not change.

When the graph $G^Q$ is *finite*, the procedure along the boundary is different. Consider a boundary rhombus pairs of $\mathcal{R}^\theta$, the following notation will be used throughout the paper and is illustrated in Figure 7 below. Let $v^\ell, v^c, v^r$ be the vertices of $G$ in cw order and let $f^c$ be the vertex of $G^*$; note that $v^c$ is the additional vertex on the edge $v^\ell v^r$ of the original graph. Denote by $w^\ell, w^r$ the white vertices of the double graph $G^D$ and by $w^\ell, b^\ell, w^c, b^r, w^r$ the vertices of $G^Q$. Then, taking the same convention as in the infinite case for the embedding gives Figure 7 (left); but it turns out that the appropriate embedding to obtain Theorem 4.2 is that of Figure 7 (center), see also Figure 6 (right), where the boundary quadrangle edge $b^\ell w^c$ has length-0 and is "replaced" by the external edge $w^\ell b^\ell$. This change of embedding preserves the combinatorics of the graph; it has the effect of exchanging the colors of the bipartite coloring of $G^Q$ in the left rhombus of the rhombus pair.

We will often be using the fact that vertices/edges of the boundary rhombus pairs of $\mathcal{R}^\theta$ encode: boundary vertices/edges of $G$, boundary vertices of the restricted dual $G^*$, where a *boundary vertex* of $G^*$ is defined to be a vertex adjacent to the vertex $o$ in $\bar{G}^*$; boundary quadrangle vertices/edges of $G^Q$, where recall that boundary quadrangles are degenerate and reduced to edges in bijection with boundary edges of $G$. We will use the notation $\{v^\ell \in \mathcal{R}^\theta\}$ for the set of vertices of type $v^\ell$ belonging to boundary rhombus pairs, and similarly for other vertices or edges of $\mathcal{R}^\theta$.
<figure><img src="image_9.png" /><figcaption>Figure 7: Notation for vertices of boundary pairs of rhombi of $\mathcal{R}^\theta$. Left: isoradial embedding of $G^Q$ under the convention that external edges have length 0. Center: isoradial embedding used in this paper. Right: rhombus vectors of $G^{\diamond/2}$ and half-angle assigned to the edges $(v^c, w^\ell)$, $(v^c, w^r)$ of $G^D$.</figcaption></figure>

The embeddings of $G^D$ and $G^Q$ are both isoradial with circumcircles having common radius 1. Consider the graph obtained from the diamond graph $G^\diamond$ by cutting rhombi into four equal length-1 rhombi. Denote this graph by $G^{\diamond/2}$ in the infinite case and, in the finite case, let $G^{\diamond/2}$ be this graph where the boundary quarter rhombi crossed by no edge of $G^D$ are removed. Then $G^{\diamond/2}$ is the diamond graph of $G^D$. Note that $G^{\diamond/2}$ is nearly the diamond graph of $G^Q$: it is slightly extended along the boundary and one should think of it as having flat rhombi associated to length-0 edges of $G^Q$. We will nevertheless refer to it as the *diamond graph* of $G^D$ or $G^Q$. An example of graph $G^{\diamond/2}$ is given in Figure 6 (left and right, grey).

Consider a boundary rhombus pair of $\mathcal{R}^\theta$, and let $e^{i\bar{\alpha}^\ell}$, $e^{i\bar{\beta}^\ell}$, resp. $e^{i\bar{\alpha}^r}$, $e^{i\bar{\beta}^r}$, be the rhombus vectors of the diamond graph $G^{\diamond/2}$ assigned to the edge $(v^c, w^\ell)$, resp. $(v^c, w^r)$, of $G^D$ see Figure 7 (right). By definition we have $\frac{\bar{\beta}^\ell - \bar{\alpha}^\ell}{2}$, $\frac{\bar{\beta}^r - \bar{\alpha}^r}{2} \in (\varepsilon, \frac{\pi}{2} - \varepsilon)$, and by construction the two rhombi have the same half-angle denoted $\bar{\theta}^\theta = \frac{\bar{\beta}^\ell - \bar{\alpha}^\ell}{2} = \frac{\bar{\beta}^r - \bar{\alpha}^r}{2}$. We further impose that $\bar{\beta}^\ell = \bar{\alpha}^r + 2\pi$ or equivalently that $\frac{\bar{\alpha}^\ell - \bar{\beta}^r}{2} \in (2\varepsilon, \pi - 2\varepsilon)$.

Amongst the boundary rhombus pairs of $\mathcal{R}^\theta$ the one containing the fixed vertex $r$, i.e. the one for which $v^c = r$, plays a special role; it will be referred to as the *root-boundary rhombus pair* or simply *root-pair*, see Figure 6 where the root pair is highlighted in light grey. We denote by $\mathcal{R}^{\theta,r}$ the set $\mathcal{R}^\theta$ without the root pair. The isoradial embedding of the graph $G^{D,r}$ is obtained from $G^D$ by removing the vertex $r$ and the edges $w^\ell r, w^r r$ of the root pair. Whenever needed, we add a superscript $r$ to the notation of Figure 7 to specify vertices of the root pair.

### 2.5.3 Train-tracks

A *train-track* of a finite isoradial graph $G$, also known as a *de Bruijn line* [26, 27] or a *rapidity line* [5] is a maximal chain of edge-adjacent rhombi of the diamond graph $G^\diamond$, such that when entering a rhombus one exits along the opposite edge [48]; each train-track $\tau$ has a parallel direction $\pm 2e^{i\bar{\alpha}_\tau}$. Consider the simply connected domain $D(G)$ obtained by taking the union of the faces of $G$. Then a train-track $\tau$ enters and exits $D(G)$, and there are exactly two parallel edges of $\tau$ outside of $D(G)$, see Figure 5, (right). The boundary rhombus vectors $\{2e^{i\bar{\alpha}^\ell}, 2e^{i\bar{\beta}^r} \in \mathcal{R}^\theta\}$, $\{2e^{i\bar{\alpha}^r}, 2e^{i\bar{\beta}^\ell} \in \mathcal{R}^\theta\}$ come in parallel pairs, and all the parallel directions of the train-tracks are encoded in $\{\pm 2e^{i\bar{\alpha}^\ell}, \pm 2e^{i\bar{\beta}^\ell} \in \mathcal{R}^\theta\} = \{\pm 2e^{i\bar{\alpha}^r}, \pm 2e^{i\bar{\beta}^r} \in \mathcal{R}^\theta\}$.

### 2.5.4 Elliptic angles

Consider an *elliptic modulus* $k$, then $k' = (1 - k^2)^{\frac{1}{2}}$ is the *complementary elliptic modulus*. Suppose that $k$ is such that $(k')^2 \in (0, \infty)$. The *complete elliptic integral of the first kind*,
<figure><img src="image_9.png" /><figcaption>Figure 7: Notation for vertices of boundary pairs of rhombi of $\mathcal{R}^\theta$. Left: isoradial embedding of $G^Q$ under the convention that external edges have length 0. Center: isoradial embedding used in this paper. Right: rhombus vectors of $G^{\diamond/2}$ and half-angle assigned to the edges $(v^c, w^\ell)$, $(v^c, w^r)$ of $G^D$.</figcaption></figure>

The embeddings of $G^D$ and $G^Q$ are both isoradial with circumcircles having common radius 1. Consider the graph obtained from the diamond graph $G^\diamond$ by cutting rhombi into four equal length-1 rhombi. Denote this graph by $G^{\diamond/2}$ in the infinite case and, in the finite case, let $G^{\diamond/2}$ be this graph where the boundary quarter rhombi crossed by no edge of $G^D$ are removed. Then $G^{\diamond/2}$ is the diamond graph of $G^D$. Note that $G^{\diamond/2}$ is nearly the diamond graph of $G^Q$: it is slightly extended along the boundary and one should think of it as having flat rhombi associated to length-0 edges of $G^Q$. We will nevertheless refer to it as the *diamond graph* of $G^D$ or $G^Q$. An example of graph $G^{\diamond/2}$ is given in Figure 6 (left and right, grey).

Consider a boundary rhombus pair of $\mathcal{R}^\theta$, and let $e^{i\bar{\alpha}^\ell}$, $e^{i\bar{\beta}^\ell}$, resp. $e^{i\bar{\alpha}^r}$, $e^{i\bar{\beta}^r}$, be the rhombus vectors of the diamond graph $G^{\diamond/2}$ assigned to the edge $(v^c, w^\ell)$, resp. $(v^c, w^r)$, of $G^D$ see Figure 7 (right). By definition we have $\frac{\bar{\beta}^\ell - \bar{\alpha}^\ell}{2}$, $\frac{\bar{\beta}^r - \bar{\alpha}^r}{2} \in (\varepsilon, \frac{\pi}{2} - \varepsilon)$, and by construction the two rhombi have the same half-angle denoted $\bar{\theta}^\theta = \frac{\bar{\beta}^\ell - \bar{\alpha}^\ell}{2} = \frac{\bar{\beta}^r - \bar{\alpha}^r}{2}$. We further impose that $\bar{\beta}^\ell = \bar{\alpha}^r + 2\pi$ or equivalently that $\frac{\bar{\alpha}^\ell - \bar{\beta}^r}{2} \in (2\varepsilon, \pi - 2\varepsilon)$.

Amongst the boundary rhombus pairs of $\mathcal{R}^\theta$ the one containing the fixed vertex $r$, i.e. the one for which $v^c = r$, plays a special role; it will be referred to as the *root-boundary rhombus pair* or simply *root-pair*, see Figure 6 where the root pair is highlighted in light grey. We denote by $\mathcal{R}^{\theta,r}$ the set $\mathcal{R}^\theta$ without the root pair. The isoradial embedding of the graph $G^{D,r}$ is obtained from $G^D$ by removing the vertex $r$ and the edges $w^\ell r, w^r r$ of the root pair. Whenever needed, we add a superscript $r$ to the notation of Figure 7 to specify vertices of the root pair.

### 2.5.3 Train-tracks

A *train-track* of a finite isoradial graph $G$, also known as a *de Bruijn line* [26, 27] or a *rapidity line* [5] is a maximal chain of edge-adjacent rhombi of the diamond graph $G^\diamond$, such that when entering a rhombus one exits along the opposite edge [48]; each train-track $\tau$ has a parallel direction $\pm 2e^{i\bar{\alpha}_\tau}$. Consider the simply connected domain $D(G)$ obtained by taking the union of the faces of $G$. Then a train-track $\tau$ enters and exits $D(G)$, and there are exactly two parallel edges of $\tau$ outside of $D(G)$, see Figure 5, (right). The boundary rhombus vectors $\{2e^{i\bar{\alpha}^\ell}, 2e^{i\bar{\beta}^r} \in \mathcal{R}^\theta\}$, $\{2e^{i\bar{\alpha}^r}, 2e^{i\bar{\beta}^\ell} \in \mathcal{R}^\theta\}$ come in parallel pairs, and all the parallel directions of the train-tracks are encoded in $\{\pm 2e^{i\bar{\alpha}^\ell}, \pm 2e^{i\bar{\beta}^\ell} \in \mathcal{R}^\theta\} = \{\pm 2e^{i\bar{\alpha}^r}, \pm 2e^{i\bar{\beta}^r} \in \mathcal{R}^\theta\}$.

### 2.5.4 Elliptic angles

Consider an *elliptic modulus* $k$, then $k' = (1-k^2)^{\frac{1}{2}}$ is the *complementary elliptic modulus*. Suppose that $k$ is such that $(k')^2 \in (0, \infty)$. The *complete elliptic integral of the first kind*,
denoted $K = K(k)$ is

$$K = \int_{0}^{\frac{\pi}{2}} \frac{1}{(1 - k^2 \sin \tau)} d\tau,$$

and for later purposes we also need $K' = K(k')$. As in [11, 12], we need the following linear transformation of rhombus angles and vectors of the diamond graph $G^{\diamond}$ associated to edges:

$$\theta = \bar{\theta} \frac{2K}{\pi}, \quad \alpha = \bar{\alpha} \frac{2K}{\pi}, \quad \beta = \bar{\beta} \frac{2K}{\pi}.$$

### 2.5.5 Z-invariant Ising model and corresponding dimer models

Underlying Z-invariance is the star-triangle transformation, also known as the $Y-\Delta$ move, on isoradial graphs. Suppose that an isoradial graph $G$ has a triangle, then this triangle can be transformed into a three-legged star while preserving isoradiality. This amounts to performing a cubic flip in the underlying diamond graph $G^{\diamond}$; the embedding of the additional vertex of the triangle is given by the cubic flip, see Figure 8.

<figure><img src="image_10.png" /><figcaption>Figure 8: A star-triangle transformation of an isoradial graph.</figcaption></figure>

Z-invariance [4, 5] phrased in the context of the Ising model requires that when decomposing the partition function according to the $2^3$ possible spin configurations at the three vertices bounding the star/triangle, it only changes by an overall constant when performing a $Y-\Delta$ move, and this constant is independent of the choice of spin configuration. This yields a set of equations for the coupling constants, known as the Ising Yang-Baxter equations, see also [65]. Extending the form of the solutions to the whole of the graph naturally leads to introducing isoradial graphs: the solution is parameterized by the rhombus half-angles assigned to edges and by the elliptic modulus $k$, where $k$ is such that $(k')^2 = 1 - k^2 \in (0, \infty)$, see [6]. The Z-invariant coupling constants are explicitly given by:

$$\forall e \in E, \quad J_e = \frac{1}{2} \ln \left( \frac{1 + \operatorname{sn}(\theta_e | k)}{\operatorname{cn}(\theta_e | k)} \right), \qquad (2.6)$$

where $\operatorname{sn}$, $\operatorname{cn}$ are two of the twelve Jacobi elliptic trigonometric functions. We refer the reader to [1, 54] for more on elliptic and related functions.

Suppose that the Z-invariant coupling constants are chosen for the Ising model. Then, the dimer weight function $\mu^j$ on the Fisher graph $G^F$ arising from Fisher's correspondence is given by, for every edge $e$ of $G^F$,

$$\mu_e^j = \begin{cases} 1 & \text{if } e \text{ is an internal edge} \\ \frac{\operatorname{cn}(\theta_e | k)}{1 + \operatorname{sn}(\theta_e | k)} & \text{if } e \text{ is an external edge arising from an edge } e^* \text{ of } G^*. \end{cases} \qquad (2.7)$$

The dimer weight function $\nu^j$ on the bipartite graph $G^Q$ arising from the XOR-Ising
denoted $K = K(k)$ is

$$K = \int_{0}^{\frac{\pi}{2}} \frac{1}{(1 - k^2 \sin \tau)} d\tau,$$

and for later purposes we also need $K' = K(k')$. As in [11, 12], we need the following linear transformation of rhombus angles and vectors of the diamond graph $G^{\diamond}$ associated to edges:

$$\theta = \bar{\theta} \frac{2K}{\pi}, \quad \alpha = \bar{\alpha} \frac{2K}{\pi}, \quad \beta = \bar{\beta} \frac{2K}{\pi}.$$

### 2.5.5 Z-invariant Ising model and corresponding dimer models

Underlying Z-invariance is the star-triangle transformation, also known as the $Y-\Delta$ move, on isoradial graphs. Suppose that an isoradial graph $G$ has a triangle, then this triangle can be transformed into a three-legged star while preserving isoradiality. This amounts to performing a cubic flip in the underlying diamond graph $G^{\diamond}$; the embedding of the additional vertex of the triangle is given by the cubic flip, see Figure 8.

<figure><img src="image_10.png" /><figcaption>Figure 8: A star-triangle transformation of an isoradial graph.</figcaption></figure>

Z-invariance [4, 5] phrased in the context of the Ising model requires that when decomposing the partition function according to the $2^3$ possible spin configurations at the three vertices bounding the star/triangle, it only changes by an overall constant when performing a $Y-\Delta$ move, and this constant is independent of the choice of spin configuration. This yields a set of equations for the coupling constants, known as the Ising Yang-Baxter equations, see also [65]. Extending the form of the solutions to the whole of the graph naturally leads to introducing isoradial graphs: the solution is parameterized by the rhombus half-angles assigned to edges and by the elliptic modulus $k$, where $k$ is such that $(k')^2 = 1 - k^2 \in (0, \infty)$, see [6]. The Z-invariant coupling constants are explicitly given by:

$$\forall e \in E, \quad J_e = \frac{1}{2} \ln \left( \frac{1 + \operatorname{sn}(\theta_e | k)}{\operatorname{cn}(\theta_e | k)} \right), \qquad (2.6)$$

where $\operatorname{sn}$, $\operatorname{cn}$ are two of the twelve Jacobi elliptic trigonometric functions. We refer the reader to [1, 54] for more on elliptic and related functions.

Suppose that the Z-invariant coupling constants are chosen for the Ising model. Then, the dimer weight function $\mu^j$ on the Fisher graph $G^F$ arising from Fisher's correspondence is given by, for every edge $e$ of $G^F$,

$$\mu_e^j = \begin{cases} 1 & \text{if } e \text{ is an internal edge} \\ \frac{\operatorname{cn}(\theta_e | k)}{1 + \operatorname{sn}(\theta_e | k)} & \text{if } e \text{ is an external edge arising from an edge } e^* \text{ of } G^*. \end{cases} \qquad (2.7)$$

The dimer weight function $\nu^j$ on the bipartite graph $G^Q$ arising from the XOR-Ising
model is given by, for every edge e of $G^Q$,

$$ \nu_e^j = \begin{cases} 1 & \text{if e is an external edge} \\ 1 & \text{if e is a boundary quadrangle edge in the finite case} \\ \operatorname{sn}(\theta_e|k) & \text{if e is a quadrangle edge/non-boundary quadrangle edge} \\ \operatorname{cn}(\theta_e|k) & \text{in the infinite/finite case, parallel to an edge e of G} \\ \end{cases} \quad (2.8) $$

Recall that the dimer model on $G^Q$ corresponds to a free-fermion, zero-field 6-vertex model [62, 76, 32]. Then, the weights (2.8) also correspond to the 6-vertex version of the elliptic parameterization of the Z-invariant 8-vertex model weights, computed by Baxter [4, 6].

### 2.5.6 The Z-invariant massive Laplacian

In the paper [11], we consider an infinite, isoradial graph G and introduce conductances and masses defining the Z-invariant massive Laplacian operator or simply Z-massive Laplacian, related to Z-invariant rooted spanning forests. Recall that every edge $(v, v')$ of G is assigned two rhombus vectors $2e^{i\bar{\alpha}}$, $2e^{i\bar{\beta}}$ and a half-angle $\bar{\theta}$ of the diamond graph $G^\diamond$. Denote by $v_1, \dots, v_d$ the neighbors of a vertex v of degree d, and for every edge $(v, v_j)$ use the notation $2e^{i\bar{\alpha}_j}$, $2e^{i\bar{\alpha}_{j+1}}$, $\theta_j$ for the associated rhombus vectors and half-angle. Fix an elliptic modulus k such that $(k')^2 \in (0, \infty)$. Then, for every edge $(v, v')$ and every vertex v of G, the conductances $\rho^k$ and masses $m^k$ of [11] are defined by:

$$ \rho_{v,v'}^k = \operatorname{sc}(\theta|k), \quad m_v^k = \sum_{j=1}^d [A(\theta_j|k) - \operatorname{sc}(\theta_j|k)], $$

where $\operatorname{sc} = \frac{\operatorname{sn}}{\operatorname{cn}}$, and

$$ A(u|k) = (k')^{-1} \left( \operatorname{Dc}(u|k) + \frac{E-K}{K} u \right), $$

$\operatorname{Dc}(u|k) = \int_0^u \operatorname{dc}^2(v|k)dv$, and $E = E(k) = \int_0^{\frac{\pi}{2}} (1 - k^2 \sin^2(\tau))^{\frac{1}{2}} d\tau$, is the complete elliptic integral of the second kind. Note that $m^k \ge 0$, by Proposition 6 of [11].

The corresponding Z-massive Laplacian matrix $\Delta^{m(k)}$ has non-zero coefficients given by:

$$ \forall v, v' \in V, \quad \Delta_{v,v'}^{m(k)} = \begin{cases} -\operatorname{sc}(\theta|k) & \text{if } (v, v') \text{ is an edge of G} \\ \sum_{j=1}^d A(\theta_j|k) & \text{if } v' = v. \end{cases} \quad (2.9) $$

The dual graph $G^*$ is also isoradial, we denote by $\rho^{k,*}$, $m^{k,*}$ the associated conductances and masses, and by $\bar{\theta}^*$ the half-angle of a dual edge. We let $\Delta^{m(k),*}$ be the corresponding Z-massive Laplacian operator and refer to it as the dual Z-massive Laplacian.

The inverse of the Z-massive Laplacian $\Delta^{m(k)}$ is the Z-massive Green function; it is denoted $G^{m(k)}$. In [11], we prove the following explicit local expression for coefficients of $G^{m(k)}$: for every pair of vertices x, y of G,

$$ G^{m(k)} = \frac{k'}{4\pi} \int_{\Gamma_{x,y}} e_{(x,y)}(u|k)du, \quad (2.10) $$

where $\Gamma_{x,y}$ is a vertical contour on the torus $\mathcal{T}(k) = \mathbb{C}/(4K\mathbb{Z} + 4iK'\mathbb{Z})$, and $e(\cdot|k): V \times V \times C \to C$ is the massive exponential function defined in [11]. To compute $e_{(x,y)}(u|k)$, consider a path $x = x_1, \dots, x_n = y$ of the diamond graph $G^\diamond$ from $x$ to $y$, let $2e^{i\bar{\alpha}_j}$ be the rhombus vector corresponding to the edge $x_j x_{j+1}$, then

$$ e_{(x_j, x_{j+1})}(u|k) = i(k')^{\frac{1}{2}} \operatorname{sc}(u_{\alpha_j}|k), \quad \text{and} \quad e_{(x,y)}(u|k) = \prod_{j=1}^{n-1} e_{(x_j, x_{j+1})}(u|k), \quad (2.11) $$
model is given by, for every edge e of $G^Q$,

$$ \nu_e^j = \begin{cases} 1 & \text{if e is an external edge} \\ 1 & \text{if e is a boundary quadrangle edge in the finite case} \\ \operatorname{sn}(\theta_e|k) & \text{if e is a quadrangle edge/non-boundary quadrangle edge} \\ \operatorname{cn}(\theta_e|k) & \text{in the infinite/finite case, parallel to an edge e of G} \\ \end{cases} \quad (2.8) $$

Recall that the dimer model on $G^Q$ corresponds to a free-fermion, zero-field 6-vertex model [62, 76, 32]. Then, the weights (2.8) also correspond to the 6-vertex version of the elliptic parameterization of the Z-invariant 8-vertex model weights, computed by Baxter [4, 6].

### 2.5.6 The Z-invariant massive Laplacian

In the paper [11], we consider an infinite, isoradial graph G and introduce conductances and masses defining the Z-invariant massive Laplacian operator or simply Z-massive Laplacian, related to Z-invariant rooted spanning forests. Recall that every edge $(v, v')$ of G is assigned two rhombus vectors $2e^{i\bar{\alpha}}$, $2e^{i\bar{\beta}}$ and a half-angle $\bar{\theta}$ of the diamond graph $G^\diamond$. Denote by $v_1, \dots, v_d$ the neighbors of a vertex v of degree d, and for every edge $(v, v_j)$ use the notation $2e^{i\bar{\alpha}_j}$, $2e^{i\bar{\alpha}_{j+1}}$, $\theta_j$ for the associated rhombus vectors and half-angle. Fix an elliptic modulus k such that $(k')^2 \in (0, \infty)$. Then, for every edge $(v, v')$ and every vertex v of G, the conductances $\rho^k$ and masses $m^k$ of [11] are defined by:

$$ \rho_{v,v'}^k = \operatorname{sc}(\theta|k), \quad m_v^k = \sum_{j=1}^d [A(\theta_j|k) - \operatorname{sc}(\theta_j|k)], $$

where $\operatorname{sc} = \frac{\operatorname{sn}}{\operatorname{cn}}$, and

$$ A(u|k) = (k')^{-1} \left( \operatorname{Dc}(u|k) + \frac{E-K}{K} u \right), $$

$\operatorname{Dc}(u|k) = \int_0^u \operatorname{dc}^2(v|k)dv$, and $E = E(k) = \int_0^{\frac{\pi}{2}} (1 - k^2 \sin^2(\tau))^{\frac{1}{2}} d\tau$, is the complete elliptic integral of the second kind. Note that $m^k \ge 0$, by Proposition 6 of [11].

The corresponding Z-massive Laplacian matrix $\Delta^{m(k)}$ has non-zero coefficients given by:

$$ \forall v, v' \in V, \quad \Delta_{v,v'}^{m(k)} = \begin{cases} -\operatorname{sc}(\theta|k) & \text{if } (v, v') \text{ is an edge of G} \\ \sum_{j=1}^d A(\theta_j|k) & \text{if } v' = v. \end{cases} \quad (2.9) $$

The dual graph $G^*$ is also isoradial, we denote by $\rho^{k,*}$, $m^{k,*}$ the associated conductances and masses, and by $\bar{\theta}^*$ the half-angle of a dual edge. We let $\Delta^{m(k),*}$ be the corresponding Z-massive Laplacian operator and refer to it as the dual Z-massive Laplacian.

The inverse of the Z-massive Laplacian $\Delta^{m(k)}$ is the Z-massive Green function; it is denoted $G^{m(k)}$. In [11], we prove the following explicit local expression for coefficients of $G^{m(k)}$: for every pair of vertices x, y of G,

$$ G^{m(k)} = \frac{k'}{4\pi} \int_{\Gamma_{x,y}} e_{(x,y)}(u|k)du, \quad (2.10) $$

where $\Gamma_{x,y}$ is a vertical contour on the torus $\mathbb{T}(k) = \mathbb{C}/(4K\mathbb{Z} + 4iK'\mathbb{Z})$, and $e(\cdot|k): V \times V \times C \to C$ is the massive exponential function defined in [11]. To compute $e_{(x,y)}(u|k)$, consider a path $x = x_1, \dots, x_n = y$ of the diamond graph $G^\diamond$ from $x$ to $y$, let $2e^{i\bar{\alpha}_j}$ be the rhombus vector corresponding to the edge $x_j x_{j+1}$, then

$$ e_{(x_j, x_{j+1})}(u|k) = i(k')^{\frac{1}{2}} \operatorname{sc}(u_{\alpha_j}|k), \quad \text{and} \quad e_{(x,y)}(u|k) = \prod_{j=1}^{n-1} e_{(x_j, x_{j+1})}(u|k), \quad (2.11) $$
where $u_{\alpha_j} := \frac{u-\alpha_j}{2}$.

Up to an explicit multiplicative constant $G^{m(k)}(x, y)$ is the expected number of visits to $y$ of the associated massive random walk on the infinite graph $G$ started at $x$, see for example Appendix D4 of [11].

From now on, we consider $k$ such that $(k')^2 \in (0, \infty)$ as fixed and omit all reference to $k$ in the notation.

### 2.5.7 Bipartite dimer models on isoradial graphs

When considering a dimer model on a bipartite isoradial graph $G = (W \cup B, E)$ with weight function $\nu$, instead of a Kasteleyn orientation, one can multiply edge-weights by a complex phase [53, 45]. This defines the complex, bipartite Kasteleyn matrix, denoted $\tilde{K}$ in the context of this section, whose non-zero coefficients are given by:

$$ \forall \text{ edge } w, b \text{ of } G, \quad \tilde{K}_{w,b} = e^{i \frac{\bar{\alpha}_e + \bar{\beta}_e}{2}} \nu_{wb}, $$

where $e^{i\bar{\alpha}_e}, e^{i\bar{\beta}_e}$ are the rhombus vectors of $G^\diamond$ associated to the edge $e = (w, b)$. The real and complex bipartite Kasteleyn matrices satisfy the alternating cycle condition around every inner face of $G$ and are gauge equivalent, see Section A.2 of Appendix A. The results of [41, 42, 44] recalled in Section 2.2 also hold with the complex, bipartite Kasteleyn matrix. In the sequel the graph $G$ will be the double graph $G^D$ or the bipartite graph $G^Q$.

## 3 $Z^u$-Dirac and $Z$-massive Laplacian operators

We let $G$ be an infinite/finite isoradial graph; its dual graph is $G^*/\bar{G}^*$, and in the finite case $G^*$ is its restricted dual. We consider the isoradial embedding of the double graph $G^D = (V^D, E^D)$ given in Section 2.5.2, see also Figure 6. Fix an elliptic parameter $k$. Recall the definition of the torus $T(k) = C/(4KZ+4iK'Z)$, and of the subset $C = C(k) := \mathbb{R}/4KZ$ of $T(k)$.

In Section 3.1 we introduce a family of bipartite Kasteleyn matrices/operators $(K(u))_{u \in C}$ on the double graph $G^D$, referred to as the $Z^u$-Dirac operators. Fixing $u \in C$, a function $F \in C^B$ is said to be $Z^u$-holomorphic if, $K(u)F = 0$. When $k = 0$, the dependence in $u$ disappears and we recover the discrete Dirac operator $\bar{\partial}$ introduced in [45], see Remark 3.4. As a consequence of Theorem 3.6 of Section 3.3, we have that if $F$ is a $Z^u$-holomorphic function, then $F|_V$ is massive harmonic on $G$ and $F|_{V^*}$ is massive harmonic on $G^*$, for the $Z$-massive Laplacian $\Delta^m$ and its dual $\Delta^{m,*}$ of [11]; explaining the part Dirac of the terminology. In the finite case, we moreover introduce the operator $K^\theta(u)$ with specific boundary conditions arising from the forthcoming Theorem 4.2 related to the Ising model, that are different from the natural dimer ones.

In Section 3.4 we restrict to the finite case. Theorem 3.8 proves that, for every $u \in C$, the determinant of the $Z^u$-Dirac operator $K(u)$ is equal, up to an explicit multiplicative constant, to the determinant of the dual massive Laplacian $\Delta^{m,*}$; we show a similar result for the operator $K^\theta(u)$ and the massive Laplacian $\Delta^{m,\theta}(u)$, where $\Delta^{m,\theta}(u)$ has specific boundary conditions and depends on $u$ along the boundary only. Interpreting these determinants as partition functions, this proves that the weighted sum of pairs of dual directed spanning trees is equal, up to an explicit constant, to the weighted sum of rooted spanning forests. In the case $k = 0$, pairs of directed spanning trees become undirected and rooted spanning forests are un-rooted so that this theorem is a consequence of Temperley's bijection [71] and of the matrix-tree theorem [50]. For $k \neq 0$, this result is non-trivial; the proof uses gauge equivalences on bipartite Kasteleyn matrices and on weighted adjacency matrices of directed graphs (digraphs), see Appendix A.
where $u_{\alpha_j} := \frac{u-\alpha_j}{2}$.

Up to an explicit multiplicative constant $G^{m(k)}(x, y)$ is the expected number of visits to y of the associated massive random walk on the infinite graph G started at x, see for example Appendix D4 of [11].

From now on, we consider $k$ such that $(k')^2 \in (0, \infty)$ as fixed and omit all reference to $k$ in the notation.

### 2.5.7 Bipartite dimer models on isoradial graphs

When considering a dimer model on a bipartite isoradial graph $G = (W \cup B, E)$ with weight function $\nu$, instead of a Kasteleyn orientation, one can multiply edge-weights by a complex phase [53, 45]. This defines the complex, bipartite Kasteleyn matrix, denoted $\tilde{K}$ in the context of this section, whose non-zero coefficients are given by:

$$ \forall \text{ edge } w, b \text{ of } G, \quad \tilde{K}_{w,b} = e^{i \frac{\bar{\alpha}_e + \bar{\beta}_e}{2}} \nu_{wb}, $$

where $e^{i\bar{\alpha}_e}, e^{i\bar{\beta}_e}$ are the rhombus vectors of $G^{\diamond}$ associated to the edge $e = (w, b)$. The real and complex bipartite Kasteleyn matrices satisfy the alternating cycle condition around every inner face of $G$ and are gauge equivalent, see Section A.2 of Appendix A. The results of [41, 42, 44] recalled in Section 2.2 also hold with the complex, bipartite Kasteleyn matrix. In the sequel the graph $G$ will be the double graph $G^D$ or the bipartite graph $G^Q$.

## 3 $Z^u$-Dirac and $Z$-massive Laplacian operators

We let $G$ be an infinite/finite isoradial graph; its dual graph is $G^*/\bar{G}^*$, and in the finite case $G^*$ is its restricted dual. We consider the isoradial embedding of the double graph $G^D = (V^D, E^D)$ given in Section 2.5.2, see also Figure 6. Fix an elliptic parameter $k$. Recall the definition of the torus $T(k) = C/(4KZ+4iK'Z)$, and of the subset $C = C(k) := \mathbb{R}/4KZ$ of $T(k)$.

In Section 3.1 we introduce a family of bipartite Kasteleyn matrices/operators $(K(u))_{u \in C}$ on the double graph $G^D$, referred to as the $Z^u$-Dirac operators. Fixing $u \in C$, a function $F \in C^B$ is said to be $Z^u$-holomorphic if, $K(u)F = 0$. When $k = 0$, the dependence in $u$ disappears and we recover the discrete Dirac operator $\bar{\partial}$ introduced in [45], see Remark 3.4. As a consequence of Theorem 3.6 of Section 3.3, we have that if $F$ is a $Z^u$-holomorphic function, then $F|_{\bar{\partial}}$ is massive harmonic on $G$ and $F|_{\bar{\partial}^*}$ is massive harmonic on $G^*$, for the $Z$-massive Laplacian $\Delta^m$ and its dual $\Delta^{m,*}$ of [11]; explaining the part Dirac of the terminology. In the finite case, we moreover introduce the operator $K^{\theta}(u)$ with specific boundary conditions arising from the forthcoming Theorem 4.2 related to the Ising model, that are different from the natural dimer ones.

In Section 3.4 we restrict to the finite case. Theorem 3.8 proves that, for every $u \in C$, the determinant of the $Z^u$-Dirac operator $K(u)$ is equal, up to an explicit multiplicative constant, to the determinant of the dual massive Laplacian $\Delta^{m,*}$; we show a similar result for the operator $K^{\theta}(u)$ and the massive Laplacian $\Delta^{m,\theta}(u)$, where $\Delta^{m,\theta}(u)$ has specific boundary conditions and depends on $u$ along the boundary only. Interpreting these determinants as partition functions, this proves that the weighted sum of pairs of dual directed spanning trees is equal, up to an explicit constant, to the weighted sum of rooted spanning forests. In the case $k = 0$, pairs of directed spanning trees become undirected and rooted spanning forests are un-rooted so that this theorem is a consequence of Temperley's bijection [71] and of the matrix-tree theorem [50]. For $k \neq 0$, this result is non-trivial; the proof uses gauge equivalences on bipartite Kasteleyn matrices and on weighted adjacency matrices of directed graphs (digraphs), see Appendix A.
For every $u \in C$, the operator $K(u)$ is gauge equivalent to an operator $K^g(u)$ associated to a model of directed spanning trees. In Proposition 3.13 of Section 3.5, we prove that this model of directed spanning trees is Z-invariant, thus explaining the part "$Z^u$" in the terminology *Z$^u$-Dirac operator*.

Using Theorem 3.6, in Corollary 3.15 of Section 3.6, we express the inverse $Z^u$-Dirac operator using the Z-massive Green function $G^m$ and its dual $G^{m,*}$ of [11] in the infinite case. This proves in Theorem 3.19 an explicit local expression for a Gibbs measure for the dimer model on $G^D$ with operator $K(u)$, generalizing to the full Z-invariant case the results of [45] proved in the case $k = 0$. In Corollary 3.17, we explicitly express the inverse of the $Z^u$-Dirac operator $K^\theta(u)$ as a function of the finite versions of the Z-massive Green functions. Theorem 3.19 is a planar, directed version of the transfer-impedance theorem of [14], where probabilities of pairs of dual directed spanning trees are computed using the Green functions of massive non-directed random walks. Apart from the locality property which is specific to Z-invariance, a similar result is obtained by Chhita [20] in the case of the square lattice with a specific choice of weights. Sun [70] expresses probabilities of directed spanning trees using the Green function of directed random walks, and Kenyon [46] proves that probabilities of rooted spanning forests are determinantal, without connecting them to directed spanning trees. It might be that the techniques of this paper, in particular gauge transformations on weighted adjacency matrices of digraphs, extend in some respect and allow to relate probabilities of pairs of dual directed spanning trees to massive non-directed Green functions.

## 3.1 Family of $Z^u$-Dirac operators ($K(u)_{u \in C}$

We will be using Section 2.3.1 on the double graph and Temperley's bijection. Recall that vertices of the double graph $G^D$ are partitioned as $V^D = B \cup W$, where $B = V \cup V^*$ and $W \leftrightarrow E$. The diamond graph of $G^D$ is $G^{\diamond/2}$, see Section 2.5.2 and Figure 6. Let $u \in C$ be fixed; we now define the $Z^u$-Dirac operator, first in the infinite case, then in the finite case.

**Infinite case.** Consider the weight function $c(u)$ on edges of $G^D$ defined by, $\forall w \in W, \forall x \in V \cup V^*$ such that $wx$ is an edge of $G^D$,

$$c(u)_{wx} = \begin{cases} f(u_{\alpha_e}, u_{\beta_e}) & \text{if } x \in V \\ f((u_{\alpha_e})^*, (u_{\beta_e})^*) & \text{if } x \in V^*, \end{cases} \quad (3.1)$$

where $e^{i\bar{\alpha}_e}, e^{i\bar{\beta}_e}$ are the rhombus vectors of $G^{\diamond/2}$ associated to the edge $e = (w, x)$; $u_\alpha := \frac{u-\alpha}{2}, u^* := K - u$, and

$$f(u_\alpha, u_\beta) = [\mathrm{sc}(u_\alpha - u_\beta) \mathrm{dn}(u_\alpha) \mathrm{dn}(u_\beta)]^{\frac{1}{2}}.$$

**Remark 3.1.** The function $dn$ is periodic in two directions and naturally defined on the torus $C/(2KZ + 4iK'Z)$ [1]. Since $dn(u + 2iK') = -dn(u)$, and since the function $c(u)$ involves products of two $dn$'s and half arguments, it is defined on the torus $T(k)$. This argument is similar to that used in [11] to define the domain of the massive exponential function.

We restrict to $u \in C$ because the weight function $c(u)$ is then positive; indeed the function $dn$ is, and $\alpha_e, \beta_e$ are such that $\frac{\beta_e - \alpha_e}{2} \in (\varepsilon, \frac{\pi}{2} - \varepsilon)$. Also, on $C$ the (pure imaginary) poles of $c(u)$ are avoided and the weights are thus finite. Results in the sequel which use elliptic trigonometric identities actually hold for all $u \in T(k)$; it is when considering the corresponding dimer model that we use positivity of the weights.

Let $K(u)$ be the complex, bipartite Kasteleyn matrix defined in Section 2.5.7 corresponding to the weight function $c(u)$, with rows indexed by white vertices of $G^D$. Non-zero
For every $u \in C$, the operator $K(u)$ is gauge equivalent to an operator $K^g(u)$ associated to a model of directed spanning trees. In Proposition 3.13 of Section 3.5, we prove that this model of directed spanning trees is Z-invariant, thus explaining the part "$Z^u$" in the terminology *Z$^u$-Dirac operator*.

Using Theorem 3.6, in Corollary 3.15 of Section 3.6, we express the inverse $Z^u$-Dirac operator using the Z-massive Green function $G^m$ and its dual $G^{m,*}$ of [11] in the infinite case. This proves in Theorem 3.19 an explicit local expression for a Gibbs measure for the dimer model on $G^D$ with operator $K(u)$, generalizing to the full Z-invariant case the results of [45] proved in the case $k = 0$. In Corollary 3.17, we explicitly express the inverse of the $Z^u$-Dirac operator $K^\theta(u)$ as a function of the finite versions of the Z-massive Green functions. Theorem 3.19 is a planar, directed version of the transfer-impedance theorem of [14], where probabilities of pairs of dual directed spanning trees are computed using the Green functions of massive non-directed random walks. Apart from the locality property which is specific to Z-invariance, a similar result is obtained by Chhita [20] in the case of the square lattice with a specific choice of weights. Sun [70] expresses probabilities of directed spanning trees using the Green function of directed random walks, and Kenyon [46] proves that probabilities of rooted spanning forests are determinantal, without connecting them to directed spanning trees. It might be that the techniques of this paper, in particular gauge transformations on weighted adjacency matrices of digraphs, extend in some respect and allow to relate probabilities of pairs of dual directed spanning trees to massive non-directed Green functions.

## 3.1 Family of $Z^u$-Dirac operators ($K(u)_{u \in C}$

We will be using Section 2.3.1 on the double graph and Temperley's bijection. Recall that vertices of the double graph $G^D$ are partitioned as $V^D = B \cup W$, where $B = V \cup V^*$ and $W \leftrightarrow E$. The diamond graph of $G^D$ is $G^{\diamond/2}$, see Section 2.5.2 and Figure 6. Let $u \in C$ be fixed; we now define the $Z^u$-Dirac operator, first in the infinite case, then in the finite case.

**Infinite case.** Consider the weight function $c(u)$ on edges of $G^D$ defined by, $\forall w \in W, \forall x \in V \cup V^*$ such that $wx$ is an edge of $G^D$,

$$c(u)_{wx} = \begin{cases} f(u_{\alpha_e}, u_{\beta_e}) & \text{if } x \in V \\ f((u_{\alpha_e})^*, (u_{\beta_e})^*) & \text{if } x \in V^*, \end{cases} \quad (3.1)$$

where $e^{i\bar{\alpha}_e}, e^{i\bar{\beta}_e}$ are the rhombus vectors of $G^{\diamond/2}$ associated to the edge $e = (w, x)$; $u_\alpha := \frac{u-\alpha}{2}, u^* := K - u$, and

$$f(u_\alpha, u_\beta) = [\mathrm{sc}(u_\alpha - u_\beta) \mathrm{dn}(u_\alpha) \mathrm{dn}(u_\beta)]^{\frac{1}{2}}.$$

**Remark 3.1.** The function $dn$ is periodic in two directions and naturally defined on the torus $C/(2KZ + 4iK'Z)$ [1]. Since $dn(u + 2iK') = -dn(u)$, and since the function $c(u)$ involves products of two $dn$'s and half arguments, it is defined on the torus $T(k)$. This argument is similar to that used in [11] to define the domain of the massive exponential function.

We restrict to $u \in C$ because the weight function $c(u)$ is then positive; indeed the function $dn$ is, and $\alpha_e, \beta_e$ are such that $\frac{\beta_e - \alpha_e}{2} \in (\varepsilon, \frac{\pi}{2} - \varepsilon)$. Also, on $C$ the (pure imaginary) poles of $c(u)$ are avoided and the weights are thus finite. Results in the sequel which use elliptic trigonometric identities actually hold for all $u \in T(k)$; it is when considering the corresponding dimer model that we use positivity of the weights.

Let $K(u)$ be the complex, bipartite Kasteleyn matrix defined in Section 2.5.7 corresponding to the weight function $c(u)$, with rows indexed by white vertices of $G^D$. Non-zero
coefficients of $K(u)$ are given by

$$ \forall \text{ edge } wx \text{ of } G^D, \quad K(u)_{w,x} = e^{i \frac{\bar{\alpha}_e + \bar{\beta}_e}{2}} c(u)_{wx}. \qquad (3.2) $$

Recall that this Kasteleyn matrix can also be interpreted as an operator mapping $C^B$ to $C^W$. We refer to this matrix/operator as the $Z^u$-Dirac operator. As an example, we explicitly compute $K(u)$ around a white vertex $w$ of $G^D$.

**Example 3.2.** Figure 9 below sets the notation. A white vertex $w$ of $G^D$ is adjacent to the black vertices $v_1, f_1, v_2, f_2$ of $G^D$ defining a rhombus of the diamond graph $G^\diamond$, such that $v_1, v_2$ belong to $G$ and $f_1, f_2$ to $G^*$. Denote by $2e^{i\bar{\alpha}}$, $2e^{i\bar{\beta}}$ the rhombus vectors of $G^\diamond$ associated to the edge $(v_1, v_2)$, and by $\bar{\theta} = \frac{\bar{\beta}-\bar{\alpha}}{2}$ the rhombus half-angle of this edge.

<figure><img src="image_11.png" /><figcaption>Figure 9: Rhombus $v_1, f_1, v_2, f_2$ of the diamond graph $G^\diamond$ (left) and rhombi of the diamond graph $G^{\diamond/2}$ associated to the edges $wv_1, wf_1, wv_2, wf_2$, with corresponding rhombus vectors and half-angles.</figcaption></figure>

Then, we have:

$$ K(u)_{w,v_1} = -e^{i \frac{\bar{\alpha}+\bar{\beta}}{2}} [\mathrm{sc}(\theta) \mathrm{dn}(u_{\alpha+2K}) \mathrm{dn}(u_{\beta+2K})]^{\frac{1}{2}} = -e^{i \frac{\bar{\alpha}+\bar{\beta}}{2}} [(k')^2 \mathrm{sc}(\theta) \mathrm{nd}(u_{\alpha}) \mathrm{nd}(u_{\beta})]^{\frac{1}{2}} $$

$$ K(u)_{w,v_2} = e^{i \frac{\bar{\alpha}+\bar{\beta}}{2}} [\mathrm{sc}(\theta) \mathrm{dn}(u_{\alpha}) \mathrm{dn}(u_{\beta})]^{\frac{1}{2}} $$

$$ K(u)_{w,f_1} = -ie^{i \frac{\bar{\alpha}+\bar{\beta}}{2}} [(k')^2 \mathrm{sc}(\theta^*) \mathrm{nd}(u_{\beta-2K}) \mathrm{nd}(u_{\alpha})]^{\frac{1}{2}} = -ie^{i \frac{\bar{\alpha}+\bar{\beta}}{2}} [\mathrm{cs}(\theta) \mathrm{dn}(u_{\beta}) \mathrm{nd}(u_{\alpha})]^{\frac{1}{2}} $$

$$ K(u)_{w,f_2} = ie^{i \frac{\bar{\alpha}+\bar{\beta}}{2}} [(k')^2 \mathrm{sc}(\theta^*) \mathrm{nd}(u_{\beta}) \mathrm{nd}(u_{\alpha+2K})]^{\frac{1}{2}} = ie^{i \frac{\bar{\alpha}+\bar{\beta}}{2}} [\mathrm{cs}(\theta) \mathrm{nd}(u_{\beta}) \mathrm{dn}(u_{\alpha})]^{\frac{1}{2}}, $$

$$ \text{using that, } \mathrm{sc}(\theta^*) = \mathrm{sc}(K - \theta) = (k')^{-1} \mathrm{cs}(\theta), \mathrm{dn}(u - K) = k' \mathrm{nd}(u), \mathrm{dn}(u + 2K) = \mathrm{dn}(u). $$

**Definition 3.3.** A function $F \in C^B$ is said to be $Z^u$-holomorphic if

$$ K(u)F = 0. $$

With the notation of Figure 9, this is equivalent to asking that the function $F$ satisfies,

$$ \forall w \in W, \quad K(u)_{w,v_1} \cdot F_{v_1} + K(u)_{w,v_2} \cdot F_{v_2} + K(u)_{w,f_1} \cdot F_{f_1} + K(u)_{w,f_2} \cdot F_{f_2} = 0 $$

$$ \Leftrightarrow \forall w \in W, \quad \mathrm{sc}(\theta)^{\frac{1}{2}} \left( \left[ \mathrm{dn}(u_{\alpha}) \mathrm{dn}(u_{\beta}) \right]^{\frac{1}{2}} \cdot F_{v_2} - \left[ (k')^2 \mathrm{nd}(u_{\alpha}) \mathrm{nd}(u_{\beta}) \right]^{\frac{1}{2}} \cdot F_{v_1} \right) + \\ + i[\mathrm{cs}(\theta)]^{\frac{1}{2}} \left( \left[ \mathrm{nd}(u_{\beta}) \mathrm{dn}(u_{\alpha}) \right]^{\frac{1}{2}} \cdot F_{f_2} - \left[ \mathrm{dn}(u_{\beta}) \mathrm{nd}(u_{\alpha}) \right]^{\frac{1}{2}} \cdot F_{f_1} \right) = 0. $$

**Remark 3.4.** When $k = 0$, then $k' = 1$, $\mathrm{dn} \equiv 1$, $\mathrm{sc} = \mathrm{tan}$, and the $Z^u$-Dirac operator is the discrete Dirac operator $\bar{\partial}$ introduced in [45]. A $Z^u$-holomorphic function is then a holomorphic function as defined in [33, 59, 45, 17]:

$$ \forall w \in W, \quad \tan(\theta)^{\frac{1}{2}} (F_{v_2} - F_{v_1}) + i[\cot(\theta)]^{\frac{1}{2}} (F_{f_2} - F_{f_1}) = 0. $$

Finite case. We will be using Sections 2.3.1 and 2.5.2: the graph $G^{D,r}$ is obtained from $G^D$ by removing the vertex $r$ and its incident edges; the set of black vertices of $G^{D,r}$ is $B^r = V^r \cup V^*$, the set of white vertices is $W^r = W \leftrightarrow E$, and $|B^r| = |W^r|$, see Figure 6
coefficients of $K(u)$ are given by

$$ \forall \text{ edge } wx \text{ of } G^D, \quad K(u)_{w,x} = e^{i \frac{\bar{\alpha}_e + \bar{\beta}_e}{2}} c(u)_{wx}. \qquad (3.2) $$

Recall that this Kasteleyn matrix can also be interpreted as an operator mapping $C^B$ to $C^W$. We refer to this matrix/operator as the $Z^u$-Dirac operator. As an example, we explicitly compute $K(u)$ around a white vertex $w$ of $G^D$.

**Example 3.2.** Figure 9 below sets the notation. A white vertex $w$ of $G^D$ is adjacent to the black vertices $v_1, f_1, v_2, f_2$ of $G^D$ defining a rhombus of the diamond graph $G^\diamond$, such that $v_1, v_2$ belong to $G$ and $f_1, f_2$ to $G^*$. Denote by $2e^{i\bar{\alpha}}$, $2e^{i\bar{\beta}}$ the rhombus vectors of $G^\diamond$ associated to the edge $(v_1, v_2)$, and by $\bar{\theta} = \frac{\bar{\beta}-\bar{\alpha}}{2}$ the rhombus half-angle of this edge.

<figure><img src="image_11.png" /><figcaption>Figure 9: Rhombus $v_1, f_1, v_2, f_2$ of the diamond graph $G^\diamond$ (left) and rhombi of the diamond graph $G^{\diamond/2}$ associated to the edges $wv_1, wf_1, wv_2, wf_2$, with corresponding rhombus vectors and half-angles.</figcaption></figure>

Then, we have:

$$ K(u)_{w,v_1} = -e^{i \frac{\bar{\alpha}+\bar{\beta}}{2}} [\mathrm{sc}(\theta) \mathrm{dn}(u_{\alpha+2K}) \mathrm{dn}(u_{\beta+2K})]^{\frac{1}{2}} = -e^{i \frac{\bar{\alpha}+\bar{\beta}}{2}} [(k')^2 \mathrm{sc}(\theta) \mathrm{nd}(u_{\alpha}) \mathrm{nd}(u_{\beta})]^{\frac{1}{2}} $$

$$ K(u)_{w,v_2} = e^{i \frac{\bar{\alpha}+\bar{\beta}}{2}} [\mathrm{sc}(\theta) \mathrm{dn}(u_{\alpha}) \mathrm{dn}(u_{\beta})]^{\frac{1}{2}} $$

$$ K(u)_{w,f_1} = -ie^{i \frac{\bar{\alpha}+\bar{\beta}}{2}} [(k')^2 \mathrm{sc}(\theta^*) \mathrm{nd}(u_{\beta-2K}) \mathrm{nd}(u_{\alpha})]^{\frac{1}{2}} = -ie^{i \frac{\bar{\alpha}+\bar{\beta}}{2}} [\mathrm{cs}(\theta) \mathrm{dn}(u_{\beta}) \mathrm{nd}(u_{\alpha})]^{\frac{1}{2}} $$

$$ K(u)_{w,f_2} = ie^{i \frac{\bar{\alpha}+\bar{\beta}}{2}} [(k')^2 \mathrm{sc}(\theta^*) \mathrm{nd}(u_{\beta}) \mathrm{nd}(u_{\alpha+2K})]^{\frac{1}{2}} = ie^{i \frac{\bar{\alpha}+\bar{\beta}}{2}} [\mathrm{cs}(\theta) \mathrm{nd}(u_{\beta}) \mathrm{dn}(u_{\alpha})]^{\frac{1}{2}}, $$

$$ \text{using that, } \mathrm{sc}(\theta^*) = \mathrm{sc}(K - \theta) = (k')^{-1} \mathrm{cs}(\theta), \mathrm{dn}(u - K) = k' \mathrm{nd}(u), \mathrm{dn}(u + 2K) = \mathrm{dn}(u). $$

**Definition 3.3.** A function $F \in C^B$ is said to be $Z^u$-holomorphic if

$$ K(u)F = 0. $$

With the notation of Figure 9, this is equivalent to asking that the function $F$ satisfies,

$$ \forall w \in W, \quad K(u)_{w,v_1} \cdot F_{v_1} + K(u)_{w,v_2} \cdot F_{v_2} + K(u)_{w,f_1} \cdot F_{f_1} + K(u)_{w,f_2} \cdot F_{f_2} = 0 $$

$$ \Leftrightarrow \forall w \in W, \quad \mathrm{sc}(\theta)^{\frac{1}{2}} \left( \left[ \mathrm{dn}(u_{\alpha}) \mathrm{dn}(u_{\beta}) \right]^{\frac{1}{2}} \cdot F_{v_2} - \left[ (k')^2 \mathrm{nd}(u_{\alpha}) \mathrm{nd}(u_{\beta}) \right]^{\frac{1}{2}} \cdot F_{v_1} \right) + \\ + i[\mathrm{cs}(\theta)]^{\frac{1}{2}} \left( \left[ \mathrm{nd}(u_{\beta}) \mathrm{dn}(u_{\alpha}) \right]^{\frac{1}{2}} \cdot F_{f_2} - \left[ \mathrm{dn}(u_{\beta}) \mathrm{nd}(u_{\alpha}) \right]^{\frac{1}{2}} \cdot F_{f_1} \right) = 0. $$

**Remark 3.4.** When $k = 0$, then $k' = 1$, $\mathrm{dn} \equiv 1$, $\mathrm{sc} = \mathrm{tan}$, and the $Z^u$-Dirac operator is the discrete Dirac operator $\bar{\partial}$ introduced in [45]. A $Z^u$-holomorphic function is then a holomorphic function as defined in [33, 59, 45, 17]:

$$ \forall w \in W, \quad \tan(\theta)^{\frac{1}{2}} (F_{v_2} - F_{v_1}) + i[\cot(\theta)]^{\frac{1}{2}} (F_{f_2} - F_{f_1}) = 0. $$

Finite case. We will be using Sections 2.3.1 and 2.5.2: the graph $G^{D,r}$ is obtained from $G^D$ by removing the vertex $r$ and its incident edges; the set of black vertices of $G^{D,r}$ is $B^r = V^r \cup V^*$, the set of white vertices is $W^r = W \leftrightarrow E$, and $|B^r| = |W^r|$, see Figure 6
(left) for an example. The set of boundary rhombus pairs of the diamond graph $G^{\diamond}$ is $\mathcal{R}^{\partial}$; $\mathcal{R}^{\partial,r}$ is the set $\mathcal{R}^{\partial}$ without the root pair containing the vertex r. For boundary rhombus pairs of $\mathcal{R}^{\partial}$, we use the notation of Figure 7, which we recall here for convenience of the reader.

<figure><img src="image_12.png" /><figcaption>Figure 10: Notation for vertices of $G^D$ in a boundary rhombus pair of $\mathcal{R}^{\partial}$, and rhombus vectors $e^{i\bar{\alpha}\ell}$, $e^{i\bar{\beta}\ell}$, resp. $e^{i\bar{\alpha}r}$, $e^{i\bar{\beta}r}$, of $G^{\diamond/2}$ assigned to the edge $(v^c, w^l)$, resp. $(v^c, w^r)$. When this pair is not the root one, the marked black edge is where the operator $K^{\partial}(u)$ differs from $K(u)$.</figcaption></figure>

We introduce two versions $K(u)$ and $K^{\partial}(u)$ of the $Z^u$-Dirac operator, corresponding to different boundary conditions; both operators map $C^{B^r}$ to $C^{W^r}$.

The operator $K(u)$ is the finite version of the operator $K(u)$ defined in the infinite case, thus justifying the notation; it is the complex, bipartite Kasteleyn matrix corresponding to the weight function $c(u)$ of (3.1) restricted to $G^{D,r}$, and we do not repeat the definition here.

The operator $K^{\partial}(u)$ has specific boundary conditions arising from the forthcoming Theorem 4.2 related to the Ising model. It is the complex, bipartite Kasteleyn matrix associated to the weight function $c^{\partial}(u)$ differing from $c(u)$ along edges $(w^{\ell}, v^c)$ of the boundary rhombus pairs of $\mathcal{R}^{\partial,r}$. For every edge $wx$ of $G^D$, we have

$$c^{\partial}(u)_{wx} = \begin{cases} c(u)_{wx} & \text{if } (w, x) \notin \{(w^{\ell}, v^c) \in \mathcal{R}^{\partial,r}\} \\ c(u)_{w^{\ell}v^c} \frac{\mathrm{cd}(u_{\beta^r})}{\mathrm{cd}(u_{\alpha^\ell})} & \text{if } (w, x) \in \{(w^{\ell}, v^c) \in \mathcal{R}^{\partial,r}\}. \end{cases} \quad (3.3)$$

In order for $c^{\partial}(u)$ to be finite, we restrict the domain of $u$ to:

$$\begin{aligned} C' &:= C \setminus \{(\alpha^{\ell} + 2K)[4K] : e^{i\bar{\alpha}\ell} \in \mathcal{R}^{\partial}\} \\ &= C \setminus \{(\beta^{r} + 2K)[4K] : e^{i\bar{\beta}r} \in \mathcal{R}^{\partial}\}, \end{aligned} \quad (3.4)$$

where the second equality is a consequence of properties of train-tracks, see Section 2.5.3. As a consequence, the weight function $c^{\partial}(u)$ is everywhere non-zero.

The coefficients of $K^{\partial}(u)$ differing from those of $K(u)$ are those corresponding to edges $(w^{\ell}, v^c)$ of $\mathcal{R}^{\partial,r}$, and we have:

$$K^{\partial}(u)_{w^{\ell}, v^c} = K(u)_{w^{\ell}, v^c} \frac{\mathrm{cd}(u_{\beta^r})}{\mathrm{cd}(u_{\alpha^\ell})} = -e^{i \frac{\bar{\alpha}^{\ell} + \bar{\beta}^{\ell}}{2}} [\mathrm{sc}(\theta^{\partial}) \mathrm{dn}(u_{\alpha^{\ell}+2K}) \mathrm{dn}(u_{\beta^{\ell}+2K})]^{\frac{1}{2}} \frac{\mathrm{cd}(u_{\beta^r})}{\mathrm{cd}(u_{\alpha^{\ell}})}. \quad (3.5)$$

**Remark 3.5.** The weight $c^{\partial}(u)_{w^{\ell}v^c}$ might be negative. Using the physics terminology, see also Remark 5.5, this means that the corresponding dimer model has vortices on boundary faces $v^c, w^l, f^c, w^r$ where this is the case.
(left) for an example. The set of boundary rhombus pairs of the diamond graph $G^{\diamond}$ is $\mathcal{R}^{\partial}$; $\mathcal{R}^{\partial,r}$ is the set $\mathcal{R}^{\partial}$ without the root pair containing the vertex r. For boundary rhombus pairs of $\mathcal{R}^{\partial}$, we use the notation of Figure 7, which we recall here for convenience of the reader.

<figure><img src="image_12.png" /><figcaption>Figure 10: Notation for vertices of $G^D$ in a boundary rhombus pair of $\mathcal{R}^{\partial}$, and rhombus vectors $e^{i\bar{\alpha}\ell}$, $e^{i\bar{\beta}\ell}$, resp. $e^{i\bar{\alpha}r}$, $e^{i\bar{\beta}r}$, of $G^{\diamond/2}$ assigned to the edge $(v^c, w^l)$, resp. $(v^c, w^r)$. When this pair is not the root one, the marked black edge is where the operator $K^{\partial}(u)$ differs from $K(u)$.</figcaption></figure>

We introduce two versions $K(u)$ and $K^{\partial}(u)$ of the $Z^u$-Dirac operator, corresponding to different boundary conditions; both operators map $C^{B^r}$ to $C^{W^r}$.

The operator $K(u)$ is the finite version of the operator $K(u)$ defined in the infinite case, thus justifying the notation; it is the complex, bipartite Kasteleyn matrix corresponding to the weight function $c(u)$ of (3.1) restricted to $G^{D,r}$, and we do not repeat the definition here.

The operator $K^{\partial}(u)$ has specific boundary conditions arising from the forthcoming Theorem 4.2 related to the Ising model. It is the complex, bipartite Kasteleyn matrix associated to the weight function $c^{\partial}(u)$ differing from $c(u)$ along edges $(w^{\ell}, v^c)$ of the boundary rhombus pairs of $\mathcal{R}^{\partial,r}$. For every edge $wx$ of $G^D$, we have

$$c^{\partial}(u)_{wx} = \begin{cases} c(u)_{wx} & \text{if } (w, x) \notin \{(w^{\ell}, v^c) \in \mathcal{R}^{\partial,r}\} \\ c(u)_{w^{\ell}v^c} \frac{\text{cd}(u_{\beta^r})}{\text{cd}(u_{\alpha^\ell})} & \text{if } (w, x) \in \{(w^{\ell}, v^c) \in \mathcal{R}^{\partial,r}\}. \end{cases} \quad (3.3)$$

In order for $c^{\partial}(u)$ to be finite, we restrict the domain of $u$ to:

$$\begin{aligned} C' &:= C \setminus \{(\alpha^{\ell} + 2K)[4K] : e^{i\bar{\alpha}\ell} \in \mathcal{R}^{\partial}\} \\ &= C \setminus \{(\beta^{r} + 2K)[4K] : e^{i\bar{\beta}r} \in \mathcal{R}^{\partial}\}, \end{aligned} \quad (3.4)$$

where the second equality is a consequence of properties of train-tracks, see Section 2.5.3. As a consequence, the weight function $c^{\partial}(u)$ is everywhere non-zero.

The coefficients of $K^{\partial}(u)$ differing from those of $K(u)$ are those corresponding to edges $(w^{\ell}, v^c)$ of $\mathcal{R}^{\partial,r}$, and we have:

$$K^{\partial}(u)_{w^{\ell}, v^c} = K(u)_{w^{\ell}, v^c} \frac{\text{cd}(u_{\beta^r})}{\text{cd}(u_{\alpha^\ell})} = -e^{i \frac{\bar{\alpha}^{\ell} + \bar{\beta}^{\ell}}{2}} [\text{sc}(\theta^{\partial}) \text{dn}(u_{\alpha^{\ell}+2K}) \text{dn}(u_{\beta^{\ell}+2K})]^{\frac{1}{2}} \frac{\text{cd}(u_{\beta^r})}{\text{cd}(u_{\alpha^{\ell}})}. \quad (3.5)$$

**Remark 3.5.** The weight $c^{\partial}(u)_{w^{\ell}v^c}$ might be negative. Using the physics terminology, see also Remark 5.5, this means that the corresponding dimer model has vortices on boundary faces $v^c$, $w^l$, $f^c$, $w^r$ where this is the case.
## 3.2 Z-massive Laplacian operators

In the infinite case, we consider the Z-massive Laplacian associated to the isoradial graph G introduced in [11], whose definition is recalled in Section 2.5.6. We also consider the dual Z-massive Laplacian $\Delta^{m,*}$ of the dual isoradial graph G*.

The purpose of this section is to define the finite versions of the Z-massive Laplacian and dual Z-massive Laplacian that are used in this paper; we introduce two operators. The first is the finite version of the dual Z-massive Laplacian, denoted $\Delta^{m,*}$ as in the infinite case. It acts on $C^{V*}$, where recall that $V^*$ is the vertex set of the restricted dual G*, whose boundary vertices are defined to be those adjacent to the vertex o in $\bar{G}^*$. For every vertex f of G* of degree d in $\bar{G}^*$, denote its neighbors by $f_1, \dots, f_d$. Let $\underline{d}$ be the degree of f in G*, then if f is a boundary vertex of G* we have $\underline{d} \neq d$, and we label the vertices so that the first $\underline{d}$ ones are common to G* and $\bar{G}^*$. For every $j \in \{1, \dots, \underline{d}\}$, $(f, f_j)$ is an edge of G* and we let $\bar{\theta}_j^*$ denote its half-angle in $G^\circ$; for every $j \in \{\underline{d}+1, \dots, d\}$, $(f, f_j) = (f, o)$ corresponds to an edge $fw_j$ of the double graph $G^D$, where $w_j$ is the unique vertex of $G^D$ which belongs to the edge $fo$, we let $\bar{\theta}_j^*$ denote the half-angle of the edge $(f, w_j)$ in $G^{\circ/2}$. Then, non-zero coefficients of the finite version of the dual Z-massive Laplacian $\Delta^{m,*}$ are given by,

$$ \Delta_{f,f'}^{m,*} = \begin{cases} -\mathrm{sc}(\theta^*) & \text{if } (f, f') \text{ is an edge of } G^* \\ \sum_{j=1}^{d} A(\theta_j^*) & \text{if } f' = f \text{ is a vertex of } G^*. \end{cases} $$

The corresponding conductances $\rho^*$ and masses $m^*$ are as defined in Section 2.4.

The second operator acts on $C^V$; it is denoted $\Delta^{m,\partial}(u)$ and we restrict to $u \in C'$. It is not the natural finite version of the operator $\Delta^m$ but has boundary conditions inherited from those of the $Z^u$-Dirac operator $K^\partial(u)$, see the proof of Theorem 3.6 below. We use the notation of the infinite Z-massive Laplacian operator $\Delta^m$, see Section 2.5.6, and also the notation of Figure 10 for the boundary rhombus pairs of $\mathbb{R}^\partial$.

In the following, the two boundary vertices $v^r, v^\ell$ of $\mathbb{R}^\partial$ enter the same framework, and we denote such a vertex by $v^\partial$. For every $v^\partial$ of $\mathbb{R}^\partial$ of degree d in G, let $v_1, \dots, v_d$ denote its neighbors in cclw order with $v_1$ on the boundary of G on the left of $v^\partial$, see Figure 12 (right). Note that if $v^\partial = v^{\partial,r}$ belongs to the root pair, then $v^c = r$ is considered as a neighbor of $v^{\partial,r}$. The coefficients of $\Delta^{m,\partial}(u)$ differing from those of $\Delta^m$ are those corresponding to edges $(v^c, v^\ell)$ of $\mathbb{R}^{\partial,r}$, and to vertices $v^r, v^\ell, v^c$ of $\mathbb{R}^\partial$ such that $v^c \neq r$:

$$ \Delta^{m,\partial}(u)_{v,v'} = \begin{cases} -\mathrm{sc}(\theta^\partial) \frac{\mathrm{cd}(u_{\beta^r})}{\mathrm{cd}(u_{\alpha^\ell})} & \text{if } (v,v') \in \{(v^c, v^\ell) \in \mathbb{R}^{\partial,r}\} \\ k' \mathrm{sc}(\theta^\partial) \frac{\mathrm{nd}(u_{\beta^\ell}) \mathrm{nd}(u_{\beta^r})}{\mathrm{cn}(u_{\alpha^\ell}) \mathrm{cn}(u_{\alpha^\ell})} (\mathrm{cn}(u_{\beta^r}) + \mathrm{cn}(u_{\alpha^\ell})) & \text{if } v' = v = v^c \in \{v^c \in \mathbb{R}^{\partial,r}\} \\ k' \sum_{j=1}^d \mathrm{sc}(\theta_j) \mathrm{nd}(u_{\alpha_j}) \mathrm{nd}(u_{\alpha_{j+1}}) & \text{if } v' = v = v^\partial \in \{v^r, v^\ell \in \mathbb{R}^\partial\}. \end{cases} \quad (3.6) $$

The corresponding conductances and masses are denoted $\rho^\partial(u)$ and $m^\partial(u)$. Note that the diagonal term $\Delta^{m,\partial}(u)_{v^\partial, v^\partial}$ has the same form as that of the natural finite version of the Z-massive Laplacian.

In order to state Theorem 3.6 relating the $Z^u$-Dirac and the Z-massive Laplacian operators, in the finite case we need to introduce the matrix $Q(u)$. It has rows indexed by black vertices of $G^r$ and columns by those of the restricted dual G*. The only non-zero coefficients of $Q(u)$ are for rows corresponding to vertices $v^c$ of boundary rhombus pairs of $\mathbb{R}^{\partial,r}$. For such a vertex $v^c$, there is one non-zero coefficient, corresponding to the column $f^c$, given by

$$ Q(u)_{v^c, f^c} = -i \frac{\mathrm{nd}(u_{\beta^\ell})}{\mathrm{cd}(u_{\alpha^\ell})} (\mathrm{cd}(u_{\beta^r}) - \mathrm{cd}(u_{\alpha^\ell})). $$