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

Consider an elliptic parameter $k$; we introduce a family of $Z^u$-Dirac operators $(\mathsf{K}(u))_{u\in\Re(\mathbb{T}(k))}$, relate them to the $Z$-massive Laplacian of [BdTR17b], and extend to the full $Z$-invariant case the results of Kenyon [Ken02] on discrete holomorphic and harmonic functions, which correspond to the case $k=0$. We prove, in a direct statistical mechanics way, 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 ${\mathsf{G}}^{\scriptscriptstyle{\mathrm{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 ${\mathsf{G}}^{\scriptscriptstyle{\mathrm{Q}}}$ arising from the XOR-Ising model, and we prove matrix identities between the Kasteleyn matrix of ${\mathsf{G}}^{\scriptscriptstyle{\mathrm{Q}}}$ and the $Z^u$-Dirac operator, that allow to reach inverse matrices as well as determinants.


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 * .
• Messikh [Mes06] 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 [BDC12] using the FK-Ising observable of [Smi10] away from the critical point.
• In Smirnov and Chelkak-Smirnov's proof of conformal invariance of the critical Z-invariant Ising model [Smi06,Smi10,CS12], the key discrete tools are observables -spin or FK (see also [KC71]) -that are holomorphic. Discrete holomorphic functions in turn are naturally related to harmonic functions [Duf68,Mer01,Ken02,CS11]; the paper [MS10] extends some of the above to the massive case.
• 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 [dT13, dT16]. An abstract proof of this identity is given in the toroidal Z-invariant case in [BdT11,BdTR17a].
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 [Ken02]; harmonic functions are in the kernel of the Laplacian operator. The spanning tree partition function is equal to the determinant of the Laplacian operator [Kir47]. 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 [Ken02], 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 [BdTR17b] and the study of the corresponding dimer model on the double graph G D .
To study the Ising model, we use Fisher's correspondence [Fis66] relating the low (or high) temperature expansion of the model [KW41a,KW41b] 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 [TF61, Ken97, CKP01, KOS06,BdT10]. 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 [CS12], see [CCK15], and the fermionic spinor observable of [KC71], see [Dub11]. 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 [BdTR17b]. 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 38 and also Theorem 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 [Wil11] constructed from two independent Z-invariant Ising models [Dub11,BdT14]. 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 [Dub11], 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; 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 [Ons44,Ken99], fully developed by Baxter [Bax78,Bax86,Bax89], see also [PAY06,AP87,AP02]. A Zinvariant 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 ∈ (0, ∞), and the coupling constants J are [Bax89]: ∀ e ∈ E, J e = 1 2 ln 1 + sn(θ e |k) cn(θ e |k) , where sn, cn are two of the Jacobi elliptic trigonometric functions, andθ e = θ e π 2K is an angle associated to the edge e in the isoradial embedding. When k = 0, i.e. k = 1, the elliptic functions sn, cn, are the trigonometric functions sin, cos, and the Ising model is critical [Li12,CD13,Lis14b]. As (k ) 2 varies from 0 to ∞, the coupling constants range from ∞ to 0 [BdTR17a] thus covering the whole range of inverse temperatures. In the paper [BdTR17b], we introduce Z-invariant rooted spanning forests with associated operator the massive Laplacian ∆ m ; when k = 0, we recover the critical Laplacian of [Ken02]. We also prove an explicit local expression for its inverse, the Z-massive Green function G m , using the discrete massive exponential function [BdTR17b]. 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. We introduce a family of Z u -Dirac operators (K(u)) u∈ (T(k)) on the double graph G D = (W ∪ B, E D ) associated to pairs of dual directed spanning trees, extending to the full Zinvariant case the Dirac operator∂ of [Ken02], corresponding to k = 0. In the finite case, we introduce a family of operators (K ∂ (u)) u∈ (T(k)) , with boundary conditions tuned for the Ising model. 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 [Ken02,Tem74,Ken99]. We prove, see also Theorem 18: Theorem 1.
• Infinite case. Let u ∈ (T(k)), then the Z u -Dirac operator K(u), the Z-massive Laplacian ∆ m and the dual ∆ m, * of [BdTR17b] satisfy the following identity: • Finite case. Let u ∈ (T(k)) , then the Z u -Dirac operators K(u), K ∂ (u), the Z-massive Laplacian ∆ m,∂ (u) and the dual ∆ m, * satisfy the following identity: A function F ∈ C B is said to be Z u -holomorphic if K(u)F = 0. As a consequence of Theorem 1, if F is Z u -holomorphic, then F |V is Z-massive harmonic on G and F |V * is Zmassive harmonic on G * , thus explaining the part "Dirac" in "Z u -Dirac operator".
In the infinite case, Theorem 1 yields the following relations for inverse operators, see also Corollary 27; the statement in the finite case is given in Corollary 29.
Corollary 2 (Infinite case). For every u ∈ (T(k)), consider the operator K(u) −1 mapping C W to C B whose coefficients are defined by, for everyv,f, w as in Figure 14, 1 2 G m, * f,f 2 + function of a massive, non-directed random walk. Apart from the locality property which is specific, a similar result is obtained by Chhita [Chh12] in the case of Z 2 with a specific choice of weights.
In Theorem 20 and Corollary 22, we restrict to the finite case and prove relations on determinants; we show, det ∆ m,∂ (u).
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 [Tem74], but the correspondence does not extend when k = 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 25, we prove that the latter is Z-invariant, thus explaining the part "Z u " of the terminology "Z u -Dirac operator".
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 2 yields an expression using the Z-massive Green function of [BdTR17b]. We obtain, see also Corollary 48 and Corollary 50 for the finite case, Corollary 6 (Infinite case). For everyw, b of G Q as in Figure 22, This proves, in an alternative way, an explicit local expression for a Gibbs measure of the dimer model on the graph G Q [BdTR17a], where the locality property is seen as directly inherited from that of the Z-massive Green function.
Using Theorem 4 and additional combinatorial arguments, we prove in Theorem 36 that the determinants of K Q and of the Z u -Dirac operator are equal, up to an explicit constant. By [Dub11], the determinant of K Q is equal up to a constant, to the squared partition function of the Ising model. Combining this with Theorem 3 gives, see also Corollary 38 for the explicit value of C(u), Theorem 7. For every u ∈ (T(k)) , [Z + Ising (G, J)] 2 = C(u)| det ∆ m,∂ (u)|.
When k = 0, we essentially recover the result of [dT16]; see Section 6.1.
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 matrixK Q . Building on the work of Dubédat [Dub11] and proving additional matrix relations, we express coefficients of the inverse operator (K F ) −1 using coefficients of the inverse operator (K Q ) −1 . Partitioning vertices of G F as A ∪ B as in [Dub11], we obtain, see also Theorem 55, Theorem 8 (Finite and infinite cases). Using the notation of Figure 26, there are four cases to consider: 1. For everyā ∈ A and every b ∈ B such that, when the graph G F is moreover finite, b is not a boundary vertex: 2. When the graph G F is finite, for everyā ∈ A and every boundary vertex b of B, we have 3. For everyā, a ∈ A, where κā ,a = 0 ifā and a do not belong to the same decoration, and to ± 1 4 if they do. 4. For everyb, b ∈ B, where (K F ) −1 a 1 ,b , (K F ) −1 a 2 ,b 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 −2Je on edges, and is related to a bipartite dimer model with some negative weights, see also Remark 56. As mentioned in [Dub11], bozonisation identities somehow prove the existence of such linear relations, but working them out requires more work, which is the subject of the above theorem.
Next we restrict to the Z-invariant case. The coefficient (K F ) −1 a,a is equal, up to an additive constant, to the coefficient (K Q ) −1 w,b , and is thus expressed using the inverse Z u -Dirac operator using Corollary 6 in the infinite case, and Corollary 50 in the finite case. The same holds for (K F ) −1 a,b when b is a boundary vertex. The coefficient (K F ) −1 b,b is a simple linear combination of two coefficients (K F ) −1 a 1 ,b , (K F ) −1 a 2 ,b , so we are left with expressing the coefficient (K F ) −1 a,b . Choosing a specific value of u in Corollary 5, and using Corollary 2 gives, see also Corollary 59 for the finite case, Corollary 9 (Infinite case). Let u f = α f +β f 2 + K. Then, sn(θ f ) nd K−θ f 2 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 Ising , defined by: where Z Ising (G, J) = σ∈{−1,1} V exp e=vv ∈E J e σ v σ v is the normalizing constant known as the Ising partition function.
From now on, we suppose that the planar graph G is embedded and simply connected. 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 + Ising and Z + Ising (G, J) the corresponding Boltzmann measure and partition function 1 .
Denote byḠ * = (V * ,Ē * ) the dual graph of G, and by o the vertex ofḠ * corresponding to the unbounded face of G. Consider also the restricted dual graph G * = (V * , E * ) obtained fromḠ * 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 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 [KW41a,KW41b]: 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 [CS12] 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 35. 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Ḡ * 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 * .  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.

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 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, simply connected 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 ν is assigned to edges of G.

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 dimer , is given by e∈M ν e is the normalizing constant known as the dimer partition function.
The main tool used to study the dimer model is the Kasteleyn matrix [Kas61, Kas67,TF61], 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 [Kas67], 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 ν 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, Note that the matrix K is skew symmetric.
When the graph G is bipartite, the set of vertices can be split into V = W ∪ 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 byK. It has rows indexed by white vertices of G and column by black ones. Non-zero coefficients are defined by: Note that the bipartite Kasteleyn matrix can also be defined as minus the transpose of the above matrixK; 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 [Kas61,Kas67,TF61] and for the dimer Boltzmann measure [Ken97]. Here are their statements.
Theorem 10 ([Kas61, Kas67,TF61]). The dimer partition function of the graph G with weight function ν is equal to: When the graph G is moreover bipartite, we have: Theorem 11 ( [Ken97]). The probability of occurrence of a subset E = {e 1 = x 1 y 1 , . . . , e l = x l y l } of edges of G, chosen with respect to the dimer Boltzmann measure P dimer is equal to where (K −1 ) E is the sub-matrix of the inverse Kasteleyn matrix K −1 whose rows and columns are indexed by vertices x 1 , y 1 , . . . , x l , y l .
When the graph G is moreover bipartite, the subset of edges E is written as E = {e 1 = w 1 b 1 , . . . , e l = w l b l }, and we also have, where (K −1 ) E is the sub-matrix of the inverse bipartite Kasteleyn matrixK −1 whose rows are indexed by black vertices b 1 , . . . , b l and columns by white vertices w 1 , . . . , w l .

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 : When G is bipartite, the bipartite Kasteleyn matrixK is an operator mapping C W to C B : Explicit expressions for Gibbs measures involve inverse Kasteleyn matrices/operators. An inverse Kasteleyn operator L is asked to satisfy the following conditions: • L x,y → 0 as the distance between x and y tends to infinity. Uniqueness is established when the graph G is Z 2 -periodic [She05,BdT10]. 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 [Coo14].
Consider the σ-field generated by cylinder sets of M(G). In all of the above cases, there is an explicit expression for a Gibbs measure P dimer on (M(G), F) whose probabilities on cylinder sets is given by the formulas of Theorem 11 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≥1 , where G n = G/nZ 2 . We refer to the original papers for an exact statement, see also Theorem 31 which has the same form.

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 [Tem74,BP93,KPW00]. Furthermore, we consider two dimer representations of the Ising model. The first is related to the LTE of the Ising model [KW41a,KW41b], while the second arises from the XOR-Ising model, built from two independent copies of the Ising model [Dub11,BdT14]. 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Ḡ * 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 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.

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Ḡ * . 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Ḡ * 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, Suppose again that G is finite, fix a vertex r of G amongst the additional vertices on boundary edges, and let V r = V \ {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 ∪ V * and its set of white vertices is W r = W ↔ E, see Figure 2 (right) for an example. Note that G D,r has the same number of black and white vertices: 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Ḡ * , known as the dual spanning tree.
Consider the fixed boundary vertex r of G as above. Denote by T r (G) the set of r-dST of G, by T o (Ḡ * ) the set of o-dST ofḠ * , and by T r,o (G,Ḡ * ) the set of pairs of dual directed spanning trees (dST-pairs) of G andḠ * 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 [Tem74], extended by [BP93] to general non-directed graphs and by [KPW00] to the directed case, proves a weight preserving bijection between dimer configurations of the double graph G D,r and dST-pairs of T r,o (G,Ḡ * ). It relies on the following bijection between edges of G D,r and directed edges of G andḠ * . Let w ∈ W r = W , x ∈ V r ∪ V * such that wx is an edge of G D,r , then (2) Note that there are no directed edges of G exiting the vertex r, and no directed edges of 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Ḡ * ; Temperley's bijection states that subsets defining dimer configurations are in correspondence with subsets defining dST-pairs of T r,o (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 andc a weight function on directed edges of G,Ḡ * . The relation between c andc which makes Temperley's bijection weight preserving naturally arises from the bijection between edges of G D,r and directed edges of G,Ḡ * . Let w ∈ W r , x ∈ V r ∪ V * , such that wx is an edge of G D,r . Using the notation of (2), we have andc r,v = 0 for every vertex v ∈ V r adjacent to r,c o,f = 0 for every vertex f ∈ V * adjacent to o.
Model on pairs of dual directed spanning trees. Suppose that directed edges of G,Ḡ * are assigned the weight functionc. Consider the Boltzmann measure on dST-pairs, denoted P r,o dST-pairs , defined by . There is also a natural correspondence between the dST-pairs Boltzmann measure P r,o dST-pairs and the dimer Boltzmann measure P D dimer on G D,r with weight function c. Note that ifc ≡ 1 on edges ofḠ * , resp. on edges of G, then Z r,o dST-pairs ((G,Ḡ * ),c) is equal to the partition function Z r dST (G,c) of r-directed spanning trees of G, resp. Z o dST (Ḡ * ,c) of o-directed spanning trees ofḠ * .

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 [Fis66,Dub11], 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Ḡ * and do the same procedure as in the infinite case. Then, remove the decoration of the vertex o as well as all edges ofḠ * 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 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.
Mapping between LTE polygon configurations and dimers. Suppose that G is finite. Fisher [Fis66] 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 [Fis66], 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, the dimer weight function µ J on edges of G F is defined to be, see Equation (1): if e is an internal edge e −2Je if e is an external edge arising from a dual edge e * of G * .
Let P F dimer and Z dimer (G F , µ J ) be the corresponding dimer Boltzmann measure and partition function. Then, as a consequence of Fisher's correspondence we have,

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 [dT07]. 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∪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.
We consider the dimer model on the bipartite graph G Q arising from the XOR-Ising model [Wil11], also known as the polarization of the Ising model [KB79], 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: approach, and [BdT14] based on results of [Nie84,WL75] for the second one. The dimer weight function ν J on G Q is defined by, for every edge e of G Q , if e is a boundary quadrangle edge in the finite case tanh(2J e ) if e is a quadrangle edge/non-boundary quadrangle edge in the infinite/finite case, "parallel" to an edge e of G cosh −1 (2J e ) if e is a quadrangle edge, "parallel" to a dual edge e * of an edge e of G.
When the graph G Q is finite, we let P Q dimer and Z dimer (G Q , ν J ) be the corresponding dimer Boltzmann measure and partition function. As a consequence of [Dub11], see also Corollary 54, we have: Combining this with Equation (4) for the Ising partition function, and denoting by E ∂ the set of boundary edges of the graph G, we obtain where in the last line we used that |E| = |E * | + |E ∂ | and Euler's formula: |E| = |V| + |V * | − 1.

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 r 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∈V , and that directed edges have positive conductances, denoted ρ, meaning that every directed edge (x, x ) has conductance ρ 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 T rooted at a vertex of G, denoted x T . Let F(G) denoted the set of rdSF of the graph G.
Suppose that G is finite and consider the Boltzmann measure on rdSF, denoted P rdSF , defined by: x is the rdSF partition function. Whenever conductances are symmetric, i.e., ρ x,x = ρ x ,x , we will remove the "d" in rdSF.
As a consequence of the directed version of Kirchhoff's matrix-tree theorem [Kir47,Tut48], the rdSF partition function is computed using the massive Laplacian operator/matrix as follows. The massive Laplacian operator ∆ m : C V → C V is defined by: The operator ∆ m is represented by a matrix, also denoted ∆ m , whose non-zero coefficients are given by: Consider the graph G † constructed from G by adding a cemetery vertex † and an edge (x, †) for every vertex x such that m x = 0. Define the modified weight function ρ m on (directed) edges of G † by, There is a natural weight-preserving bijection between T † (G † ) and F(G): a †-directed spanning tree of G † corresponds to the rooted directed spanning forest of G obtained by replacing every edge (x, †) of the dST by a root of the rdSF.
Denote by ∆ † the (non-massive) Laplacian matrix of G † with weight function ρ m on the edges. Then, ∆ m is the Laplacian matrix ∆ † from which one has removed the row and column corresponding to the cemetery † and thus, by Kirchhoff's matrix-tree theorem [Kir47,Tut48], the determinant of ∆ m counts ρ m weighted †-dST of G † . Using the bijection between †-dST of G † and rdSF of G, we thus have, When there is at least one vertex with positive mass, the massive Green function, denoted G m , is the inverse of the massive Laplacian ∆ m . In the chore of the paper, graphs are written with the letter G with or without superscripts, so that 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 ρ and masses m, see for example Appendix D of [BdTR17b], where a number of facts are recalled.

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

Isoradial graphs, diamond graphs and angles
Isoradial graphs naturally appear when considering a discrete version of the Cauchy-Riemann equations, see [Duf68] and [Mer01,Ken02,CS11]; they also arise in Z-invariant models when solving the corresponding Yang-Baxter equations [Bax89,CS06]; the name isoradial comes from the paper [Ken02].
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 [CS12]: 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).
In the infinite case, the diamond graph, denoted G , is constructed from an isoradial graph G and its dual G * as follows: its vertex set is V ∪ 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 are side-length-2 rhombi.
There is a bijection between rhombi of G 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θ e ∈ (0, π 2 ) defined to be the half-angle of the corresponding rhombus at the edge e. We furthermore ask thatθ e ∈ (ε, π 2 −ε), for some ε > 0. The rhombus angle of the dual edge e * isθ e * = π 2 −θ e :=θ * e . To a directed edge e = (v, v ) of G we further assign two rhombus vectors 2e iᾱe , 2e iβe of G , such that 2e iᾱe is on the right of the edge (v, v ), see Figure 5 (right). The anglesᾱ e andβ e are defined so thatβ e−ᾱe 2 =θ e . Whenever no confusion occurs, we remove the subscript e from the notation. In the finite case, the diamond graph, also denoted G , 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 come in pairs, both having the same rhombus half-angle; let us denote by R ∂ the set of boundary rhombus pairs, an instance is highlighted in Figure 5 (right, light grey).

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 , see Figure 6 (left). r r 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 /2 (grey lines), a boundary rhombus pair of R ∂ and the root pair (light grey).
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 , 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 R ∂ , the following notation will be used throughout the paper and is illustrated in Figure 7 below. Let v , 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 v r of the original graph. Denote by w , w r the white vertices of the double graph G D and by w , b , 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 35 is that of Figure 7 (center), see also Figure 6 (right), where the boundary quadrangle edge b w c has length-0 and is "replaced" by the external edge w b . 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 R ∂ 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Ḡ * ; 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 ∈ R ∂ } for the set of vertices of type v belonging to boundary rhombus pairs, and similarly for other vertices or edges of 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 by cutting rhombi into four equal length-1 rhombi. Denote this graph by G /2 in the infinite case and, in the finite case, let G /2 be this graph where the boundary quarter rhombi crossed by no edge of G D are removed. Then G /2 is the diamond graph of G D . Note that G /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 /2 is given in Figure 6 (left and right, grey).
Amongst the boundary rhombus pairs of R ∂ 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 R ∂,r the set R ∂ 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 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.

Train-tracks
A train-track of a finite isoradial graph G, also known as a de Bruijn line [dB81a, dB81b] or a rapidity line [Bax86] is a maximal chain of edge-adjacent rhombi of the diamond graph G , such that when entering a rhombus one exits along the opposite edge [KS05]; each train-track τ has a parallel direction ±2e iᾱτ . Consider the simply connected domain D(G) obtained by taking the union of the faces of G. Then a train-track τ enters and exits D(G), and there are exactly two parallel edges of τ outside of D(G), see Figure 5, (right). The boundary rhombus vectors {2e iᾱ , 2e iβ r ∈ R ∂ }, {2e iᾱ r , 2e iβ ∈ R ∂ } come in parallel pairs, and all the parallel directions of the train-tracks are encoded in {±2e iᾱ , ±2e iβ ∈ R ∂ } = {±2e iᾱ r , ±2e iβ r ∈ R ∂ }.

Elliptic angles
Consider an elliptic modulus k, then k = (1 − k 2 ) 1 2 is the complementary elliptic modulus. Suppose that k is such that (k ) 2 ∈ (0, ∞). The complete elliptic integral of the first kind, and for later purposes we also need K = K(k ). As in [BdTR17b,BdTR17a], we need the following linear transformation of rhombus angles and vectors of the diamond graph G associated to edges:

Z-invariant Ising model and corresponding dimer models
Underlying Z-invariance is the star-triangle transformation, also known as the Y-∆ 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 ; the embedding of the additional vertex of the triangle is given by the cubic flip, see Figure 8. Bax86] 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-∆ 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 [PAY06]. Extending the form of the solutions to the whole of the graph naturally leads to introducing isoradial graphs: the solution is parametrized by the rhombus half-angles assigned to edges and by the elliptic modulus k, where k is such that (k ) 2 = 1 − k 2 ∈ (0, ∞), see [Bax89]. The Z-invariant coupling constants are explicitly given by: where sn, cn are two of the twelve Jacobi elliptic trigonometric functions. We refer the reader to [AS64,Law89] 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 µ J on the Fisher graph G F arising from Fisher's correspondence is given by, for every edge e of G F , if e is an external edge arising from an edge e * of G * .
The dimer weight function ν J on the bipartite graph G Q arising from the XOR-Ising model is given by, for every edge e of G Q , if e is a boundary quadrangle edge in the finite case sn(θ e |k) if e is a quadrangle edge/non-boundary quadrangle edge in the infinite/finite case, parallel to an edge e of G cn(θ e |k) if e is parallel to a dual edge e * of an edge e of G. (8)

The Z-invariant massive Laplacian
In the paper [BdTR17b], 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ᾱ , 2e iβ and a half-angleθ of the diamond graph G . Denote by v 1 , . . . , v d the neighbors of a vertex v of degree d, and for every edge (v, v j ) use the notation 2e iᾱ j , 2e iᾱ j+1 ,θ j for the associated rhombus vectors and half-angle. Fix an elliptic modulus k such that (k ) 2 ∈ (0, ∞). Then, for every edge (v, v ) and every vertex v of G, the conductances ρ k and masses m k of [BdTR17b] are defined by: where sc = sn cn , and 1 2 dτ, is the complete elliptic integral of the second kind. Note that m k ≥ 0, by Proposition 6 of [BdTR17b].
The corresponding Z-massive Laplacian matrix ∆ m(k) has non-zero coefficients given by: The dual graph G * is also isoradial, we denote by ρ k, * , m k, * the associated conductances and masses, and byθ * the half-angle of a dual edge. We let ∆ 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 ∆ m(k) is the Z-massive Green function; it is denoted G m(k) . In [BdTR17b], we prove the following explicit local expression for coefficients of G m(k) : for every pair of vertices x, y of G, where Γ x,y is a vertical contour on the torus T(k) := C/(4KZ + 4iK Z), and e( · |k) : V × V × C → C is the massive exponential function defined in [BdTR17b]. To compute e (x,y) (u|k), consider a path x = x 1 , . . . , x n = y of the diamond graph G from x to y, let 2e iᾱ j be the rhombus vector corresponding to the edge x j x j+1 , then where u α j := u−α 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 [BdTR17b].
From now on, we consider k such that (k ) 2 ∈ (0, ∞) as fixed and omit all reference to k in the notation.

Bipartite dimer models on isoradial graphs
When considering a dimer model on a bipartite isoradial graph G = (W ∪ B, E) with weight function ν, instead of a Kasteleyn orientation, one can multiply edge-weights by a complex phase [Kup98,Ken02]. This defines the complex, bipartite Kasteleyn matrix, denotedK in the context of this section, whose non-zero coefficients are given by: where e iᾱe , e iβe are the rhombus vectors of G 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 [Kas61,Kas67,Ken97] 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 .

Z u -Dirac and Z-massive Laplacian operators
We let G be an infinite/finite isoradial graph; its dual graph is 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 k and recall the definition of the torus T(k) = C/(4KZ + 4iK Z).
In Section 3.1 we introduce a family of bipartite Kasteleyn matrices/operators (K(u)) u∈ (T(k)) on the double graph G D , referred to as the Z u -Dirac operators. Fixing u ∈ (T(k)), a function F ∈ 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∂ introduced in [Ken02], see Remark 16. As a consequence of Theorem 18 of Section 3.3, we have that if F is a Z uholomorphic function, then F |V is massive harmonic on G and F |V * is massive harmonic on G * , for the Z-massive Laplacian ∆ m and its dual ∆ m, * of [BdTR17b]; explaining the part Dirac of the terminology. In the finite case, we moreover introduce the operator K ∂ (u) with specific boundary conditions arising from the forthcoming Theorem 35 related to the Ising model, that are different from the natural dimer ones.
In Section 3.4 we restrict to the finite case. Theorem 20 proves that, for every u ∈ (T(k)), 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 ∆ m, * ; we show a similar result for the operator K ∂ (u) and the massive Laplacian ∆ m,∂ (u), where ∆ m,∂ (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 [Tem74] and of the matrix-tree theorem [Kir47]. For k = 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 ∈ (T(k)), the operator K(u) is gauge equivalent to an operator K g (u) associated to a model of directed spanning trees. In Proposition 25 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 18, in Corollary 27 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 [BdTR17b] in the infinite case. This proves in Theorem 31 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 [Ken02] proved in the case k = 0. In Corollary 29, we explicitly express the inverse of the Z u -Dirac operator K ∂ (u) as a function of the finite versions of the Z-massive Green functions. Theorem 31 is a planar, directed version of the transfer-impedance theorem of [BP93], 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 [Chh12] in the case of the square lattice with a specific choice of weights. Sun [Sun16] expresses probabilities of directed spanning trees using the Green function of directed random walks, and Kenyon [Ken17] 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 nondirected Green functions.

Family of Z u -Dirac operators (K(u)) u∈ (T(k))
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 The diamond graph of G D is G /2 , see Section 2.5.2 and Figure 6. Let u ∈ (T(k)) 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 where e iᾱe , e iβe are the rhombus vectors of G /2 associated to the edge e = (w, x); u α := u−α 2 , 2 . Remark 13. The function dn is periodic in two directions and naturally defined on the torus C/(2KZ + 4iK Z) [AS64]. 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 [BdTR17b] to define the domain of the massive exponential function.
We restrict to u ∈ (T(k)) because the weight function c(u) is then positive; indeed the function dn is, and α e , β e are such thatβ e−ᾱe 2 ∈ (ε, π 2 − ε). Also, on (T(k)) 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 ∈ 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 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 14. 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 , such that v 1 , v 2 belong to G and f 1 , f 2 to G * . Denote by 2e iᾱ , 2e iβ the rhombus vectors of G associated to the edge (v 1 , v 2 ), and byθ =β −ᾱ 2 the rhombus half-angle of this edge. Figure 9: Rhombus v 1 , f 1 , v 2 , f 2 of the diamond graph G (left) and rhombi of the diamond graph G /2 associated to the edges wv 1 , wf 1 , wv 2 , wf 2 , with corresponding rhombus vectors and half-angles.
With the notation of Figure 9, this is equivalent to asking that the function F satisfies, Remark 16. When k = 0, then k = 1, dn ≡ 1, sc = tan, and the Z u -Dirac operator is the discrete Dirac operator∂ introduced in [Ken02]. A Z u -holomorphic function is then a holomorphic function as defined in [Duf68,Mer01,Ken02,CS11]: 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 ∪ V * , the set of white vertices is W r = W ↔ E, and |B r | = |W r |, see Figure 6 (left) for an example. The set of boundary rhombus pairs of the diamond graph G is R ∂ ; R ∂,r is the set R ∂ without the root pair containing the vertex r. For boundary rhombus pairs of R ∂ , we use the notation of Figure 7, which we recall here for convenience of the reader. v e iβ e iᾱ r e iβ r θ ∂ Figure 10: Notation for vertices of G D in a boundary rhombus pair of R ∂ , and rhombus vectors e iᾱ , e iβ , resp. e iᾱ r , e iβ r , of G /2 assigned to the edge (v c , w ), resp. (v c , w r ). When this pair is not the root one, the marked black edge is where the operator K ∂ (u) differs from K(u).
We introduce two versions K(u) and K ∂ (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 (12) restricted to G D,r , and we do not repeat the definition here.
The operator K ∂ (u) has specific boundary conditions arising from the forthcoming Theorem 35 related to the Ising model. It is the complex, bipartite Kasteleyn matrix associated to the weight function c ∂ (u) differing from c(u) along edges (w , v c ) of the boundary rhombus pairs of R ∂,r . For every edge wx of G D , we have In order for c ∂ (u) to be finite, we restrict the domain of u to: where the second equality is a consequence of properties of train-tracks, see Section 2.5.3. As a consequence, the weight function c ∂ (u) is everywhere non-zero.
The coefficients of K ∂ (u) differing from those of K(u) are those corresponding to edges (w , v c ) of R ∂,r , and we have: Remark 17. The weight c ∂ (u) w v c might be negative. Using the physics terminology, see also Remark 56, this means that the corresponding dimer model has vortices on boundary faces v c , w , f c , w r where this is the case.

Z-massive Laplacian operators
In the infinite case, we consider the Z-massive Laplacian associated to the isoradial graph G introduced in [BdTR17b], whose definition is recalled in Section 2.5.6. We also consider the dual Z-massive Laplacian ∆ 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 ∆ 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Ḡ * . For every vertex f of G * of degree d inḠ * , denote its neighbors by f 1 , . . . , f d . Letd be the degree of f in G * , then if f is a boundary vertex of G * we haved = d, and we label the vertices so that the firstd ones are common to G * andḠ * . For every j ∈ {1, . . . ,d}, (f, f j ) is an edge of G * and we letθ * j denote its half-angle in G ; for every j ∈ {d + 1, . . . , d}, (f, f j ) = (f, o) corresponds to an edge f w j of the double graph G D , where w j is the unique vertex of G D which belongs to the edge f o, we letθ * j denote the half-angle of the edge (f, w j ) in G /2 . Then, non-zero coefficients of the finite version of the dual Z-massive Laplacian ∆ m, * are given by, The corresponding conductances ρ * and masses m * are as defined in Section 2.4.
The second operator acts on C V r ; it is denoted ∆ m,∂ (u) and we restrict to u ∈ (T(k)) . It is not the natural finite version of the operator ∆ m but has boundary conditions inherited from those of the Z u -Dirac operator K ∂ (u), see the proof of Theorem 18 below. We use the notation of the infinite Z-massive Laplacian operator ∆ m , see Section 2.5.6, and also the notation of Figure 10 for the boundary rhombus pairs of R ∂ .
In the following, the two boundary vertices v r , v of R ∂ enter the same framework, and we denote such a vertex by v ∂ . For every v ∂ of R ∂ of degree d in G, let v 1 , . . . , v d denote its neighbors in cclw order with v 1 on the boundary of G on the left of v ∂ , see Figure 12 (right). Note that if v ∂ = v ∂,r belongs to the root pair, then v c = r is considered as a neighbor of v ∂,r . The coefficients of ∆ m,∂ (u) differing from those of ∆ m are those corresponding to edges (17) The corresponding conductances and masses are denoted ρ ∂ (u) and m ∂ (u). Note that the diagonal term ∆ m,∂ (u) v ∂ ,v ∂ has the same form as that of the natural finite version of the Z-massive Laplacian.
In order to state Theorem 18 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 R ∂,r . For such a vertex v c , there is one non-zero coefficient, corresponding to the column f c , given by

Relating the Z u -Dirac and the Z-massive Laplacian operators
Theorem 18 below relates the Z u -Dirac operator, the Z-massive Laplacian and the dual Zmassive Laplacian in the infinite and finite cases. In the infinite case, this extends to the full Z-invariant case the results of Section 6 of [Ken02] which correspond to k = 0.
• Infinite case. Let u ∈ (T(k)), then the Z u -Dirac operator K(u), the Z-massive Laplacian ∆ m and the dual ∆ m, * satisfy the following identity: • Finite case. Let u ∈ (T(k)) , then the Z u -Dirac operators K(u), K ∂ (u), the Z-massive Laplacian ∆ m,∂ (u) and the dual ∆ m, * satisfy the following identity: Remark 19.
• As a consequence of Theorem 18 in the infinite case, we have that • From the proof of Theorem 18 it follows that, in the finite case, we also have the following matrix relation: where ∆ m (u) is the natural finite version of the Z-massive Laplacian, i.e., where diagonal coefficients of boundary vertices of ∆ m (u) are defined as in the last line of (17).
Proof. As in the critical case [Ken02], the proof consists in showing that matrix coefficients on both sides are equal. We separate the infinite case together with the part of the finite case which enters the infinite framework, from the part of the finite case which is specific. To simplify notation, we omit the argument u from the operators.
When two black vertices x, y of G D , resp. G D,r , are at distance more than two, the corre- x,y , is trivially equal to 0, we thus suppose that x, y are at distance 2 or 0.
Infinite case and part of the finite case.
with v a vertex of G and f a vertex of G * , at distance two. Denote by w, w the two white vertices of the quadrangle of G D /G D,r containing v and f ; let e iᾱ , e iβ , resp. e iᾱ , e iβ , be the rhombus vectors of G /2 associated to the edge (w, v), resp. (w , v), and letθ, resp.θ , be the corresponding halfangle, see Figure 11. In the finite case, we moreover suppose that (v, f ) = (v c , f c ) for all boundary rhombus pairs of R ∂,r , then we have K ∂ = K for all coefficients involved. Let us prove that e iᾱ e iβ θθ Figure 11: Notation for computing Next, let (x, y) = (v 1 , v 2 ) be an edge of the graph G/G r corresponding to a path v 1 , w, v 2 of G D /G D,r . In the finite case, we moreover suppose that (v 1 , v 2 ) = (v c , v ) for all boundary rhombus pairs of R ∂,r , then we have K ∂ = K for all coefficients involved. Using the notation of Figure 9, we have We now handle the case where x, y are at distance 0, i.e., when x = y and we suppose that Remaining part of the finite case.
is not skew-symmetric, we have to consider the coefficients (v c , f ) and (f, v c ) separately. For the coefficient (f c , v c ), using the notation of Figure 10, we have This is equal to 0 as in the finite non-boundary case, because K ∂ is equal to K on the edges (w , f c ) and for some boundary rhombus pair of R ∂,r . We have, by the finite non-boundary computation.
When x, y are at distance 0, and x = y = v is a vertex of G, we are left with the case where v = v c for some boundary rhombus pair of R ∂,r . We have, by the finite non-boundary computation

Determinants of the Z u -Dirac and Z-massive Laplacian operators
We restrict to the finite case and consider the Z u -Dirac operators K(u) and K ∂ (u) of Section 3.1 (finite case), the Z-massive Laplacian ∆ m,∂ (u), and the dual Z-massive Laplacian ∆ m, * of Section 3.2 (finite case). Theorem 20 below proves that the determinants of K(u) and ∆ m, * are equal up to an explicit constant depending on u. Corollary 22 establishes a similar result for the determinants of K ∂ (u) and ∆ m,∂ (u), thus implying identities between partition functions.

Determinants as partition functions
Returning to Section 2.3.1 on Temperley's bijection, in particular to Equation (3), letc(u), be the weight function on directed edges of G,Ḡ * corresponding to the weight function c(u) of Equation (12). Then, the partition function of pairs of dual r-rooted and o-rooted directed spanning trees of G,Ḡ * with weight functionc(u), is equal to the partition function of the dimer model on the graph G D,r with weight function c(u) [Tem74,KPW00], which is equal to Kas61]. We proceed in a similar way with dimers weighted by c ∂ (u) and thus have, Returning to Section 2.4, we have that det ∆ m, * counts weighted rooted spanning forests of G * with conductances ρ * and masses m * , where conductances are symmetric. In a similar way, det ∆ m,∂ (u) counts weighted rooted directed spanning forests of G r with conductances ρ ∂ (u), m ∂ (u), where the dependence in u is along the boundary only, and the conductances are symmetric away from the boundary. That is,

Statements
Recall that every edge e = (w, x) of G D,r is assigned two rhombus vectors e iᾱe , e iβe of the diamond graph G /2 and a half-angleθ e . Moreover, every white vertex w of G D,r is in the center of a rhombus of the diamond graph G ; we letθ w be the half-angle of this rhombus at one of the two primal vertices.
Theorem 20. Let M 0 be a dimer configuration of G D,r . Then, for every u ∈ (T(k)), we have Comments on this theorem are given in the introduction to Section 3 especially how, in the critical case k = 0, it is an easy consequence of Temperley's bijection and the matrix-tree theorem, and how the argument does not directly extend to the non-critical case.
As one expects, the quantity e=wx∈M 0 [dn(u αe ) dn(u βe )] 1 2 is independent of the choice of perfect matching M 0 . To see this, it suffices to check that the alternating product around every inner quadrangle face of G D,r is equal to 1; this is proved in Lemma 23 below.
The following is an immediate corollary of Theorem 20.
Corollary 22. Let M 0 be a dimer configuration of G D,r . Then, for every u ∈ (T(k)) we have, Proof. From Theorem 18, we have: Replacing det K(u) in the above by Equation (20) and recalling that |W r | = |B r | = |V r | + |V * | and that |V r | = |V| − 1, yields the result. Note that an argument similar to that used for proving Theorem 20 would give an alternative direct proof of Corollary 22.

Proof of Theorem 20
The proof of Theorem 20 is postponed until the end of this section. It is a consequence of four intermediate results, see Equations (21), (24), (26), (27) below, which we now establish. We will use Section 2.3.1 on Temperley's bijection and Appendix A on gauge transformations.
From pairs of directed spanning trees to directed spanning trees. We start by defining a gauge transformation K g (u) of the matrix K(u). Being gauge equivalent, the determinants of K g (u) and K(u) are equal up to an explicit constant. We then use Temperley's bijection to deduce that the determinant of K g (u) counts weighted o-directed spanning trees of the dual graphḠ * .
Let g(u) be the weight function on white-to-black edges of G D,r defined by, for every edge wx of G D,r , Consider the matrix K g (u) obtained from K(u) by multiplying edge-weights by g(u). Returning to the definition of K(u), see (12) and (13), we obtain that non-zero coefficients of K g (u) are given by, Lemma 23. Consider a perfect matching M 0 of G D,r . The bipartite, weighted adjacency matrices K g (u) and K(u) are gauge equivalent and we have, Proof. Details on gauge equivalences for bipartite, weighted adjacency matrices are given in Section A.2 of Appendix A. To prove that K and K g are gauge equivalent, it suffices to show that the alternating products of K and K g around inner face-cycles of G D,r are equal. Inner face-cycles of G D,r are quadrangles; using the notation of Figure 11 we have, By definition of g we have, . The equality between determinants comes from [Kup98], see also Corollary 67, and from the fact that, wx∈M 0 sc(θ w ) = w∈W sc(θ w ), since a dimer configuration covers all white vertices W r of G D,r , and W r = W . Now, the matrix K g (u) is a complex, bipartite Kasteleyn matrix of the graph G D,r , where edges are assigned the positive weight function c g (u) given by: Returning to Section 2.3.1, the weight function c g (u) determines a weight functionc g (u) on directed edges of G,Ḡ * , see (3). Then,c g (u) v,v = 1, for all directed edges (v, v ) of G such that v is a vertex of G r . Let us now express the weight functionc g (u) on directed edges of G * using the rhombus vectors and half-angles assigned to those edges. Recall that an edge (f, f ) of the restricted dual G * is assigned two rhombus vectors 2e iᾱ , 2e iβ and a half-anglē θ * of G . An edge (f, o) ofḠ * corresponds to an edge f w of G D , and one assigns to (f, o) the rhombus vectors e iᾱ , e iβ and half-angleθ * of G /2 associated to the edge (f, w). Then, for every directed edge (f, f ) ofḠ * such that f is a vertex of G * we have, using the notation of (3):c Since edges of the primal graph have weight 1, by the KPW-Temperley bijection [Tem74,KPW00], the determinant of K g (u) counts weighted o-directed spanning trees of T o (Ḡ * ), where directed edges ofḠ * are assigned the weight function γ * (u) given by, for every directed That is, Matrix-tree theorem. An alternative way of computing the partition function (24) is to use the directed version of the matrix-tree theorem [Kir47,Tut48]. Let ∆ * (u) be the (nonmassive) Laplacian matrix of the graphḠ * with conductances γ * (u) on directed edges. The non-zero coefficients of the matrix ∆ * (u) are given by: Let ∆ * ,o (u) be the matrix obtained from ∆ * (u) by removing the row and column corresponding to the vertex o. Then, From o-directed spanning trees to rooted spanning forests. The last step consists in going from o-directed spanning trees ofḠ * counted by det ∆ * ,o (u) to rooted spanning forests of G * counted by det ∆ m, * using gauge equivalences on weighted adjacency matrices of digraphs, see Section A.1 of Appendix A.
Lemma 24. The weighted adjacency matrices (k ) − 1 2 ∆ * ,o (u) and ∆ m, * are gauge equivalent and we have, Proof. Both matrices (k ) − 1 2 ∆ * ,o (u) and ∆ m, * have the same associated digraph which is the restricted dual G * where a loop is added at every vertex, and each undirected edge is replaced by the two possible directed edges. The graph G * being connected, the associated digraph is strongly connected. We use Lemma 64 to prove gauge equivalence of the matrices.
For every vertex f of G * , define q f,f = 1. Next, consider two distinct vertices f 1 , f 2 of G * and a simple di-path γ from f 1 to f 2 . Set, The function q is well defined because it is the exponential function of [BdTR17b] evaluated at u − 2K − 2iK , see also Equation (11). Indeed, Fix a vertex f 0 of G * , and define D f 0 to be the diagonal matrix whose diagonal coefficient Recall the definition of the Laplacian matrix ∆ * ,o (u), see (25) and (23). For the diagonal coefficient corresponding to a vertex f of G * , we have For an edge (f, f ) of G * , we have By Lemma 64, Equation (28)  By Lemma 23 and Equation (24), | det K(u)| is equal up to a constant to | det K g (u)| which counts o-directed spanning trees of T o (Ḡ * ) with conductances γ * (u) on directed edges ofḠ * given by (23); that is: where C(u) is given in Lemma 23. By Temperley's bijection again, there is a one-to-one correspondence between dimer configurations of G D,r and o-directed spanning trees of T o (Ḡ * ) (such that the primal tree is r-rooted). We prove Z-invariance of this o-rooted directed spanning tree model onḠ * ; using the above, the decomposition of the partition function has a direct interpretation in terms of the dimer model on G D,r .
Since duality preserves isoradiality, we actually show Z-invariance of the r-directed spanning tree model on G, where directed edges of G are assigned conductances γ(u) given by, for every directed edge (v, v ) of G, where 2e iᾱ , 2e iβ ,θ are the rhombus vectors and half-angle of G associated to the edge (v, v ); γ(u) is the primal version of the conductances γ * (u) of (23).
Proposition 25. Consider a finite isoradial graph G, and let u ∈ (T(k)). Then, the model of r-directed spanning trees on G, with weight function γ(u) on the edges is Z-invariant.
• By Theorem 20, we have: and in the paper [BdTR17b], the model of rooted spanning forests with these weights is shown to be Z-invariant. But, since the proof of Theorem 20 does not provide a bijection between directed spanning trees and rooted spanning forests, the two decompositions of the partition functions are not directly comparable; they should nevertheless be compatible. Note that the computations of the proof of Proposition 25 are reminiscent but much simpler than those of [BdTR17b] (the latter have been removed from the published version).
• The critical spanning tree model of [Ken02] with conductances tan(θ) is Z-invariant [Ken99]. The result of [BdTR17b] extends this to rooted spanning forest while Proposition 25 extends it to directed spanning trees.
Proof. Let G Y and G ∆ be two finite isoradial graphs differing by a star-triangle transformation, and let r be a fixed root on the boundary of the graph, outside of the hexagon defining the star/triangle. Let γ ∆ (u), resp. γ Y (u), be the weight function on r-dST of G ∆ , resp. G Y . We use the notation of Figure 13, and write γ ∆ i,j /γ Y i,j for the weight of the edge (v i , v j ). Using the identities dn(u + K) = k nd(u), sc(θ * ) = sc(K − θ) = (k ) −1 cs(θ), we have the following: for every j ∈ {1, 2, 3}, with cyclic notation for indices: Let Z(G Y | i ), resp. Z(G ∆ | i ), be the partition function of G Y , resp. G ∆ , restricted to an outer configuration T belonging to the set defined in i, i ∈ {I, II, III}, divided by the contribution of the configuration T. Then, proving Z-invariance amounts to showing that there exists a constant C, such that Let us prove that this is indeed the case, with C = (k ) 3 2 3 j=1 sc(θ j ).
The proof that Z(G Y | III ) = C Z(G ∆ | III ) is concluded by using Equation (29) again.
3.6 Inverse Z u -Dirac operator and dimer model on the double graph Using Theorem 18, in Corollaries 27 and 29, we express the inverse Z u -Dirac operators K(u) −1 and K ∂ (u) −1 using the Z-massive and dual Z-massive Green functions. In the infinite case, this allows to prove in Theorem 31 an explicit local expression for a Gibbs measure for the dimer model on the double graph with weight function c(u).

Inverse Z u -Dirac operator and Z-massive Green functions
Infinite case. Consider an infinite isoradial graph G, and the Z u -Dirac operator K(u). Consider also the Z-massive Laplacian ∆ m on G, the dual Z-massive Laplacian ∆ m, * of the dual G * [BdTR17b]. When k = 0, let G m and G m, * be the Z-massive and dual massive Green functions of G and G * of [BdTR17b], whose definition is recalled in Section 2.5.6. When k = 0, the mass is 0 and we let G 0 and G 0, * be the Green and dual Green functions of [Ken02].
Notation for coefficients of Corollary 27. Letv, resp.f , be a vertex of G, resp. G * . Let w be a white vertex of G D , its neighbors in G D are v 1 , f 1 , v 2 , f 2 , and let e iᾱ f , e iβ f ,θ f be the rhombus vectors and half-angle of G /2 associated to the edge (w, v 2 ), where the subscript "f" stands for "final", see Figure 14.
vf Figure 14: Notation for coefficients of Corollary 27 and 29.
As a consequence of Theorem 18 (infinite case) we obtain, Corollary 27. For every u ∈ (T(k)), consider the operator K(u) −1 mapping C W to C B defined by: • Coefficients. For everyv,f, w as in the notation above, Then K(u) −1 is an inverse of the Kasteleyn operator K(u). When the graph G D is moreover Z 2 -periodic, it is the unique inverse decreasing to zero at infinity.
• When k = 0, then dn ≡ 1 and we recover Corollary 7.2 of [Ken02]. When the graph is Z 2 and the weights are specific, Chhita [Chh12] has a result in the same flavor, relating the inverse Kasteleyn operator of pairs of dual directed spanning trees to a massive Green function.
• The operator K(u) −1 is local. Indeed the expression for K(u) −1 x,w , withx =v ∈ V or x =f ∈ V * , only depends on: two paths of the diamond graph G fromx to the two neighbors x 1 , x 2 of w in G or G * , where the two paths are those used in computing the massive exponential function of coefficients of the massive Green function of [BdTR17b]; and from the rhombus vectors and half-angles of G /2 associated to the edges (w, x 1 ), (w, x 2 ).
Proof. By definition, see Section 2.2.2, to show that K(u) −1 is an inverse, we need to prove that K(u) −1 K(u) = Id, and that K(u) −1 x,w → 0 as |w − x| → ∞. For the second part, when k = 0 we prove in [BdTR17b] that the Green function G m and G m, * decrease exponentially fast to 0 at infinity; since the function dn is uniformly bounded, the same holds for K(u) −1 . When k = 0, the Green function explodes like − log |x − x | at infinity, but the difference converges polynomially fast to 0 as is proved in [Ken02].
For the first part, we apply K(u) to the right of the definition of K(u) −1 in matrix form and use Theorem 18. Note that this step also uses associativity of the infinite matrix product, which holds in this case.
Finite case. Consider a finite isoradial graph G and the Z u -Dirac operators K(u) and K ∂ (u) of Section 3.1 (finite case). Consider also the Z-massive Laplacian ∆ m,∂ (u) and dual Laplacian ∆ m, * of Section 3.2 (finite case). For the purpose of handling the Ising model, our goal is to obtain an explicit expression for the inverse of the operator K ∂ (u). In order to ensure that K ∂ (u) and ∆ m,∂ (u) are invertible, we restrict u to: Indeed, by the forthcoming Remark 37 this ensures that det K ∂ (u) = 0 and by Corollary 22 that det ∆ m,∂ (u) = 0 (since nd is positive when u is real). Denote by G m, * the inverse of the (finite) matrix ∆ m, * and by G m,∂ (u) the inverse of ∆ m,∂ (u). Note that when some of the conductances (ρ ∂ v c ,v (u)) of the boundary rhombus pairs of R ∂,r are negative, there is a twist in defining the associated random walk and in giving the random walk interpretation of the Green function, but this might at most happen for half of the boundary edges.
Notation for coefficients of Corollary 29. Let w be a white vertex of G D,r ; if w / ∈ {w , w r ∈ R ∂ }, then we use the notation of the infinite case, see Figure 14; if w ∈ {w , w r ∈ R ∂ }, then the vertex f 2 is absent. As an immediate consequence of Theorem 18 (finite case) we obtain Corollary 29. For every u ∈ (T(k)) , the inverse matrix K ∂ (u) −1 has the following explicit expression.
• Coefficients. For every vertexv of G r ,f of G * , and every white vertex w of G D,r , using the notation of Figure 14, where I {w / ∈{w ,w r ∈R ∂ }} is equal to 0 if w is a boundary vertex of G D,r and 1 otherwise. In the sum over (v c , f c ) ∈ R ∂,r , we use the notation of Figure 10 for vertices and rhombus vectors of the boundary rhombus pairs of R ∂,r ; K ∂ (u) −1 v c ,w in the formula for K ∂ (u) −1 f,w is given by the first formula withv = v c .
Example 30. Let us compute the explicit values of Corollaries 27 and 29 in the case where u = u f := α f +β f 2 + K. This will be used again in Sections 4.3 and 5.3. In the finite case, we will also need to evaluate it at u = v f : Using that dn(u ± K) = k nd(u), dn(u ± 2K) = dn(u), dn(−u) = dn(u) and the above relations, we obtain, for u ∈ {u f , v f }, an thus, with the notation of Corollaries 27 and 29: • Infinite case.

Dimer model on an infinite isoradial double graph G D
Suppose that the isoradial graph G is infinite; when it is moreover Z 2 -periodic, we consider the natural exhaustion (G D n ) n≥1 of G D by toroidal graphs, where G D n = G D /nZ 2 . Let F denote the σ-field generated by cylinder sets of G D . Using arguments of [CKP01, KOS06,dT07], we obtain an explicit, local expression for a Gibbs measure of the dimer model on G D with weight function c(u). Since the proof closely follows that done in the papers [dT07, BdT11, BdTR17a], we do not repeat it here. The key requirements are that the operator K(u) −1 is local and unique in the Z 2 -periodic case.
Theorem 31. For every u ∈ (T(k)), there exists a unique probability measure on (M(G D ), F), denoted P D,u dimer , such that the probability of occurrence of a subset of edges E = {w 1 x 1 , · · · , w l x l } in a dimer configuration of G D is given by: where (K(u) −1 ) E is the sub-matrix of K(u) −1 given by Corollary 27, whose rows are indexed by x 1 , . . . , x l and columns by w 1 , . . . , w l . The measure P D,u dimer is a Gibbs measure. Moreover, when the graph G D is Z 2 -periodic, the probability measure P D,u dimer is obtained as weak limit of the Boltzmann measures on the toroidal exhaustion (G D n ) n≥1 . Remark 32.
• As mentioned in the introduction to Section 3, this theorem is a directed version of the transfer impedance theorem of [BP93]. A result in the same flavor, i.e., computing probabilities of pairs of directed spanning trees using the massive Green functions of massive non-directed random walks, is obtained by [Chh12] in the case where G = Z 2 with specific weights.
• A version of this theorem in the finite case can be obtained using Remark 19.
Example 33. As an example of application we express the probability of single edges occurring in dimer configurations of G D chosen with respect to the measure P D,u dimer , using the Z-massive and dual Z-massive Green functions of [BdTR17b]. We use the notation of Figure 14 and omit the subscript "f" since there is no confusion possible between the initial and final vertices. Details of computations are given in Appendix B.1.

Kasteleyn operator of the graph G Q and Z u -Dirac operator
Let G be an isoradial graph, infinite or finite. We consider the isoradial embedding of the bipartite graph G Q = (V Q , E Q ) given in Section 2.5.2, and the weight function ν J of Equation (8) arising from the Z-invariant Ising model. Let K Q be the associated complex, bipartite Kasteleyn matrix defined in Section 2.5.7, with rows indexed by black vertices. In the finite case we moreover consider the diagonal matrix D Q,B , resp. D Q,W , whose rows/columns are indexed by black/white vertices of G Q , and whose diagonal coefficients are: Let K Q be the modified, complex, bipartite Kasteleyn matrix defined by that is, K Q is obtained from K Q by multiplying the weight of the edges b w , b r w r of all boundary rhombus pairs of R ∂ by sn(θ ∂ ).
Consider also the isoradial embedding of the double graph G D of Section 2.5.2. For every u ∈ (T(k)), let K(u) be the Z u -Dirac operator and, in the finite case, for every u ∈ (T(k)) , let K ∂ (u) be the Z u -Dirac operator with specific boundary conditions, as defined in Equation (12).
The main result of this section, and actually the key result of this paper, is Theorem 35 of Section 4.1 proving, for every u ∈ (T(k)), explicit linear relations between the matrices K Q and K(u) in the infinite case, and between K Q and K ∂ (u) in the finite case.
In Section 4.2 we restrict to the finite case; using Theorem 35 and a combinatorial argument, we prove in Theorem 36 that the determinant of the Kasteleyn matrix K Q is equal, up to an explicit multiplicative constant depending on u, to the determinant of the Z u -Dirac operator K ∂ (u). Combining this with Theorem 20 proves that the determinant of K Q is equal, up to an explicit constant, to the determinant of the Z-massive Laplacian ∆ m,∂ (u). Interpreting these determinants as partition functions, Theorem 36 proves that the squared partition function of the Z-invariant Ising model with + boundary conditions is equal, up to an explicit constant, to the partition function of weighted rooted directed spanning forests counted by ∆ m,∂ (u), where the dependence in u is along the boundary only. This generalizes to the full Z-invariant case the results of [dT13, dT16] proved in the Z-invariant critical case, and to the case of simply connected domains the result of [BdTR17a] proved in the toroidal case. The proof we provide here has a slight combinatorial flavor but is mainly based on matrix relations, so quite different from [dT13, dT16]. Note that the combinatorics argument of [dT16] can be generalized to the full Z-invariant case and would give an alternative proof. Note also that the boundary trick of Chelkak and Smirnov [CS12] allows us to remove dual trees along the boundary which we could not do in [dT16].
Using Theorem 35, Corollaries 45 and 46 of Section 4.3 prove linear relations between the inverse operator (K Q ) −1 and the inverse Z u -Dirac operator. Choosing specific values of u allows us to express the dimer measure of the graph G Q using the inverse Z u -Dirac operator and the Z-massive Green functions, see Corollaries 48 and 50. This also provides an alternative direct way of finding a local formula for (K Q ) −1 [BdTR17a], explicitly relating it to the Z-massive Green functions.

Relating the Kasteleyn operator K Q and the Z u -Dirac operator
The main result of this section is Theorem 35 proving an explicit relation between the matrices K Q and K(u) in the infinite case, and between K Q and K ∂ (u) in the finite case. In order to state this theorem we need to introduce two additional matrices S(u) and T (u). Both of them are "rectangular" with "twice" more rows than columns.
The matrix S(u) has rows indexed by black vertices of G Q and columns by white vertices of G D , resp. of G D,r , in the infinite case, resp. finite case. If G Q is infinite, let b be a black vertex; if G Q is finite, let b be a black vertex of an inner quadrangle. Let w be the white vertex of G D corresponding to the quadrangle to which b belongs. Then, the only non-zero coefficient of the row corresponding to b is: with the following notation, see Figure 15 (left): w is the white vertex of G Q such that bw is parallel to an edge of G; e iᾱ , e iβ ,θ are the rhombus vectors and half-angle of G /2 assigned to the edge (b, w).
Suppose that G Q is finite, and let b be a black vertex of a boundary quadrangle of G Q . Then b ∈ {b r , b } for some boundary rhombus pair of R ∂ , see Figure 15 (right), and the non-zero coefficient of the row corresponding to b is: where for b r the definition is coherent with that of the non-boundary case, since the rhombus vectors of G /2 assigned to the edge (b r , w r ) are e iᾱ r , e iβ r . For b , the rhombus vectors assigned to the edge (b , w ) are e iᾱ , e iβ so that there is a change of definition; this comes from our choice of embedding of G Q which exchanges the bipartite coloring of vertices in the left rhombus of the pair. The matrix T (u) has rows indexed by white vertices of G Q and columns by black vertices of G D , resp. of G D,r , in the infinite case, resp. finite case. If G Q is infinite, let w be a white vertex; if G Q is finite, let w be a white vertex such that w = w c for all boundary rhombus pairs of R ∂ . The vertex w is on a rhombus edge vf of the diamond graph G , where v is a vertex of G and f a vertex of G * , see Figure 15 (left). Then the row of T (u) corresponding to w has two non-zero coefficients defined by, where e iβ is the rhombus vector (w, v) and e i(β+π) is the rhombus vector (w, f ).
Suppose that G Q is finite, and let w = w c for some boundary rhombus pair of R ∂ , then w c is on a rhombus edge v c f c of G , see Figure 15 (right). As long as the rhombus pair is not the root one, i.e., the one where v c = r, the row of T (u) corresponding to w c has non-zero coefficients given by, When the boundary rhombus pair is the root pair, then only the term t(u) w c ,f c is defined. Note that the definition is specific for v c , and coherent with the non-boundary case for f c since the rhombus vector (w c , f c ) is e iβ .
We are now ready to state our main result.
• Infinite case. Let u ∈ (T(k)), 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: • Finite case. Let u ∈ (T(k)) , 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: Proof. In the whole of the proof, we omit the argument u from the matrices.
Infinite case and finite non-boundary case. Figure 16 below sets the notation. Let b be a black vertex of G Q , then b belongs to a quadrangle corresponding to a vertex w of G D . If G Q is finite, suppose further that the quadrangle is not a boundary one, or equivalently that w is not a boundary vertex of G D . Let v 1 , f 1 , v 2 , f 2 be the four black vertices of G D incident to w. Denote by w 1 , w 2 , w 3 the three white vertices of G Q incident to b, and let e iᾱ , e iβ be the rhombus vectors of the edge (b, w 1 ).
The coefficient [K Q T ] b,x of the LHS of (37) and (38) is non-zero only when x ∈ {v 1 , v 2 , f 1 , f 2 }, and Replacing coefficients of K Q yields for the LHS: Note that in the finite case, we have K ∂ = K for coefficients involved. The coefficient [S K] b,x of the RHS of (37) and (38) is also non-zero when x ∈ {v 1 , v 2 , f 1 , f 2 } and we have [S K] b,x = s b,w K w,x . To compute these terms, we first express the part cn(u β )[nd(u α ) nd(u β )] 1 2 := ( ) in s b,w , see (33), using the angles and parameters involved in the four coefficients of K.
where in the first line we used that u α = K + u α+2K , in the third that u β = −(u β−2K ) * , u α = K − (u α ) * , and in the fourth that u β = K − (u β ) * , u α = 2K − (u α+2K ) * . Replacing coefficients of K by their definition gives for the RHS of (37), x is then straightforward when x = v 2 and x = f 2 . When x = v 1 , this is a consequence of the identity, see [Law89, chap.2, ex.32 (i)], We are left with proving the case x = f 1 . Multiplying the identity (39) by nd(u − v), using The proof is concluded by evaluating the above at Finite boundary case. The notation used are those of Figure 15 (right). Let b ∈ {b , b r } be a black vertex of some boundary rhombus pair of R ∂ . Suppose first that this rhombus pair is not the root one. Then, if b = b r , resp. b = b , the coefficients of the LHS and RHS of (38) are non-zero when x ∈ {v c , v r , f c }, resp. x ∈ {v c , v , f c }. We need to prove: In both cases, the last two equalities are as in the full plane with the appropriate change of notation. We thus need to check the first equality of each case. Returning to the definition of K Q and boundary values of the matrices S and T defined in Equations (34) and(36) we have, In a similar way, Note that we always have s b r ,w r , s b ,w = 0 on (T(k)) , see (15), so that it makes sense to divide by these quantities.
The last case we need to consider is if b belongs to the root boundary rhombus pair of R ∂ . But then, of the three equations above, the first one is absent since we have v c = r, so we are left with the last two which are as in the full plane case. This ends the proof of (38) and thus finishes the proof of Theorem 35.

Determinants of the Kasteleyn matrix K Q and of the Z u -Dirac operator
We restrict to the finite case. Theorem 36 proves that the determinants of the matrices K Q and K ∂ (u) are equal up to an explicit multiplicative constant depending on u. By [Dub11], see also Equation (5) In the whole of this section, we restrict the domain of u to Note that by Section 2.5.3 on train-tracks, this amounts to removing all the parallel directions of the train-tracks of the isoradial graph G.

Results
Every edge e = (w, x) of G D,r is assigned two rhombus vectors e iᾱe , e iβe of the diamond graph G /2 and a half-angleθ e . We partition the set of white vertices W r = W of G D,r as W ∂ ∪ W • , where W ∂ consists of boundary vertices of W and W • of inner ones. Every white vertex w of G D,r is in the center of a rhombus of the diamond graph G ; we letθ w be the half-angle of this rhombus at one of the two primal vertices.
In the statement below, a specific role is played by the boundary rhombus pair of R ∂ containing the root r. We use the notation of Figure 15 (right) and add a superscript r to specify vertices/angles of this root pair.
Theorem 36. Let M 1 be a dimer configuration of G D,r . Then, for every u ∈ (T(k)) , we have whereC(u) is equal to: Remark 37.
• Since for u ∈ (T(k)) , sc(u αe ) sc(u βe ) = 0, and since det K Q = 0, we have that the matrix K ∂ (u) is invertible for these values of u.
• The quantity wx∈M 1 :w∈W | sc(u αe ) sc(u βe )| 1 2 is independent of the choice of perfect matching M 1 . Similarly to Remark 21, this consists in showing that the alternating product around every inner face of G D,r is equal to 1; the proof is similar to that of Lemma 23.
• A surprising fact is that the LHS is independent of u while the RHS does not seem to be, also the RHS seems to depend on the choice of root vertex r. It is not straightforward to see why this indeed not the case. An alternative way of proving this theorem is to extend to the full Z-invariant case the combinatorial argument of [dT16]; one then better sees the parameter u and the choice of root vertex r appearing.
Combining Theorems 36 and 20, we deduce that the squared partition function of the Zinvariant Ising model with + boundary conditions is equal, up to an explicit constant, to the determinant of the massive Laplacian ∆ m,∂ (u), i.e., to the partition function of rooted directed spanning forests. This generalizes to the full Z-invariant case the result of [dT13, dT16] and to simply connected domains the result of [BdTR17a] proved for the characteristic polynomial in the toroidal case in an abstract way.
Corollary 38. For every u ∈ (T(k)) , we have sn(θ ∂,r ) cd(u β r,r ) sn(u β r,r ) cd(u α r,r ) sn(u α r,r ) Proof. Equality 2. is obtained by writing [Z + 1 2 , using Equality 1. evaluated once at u and once at u+2K, and using the identities | sc(u−K)| = Let us prove Equality 1. From Equation (5) we have, for all coupling constants J, Returning to the definition of the Z-invariant weights (6) and (8) gives, where in the second equality we use the notation introduced before the statement of Theorem 36. Since Z dimer (G Q , ν J ) = | det K Q |, Theorem 36 yields The proof is concluded by using Corollary 22, choosing M 1 as perfect matching of G D,r and using that |V| − 1 = −|V * | + |E|.

Proof of Theorem 36
Let us prove that for the modified Kasteleyn matrix K Q , we have To obtain the result for K Q , we use Relation (32) which implies that The proof has three main steps: the first consists in using a partition of black/white vertices of G Q and Theorem 35 for comparing det K Q and det K ∂ (u); in the second step, we specify this partition so that the computation required by the first step become tractable; finally, we perform these computations in the third step.
Partition of the vertices of G Q and Theorem 35. We partition black vertices of G Q into two subsets: B = B 1 ∪ B 2 , where B 1 has one black vertex per quadrangle of G Q and B 2 = B \ B 1 . Since boundary quadrangles of G Q are reduced to edges, B 1 contains all black vertices of boundary quadrangles, and half of the black vertices of inner quadrangles. As a consequence, B 1 has a natural partition as B ∂ 1 ∪ B • 1 , where B ∂ 1 , resp. B • 1 , consists of boundary quadrangle, resp. inner quadrangles, black vertices. Then B 2 has no boundary quadrangle black vertices. In a similar way, we partition white vertices of G Q : Recalling that there is a natural bijection between quadrangles of G Q and white vertices of G D,r , we have the following natural bijections: Here are some notation for sub-matrices of the matrices S(u) and T (u) defined in Section 4.1, and for sub-matrices of the matrix K Q .
Using Theorem 35, we obtain the following lemma.
Lemma 39. Consider a partition of the black and white vertices of G Q as above. Then, for every u ∈ (T(k)) , Proof. Note that S 2 (u) and S • 1 (u) are invertible because we choose u ∈ (T(k)) . By Theorem 35, we have the following identity: . We extractR(u) W 2 W • from the first equation. PluggingR(u) W 2 W • in the second equation gives R(u); taking the determinant ends the proof.
Combinatorial partition of the vertices of G Q . Let us specify the partition of the black/white vertices of G Q . It is constructed from a well chosen perfect matching M 1 of G D,r , which we now define. Recall that by Temperley's bijection, perfect matchings of G D,r are in bijection with pairs of dual directed spanning trees of T r,o (G,Ḡ * ), see Section 2.3.1. Consider a spanning tree of the restricted dual G * , and root it at the vertex f c,r incident to the root vertex v c = r in the diamond graph G . To this spanning tree, add the edge (f c,r , o) which is the dual of the edge rv ,r of G. This defines an o-directed spanning tree T * ofḠ * . Consider the r-directed spanning tree T of G which is the dual of T * , rooted at the vertex r. Then (T, T * ) is a pair of dual spanning trees of T r,o (G,Ḡ * ), and we let M 1 be the corresponding perfect matching of G D,r , see Figure 17. We now define the partition of black/white vertices of G Q arising from M 1 ; this amounts to specifying the partition for vertices of inner quadrangle since those of boundary quadrangles define B ∂ 1 /W ∂ 1 . For j ∈ {1, 2}, b j /w j denotes a black/white vertex of B j /W j . Consider an inner quadrangle of G Q corresponding to a white vertex w of G D,r . Then, exactly one edge of the quadrangle is crossed by an edge wx of the perfect matching M 1 . The partition is defined as follows: if x = v ∈ V r , then w 1 is the white vertex on the right of the edge (w, v), and b 1 is the black vertex on the left; if x = f ∈ V * , then w 1 is the white vertex on the left of the edge (w, f ), and b 1 is the black vertex on the right. An example is provided in Figure 17 (right) and Figure 19.
Let us prove a combinatorial lemma. Note that because u ∈ (T(k)) , the fact that a coefficient of T 1 (u) or R(u) is non-zero is independent of u. As in Section A.2 of Appendix A, to the matrix T 1 (u) corresponds a bipartite graph G(T 1 ) = (W 1 ∪ B r , E(T 1 )), where there is an edge w 1 x iff t(u) w 1 ,x = 0, with w 1 ∈ W 1 , x ∈ B r = V r ∪ V * . In a similar way, to the matrix R(u) of Lemma 39 corresponds a bipartite graph G(R) = (W • ∪ W 2 , E(R)) where there is an edge ww 2 iff r(u) w,w 2 = 0, with w ∈ W • , w 2 ∈ W 2 . The matrix T 1 (u), resp. R(u), is then a bipartite, weighted adjacency matrix of the graph G(T 1 ), resp. G(R). r Figure 18: Example of graph G(T 1 ) (red edges) and of graph G(R) (blue edges).
Lemma 40. Consider a partition of the black/white vertices of G Q arising from the perfect matching M 1 of G D,r . Then, the graph G(T 1 ) is a spanning tree on the vertex set W 1 ∪ B r , and the graph G(R) is a union of trees on the vertex set W • ∪ W 2 , spanning all vertices of Proof. Let us prove that G(T 1 ) is a spanning tree. Every vertex w 1 of W 1 has degree 2 in G(T 1 ): it is adjacent to a vertex x of B r such that wx is an edge of the perfect matching M 1 , and to a vertex x such that xx is an edge of the diamond graph G with w 1 in its middle. Using the natural bijection W 1 ↔ W , that is, identifying every vertex w 1 with the corresponding vertex w of W r does not change the combinatorics of the graph G(T 1 ). With this identification, the graph G(T 1 ) contains all edges of the perfect matching M 1 , and the second edge wx incident to w is such that (w, x ) is on the right, resp. left, of (x, w) if x ∈ V * , resp. x ∈ V r . Then, by Proposition 7.3. of [dT16], the graph G(T 1 ) is a spanning tree. An example of G(T 1 ) is pictured with red edges in Figure 18.
Let us prove that G(R) is a union of trees spanning all vertices of W • ∪ W 2 . By definition, a white vertex w of W • is adjacent to all white vertices of W 2 that are adjacent to the black vertices b 1 and b 2 of the quadrangle w corresponds to; that is, w is adjacent to the vertex w 2 of the quadrangle of w and maybe to another vertex of W 2 . As a consequence G(R) is spanning all vertices of W • ∪ W 2 . It cannot contain a cycle for otherwise it would mean that there is a vertex of B r which does not belong to M 1 , which is a contradiction with it being a perfect matching. An example of G(R) is pictured with blue edges in Figure 18.
Corollary 41. Consider a partition of the black/white vertices of G Q arising from the perfect matching M 1 of G D,r .
• Using the identification W 1 ↔ W r = W , we have, Proof. Writing the determinant as a sum over permutations, we have that non-zero terms in the expansion of det T 1 (u), resp. det R(u), correspond to perfect matchings of the bipartite graph G(T 1 ), resp. G(R). Since these two graphs are trees or union of trees spanning all vertices, they have at most one perfect matching; indeed if they had more, the union of two different ones would yield a cycle which is in contradiction with being a tree. Using the identification W 1 ↔ W , The graph G(T 1 ) has one perfect matching given by edges of M 1 (pictured in thick red lines in Figure 18), while the graph G(R) has one perfect matching given by the natural identification W • ↔ W 2 (pictured in thick blue lines in Figure 18). Since, the contribution of a perfect matching to the determinant is the product of the edge-weights (up to a sign), this ends the proof of the corollary.
Since S ∂ 1 (u), S • 1 (u), S 2 (u) are diagonal matrices, their determinant is the product of the diagonal terms. Combining Lemma 39 and Corollary 41 we thus obtain the following.
Corollary 42. Consider a partition of the black/white vertices of G Q arising from the perfect matching M 1 of G D,r . Then, for every u ∈ (T(k)) , where r (u) w, Computation of (I)(II) in Identity (41). For every u ∈ (T(k)) , define the weight function η(u) on edges of G D,r as follows.
Lemma 43. The product (I)(II) is equal to: Proof. In the whole of the proof, we simply denoteθ w byθ, and omit the argument u from matrix coefficients.
We first handle Part (II) involving inner vertices. Let w ∈ W • and wx be an edge of the perfect matching M 1 , with x = v or f . By definition of the partition of black/white vertices arising from M 1 , we have b 1 , b 2 , w 1 , w 2 as in Figure 19. Let v be the primal vertex such that the edge v v crosses the quadrangle, and let 2e iᾱ , 2e iβ be the two rhombus vectors of G associated to the edge (v , v). Figure 19: Notation for computing r w,w 2 tw 1 ,x , when w ∈ W • and x = v (left), x = f (right).
We compute the term r w,w 2 when x = v. Using the notation of Figure 19, the rhombus vectors of G /2 assigned to the edge (w, v) are e iᾱe = e iᾱ , e iβe = e iβ . Returning to the definition of the matrices K Q and S, see (33), we have As a consequence, By definition of T , see (35), we have t w 1 ,v = e −iβ 2 cn(u β ) = e −iβ e 2 cn(u βe ), and we conclude that We now turn to the term r w,w 2 in the case where x = f . The rhombus vectors of G /2 assigned to the edge (w, f ) are e iᾱe = e i(β−π) , e iβe = e iᾱ . Moreover, referring to Figure 19, we see that taking x = f has the effect of exchanging b 1 and b 2 and leaving w 1 , w 2 fixed. The quantity r w,w 2 being skew-symmetric in b 1 , b 2 , we have that r w,w 2 is equal to the opposite of (43). As a consequence, using that e iᾱe = e i(β−π) , e iβe = e iᾱ , and that dn(u − K) = k nd(u). By definition of T , we have t w 1 ,f = e −iβ −π 2 cd(u β−2K ) = e −iᾱ e 2 cd(u αe ), and we deduce that Summarizing, we have proved that, for every edge wx ∈ M 1 such that w ∈ W • , and thus, We now compute Part (I) involving boundary vertices of W ∂ . We will be using the notation of Figure 15 for vertices, rhombus vectors and angles of boundary rhombus pairs of R ∂ , and add a superscript r when the pair is the root pair, i.e., the one where v c = r. With our choice of perfect matching M 1 , the boundary contribution (I) of (41) can be rewritten as, see also Figure 17 (right): Suppose first that w r , w ∈ W ∂ \ {w r,r , w ,r }. Recalling the definition of S and T along the boundary, see (34) and (36), we have On the other hand, by definition of η, the product of weights of the edges w r v r , w v c of the perfect matching M 1 is equal to, writing u α +2K = u α −K, using elliptic trigonometric identities and the fact thatᾱ r =β [2π].
so that, Let us now consider the term s b r,r ,w r,r s b ,r ,w ,r t w r,r ,v r,r t w c,r ,f c,r . For the purpose of this computation, it is useful to imagine that the vertex v c = r is present and that t w c,r ,v c,r , η w ,r v c,r are defined as for the other pairs of rhombi. We then have, omitting to write the superscript r, The product of the first two terms is equal to | cn(θ ∂ )| by the above computation. Then, returning to the definition of the weight function η, we have Returning to the definition of the matrix T along the boundary, we have using thatβ =ᾱ r [2π]. Putting the three computations together, and writing the superscript r again, we deduce that s b r,r ,w r,r s b ,r ,w ,r t w r,r ,v r,r t w c,r ,f r 1 η w r,r v r,r η w ,r f r = sn(θ ∂,r ) cn(θ ∂,r ) cd(u β r,r ) sn(u α r,r ) , and thus (I) = sn(θ ∂,r ) cn(θ ∂,r ) cd(u β r,r ) sn(u α r,r ) using that W ∂ = {w , w r ∈ R ∂ }. Combining (44) and (45) allows to conclude the proof of Lemma 43.
The next lemma proves a simplified expression for the product of the weights η wx in Lemma 43.
Lemma 44. For every u ∈ (T(k)) , we have the following identity, where the weight function η is defined in (42).
Proof. We have the following identities: As a consequence, for every edge wx of G D,r , the weight function η wx can be rewritten as: Now consider a vertex w of G D,r and, using the notation of Figure 16, the corresponding rhombus v 1 , f 1 , v 2 , f 2 of the diamond graph G . Introduce the following notation for the rhombus vectors: that is, the notation α, resp. β, is for vectors on the right, resp. left, of the primal edge of the rhombus. With this notation, the weight η wx can be written as η wx = sn(u α 1 (w) ) sn(u α 2 (w) ) sn(u β 1 (w) ) sn(u β 2 (w) ) 1 2 , and this, independently of whether x = v or f . Let us prove that wx∈M 1 :w∈W η wx = w∈W sn(u α 1 (w) ) sn(u α 2 (w) ) sn(u β 1 (w) ) sn(u β 2 (w) ) 1 2 = 1.
Because of (46), the product over white vertices W can be seen as a product over rhombus vectors of G . Then, every inner rhombus vector of G occurs twice exactly and contributes once to the numerator and once to the denominator, so that the contributions cancel. Boundary rhombus vectors of G occur once but, referring to Section 2.5.3 on train-tracks, we know that they come in parallel pairs and contribute once to the numerator and once to the denominator; the contributions thus also compensate ending the proof of this lemma.
Putting together Corollary 42, Lemmas 43 and 44 ends the proof of Theorem 36.

Dimer model on the graph G Q and inverse Z u -Dirac operator
Using Theorem 35, in Corollaries 45 and 46, we prove linear relations satisfied by the inverse Kasteleyn operator (K Q ) −1 and the inverse of the Z u -Dirac operators K(u) and K ∂ (u). Section 4.3.2 is about applications of these results to the dimer model on G Q . In particular, when the graph G Q is infinite, we prove an alternative way of obtaining a local formula for the inverse [BdTR17a] which is seen as directly related to the Z-massive Green functions.

Inverse Kasteleyn operator (K Q ) −1 and inverse Z u -Dirac operator
Infinite case. In the paper [BdTR17a] we prove an explicit local expression for an inverse (K Q ) −1 of the operator K Q , which decreases to 0 exponentially fast in the distance when k = 0, and as the inverse distance when k = 0. When k = 0, the local expression is actually computed in [Ken02]. When the graph G Q is Z 2 -periodic, the operator (K Q ) −1 is the unique inverse decreasing to 0 at infinity. In Corollary 45 below, we use the existence and uniqueness of this inverse operator but not the explicit expression; we also need the following notation.
Notation for coefficients of Corollary 45. Letw be a white vertex of G Q andv,f be its adjacent vertices in the diamond graph G /2 , such thatv ∈ V andf ∈ V * . Denote by e iβ i the rhombus vector corresponding to the edge (w,v). Let w be a white vertex of G D and b, b be the black vertices of G Q of the corresponding quadrangle. To the vertex b, we assign the rhombus vectors e iᾱ f , e iβ f of G /2 of the edge (b, w), where the vertex w is such that the edge bw is parallel to an edge of G. Then, the rhombus vectors assigned to the vertex b are e iᾱ f +π , e iβ f +π , see Figure 20; the subscripts "i" and "f" stand for "initial" and "final".
As a consequence of Theorem 35 (infinite case), we obtain the following. Corollary 45. For every u ∈ (T(k)), as long as they are unique (which happens for sure in the Z 2 -periodic case), the inverse operators (K Q ) −1 , K(u) −1 , and the matrices S(u), T (u) satisfy the following identity.
• Coefficients. For everyw,v,f , and every w, b, b as in the notation above, we have: Or equivalently, Proof. The matrix form is obtained by left multiplying by (K Q ) −1 and right multiplying by K(u) −1 Equation (37) of Theorem 35. We are allowed to do so because, coefficients of the inverses decrease to 0 at infinity and the other matrices involved only have finitely many non-zero terms per row and column, implying associativity of the infinite matrix products. We also use uniqueness of the right (or left) inverse. Indeed, together with the fact that the products K Q (K Q ) −1 K Q and K(u)K(u) −1 K(u) are associative, this implies that they each are inverses on both sides [Coo14].
For coefficients, we return to the definition of the matrix S(u), see (33), and obtain using that cn(u − K) = k sd(u − K), nd(u − K) = (k ) −1 dn(u) in the penultimate line.
Returning to the definition of coefficients of the matrix T , see (35), we have using that cd(u − K) = sn(u) in the last equality. This ends the proof of the formula for coefficients.
Finite case. We restrict to u ∈ (T(k)) so that the Z u -Dirac operator K ∂ (u) is invertible.
Notation for coefficients of Corollary 46. As in the notation for coefficients of Corollary 45, we add a subscript i/f for rhombus vectors and half-angles of initial/final vertices. If w ∈ {w c ∈ R ∂ }, or w ∈ {w , w r ∈ R ∂ }, we thus have the notation of Figure 21. Figure 21: Notation for boundary coefficients of Corollary 46.
As a consequence of Theorem 35 (finite case), we obtain the following.
• Coefficients. We have two cases to consider.
1. For everyw,v,f ; for every w, b, b such that w / ∈ {w , w r ∈ R ∂ }, using the notation of Figure 20 and 21 (left), we have 2. For everyw,v,f , for every w such that w ∈ {w , w r } for one of the rhombus pairs of R ∂ , then using the notation of Figures 20 and 21, we have • Finite case. Point 1.
For everyw,v,f ; for every w, b, b such that w / ∈ {w , w r ∈ R ∂ }, with the notation of Figures 20 and 21 (left), we have

The dimer model on an isoradial graph G Q and the Z u -Dirac operator
Infinite case. In the paper [BdTR17a], we prove an explicit local expression for an inverse (K Q ) −1 ; as a byproduct we obtain a local formula for a Gibbs measure P Q dimer on (M(G Q ), F), involving the operators K Q and (K Q ) −1 ; we refer to the paper [BdTR17a] for the explicit formula for (K Q ) −1 and to Section 2.2.2 for the explicit formula of the Gibbs measure P Q dimer . Using Corollary 45 and Corollary 27, we provide an alternative direct way of finding the local formula for the inverse operator (K Q ) −1 of [BdTR17a], where the locality property is directly seen as inherited from that of the Z-massive and dual massive Green functions: we first express coefficients of (K Q ) −1 using the inverse Z u -Dirac operator for appropriate values of u, and then express the latter using the Z-massive and dual massive Green function of [BdTR17b]. Note that it is not immediate to see equality between the formulas of [BdTR17a] and Corollay 48; it probably requires to use elliptic trigonometric identities. The approach we propose here also extends to the finite case.
Notation for Corollary 48. Letw be a white vertex of G Q and b be a black one.
Consider v,f, e iβ i , and w, b, w, e iᾱ f , e iβ f as in Figure 20. The quadrangle of the vertex b corresponds to a rhombus v 1 , f 1 , v 2 , f 2 of the diamond graph G , where vertices are labeled so that the edge (v 1 , v 2 ) of G is parallel to the edge (b, w) of G Q , and f 1 is on the right of (v 1 , v 2 ), see Figure 22.w Figure 22: Notation for Corollary 48.
Corollary 48. For every white vertexw and every black vertex b of G Q , using the notation of Figure 22, we have Proof. We set u = β f in Corollary 45. Then, u α f = θ f , u β f = 0, and we have cn(u β f ) = 1, We now set u = β f in Corollary 27. Using that dn(−u) = dn(u), dn(K ± u) = k nd(u), dn(K) = k , sc(θ * f ) = (k ) −1 cs(θ f ), we obtain: Plugging this into the first equality of the corollary yields the second and concludes the proof.
Example 49. As an example of application we express the probability of single edges occurring in dimer configurations of G Q chosen with respect to the measure P Q dimer , as a function of single edge probabilities of the dimer model on G D with Gibbs measure P D,β dimer . We then use Example 33 evaluated at u = β to obtain explicit expressions using the function H; details of computations are given in Appendix B.2. Using the notation of Figure 16, we have We recover the results of the computations of Theorem 37 of [BdTR17a] after using addition formulas.
Finite case. An explicit expression for the Boltzmann measure P Q dimer as a function of the matrix K Q and its inverse (K Q ) −1 is given in Section 2.2. We now express the inverse Kasteleyn matrix (K Q ) −1 as a function of the Z u -Dirac operator for some choices of u. Note that we cannot proceed as in the infinite case using u = β f since then the matrix K ∂ (β f ) is not invertible.
1. For every white vertexw of G Q and every black vertex b / ∈ {b , b r ∈ R ∂ }, using the notation of Figures 20 and 21 (left), Then, Remark 51. Example 30 expresses coefficients of K(u f ) −1 , K(v f ) −1 using the Z-massive and dual massive Green function. Plugging this into Corollary 50 allows to write coefficients of (K Q ) −1 using the Green functions G m, Proof. Considerw,v,f and w, b, b such that w / ∈ {w , w r ∈ R ∂ } as in the statement of Corollary 46. Then, from Example 47, we have the following linear system of equations: Solving for (K Q ) −1 w,b gives, Now, using [Law89,2.4.8] and the fact that dn(2K − u) = dn(u), cn(2K − u) = − cn(u), we have for every u ∈ T(k), cn 2 (u) + cn 2 (K − u) = cn(2u)+dn(2u) 1+dn(2u) Putting everything together ends the proof of Point Suppose that edges of G F are oriented according to a Kasteleyn orientation, and denote by K F the corresponding Kasteleyn matrix, as defined in Section 2.2. Following Dubédat [Dub11], we partition vertices of G F as V F = A∪B, where A consists of vertices incident to four internal edges of G F , and B consists of those incident to either two internal and one external edges or to two internal edges (this possibility only occurs in the finite case). Vertices of type A, resp. B, will be denoted by a, resp. b, with or without sub/super-scripts. Up to a reordering of the rows and columns, the matrix K F can be written in block form as Recall that in the finite case, boundary quadrangles of G Q are degenerate and consist of edges in bijection with boundary edges of G. Boundary B-vertices of G F are defined to be B-vertices incident to two internal edges only; they are in natural bijection with black, resp. white, vertices of boundary quadrangles of G Q , and also with boundary white vertices of the double graph G D .
In [Dub11] Dubédat shows how, in the case where G is Z 2 , a Kasteleyn orientation on G F induces a Kasteleyn orientation on G Q ; this generalizes to the case where G is planar: using the notation of Figure 23 below, define then it is straightforward to check that the orientation so defined on G Q is Kasteleyn. Note that when G is finite, the case ε b,w is not present when b belongs to a boundary quadrangle. Wrapping up and using the notation of Figure 23, we have: for the Kasteleyn matrix K F , The first contribution of this section is Theorem 55 of Section 5.2 expressing the inverse Kasteleyn operator (K F ) −1 using the inverse bipartite Kasteleyn operator (K Q ) −1 ; proving that the contour Ising Boltzmann/Gibbs measures can be computed from the bipartite dimer model on G Q ; note that this result is not restricted to the Z-invariant case. The proof of Theorem 55 builds on matrix relations of [Dub11]; this is the subject of Section 5.1. In Section 5.3 we restrict to the Z-invariant case and obtain Corollary 59, one of the main results of this paper, expressing the inverse Kasteleyn operator (K F ) −1 using the inverse Z u -Dirac operator and also using the Z-massive and dual Green functions. This shows that the contour Ising Boltzmann/Gibbs measures can be computed using information from random walks only (with specific boundary conditions in the finite case). As written in the introduction to this paper, this has implications for other observables of the Ising model, as the spin-Ising observable of [CS12] or the fermionic spinor observable of [KC71]. Note that in the infinite case, this also gives an alternative direct way of finding the local formula for (K F ) −1 of [BdTR17a], explicitly relating it to the Green functions.
5.1 Relating Kasteleyn matrices of the Fisher graph G F and the bipartite graph G Q Dubédat [Dub11] establishes a matrix relation between the matrix K F and a block triangular matrix containingK Q as one of the diagonal blocks. Using this matrix relation, he proves that the squared dimer partition function of G F is equal, up to a constant, to the dimer partition function of G Q in the finite case, and that the characteristic polynomials of the two models are equal in the case of infinite Z 2 -periodic graphs. By adding defects to Ising coupling constants, this allows him to prove bozonisation identities.
We attribute the forthcoming Proposition 53, consisting of two matrix relations, to Dubédat. The first is the actual identity of [Dub11]; it is appropriate for comparing determinants of the matrices K F andK Q (related matrix relations can also be found in [CCK15]). The second proves an identity between K F and a block diagonal matrix containingK Q in both diagonal blocks; it is not present in the paper [Dub11] but does not require much more work; it is useful for comparing matrix inverses. For convenience of the reader we provide a proof because: we write weights in a different way, directly write the proof for all planar graphs (and not Z 2 ), handle the boundary conditions very carefully, and the second identity needs an additional argument.
In order to state Proposition 53, we need to introduce the following matrices, all of which are "square".
The matrix I W,A has rows indexed by white vertices of G Q and columns by A-vertices of G F . It is the identity matrix associated to the following bijection between W and A. Using the notation of Figure 23, a vertex w of G Q is incident to a unique external edge wb of G Q ; the edge wb is naturally "parallel" to a path fromb to b in the external cycle of the closest decoration of G F ; then the vertex a in bijection with w is the unique vertex of type A in the path fromb to b.
The matrix X has rows indexed by B-vertices of G F and columns by black vertices of G Q . It is block diagonal, with blocks of size 2 × 2 corresponding to edges of G * . For each such edge, the rows are indexed by b and b , the two corresponding adjacent B-vertices, and the columns are indexed by the two black vertices b, b of the quadrangle of G Q traversed by the edge bb , with b closest to b, see Figure 23. The non-zero coefficients of the row corresponding to the vertex b are, In the finite case, the matrix X also has size 1 blocks corresponding to boundary B-vertices of G F . For such a vertex b, let b be the closest black vertex of G Q . Then, the only non-zero coefficient of the row corresponding to b is: The matrix M has rows indexed by B-vertices and columns by A-vertices of G F . It is block diagonal, with blocks corresponding to decorations, each block having per size the number of B-vertices times the number of A-vertices of the decoration. The matrix M is the matrix K F B,A with some signs reversed. That is, for a B-vertex b, denote by a, a its two neighbors of type A so that in cclw order around the triangle we have a, a , b, see Figure 23. Then the non-zero coefficients of the row corresponding to the vertex b are, The matrix M has rows indexed by A-vertices and columns by B-vertices of G F . It is defined as, Indeed, it is block diagonal, with blocks corresponding to decorations; for each decoration, the block is a directed adjacency matrix of the bipartite graph consisting of the outer cycle of the decoration. The orientation on this cycle is Kasteleyn because it is on the whole graph and this cycle contains no vertex in its interior [Kas67]. As a consequence, the determinant of this block is equal to ± the number of dimer configurations of this cycle, that is ±2, and the block is thus invertible.
In the sequel, it would have been tempting to sometimes use the inverse of the matrix K F A,A , but this matrix is not always invertible. Indeed, it is block diagonal with blocks corresponding to decorations, and when the decoration is associated to a dual vertex of odd degree, then the corresponding block of K F A,A is not invertible.
[Dub11] The Kasteleyn matrix K F of the Fisher graph G F and the bipartite Kasteleyn matrixK Q of the graph G Q are related by the following identities: Proof. The first identity is an easy consequence of the second, so let us prove the second; unless specified, the arguments hold in the infinite and finite cases. We need to show the four identities below.
Note that even in the infinite case, they all make sense since matrices involved have finitely many non-zero coefficients per rows and columns.
Identity (52) is immediate by definition of M . We now prove (53) and (54) and then show that (55) follows.
Proof of (53). Let us show that K F A,B M = −K F A,A . Consider an A-vertex a of G F , and let a 1 , a 2 , b 1 , b 2 be its four neighbors in G F with the notation of Figure 24. Then the coefficient (K F A,B M ) a,a is a priori non zero when a ∈ {a, a 1 , a 2 }. Returning to the definition of K F and M , we have When a ∈ {a 1 , a 2 }, using moreover that the orientation around the triangles a, b 1 , a 1 and a, a 2 , b is Kasteleyn, we have thus ending the proof of (53).
where recall m b ,a 3 = −ε b ,a 3 and m b ,a 4 = ε b ,a 4 .
The coefficient (XK Q I W,A ) b,a of the RHS of (54) is non-zero for the same choices of a. Re-turning to the definition of X andK Q we have, Finite case. The proof is as in the infinite case as long as b is not a boundary B-vertex of G F , so let b be a boundary B-vertex and refer to Figure 25 (right) for notation. The coefficient Similarly to the infinite case computation, we have: The coefficient (XK Q I W,A ) b,a of the RHS is non-zero for the same choices of a. Returning to the definition of X andK Q (boundary case) we have: This ends the proof of (54).
Proof of (55). From Remark 52 the matrix K A,B is invertible, thus from (53) we have Plugging this into the LHS of (54) gives that it is equal to: Returning to the definition of M (or to (52)), we have that the LHS of(55) is Using that the matrix K F is skew-symmetric, we deduce that LHS (55) = −[LHS (54)] t . The same clearly holds for the RHS of the two equations; they are thus equivalent and we have proved (54).
Corollary 54 ( [Dub11]). Suppose that the graph G is finite. Then, Proof. By the first identity of Proposition 53, we have | det The determinant of X is computed by calculating that of its blocks. Let bb be an edge of G F corresponding to an edge e * of G * and let b, b be the black vertices of the quadrangle of G Q traversed by the edge bb . Then by definition, the corresponding block of X is: .
Its determinant is equal to 1 + |K F b,b | 2 = 1 + e −4Je , thus ending the proof of the corollary.

Relating inverse Kasteleyn matrices of G F and G Q
Consider the inverse Kasteleyn operators (K F ) −1 and (K Q ) −1 . When the graph G is infinite and Z 2 -periodic, these inverses denote the unique ones decreasing to 0 at infinity [KOS06,BdT10]. When the graph G is infinite and isoradial (not necessarily Z 2 -periodic) and the corresponding dimer weights on G Q and G F are Z-invariant, then (K Q ) −1 and (K F ) −1 are the operators decreasing to 0 at infinity with local expression given in [Ken02,BdT11,BdTR17a].
From Proposition 53 and proving additional relations (not present in the paper [Dub11]) we show, in Theorem 55 below, identities relating the inverse operator (K F ) −1 to the inverse operator (K Q ) −1 . Using Section 2.2.2, Theorem 55 allows to express the dimer Boltzmann measure (finite case) and the Gibbs measure (infinite case) of the Fisher graph G F , denoted P F dimer , using coefficients of the matrix K F and of the inverse bipartite Kasteleyn operator (K Q ) −1 , see also Example 57. To state Theorem 55, we need to define two additional matrices D B,A and κ.
The matrix D B,A has rows indexed by black vertices of G Q and columns by A-vertices of G F . It is a diagonal matrix associated to the following bijection between B and A. Using the notation of Figure 23, a vertex b of G Q belongs to a unique external edge bw of G Q ; then the vertex a in bijection with b is the vertex in bijection with w in the construction of the matrix I W,A . The corresponding diagonal coefficient d b,a is, The matrix κ has rows and columns indexed by A-vertices of G F . It is block diagonal with blocks corresponding to decorations, each block having per size the number of A-vertices of the decoration. Given two vertices a, a of a decoration of G F , we have, where n(a, a ) is the number of edges oriented cw in the cclw path going from a to a in the A-cycle of the decoration.
Notation for coefficients (K F ) −1 u,v of Theorem 55. Ifū =ā is an A-vertex, thenw is the white vertex of G Q corresponding toā in the bijection defining I W,A . Ifū =b is a B-vertex, thenā 1 andā 2 are the two A-vertices belonging to the same triangle, withb,ā 1 ,ā 2 in cclw order around the triangle. Note that this definition also holds ifb is a boundary vertex.
If v = b is a B-vertex, we let b be the closest black vertex of G Q . When moreover b is not a boundary vertex, we let b be the B-vertex such that bb defines an edge e * of G * ; we let b be the black vertex of G Q closest to b (b and b are the black vertices of the quadrangle of G Q traversed by the edge bb ); the coupling constant of the edge e, dual of e * , is denoted J f e , where "f" stands for "final". If v = a is an A-vertex, we let b and b be as defined in the matrix D B,A , see Figure 26. Theorem 55. As long as they are unique (which happens for sure in the finite and Z 2periodic cases), the inverse Kasteleyn operator (K F ) −1 can be expressed using the inverse bipartite Kasteleyn operator (K Q ) −1 as follows.
• Coefficients. We have four cases to consider and use the notation of Figure 26.
1. For everyā ∈ A and every b ∈ B such that, when the graph G F is moreover finite, b is not a boundary vertex: 2. When the graph G F is finite, for everyā ∈ A and every boundary vertex b of B, we have 3. For everyā, a ∈ A, 4. For everyb, b ∈ B, where (K F ) −1 a 1 ,b , (K F ) −1 a 2 ,b are given by (59). Remark 56.
• When proving the local formula for (K F ) −1 a,a [BdT11,BdTR17a] in the Z-invariant case, we obtained a formula of the form (61) -with the constant κā ,a -without explicitly relating it to a coefficient of (K Q ) −1 . It is quite remarkable that this formula holds in the full planar case (without assuming Z-invariance), in the finite and infinite cases.
• In the finite case, we do not need positivity of the coupling constants J. In particular if the coupling constants are all negative, Theorem 55 expresses probabilities of the dimer model on the non-bipartite graph G F with positive weights e −2Je > 1 on external edges, as a function of the inverse Kasteleyn operator of a "dimer model" on the bipartite graph G Q with negative weights (tanh(J e ) < 0) on quadrangle edges parallel to edges of G. Having a negative weight for an edge amounts to reversing its orientation. From a physics point of view, this amounts to adding defects or creating vortices at each face where this new orientation is not Kasteleyn. The physics paper [NO17] considers bipartite models with negative weights and somehow describes the non-Harnacity of the associated spectral curves.
• Coefficients of the inverse Kasteleyn operator (K F ) −1 not only allow to express the dimer Boltzmann/Gibbs measure P F dimer , but are also related to important observables of the Ising model. By [CCK15], the coefficient (K F ) −1 b,b is essentially the spin-observable of [CS12] when fixing one vertex to be on the boundary of the domain; by [Dub11] the coefficient (K F ) −1 a,a is the fermionic spinor correlator of [KC71] and by [NK85], it is related to the FK-Ising observable [Smi06, Smi10, CS12] when fixing one vertex on the boundary of the domain, taking appropriate boundary conditions, up to normalization.
εb ,ā 1 (K F ) −1 a 1 ,a + εb ,ā 2 tanh(2J e )(K F ) −1 a 2 ,a + εb ,b εb ,ā 3 cosh −1 (2J e )(K F ) −1 a 3 ,a = 0. For an edge bb of G F corresponding to an edge e * of G * , let b, b be the black vertices of the quadrangle of G Q traversed by the edge bb . Then the inverse of the corresponding block is, In the finite case, the matrix X also has a size 1, identity block for each boundary B-vertex b of G F and its closest black vertices b of G Q . This ends the proof of formulas (59), (61), (62) for coefficients and we now turn to the proof of (58).
The expressions for ( are a direct consequence of Proposition 53. This is not the case of (K F ) −1 A,A which requires proving additional identities. From Proposition 53 we know that, where in the second and third equalities we used skew-symmetry of K F , and in the third the definition of M . We thus need to prove that so let us prove (64). We will be using the notation of Figure 28 below. We need to introduce an additional matrix, the matrix D A,B , which has rows indexed by A-vertices and columns by B-vertices of G F . It is diagonal: to an A-vertex a corresponds the unique B-vertex b such that b comes before a in the cw ordering of the triangle containing a and b. The diagonal coefficient d a,b is: For example to a 2 of Figure 28 corresponds the vertex b, and the coefficient d Consider a B-vertex b of G F . Then, the coefficient (κ A,A K F A,B ) a,b of the LHS of (65) is a priori non-zero for all A-vertices a belonging to the same decoration as b. We have, Returning to the definition of the matrix Λ, as long as a = a 2 , we have that κ a,a 1 = κ a,a 2 ε a 2 ,a 1 implying, since the orientation around the triangle a 1 , a 2 , b is Kasteleyn and using the definition of D A,B . When a = a 2 , then κ a 2 ,a 1 = − 1 4 ε a 2 ,a 1 , κ a 2 ,a 2 = 1 4 ; thus using again the Kasteleyn orientation around the triangle and the definition of D A,B , thus ending the proof of (65). Plugging (65) into (64) leaves us with showing the equivalent Infinite case. Let b be a vertex of G Q , then the coefficient ( (66) is non-zero whenb ∈ {b, b }. Recalling the computation of X −1 given in (63), we have = −e −2J cosh −1 (2J) = −1 + tanh(2J), using (56).
and hence Equation (66) is proved in the infinite case.
Finite, boundary case. Consider a black vertex b of G Q . As long as b is not a boundary vertex of G Q , the argument is as in the infinite case, so we suppose that b is a boundary vertex. The coefficient (2X −1 K F B,B ) b,b of the LHS of (66) is zero for allb ∈ B. The coefficient (D B,A K F A,B + 2K Q I W,A D A,B ) b,b of the RHS of (66) is a priori non-zero whenb ∈ {b, b 1 }, and we have where in the computation forb = b we have used that the boundary coefficientK Q b,w 2 = ε b,a 2 . The computation forb = b 1 is as in the infinite case. This ends the proof of the finite, boundary case of (66) and the proof of Theorem 55.
Example 57. As an example we give the probability of single edges occurring in dimer configurations of G F chosen with respect to the Boltzmann measure P F dimer in the finite case, or the Gibbs measure P F dimer in the infinite case, as a function of edge probabilities of the dimer model on G Q . Details of computations are given in Appendix B.3. Using the notation of Figure 23, we have

In the Z-invariant case
We restrict to the case where the graph G is isoradial, finite or infinite, with Ising coupling constants J given by (6), dimer weights on G F by (7) and dimer weights on G Q by (8). The main result of this section is Corollary 59 relating the inverse Kasteleyn operator (K F ) −1 to the inverse Z u -Dirac operator using Theorem 55, Example 47, Corollaries 50 and 48, and to the Z-massive Green functions. This proves one of the main results of this paper, namely that the contour Ising Boltzmann/Gibbs measures can be computed from the inverse Z u -Dirac operator, and also from the Z-massive Green functions. As a byproduct, in the infinite this also gives a direct alternative way of proving the local formula of [BdTR17a] for (K F ) −1 , where the locality is seen as directly inherited from the Green function. Note that as for the local expression of (K Q ) −1 , it is not immediate to see equality with the expression of [BdTR17a].
We first relate the real and complex bipartite Kasteleyn matricesK Q and K Q of the graph G Q ; we use Appendix A.2. Define the following function q on pairs of vertices of G Q , inductively on edges. For every edge bw of G Q , let For every pair of vertices x, y of G Q , let q x,y = (x ,y )∈γx,y q x ,y , where γ x,y is an edge-path of G Q from x to y. Then, since the matricesK Q and K Q satisfy the alternating product condition around every face/inner face of G Q if the graph is infinite/finite [Kup98,Ken02], the function q is well defined.
Consider a fixed vertex x 0 of G Q , and define the diagonal matrices D x 0 ,B , D x 0 ,W on black, resp. white vertices, of G Q by: By [Kup98,Ken02], see also Appendix A.2, we havẽ thus implying the following lemma: Lemma 58 ( [Kup98,Ken02]). Consider the function q as defined above. Then, coefficients of the inverse of the matricesK Q and K Q are related by the following, for every white vertex w and every black vertex b of G Q , Note that the above matrix relation holds in the finite and infinite cases because coefficients of the diagonal matrices are finite and uniformly bounded away from 0.
Cases 2. and 3. of Theorem 55 directly relate coefficients of (K F ) −1 to (K Q ) −1 w,b . Using Lemma 58, Corollary 48 (infinite case), Corollary 50 and Remark 51 (finite case), these coefficients of (K F ) −1 are easily expressed using the inverse Z u -Dirac operator and the Zmassive Green functions.
Case 4. of Theorem 55 uses Case 1. so we are left with considering Case 1. The notation used are summarized in Figure 29 below.
Figure 29: Left: notation for initial vertices. Right: notation for final vertices.
• For everyā ∈ A and every b ∈ B such that, when the graph G F is moreover finite, b is not a boundary vertex: • Moreover, as a function of the inverse Z u f -Dirac operator, where u f = α f +β f 2 + K, we have: Finite case. Ifā corresponds to a vertexw c for some rhombus pair of R ∂ , we then use the notation of Figure 29 (2nd quadrant on the left): f,w are expressed using the Green functions G m,∂ (u) and G m, * in Example 30.
Proof. Let us prove the first point. We first compare (K Q ) −1 w,b and (K Q ) −1 w,b of Theorem 55. Using the notation of Figure 29 (right), we have that is because, omitting the subscript "f", We thus have, We are left with computing the terms involving the coupling constants J f e in the Z-invariant case. By definition, where in the third equality we used [Law89, chap.2, ex.14 (iii)] and in the fourth that sn(v + K) = cd(v). From this and Identity (48), we obtain, Putting together Equation (67), (68) and (69) ends the proof of the first point; let us now prove the second.
In the infinite case, by Example 47, we have thus proving the first line. The second line is obtained by using Example 30 to express K(u f ) −1 v,w and K(u f ) −1 f,w using the Z-massive Green functions G m and G m, * . In the finite case, by Example 47, we have thus concluding the proof.

Examples
In this section we specify some of our results to two cases of interest: the critical Z-invariant Ising model, and the full Z-invariant Ising model when the underlying isoradial graph G is Z 2 with the regular embedding.

Z-invariant critical case
The Z-invariant Ising model is critical when k = 0 (k = 1) [Li12,CD13,Lis14b]. In this case, the elliptic functions sn, cn are the trigonometric functions sin, cos and dn ≡ 1.
Returning to Section 3.1, the finite Z u -Dirac operator K ∂ (u) with boundary conditions arising from the Ising model is: Away from the boundary, we recover the Dirac operator of [Ken02].
Dimer model on the bipartite graph G Q , finite case. Corollary 50 expresses coefficients of (K Q ) −1 as a function of the inverse Z u -Dirac operator. When b does not belong to a boundary quadrangle (Case 1.), we have f,w , and coefficients of t(u) are given by (35) and (36). The mass of the Z-massive Laplacian is equal to 0 away from the bound-ary [BdTR17b]; then Example 30 expresses coefficients of K ∂ (u f ) −1 using the Z-Green function G 0,∂ (u f ) and dual Z-Green function G 0, * ; we have Dimer model on the Fisher graph G F , finite case. Corollary 59 (finite case) simply becomes: Partition function of the Ising model with + boundary conditions. Corollary 38, expressing the squared Z-invariant Ising partition function, holds for every u ∈ (T(k)) . When k = 0, a nice expression is obtained by setting u = iu and taking the limit u → −∞.
We essentially recover the main result of [dT16] proving that the squared critical Z-invariant Ising model partition function is equal, up to an explicit constant, to the partition function of spanning trees with specific boundary conditions. The difference is that we here consider + boundary conditions instead of free ones, and more importantly, we use the boundary trick of Chelkak and Smirnov [CS12] allowing us to remove all contributions from dual spanning trees, which we could not do in [dT16].
6.2 Full Z-invariant case when the isoradial graph G is Z 2 The goal of this example is to relate our results to the papers [Mes06,BDC12]. When the isoradial graph G = Z 2 in its regular embedding, all rhombus half-anglesθ e are equal to π 4 . Using that sn K 2 = 1 + k .
In the infinite case, by Theorem 55 and Corollary 48, the coefficient (K F ) −1 a,a of the inverse Kasteleyn operator (K F ) −1 is equal to, as long asā and a do not belong to the same decoration, A similar expression holds in the finite case, see Corollary 50 and Example 30. According to [Dub11,NK85], this coefficient is the fermionic spinor observable of [KC71] and, up to normalization, the FK-spin observable of [Smi10,CS12]. Now, the proof of [BDC12] for showing the occurrence of large deviation estimates of the massive random walk in the correlation length of the spin correlations [Mes06] consists in proving that, in the super-critical regime, spin correlations can be approximated by the FK-spin observable, and then using massive harmonicity of the latter to relate it to the massive random walk. Our explicit expression for |(K F ) −1 a,a | in the finite case gives a direct explanation of the occurrence of these large deviation estimates, and up to handling boundary terms, should give a rather direct proof valid in the whole super-critical Z-invariant case.

A Gauge equivalence revisited
We consider gauge equivalence of weighted adjacency matrices of digraphs and rephrase gauge equivalence of bipartite weighted adjacency matrices as defined in [Kup98,KOS06] in this context.

A.1 Definitions
In the whole of this section, we consider square matrices of size n×n; let M be such a matrix. We associate two graphs to M , a non-directed one G(M ) = (V(M ), E(M )) and a directed one D(M ) = (V(M ), A(M )), both having the same vertex set of cardinality n in bijection with rows/columns of M . An edge x j x is in E(M ) iff M x j ,x = 0 or M x ,x j = 0. A directed edge (or simply an edge) (x j , x ) is in A(M ) iff the coefficient M x j ,x = 0. The matrix M is a weighted adjacency matrix of the digraph D(M ). Whenever no confusion occurs, we remove the argument M of the graphs.
Let us recall a few definitions. A di-path of a digraph D = (V, A) is a sequence (x 0 , . . . , x n ) of vertices such that for every j ∈ {0, . . . , n − 1}, (x j , x j+1 ) is an edge of A. A simple di-path is a path with pairwise disjoint vertices. A di-cycle is a di-path such that the first and last vertices are the same. A simple di-cycle is a cycle whose only common vertices are the first and the last. Note that loops and length-two di-cycles are simple. A digraph is strongly connected if any two pairs of vertices are joined by a di-path. The above definitions are easily adapted in the case of non-directed graphs. Remark 61.
• Having the same associated digraph is equivalent to asking that M x,y = 0 ⇔ N x,y = 0.
• Since loops are simple di-cycles, if M and N are gauge equivalent, they have equal diagonal coefficients.
• Definition 60 holds if and only if the product condition holds for every di-cycle of D.
Lemma 62. Let M, N be two gauge equivalent matrices, then Proof. This is proved by writing the determinant as a sum over permutations, and doing the cyclic decomposition of permutations.
Let us suppose that M and N are gauge equivalent and that the associated digraph D is strongly connected. Define the function q ∈ C V×V on pairs of vertices of V taking values in C as follows. For vertices x, y such that (x, y) is an edge of D, set q x,y = N x,y M x,y .
For vertices x, y of D, since the digraph is strongly connected, there exists a di-path γ from x to y; set q x,y = e=(x ,y )∈γ q x ,y .
Note that if y = x, then γ is a di-cycle and we have q x,x = 1.
Remark 63. The function q is well defined, i.e., independent of the choice of path from x to y. If y = x, then q x,x = 1 independently of the choice of di-cycle from x to x. If y = x, consider two di-paths γ 1 , γ 2 from x to y. Since the digraph D is strongly connected, there exists a simple di-pathγ from y to x. Then, γ 1 (resp. γ 2 ) followed byγ is a di-cycle and by definition of gauge equivalence we have, e=(x ,y )∈γ 1 q x ,y e=(x ,y )∈γ q x ,y = 1 = e=(x ,y )∈γ 2 q x ,y e=(x ,y )∈γ q x ,y , implying that e=(x ,y )∈γ 1 q x ,y = e=(x ,y )∈γ 2 q x ,y , and the function q is thus well defined.
The following lemma proves that if the associated digraph is strongly connected, gauge equivalence amounts to having the two matrices related through a diagonal matrix. Proof. Suppose that M and N are gauge equivalent. Fix a vertex x 0 of D, and define D x 0 to be the diagonal matrix whose diagonal coefficient D x 0 x,x corresponding to the vertex x is q x 0 ,x . Since the digraph D is strongly connected, the matrix D x 0 is invertible. Let us prove that M = D x 0 N (D x 0 ) −1 . Non zero coefficients of M and N correspond to edges of D; let (x, y) be such an edge. Consider a di-path γ from x 0 to x, then γ = γ ∪ {(x, y)} is a di-path from x 0 to y. Using γ to compute q x 0 ,x and γ to compute q x 0 ,y , we deduce that, Remark 65. Consider two matrices M, N having the same associated graphs G and D, such that G is connected and such that for every undirected edge e of G, the two possible oriented edges are present in D; then D is strongly connected. If moreover G is planar and embedded in a planar way, then every simple di-cycle c of length ≥ 3 is the union of face di-cycles, where edges not in c are traversed in both directions. As a consequence, in this case proving gauge equivalence for M and N is equivalent to proving that,

A.2 The bipartite case
Consider a non-directed, finite, bipartite graph G = (V = W ∪ B, E), such that |W| = |B| = n, having at least one perfect matching. Note that G being bipartite, it cannot have loops; we furthermore suppose that it has no multiple edges, i.e., that it is simple. A bipartite, weighted adjacency matrix K associated to the graph G has rows indexed by vertices of B, columns by those of W, and non-zero coefficients correspond to edges of G.
Up to a reordering of the rows and columns, we can suppose that rows of K are indexed by b 1 , . . . , b n and columns by w 1 , . . . , w n .
Instead of seeing K as a bipartite adjacency matrix, we can interpret it as an adjacency matrix of the digraph D 0 . In this interpretation, rows and columns are indexed by vertices of V 0 and diagonal coefficients correspond to edges of the perfect matching, they now represent loops.
Consider the diagonal matrix D 0 K whose j-th diagonal coefficients is the coefficient K b j ,w j corresponding to the j-th edge of M 0 . Define the matrix K 0 to be, Note that the matrix K 0 has ones on the diagonal.
Definition 66. Let G be a finite, bipartite graph, and let K, L be two associated bipartite, weighted adjacency matrices. Fix a perfect matching M 0 of G, and let D 0 be the directed graph constructed from G and M 0 . Then, K and L are said to be gauge equivalent if the matrices K 0 and L 0 , seen as weighted adjacency matrices of the digraph D 0 , are gauge equivalent.
Rephrasing Lemma 62 in the context of bipartite graphs, we obtain Corollary 67. Let K, L be two gauge equivalent bipartite, weighted adjacency matrices, then We now rephrase Definition 66 in the more usual form [Kup98,KOS06]. Consider the bipartite graph G together with the reference perfect matching M 0 . An alternating cycle of G and M 0 is a simple cycle of G whose edges alternate between edges in M 0 and edges in E \ M 0 .
Lemma 68. The matrices K and L are gauge equivalent iff for every alternating cycle c of G and M 0 of length > 2, we have Proof. By definition, K and L are gauge equivalent if K 0 and L 0 seen as adjacency matrices of the digraph G 0 are gauge equivalent, see Definition 60. Length one di-cycles of D 0 are loops corresponding to diagonal coefficients of K 0 , L 0 . The latter are all equal to 1 by definition, and thus equal. Consider a simple di-cycle c of D 0 of length m ≥ 2. Up to a relabeling of the vertices, it can be denoted as c = {x 1 , . . . , x l = x 1 }. Then, ∀ j, x j corresponds to an edge b j w j of the perfect matching M 0 . By construction of D 0 , the cycle c is in correspondence with an alternating cycle {w 1 , b 1 , w 2 , . . . , w m , b m , w 1 } of length 2m of G. By definition of K 0 we have, m j=1 K 0 x j ,x j+1 = m j=1 K b j ,w j+1 m j=1 K b j ,w j , with cyclic notation for indices. A similar equality holds for the matrices L 0 and L, thus ending the proof.
Remark 69. Definition 66 is independent of the choice of M 0 . To prove this, use the fact that if M 1 is another reference perfect matching, then the superimposition of M 0 and M 1 consists of alternating cycles of length > 2 and doubled edges.
From now on we suppose that the bipartite graph G is finite, planar and simply connected, and we let K, L be two bipartite, weighted adjacency matrices of G.
Lemma 70 ( [Kup98]). If the alternating products of the matrices K and L are equal around every inner face-cycle of G, then K and L are gauge equivalent.
Proof. This is proved by induction on the number of faces contained in an alternating cycle of length > 2.
Suppose that K and L satisfy the alternating product condition around every inner face-cycle of G. Similarly to the directed case, define the function q ∈ C V×V as follows. For every edge bw of G, For every pair of vertices x, y of G, let γ be a path from x to y and set q x,y = (x ,y )∈γ q x ,y . The function q is well defined [Kup98] and Proof. The proof can be found in [Kup98]. For the purpose of Section 5.3, we make the diagonal matrices explicit assuming the alternating product condition is satisfied. Fix a vertex x 0 of G and set, where q is given above. Then, K = D x 0 ,B L D x 0 ,W .
B Computations of probabilities of single edges.
B.1 Dimers model on an infinite isoradial double graph G D We compute the probability of single edges occurring in dimer configurations of G D chosen with respect to the measure P D,u dimer . Notation are recalled in Figure 30 below; since no confusion occurs, we omit the subscripts f from α, β, θ. Returning to the definition of K Q , this immediately gives, P Q dimer (bw 1 ) = P D,β dimer (wv 2 ), P Q dimer (bw 2 ) = P D,β dimer (wf 2 ).

B.3 Dimers on G F
We now compute single edge probabilities for dimer configurations of G F chosen with respect to the Boltzmann measure P F dimer . where we used that P Q dimer (bw ) =K Q b,w (K Q ) −1 w ,b = ε b,a tanh(2J)(K Q ) −1 w ,b , = ε a ,a − 1 4 + P Q dimer (bw ) 2 tanh(2J) , since the orientation of the triangle is Kasteleyn.
We deduce that P F dimer (aa ) = K F a,a (K F ) −1 a ,a = 1 4 − P Q dimer (bw ) 2 tanh(2J) . Settingā = a and in Formula (59) gives, where, as before, we returned to the definition of K Q and of single edge probabilities of G Q = ε b,a 2 P Q dimer (bw ) + P Q dimer (b w ) + P Q dimer (bw ) tanh(2J) = ε b,a 2 1 − P Q dimer (b w ) + P Q dimer (bw ) tanh(2J) , using that the sum of probabilities is 1 around w for P Q dimer . We deduce that P F dimer (ba ) = K F b,a (K F ) −1 a ,b = 1 2 − P Q dimer (b w ) 2 + P Q dimer (bw ) 2 tanh(2J) .
As a consequence also, we have: P F dimer (ba ) + P F dimer (aa ) = 3 4 − P Q dimer (b w ) 2 . Using that the sum of probabilities around a is 1, we deduce that P F dimer (ab) = P Q dimer (bw) 2 + P Q dimer (bw ) 2 tanh(2J) . Using that the sum of probabilities around the vertex b is 1, we obtain By symmetry P Q dimer (bw) = 1 2 = P Q dimer (b w ), implying the expressions of Example 57.