Cédric Boutillier
email: cedric.boutillier@sorbonne-universite.fr

Béatrice De Tilière

Kilian Raschel
email: raschel@math.cnrs.fr

The Z-invariant Ising model via dimers

come

Introduction

The Z-invariant Ising model, fully developed by Baxter [START_REF] Baxter | Solvable eight-vertex model on an arbitrary planar lattice[END_REF][START_REF] Baxter | Free-fermion, checkerboard and Z-invariant lattice models in statistical mechanics[END_REF][START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF], takes its roots in the work of Onsager [START_REF] Onsager | Crystal statistics. I. A two-dimensional model with an order-disorder transition[END_REF][START_REF] Wannier | The statistical problem in cooperative phenomena[END_REF], see also [START_REF] Au-Yang | Q-dependent susceptibilities in ferromagnetic quasiperiodic Z-invariant ising models[END_REF][START_REF] Perk | Yang Baxter equations[END_REF][START_REF] Martìnez | Correlation functions for the Z-invariant Ising model[END_REF][START_REF] Martìnez | Multi-spin correlation functions for the Z-invariant Ising model[END_REF][START_REF] Costa-Santos | Geometrical aspects of the Z-invariant Ising model[END_REF] for further developments in the physics community. It is defined on a planar, embedded graph G = (V, E) satisfying a geometric constraint known as isoradiality, imposing that all faces are inscribable in a circle of radius 1. In this introduction, the graph G is assumed to be infinite and locally finite. The star-triangle move (see Figure 5) preserves isoradiality; it transforms a three-legged star of the graph 1 into a triangle face. The Ising model is said to be Z-invariant if, when decomposing the partition function according to the possible spins at vertices bounding the triangle/star, the contributions only change by an overall constant. This constraint imposes that the coupling constants J = (J e ) e∈E satisfy the Ising model Yang-Baxter equations. The solution to these equations is parametrized by angles naturally assigned to edges in the isoradial embedding of the graph G, and an elliptic parameter k, with k 2 ∈ (-∞, 1): ∀ e ∈ E, J e = J(θ e |k) = 1 2 log 1 + sn(θ e |k) cn(θ e |k) ,

where sn and cn are two of the twelve Jacobi trigonometric elliptic functions. More details and precise references are to be found in Section 2.2. When k = 0, the elliptic functions sn, cn degenerate to the usual trigonometric functions sin, cos and one recovers the critical Z-invariant Ising model, whose criticality is proved in [START_REF] Li | Critical temperature of periodic Ising models[END_REF][START_REF] Cimasoni | The critical temperature for the Ising model on planar doubly periodic graphs[END_REF][START_REF] Lis | Phase transition free regions in the Ising model via the Kac-Ward operator[END_REF]. Note that the coupling constants range from 0 to ∞ as k varies, thus covering the whole range of temperatures, see Lemma

A fruitful approach for studying the planar Ising model is to use Fisher's correspondence [START_REF] Fisher | On the dimer solution of planar Ising models[END_REF] relating it to the dimer model on a decorated version G F of the graph G, see for example the book [START_REF] Mccoy | The two-dimensional Ising model[END_REF]. The dimer model on the Fisher graph arising from the critical Z-invariant Ising model was studied by two of the present authors in [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: the periodic case[END_REF][START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF]. One of the main goals of this paper is to prove a generalization to the full Z-invariant Ising model of the latter results.

Furthermore, we answer questions arising when the parameter k varies. In the same spirit, we also solve the bipartite dimer model on the graph G Q associated to two independent Zinvariant Ising models [START_REF] Dubédat | Exact bosonization of the Ising model[END_REF][START_REF] Boutillier | Height representation of XOR-Ising loops via bipartite dimers[END_REF] and related to the XOR-Ising model [START_REF] Kadanoff | Correlation functions on the critical lines of the Baxter and Ashkin-Teller models[END_REF][START_REF] Wilson | XOR-Ising loops and the Gaussian free field[END_REF]. In order to explain the main features of our results, we now describe them in more details.

The Kasteleyn matrix/operator [START_REF] Kasteleyn | The statistics of dimers on a lattice: I. the number of dimer arrangements on a quadratic lattice[END_REF][START_REF] Temperley | Dimer problem in statistical mechanics-an exact result[END_REF] is the key object used to obtain explicit expressions for quantities of interest in the dimer model, as the partition function, the Boltzmann/Gibbs measures and the free energy. It is a weighted, oriented, adjacency matrix of the dimer graph. Our first main result proves an explicit, local expression for an inverse K -1 of the Kasteleyn operator K of the dimer model on the Fisher graph G F arising from the Z-invariant Ising model; it can loosely be stated as follows, see Theorem 11 for a more precise statement.

Theorem 1. Define the operator K -1 by its coefficients:

∀ x, y ∈ V F , K -1 x,y = ik 8π Γx,y f x (u + 2K)f y (u) e (x,y) (u)du + C x,y ,
where f and e, see [START_REF] Cimasoni | The critical temperature for the Ising model on planar doubly periodic graphs[END_REF] and [START_REF] Boutillier | Height representation of XOR-Ising loops via bipartite dimers[END_REF], respectively, are elliptic functions defined on the torus T(k), whose aspect ratio depends on k. The contour of integration Γ x,y is a simple closed curve winding once vertically around T(k), which intersects the horizontal axis away from the poles of the integrand; the constant C x,y is equal to ±1/4 when x and y are close, and 0 otherwise, see [START_REF] Dubédat | Exact bosonization of the Ising model[END_REF].

Then K -1 is an inverse of the Kasteleyn operator K on G F . When k = 0, it is the unique inverse with bounded coefficients.

Remark 2.

• The expression for K -1

x,y has the remarkable feature of being local. This property is inherited from the fact that the integrand, consisting of the function f and the massive discrete exponential function, is itself local : it is defined through a path joining two vertices corresponding to x and y in the isoradial graph G. This locality property is unexpected when computing inverse operators in general.

• As for the other local expressions proved for inverse operators [START_REF] Kenyon | The Laplacian and Dirac operators on critical planar graphs[END_REF][START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF][START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF], Theorem 11 has the following interesting features: there is no periodicity assumption on the isoradial graph G, the integrand has identified poles implying that explicit computations can be performed using the residue theorem (see Appendix B), asymptotics can be obtained via a saddle-point analysis (see Theorem 13).

• The most notable feature is that Theorem 11 is a generalization to the elliptic case of Theorem 1 of [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF]. Let us explain why it is not evident that such a generalization should exist. Thinking of Z-invariance from a probabilist's point of view suggests that there should exist local expressions for probabilities. The latter are computed using the Kasteleyn operator K and its inverse, suggesting that there should exist a local expression for the inverse operator K -1 , but giving no tools for finding it. Until our recent paper [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF], local expressions for inverse operators were only proved for critical models [START_REF] Kenyon | The Laplacian and Dirac operators on critical planar graphs[END_REF][START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF], leading to the belief that not only Z-invariance but also criticality played a role in the existence of the latter. Another difficulty was that some key tools were missing. We believed that if a local expression existed in the non-critical case, it should be an elliptic version of the one of the critical case, thus requiring an elliptic version of the discrete exponential function of [START_REF] Mercat | Exponentials form a basis of discrete holomorphic functions on a compact[END_REF], which was unavailable. This was our original motivation for the paper [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF] introducing the massive discrete exponential function and the Z-invariant massive Laplacian. The question of solving the dimer representation of the full Z-invariant Ising model turned out to be more intricate than expected, but our original intuition of proving an elliptic version of the critical results turns out to be correct.

• On the topic of locality of observables for critical Z-invariant models, let us also mention the paper [START_REF] Grimmett | Universality for bond percolation in two dimensions[END_REF] by Manolescu and Grimmett, recently extended to the random cluster model [START_REF] Duminil-Copin | Universality for the random-cluster model on isoradial graphs[END_REF]. Amongst other results, the authors prove the universality of typical critical exponents and Russo-Seymour-Welsh type estimates. The core of the proof consists in iterating star-triangle moves in order to relate different lattices. This is also the intuition behind locality in Z-invariant models: if these critical exponents were somehow related to inverse operators (which could maybe be true for the q = 2 case), then one would expect local expressions for these inverses.

In Theorem 19, using the approach of [START_REF] De Tilière | Quadri-tilings of the plane[END_REF], see also [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF], we prove an explicit, local expression for a Gibbs measure on dimer configurations of the Fisher graph, involving the operator K and the inverse K -1 of Theorem 1. This allows to explicitly compute probability of edges in polygon configurations of the low or high temperature expansion of the Ising model, see Equation [START_REF] Kennelly | The equivalence of triangles and three-pointed stars in conducting networks[END_REF].

Suppose now that the isoradial graph G is Z 2 -periodic, and let G 1 = G/Z 2 be the fundamental domain. Following an idea of [START_REF] Kenyon | The Laplacian and Dirac operators on critical planar graphs[END_REF] and using the explicit expression of Theorem 1, we prove an explicit formula for the free energy of the Z-invariant Ising model, see also Corollary 21. This expression is also local in the sense that it decomposes as a sum over edges of the fundamental domain G 1 . A similar expression is obtained by Baxter [START_REF] Baxter | Solvable eight-vertex model on an arbitrary planar lattice[END_REF][START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF], see Remark 24 for a comparison between the two.

Theorem 3. The free energy F k Ising of the Z-invariant Ising model is equal to: where sc = sn cn and the function H is defined in (66)-(67).

F k Ising = -|V 1 | log 2 2 -|V 1 | K 0 2H (
It turns out that the free energy of the Ising model is closely related to that of the Z-invariant spanning forests of [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF], see also Corollary 22.

Corollary 4. One has

F k Ising = -|V 1 | log 2 2 + 1 2 F k forest .
This extends to the full Z-invariant Ising model the relation proved in the critical case [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF] between the Ising model free energy and that of critical spanning trees of [START_REF] Kenyon | The Laplacian and Dirac operators on critical planar graphs[END_REF]. Moreover, in [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF] we prove a continuous (i.e., second order ) phase transition at k = 0 for Z-invariant spanning forests, by performing an expansion of the free energy around k = 0: at k = 0, the free energy is continuous, but its derivative has a logarithmic singularity. As a consequence of Corollary 4 we deduce that the Z-invariant Ising model has a second order phase transition at k = 0 as well. This result in itself is not surprising and other techniques, such as those of [START_REF] Duminil-Copin | Connection probabilities and RSW-type bounds for the two-dimensional FK Ising model[END_REF] and the fermionic observable [START_REF] Chelkak | Universality in the 2D Ising model and conformal invariance of fermionic observables[END_REF] could certainly be used in our setting too to derive this kind of result; but what is remarkable is that this phase transition is (up to a factor 1 2 ) exactly the same as that of Z-invariant spanning forests. More details are to be found in Section 4.3.

It is interesting to note that the Z-invariant Ising model satisfies a duality relation in the sense of Kramers and Wannier [START_REF] Kramers | Statistics of the two-dimensional ferromagnet[END_REF][START_REF] Kramers | Statistics of the two-dimensional ferromagnet[END_REF]: the high temperature expansion of a Z-invariant Ising model with elliptic parameter k on an isoradial graph G, and the low temperature expansion of a Z-invariant Ising model with dual elliptic parameter k * = i k √ 1-k 2 on the dual isoradial graph G * yield the same probability measure on polygon configurations of the graph G. The elliptic parameters k and k * can be interpreted as parametrizing dual temperatures, see Section 4.2 and also [START_REF] Chan | The Ising susceptibility scaling function[END_REF][START_REF] Mccoy | The anisotropic Ising correlations as elliptic integrals: duality and differential equations[END_REF].

The next result proves a self-duality property for the Ising model free energy, see also Corollary 30. This is a consequence of Corollary 4 and of Lemma 29, proving a self-duality property for the Z-invariant massive Laplacian.

Corollary 5. The free energy of the Z-invariant Ising model on the graph G satisfies the self-duality relation

F k Ising + |V 1 | 2 log k = F k * Ising + |V 1 | 2 log k * ,
where k = √ 1 -k 2 is the complementary elliptic modulus, and k * = 1/k .

The above result extends to all isoradial graphs a self-duality relation proved by Baxter [START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF] in the case of the triangular and honeycomb lattices. Note that this relation and the assumption of uniqueness of the critical point was the argument originally used to derive the critical temperature of the Ising model on the triangular and honeycomb lattices, see also Section 4.4.

In Section 5 we consider the dimer model on the graph G Q associated to two independent Zinvariant Ising models. This dimer model is directly related to the XOR-Ising model [START_REF] Dubédat | Exact bosonization of the Ising model[END_REF][START_REF] Boutillier | Height representation of XOR-Ising loops via bipartite dimers[END_REF]. Our main result is to prove an explicit, local expression for the inverse K -1 of the Kasteleyn operator associated to this dimer model. This is a generalization, in the specific case of the bipartite graph G Q , of the local expression obtained by Kenyon [START_REF] Kenyon | The Laplacian and Dirac operators on critical planar graphs[END_REF] for all "critical" bipartite dimer models.

Theorem 6. Define the operator K -1 by its coefficients:

∀ (b, w) ∈ B Q × W Q , K -1 b,w = 1 4iπ Γ b,w f (b,w) (u)du,
where f (b,w) is an elliptic function defined on the torus T(k), defined in Section 5.2. The contour Γ b,w is a simple closed curve winding once vertically around T(k), which intersects the horizontal axis away from the poles of the integrand.

Then K -1 is an inverse operator of K. For k = 0, it is the only inverse with bounded coefficients.

We also derive asymptotics and deduce an explicit, local expression for a Gibbs measure on dimer configurations of G Q , allowing to do explicit probability computations.

Outline of the paper.

• Section 2. Definition of the Ising model, of the two corresponding dimer models and of their Z-invariant versions. Definition of the Z-invariant massive Laplacian of [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF].

• Section 3. Study of the Z-invariant Ising model on G via the dimer model on the Fisher graph G F and the corresponding Kasteleyn operator K: definition of a oneparameter family of functions in the kernel of K, statement and proof of a local formula for an inverse K -1 , explicit computation of asymptotics, specificities when the graph G is periodic (connection with the massive Laplacian), and consequences for the dimer model on G F .

• Section 4. Behavior of the model as the parameter k varies: duality in the sense of Kramers and Wannier [START_REF] Kramers | Statistics of the two-dimensional ferromagnet[END_REF][START_REF] Kramers | Statistics of the two-dimensional ferromagnet[END_REF], phase transition in k, self-duality property, connection with the modular group.

• Section 5. Study of the double Z-invariant Ising model on G via the dimer model on the bipartite graph G Q and the Kasteleyn operator K: one-parameter family of functions in the kernel of K, statement and proof of a local formula for an inverse K -1 , explicit computation of asymptotics and consequences for the dimer model on G Q .

Acknowledgments: We acknowledge support from the Agence Nationale de la Recherche (projet MAC2: ANR-10-BLAN-0123) and from the Région Centre-Val de Loire (projet MADACA).

We are grateful to the referee for his/her many insightful comments.

2 The models in question

The Ising model via dimers

In this section we define the Ising model and two of its dimer representations. The first is Fisher's correspondence [START_REF] Fisher | On the dimer solution of planar Ising models[END_REF] providing a mapping between the high or low temperature expansion of the Ising model on a graph G and the dimer model on a non-bipartite graph G F . The second is a mapping between two independent Ising models on G and the dimer model on a bipartite graph G Q [START_REF] Dubédat | Exact bosonization of the Ising model[END_REF][START_REF] Boutillier | Height representation of XOR-Ising loops via bipartite dimers[END_REF].

The Ising model

Consider a finite, planar graph G = (V, E) together with positive edge-weights J = (J e ) e∈E .

The Ising model on G with coupling constants J is defined as follows. A spin configuration σ of G is a function on vertices of G with values in {-1, 1}. The probability of occurrence of a spin configuration σ is given by the Ising Boltzmann measure, denoted P Ising :

P Ising (σ) = 1 Z Ising (G, J) exp e=xy∈E J e σ x σ y ,
where Z Ising (G, J) is the normalizing constant known as the Ising partition function.

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 high temperature expansion [START_REF] Kramers | Statistics of the two-dimensional ferromagnet[END_REF][START_REF] Kramers | Statistics of the two-dimensional ferromagnet[END_REF] of the Ising model partition function gives the following identity:

Z Ising (G, J) = 2 |V| e∈E cosh J e P∈P(G) e∈P tanh J e .

The dimer model

Consider a finite, planar graph G = (V, E) together with positive edge-weights ν = (ν e ) e∈E .

A dimer configuration M of G, also known as a perfect matching, is a subset of edges of G such that every vertex is incident to exactly one edge of M. Let M(G) denote the set of dimer configurations of the graph G. The probability of occurrence of a dimer configuration M is given by the dimer Boltzmann measure, denoted P dimer :

P dimer (M) = e∈E ν e Z dimer (G, ν) ,
where Z dimer (G, ν) is the normalizing constant, known as the dimer partition function.

Dimer representation of a single Ising model: Fisher's correspondence

Fisher's correspondence [START_REF] Fisher | On the dimer solution of planar Ising models[END_REF][START_REF] Dubédat | Exact bosonization of the Ising model[END_REF] gives a mapping between polygon configurations of a graph G and dimer configurations of a decorated version of the graph, denoted G F and called the Fisher graph. For the purpose of this paper it suffices to consider graphs with no boundary.

The decorated graph

G F = (V F , E F ) is constructed from G as follows.
Every vertex of G of degree d is replaced by a decoration containing 2d vertices: a triangle is attached to every edge incident to this vertex and these triangles are glued together in a circular way, see Figure 1.

The correspondence goes as follows. To a polygon configuration P of G one assigns 2 |V| dimer configurations of G F : edges present (resp. absent) in P are present (resp. absent) in G F ; then there are exactly two ways to fill each decoration of G F so as to have a dimer configuration, see Figure 1. From the correspondence, we know that:

G G F
Z Ising (G, J) = e∈E cosh J e Z dimer (G F , ν). (1)
Note that the above is Dubédat's version of Fisher's correspondence [START_REF] Dubédat | Exact bosonization of the Ising model[END_REF]. It is more convenient than the one used in [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: the periodic case[END_REF][START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF] because it allows to consider polygon configurations rather than complementary ones, and the Fisher graph has less vertices, thus reducing the number of cases to handle.

Dimer representation of the double Ising model

Based on results of physicists [START_REF] Kadanoff | Some critical properties of the eight-vertex model[END_REF][START_REF] Wu | Ising model with four-spin interactions[END_REF][START_REF] Fan | General lattice model of phase transitions[END_REF][START_REF] Wu | Staggered ice-rule vertex model -The Pfaffian solution[END_REF], Dubédat [START_REF] Dubédat | Exact bosonization of the Ising model[END_REF] provides a mapping between two independent Ising models, one living on the primal graph G, the other on the dual graph G * , to the dimer model on a bipartite graph G Q . Based on results of [START_REF] Nienhuis | Critical behavior of two-dimensional spin models and charge asymmetry in the Coulomb gas[END_REF][START_REF] Wu | Ising model with four-spin interactions[END_REF], two of the authors of the present paper exhibit an alternative mapping between two independent Ising models living on the same graph G (embedded on a surface of genus g) to the bipartite dimer model on

G Q [9].
Since the above mentioned mappings cannot be described shortly, we refer to the original papers and only define the bipartite graph G Q and the corresponding dimer weights. Note that dimer probabilities on the graph G Q can be interpreted as probabilities of the low temperature expansion of the XOR-Ising model [START_REF] Boutillier | Height representation of XOR-Ising loops via bipartite dimers[END_REF], also known as the polarization of the Ising model [START_REF] Kadanoff | Correlation functions on the critical lines of the Baxter and Ashkin-Teller models[END_REF][START_REF] Wilson | XOR-Ising loops and the Gaussian free field[END_REF] obtained by taking the product of the spins of the two independent Ising models.

We only consider the case where the graph G is planar and infinite. The bipartite graph

G Q = (V Q , E Q ) is obtained from G as follows.
Every edge e of G is replaced by a "rectangle", and the "rectangles" are joined in a circular way. The additional edges of the cycles are referred to as external edges. Note that in each "rectangle", two edges are "parallel" to an edge of the graph G and two are "parallel" to the dual edge of G * , see Figure 2.

Let ν = (ν e ) e∈E Q be the dimer weight function corresponding to two independent Ising models with coupling constants J. Then ν is equal to [START_REF] Dubédat | Exact bosonization of the Ising model[END_REF][START_REF] Boutillier | Height representation of XOR-Ising loops via bipartite dimers[END_REF]]

ν e =            tanh(2J e )
if e belongs to a "rectangle" and is parallel to an edge e of G, cosh(2J e ) -1 if e belongs to a "rectangle" and is parallel to the dual of an edge e of G, 1 if e is an external edge, 0 otherwise.

Z-invariant Ising model, dimer models and massive Laplacian

Although already present in the work of Kenelly [START_REF] Kennelly | The equivalence of triangles and three-pointed stars in conducting networks[END_REF], Onsager [START_REF] Onsager | Crystal statistics. I. A two-dimensional model with an order-disorder transition[END_REF] and Wannier [START_REF] Wannier | The statistical problem in cooperative phenomena[END_REF], the notion of Z-invariance has been fully developed by Baxter in the context of the integrable Figure 2: A piece of a graph G (plain grey lines) and its dual graph G * (dotted grey lines), and the corresponding bipartite graph G Q (plain black lines).

8-vertex model [START_REF] Baxter | Solvable eight-vertex model on an arbitrary planar lattice[END_REF], and then applied to the Ising model and self-dual Potts model [START_REF] Baxter | Free-fermion, checkerboard and Z-invariant lattice models in statistical mechanics[END_REF]; see also [START_REF] Perk | Yang Baxter equations[END_REF][START_REF] Au-Yang | Q-dependent susceptibilities in ferromagnetic quasiperiodic Z-invariant ising models[END_REF][START_REF] Kenyon | An introduction to the dimer model[END_REF]. Z-invariance imposes a strong locality constraint which leads to the parameters of the model satisfying a set of equations known as the Yang-Baxter equations. From the point of view of physicists it implies that transfer matrices commute, and from the point of view of probabilists it suggests that there should exist local expressions for probabilities, but it provides no tool for finding such expressions if they exist.

In Section 2.2.1 we define isoradial graphs, the associated diamond graph and star-triangle moves, all being key elements of Z-invariance. Then in Section 2.2.2 we introduce the Zinvariant Ising model [START_REF] Baxter | Solvable eight-vertex model on an arbitrary planar lattice[END_REF][START_REF] Baxter | Free-fermion, checkerboard and Z-invariant lattice models in statistical mechanics[END_REF][START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF], followed by the corresponding versions for the dimer models on G F and G Q . Finally in Section 2.2.5 we define the Z-invariant massive Laplacian and the corresponding model of spanning forests [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF].

Isoradial graphs, diamond graphs and star-triangle moves

Isoradial graphs, whose name comes from the paper [START_REF] Kenyon | The Laplacian and Dirac operators on critical planar graphs[END_REF], see also [START_REF] Duffin | Potential theory on a rhombic lattice[END_REF][START_REF] Mercat | Discrete Riemann surfaces and the Ising model[END_REF], are defined as follows. An infinite planar graph G = (V, E) is isoradial, if it can be embedded in the plane in such a way that all internal faces are inscribable in a circle, with all circles having the same radius, and such that all circumcenters are in the interior of the faces, see Figure 3 (left). This definition is easily adapted when G is finite or embedded in the torus.

From now on, we fix an embedding of the graph, take the common radius to be 1, and also denote by G the embedded graph. An isoradial embedding of the dual graph G * , with radius 1, is obtained by taking as dual vertices the circumcenters of the corresponding faces.

The diamond graph, denoted G , is constructed from an isoradial graph G and its dual G * . Vertices of G are those of G and those of G * . A dual vertex of G * is joined to all primal vertices on the boundary of the corresponding face, see Figure 3 (right). Since edges of the diamond graph G are radii of circles, they all have length 1, and can be assigned a direction ±e iα . Note that faces of G are side-length 1 rhombi. Using the diamond graph, angles can naturally be assigned to edges of the graph G as follows.

Every edge e of G is the diagonal of exactly one rhombus of G , and we let θ e be the halfangle at the vertex it has in common with e, see Figure 4. We have θ e ∈ (0, π 2 ), because circumcircles are assumed to be in the interior of the faces. From now on, we ask more and suppose that there exists ε > 0 such that θ e ∈ (ε, π 2 -ε). We further assign two rhombus vectors to the edge e, denoted by e iαe and e iβ e , see Figure 4. A train-track of G is a bi-infinite chain of edge-adjacent rhombi of G which does not turn: on entering a face, it exits along the opposite edge [START_REF] Kenyon | Rhombic embeddings of planar quad-graphs[END_REF]. Each rhombus in a train-track T has an edge parallel to a fixed unit vector ±e iα T , known as the direction of the train-track.

Train-tracks are also known as rapidity lines or simply lines in the field of integrable systems, see for example [START_REF] Baxter | Solvable eight-vertex model on an arbitrary planar lattice[END_REF].

The star-triangle move, also known as the Y -∆ transformation, underlies Z-invariance [START_REF] Baxter | Solvable eight-vertex model on an arbitrary planar lattice[END_REF][START_REF] Baxter | Free-fermion, checkerboard and Z-invariant lattice models in statistical mechanics[END_REF]. It is defined as follows: if G has a vertex of degree 3, that is a star Y , it can be replaced by a triangle ∆ by removing the vertex and connecting its three neighbors. The graph obtained in this way is still isoradial: its diamond graph is obtained by performing a cubic flip in G , that is by flipping the three rhombi of the corresponding hexagon, see Figure 5. This operation is involutive.

Figure 5: Star-triangle move on an isoradial graph G (plain lines) and cubic flip on the underlying diamond graph G (dotted lines).

Z-invariant Ising model

The Ising model defined on a graph G is said to be Z-invariant, if when decomposing the partition function according to the possible spin configurations at the three vertices of a star/triangle, it only changes by a constant when performing the Y -∆ move, this constant being independent of the choice of spins at the three vertices.

This strong constraint yields a set of equations known as the Ising model Yang-Baxter equations, see (6.4.8) of [START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF] and also [START_REF] Onsager | Crystal statistics. I. A two-dimensional model with an order-disorder transition[END_REF][START_REF] Wannier | The statistical problem in cooperative phenomena[END_REF]. The solution to these equations can be parametrized by the elliptic modulus k, where k is a complex number such that k2 ∈ (-∞, 1), and the rapidity parameters, see Equation (7.8.4) and page 478 of [START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF]. In this context it is thus natural to suppose that the graph G is isoradial. Extending the form of the coupling constants to the whole of G we obtain that they are given by, for every edge e of G,

J e = J(θ e |k) = 1 2 log 1 + sn(θ e |k) cn(θ e |k) , or equivalently sinh(2J(θ e |k)) = sc(θ e |k) 1 , ( 2
)
where k is the elliptic modulus,

θ e = θ e 2K π , K = K(k) = π 2 0 1 √ 1-k 2 sin 2 τ
dτ is the complete elliptic integral of the first kind, cn(•|k), sn(•|k) and sc(•|k) = sn(•|k) cn(•|k) are three of the twelve Jacobi trigonometric elliptic functions. More on their definition can be found in the books [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF]Chapter 16] and [START_REF] Lawden | Elliptic functions and applications[END_REF]; a short introduction is also given in the paper [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]Section 2.2]. Identities that are useful for this paper can be found in Appendix A.

For a given isoradial graph G, we thus have a one-parameter family of coupling constants (J) k , indexed by the elliptic modulus k, with k 2 ∈ (-∞, 1). For every edge e, the coupling constant J(θ e |k) is analytic in k 2 and increases from 0 to ∞ as k 2 increases from -∞ to 1, see Lemma 26; the elliptic modulus k thus parametrizes the whole range of temperatures. When k = 0, elliptic functions degenerate to trigonometric functions, and we have:

J(θ e |0) = 1 2 log 1 + sin θ e cos θ e .
The Ising model is critical at k = 0, see [START_REF] Li | Critical temperature of periodic Ising models[END_REF][START_REF] Cimasoni | The critical temperature for the Ising model on planar doubly periodic graphs[END_REF][START_REF] Lis | Phase transition free regions in the Ising model via the Kac-Ward operator[END_REF]. More on this subject is to be found in Section 4.

Corresponding dimer model on the

When k = 0 we have dn = 1 and sc = tan, which corresponds to the critical case.

Corresponding dimer model on the bipartite graph G Q

In a similar way, we compute the dimer weight function ν of the graph G Q corresponding to two independent Z-invariant Ising model. We have cosh(2J e ) = 1 2

1 + sn θ e cn θ e + cn θ e 1 + sn θ e = nc θ e , (4)
tanh(2J e ) = sinh(2J e ) cosh(2J e ) = sn θ e .

As a consequence of Section 2.1.4, the dimer weight function ν on the bipartite graph G Q is

ν e =           
sn θ e if e is parallel to an edge e of G, cn θ e if e is parallel to the dual edge of an edge e of G, 1 if e is an external edge, 0 otherwise.

(5)

The Z-invariant massive Laplacian

We will be using results on the Z-invariant massive Laplacian introduced in [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]. Let us recall its definition and the key facts required for this paper.

Following [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]Equation (1)], the massive Laplacian operator ∆ m : C V → C V is defined as follows. Let x be a vertex of G of degree n; denote by e 1 , . . . , e n edges incident to x and by θ 1 , . . . , θ n the corresponding rhombus half-angles, then

(∆ m f )(x) = n j=1 ρ(θ i |k)[f (x) -f (y)] + m 2 (x|k)f (x), (6)
where the conductances ρ and (squared) masses (m 2 ) are defined by

ρ e = ρ(θ e |k) = sc(θ e |k), (7)
(m 2 )(x) = m 2 (x|k) = n j=1 (A(θ j |k) -sc(θ j |k)), (8) with A
(u|k) = 1 k Dc(u|k) + E -K K u ,
where Dc(u|k) = u 0 dc 2 (v|k)dv, and E = E(k) is the complete elliptic integral of the second kind.

We also need the definition of the discrete k-massive exponential function or simply massive exponential function, denoted e (•,•) (•), of [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]Section 3.3]. It is a function from V × V × C to C. Consider a pair of vertices x, y of G and an edge-path x = x 1 , . . . , x n = y of the diamond graph G from x to y; let e iα j be the vector corresponding to the edge x j x j+1 . Then e (x,y) (•) is defined inductively along the edges of the path:

∀ u ∈ C, e (x j ,x j+1 ) (u) = i √ k sc u -α j 2 , e (x,y) (u) = n-1 j=1 e (x j ,x j+1 ) (u), (9)
where α j = α j 2K π . These functions are in the kernel of the massive Laplacian [START_REF] Bostan | The Ising model: from elliptic curves to modular forms and Calabi-Yau equations[END_REF], see [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]Proposition 11].

The massive Green function, denoted G m , is the inverse of the massive Laplacian operator [START_REF] Bostan | The Ising model: from elliptic curves to modular forms and Calabi-Yau equations[END_REF]. The following local formula is proved in [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]Theorem 12]:

G m (x, y) = k 4iπ Γx,y e (x,y) (u)du, (10)
where k = √ 1 -k 2 is the complementary elliptic modulus, Γ x,y is a vertical contour on the torus T(k) := C/(4KZ + 4iK Z), whose direction is given by the angle of the ray R -→ xy.

The massive Laplacian is the operator underlying the model of spanning forests, the latter being defined as follows. A spanning forest of G is a subgraph spanning all vertices of the graph, such that every connected component is a rooted tree. Denote by F(G) the set of spanning forests of G and for a rooted tree T, denote its root by x T . The spanning forest Boltzmann measure, denoted P forest , is defined by:

∀ F ∈ F(G), P forest (F) = T∈F m 2 (x T |k) e∈T ρ(θ e |k) Z forest (G, ρ, m) ,
where Z forest (G, ρ, m) is the spanning forest partition function. In [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]Theorem 41] we prove that this model is Z-invariant (thus explaining the name Z-invariant massive Laplacian). By Kirchhoff's matrix-tree theorem we have Z forest (G, ρ, m) = det(∆ m ).

3 Z-invariant Ising model via dimers on the Fisher graph G F From now on, we consider a fixed elliptic modulus k 2 ∈ (-∞, 1), so that we will remove the dependence in k from the notation.

In the whole of this section, we let G be an infinite isoradial graph and G F be the corresponding Fisher graph. We suppose that edges of G F are assigned the weight function ν of (3) arising from the Z-invariant Ising model.

We give a full description of the dimer model on the Fisher graph G F with explicit expressions having the remarkable property of being local. This extends to the Z-invariant non-critical case the results of [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: the periodic case[END_REF][START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF] obtained in the Z-invariant critical case, corresponding to k = 0. One should keep in mind that when k = 0, the "torus" T(0) is in fact an infinite cylinder with two points at infinity, and that "elliptic" functions are trigonometric series.

Prior to giving a more detailed outline, we introduce the main object involved in explicit expressions for the dimer model, namely, the Kasteleyn matrix/operator [START_REF] Kasteleyn | The statistics of dimers on a lattice: I. the number of dimer arrangements on a quadratic lattice[END_REF][START_REF] Temperley | Dimer problem in statistical mechanics-an exact result[END_REF].

Kasteleyn operator on the Fisher graph

An orientation of the edges of G F is said to be admissible if all cycles bounding faces of the graph are clockwise odd, meaning that, when following such a cycle clockwise, there is an odd number of co-oriented edges. By Kasteleyn [START_REF] Kasteleyn | Graph theory and crystal physics[END_REF], such an orientation always exists.

Suppose that edges of G F are assigned an admissible orientation, then the Kasteleyn matrix K is the corresponding weighted, oriented, adjacency matrix of G F . It has rows and columns indexed by vertices of G F and coefficients given by, for every x, y ∈ V F ,

K x,y = sgn(x, y)ν xy ,
where ν is the dimer weight function (3) and sgn(x, y) = 1 if x ∼ y and x → y, -1 if x ∼ y and y → x.

Note that K can be seen as an operator acting on C V F :

∀ f ∈ C V F , ∀ x ∈ V F , (Kf ) x = y∈V F K x,y f y .
Outline. Section 3 is structured as follows. In Section 3.2 we introduce a one-parameter family of functions in the kernel of the Kasteleyn operator K; this key result allows us to prove one of the main theorems of this paper: a local formula for an inverse K -1 of the operator K, see Theorem 11 of Section 3.3. Then in Section 3.4 we derive asymptotics of this inverse. In Section 3.5 we handle the case where the graph G is periodic. Finally in Section 3.6 we derive results for the dimer model on G F : we prove a local expression for the dimer Gibbs measure, see Theorem 19, and a local formula for the dimer and Ising free energies, see Theorem 20 and Corollary 21; we then show that up to an additive constant the Ising model free energy is equal to 1 2 of the spanning forest free energy, see Corollary 22.

Notation. Throughout this section, we use the following notation. A vertex x of G F belongs to a decoration corresponding to a unique vertex x of G. Vertices of G F corresponding to a vertex x of G are labeled as follows. Let d(x) be the degree of the vertex x in G, then the decoration consists of d(x) triangles, labeled from 1 to d(x) in counterclockwise order. For the j-th triangle, we let v j (x) be the vertex incident to an edge of G, and w j (x), w j+1 (x) be the two adjacent vertices in counterclockwise order, see Figure 6.

There is a natural way of assigning rhombus unit-vectors of G to vertices of G F : for every vertex x of G and every k ∈ {1, . . . , d(x)}, let us associate the rhombus vector e iα j (x) to w j (x), and the rhombus vectors e iα j (x) , e iα j+1 (x) to v j (x), see Figure 6; we let θ j (x) be the half-angle at the vertex x of the rhombus defined by e iα j (x) and e iα j+1 (x) , with

θ j (x) ∈ (0, π 2 ). G G F G x w j v j w j+1 v j+1
x x e iα j e iα j+1

Functions in the kernel of the Kasteleyn operator K

The definition of the one-parameter family of functions in the kernel of the Kasteleyn operator K requires two ingredients: the function f of Definition 3.1 and the massive discrete exponential function of [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF].

The function f uses the angles (α j (x)) assigned to vertices of G F , the latter being a priori defined in R/2πZ. For the function f to be well defined, we actually need them to be defined in R/4πZ, which is equivalent to a coherent choice for the determination of the square root of e iα j (x) . This construction is done iteratively, relying on our choice of Kasteleyn orientation.

Fix a vertex x 0 of G and set the value of α 1 (x 0 ) to some value, say 0. In the following, we use the index j (resp. ) to refer to vertices of G F belonging to a decoration x (resp. y) of G; with this convention, we omit the arguments x and y from the notation. For vertices in a decoration of a vertex x of G, define

α j+1 = α j + 2θ j if w j → w j+1 , α j + 2θ j + 2π if w j+1 → w j . ( 11
)
Given a directed path γ, let co(γ) be the number of co-oriented edges. Here is the rule defining angles in the decoration corresponding to a vertex y of G, neighbor of the vertex x. Let j and be such that v j is incident to v , as in Figure 7. Consider the length-three directed path w j , v j , v , w from w j to w . Then

α = α j -π if co(w j , v j , v , w ) is odd, α j + π if co(w j , v j , v , w ) is even. ( 12
) v j w j v w x y e iα j e iα
f : V F × C → C is defined by      f(w j , u) := f w j (u) = nc( u-α j 2 ), f(v j , u) := f v j (u) = K v j ,w j f w j (u) + K w j+1 ,v j f w j+1 (u) = K v j ,w j nc( u-α j 2 ) + K w j+1 ,v j nc( u-α j+1 2
).

(

) 13
Definition 3.2. The function g :

V F × V F × C → C is defined by g(x, y, u) := g x,y (u) = f x (u + 2K)f y (u) e (x,y) (u), (14)
where e (•,•) (•) is the massive exponential function of [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF], whose definition is recalled in Section 2.2.5.

Remark 8. The function g is meromorphic and biperiodic:

g (x,y) (u + 4K) = g (x,y) (u + 4iK ) = g (x,y) (u),
so that we restrict the domain of definition to T(k) := C/(4KZ + 4iK Z). Note however that taken separately, f x (• + 2K) and f y (•) are not periodic on T(k): only their product is.

The function g can also be seen as a one-parameter family of matrices (g(u)) u∈T(k) , where for every u ∈ T(k), g(u) has rows and columns indexed by vertices of G F , and g(u) x,y := g (x,y) (u).

We have the following key proposition.

Proposition 9. For every u ∈ T(k), Kg(u) = g(u)K = 0.

Proof. Note that since K is skew-symmetric, and that up to a sign, the functions x → g (x,z) (u) and x → g (z,x) (u + 2K) are equal:

g (z,x) (u + 2K) = f z (u + 4K)f x (u + 2K) e (z,x) (u + 2K) = -f x (u + 2K)f z (u) e (x,z) (u) = -g (x,z) (u),
it is enough to check the first equality, i.e., Kg(u) = 0.

Let us fix z. We need to check that for every vertex x of G F ,

d i=1 K x,x i g (z,x i ) (u) = 0,
where x 1 , . . . , x d are the d (equal to three or four) neighbors of x in G F . We distinguish two cases depending on whether the vertex x is of type w or v.

• If x = w j (x) for some j, then x has four neighbors: 6. Since all these vertices belong to the same decoration, the part f z (u + 2K) e (z,x) (u) is common to all the terms g z,x i (u). One is left with proving the following identity:

w j-1 (x) = w j-1 , w j+1 (x) = w j+1 , v j-1 (x) = v j-1 and v j (x) = v j , see Figure
(Kf)

w j = K w j ,w j-1 f w j-1 + K w j ,w j+1 f w j+1 + K w j ,v j-1 f v j-1 + K w j ,v j f v j = 0.
Using the second line in Equation ( 13) to express f v j and f v j-1 in terms of f w 's, one gets for the left-hand side of the previous equation:

(K w j ,w j-1 + K w j ,v j-1 K v j-1 ,w j-1 )f w j-1 + (K w j ,v j K v j ,w j + K w j ,v j-1 K w j ,v j-1 )f w j + (K w j ,w j+1 + K w j ,v j K w j+1 ,v j )f w j+1 .
The coefficient in front of f w j ,

K w j ,v j K v j ,w j + K w j ,v j-1 K w j ,v j-1 = -(K w j ,v j ) 2 + (K w j ,v j-1 ) 2 = -1 + 1,
is trivially equal to zero. Moreover, because of the condition on the orientation of the triangles in the Kasteleyn orientation, we have:

K w j ,w j-1 + K w j ,v j-1 K v j-1 ,w j-1 = K w j ,w j+1 + K w j ,v j K w j+1 ,v j = 0. ( 15
)
Indeed, to check this, it is enough to look at the case where the edges of a triangle are all oriented clockwise, and notice that the quantity is invariant if we simultaneously change the orientation of any pair of edges of the triangle, which is a transitive operation on all the (six) clockwise odd orientations of a triangle. So (Kf) w j is identically zero.

• If x = v j (x) = v j for some j, then x has three neighbors:

w j (x) = w j , w j+1 (x) = w j+1 and v (y) = v . Factoring out f z (u + 2K) e (z,x) (u), it is sufficient to prove that K v j ,w j f w j (u) + K v j ,w j+1 f w j+1 (u) + K v j ,v f v (u) e (x,y) (u) = 0.
Note that under inversion of the orientation of all edges around any of the vertices w j , w j+1 and v j , all the signs of three terms either stay the same, or change at the same time. To fix ideas, we can thus suppose that the edges of the triangles w j , w j+1 , v j and w , w +1 , v are all oriented clockwise, and that the edge between v j and v is oriented from v to v j , as in Figure 7. Returning to the definition of the angles mod 4π, see [START_REF] Chan | The Ising susceptibility scaling function[END_REF] and ( 12), and simplifying notation, we obtain

α = α j (x), β = α j+1 (x) = α + 2θ, α = α (y) = α -2K, β = α +1 (y) = α + 2θ = β -2K.
We have:

K v j ,w j f w j (u) + K v j ,w j+1 f w j+1 (u) = f w j (u) -f w j+1 (u) = cn( u-β 2 ) -cn( u-α 2 ) cn( u-α 2 ) cn( u-β 2 )
.

On the other hand, (13) entails that

f v (u) = f w (u) + f w +1 (u) = nc( u-α 2 ) + nc( u-β 2 ) = - 1 k sd( u-α 2 ) + sd( u-β 2 ) sd( u-α 2 ) sd( u-β 2 )
.

This has to be multiplied by K v j ,v = -sn θ 1+cn θ and by the exponential function e (x,y) (u), so that:

K v j ,v e (x,y) (u)f v (u) = sn θ 1 + cn θ 1 k sd( u-α 2 ) + sd( u-β 2 ) sd( u-α 2 ) sd( u-β 2 ) (-k ) sc( u-α 2 ) sc( u-β 2 ) = - sn θ 1 + cn θ sn( u-α 2 ) dn( u-β 2 ) + sn( u-β 2 ) dn( u-α 2 ) cn( u-α 2 ) cn( u-β 2 )
.

Proving that

K v j ,w j f w j (u) + K v j ,w j+1 f w j+1 (u) + K v j ,v e (x,y) (u)f v (u) = 0 ( 16
)
amounts to showing that

(1 + cn θ) cn( u-β 2 ) -cn( u-α 2 ) = sn θ sn( u-α 2 ) dn( u-β 2 ) + sn( u-β 2 ) dn( u-α 2 ) .
However, the addition formula (see Exercice 32 (v) in [START_REF] Lawden | Elliptic functions and applications[END_REF]Chapter 2] and also the similar relation ( 55)) reads:

cn(u + v) cn u = cn v -sn(u + v) sn u dn v. Evaluated at u = u-α 2 , v = -u-β
2 and u + v = θ (and exchanging the role of α and β for the second equation), we obtain

cn( u-β 2 ) cn θ = cn( u-α 2 ) + sn( u-β 2 ) dn( u-α 2 ) sn θ, cn( u-α 2 ) cn θ = cn( u-β 2 ) -sn( u-α 2 ) dn( u-β 2 ) sn θ.
Taking the difference of these two equations yields the result.

Remark 10. If µ is a measure on T(k) and if we define h x,y = g x,y (u)dµ(u), then by linearity of the integral, one also has Kh = hK = 0. A particular case, which will be important for what follows (see also [START_REF] Kenyon | The Laplacian and Dirac operators on critical planar graphs[END_REF][START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF]), is the case when µ is the integration along a contour on T(k).

Local expression for the inverse of the Kasteleyn operator K

We now state Theorem 11, proving an explicit, local formula for coefficients of the inverse K -1 of the Kasteleyn operator K. This formula is constructed from the function g of Definition 3.2.

Theorem 11. Consider the dimer model on the Fisher graph G F arising from the Z-invariant Ising model on the isoradial graph G, and let K be the corresponding Kasteleyn operator. Define the operator K -1 by its coefficients:

∀ x, y ∈ V F , K -1 x,y = ik 8π Γx,y g x,y (u)du + C x,y (17)
= ik 8π Γx,y f x (u + 2K)f y (u) e (x,y) (u)du + C x,y , (18)
where the contour of integration Γ x,y is a simple closed curve winding once vertically around the torus T(k) (along which the second coordinate globally increases), which intersects the horizontal axis in the angular sector (interval) s x,y of length larger than or equal to 2K (see Section 3.3.2), and the constant C x,y is given by

C x,y = - 1 4 •      1 if x = y = w j (x), (-1) n(x,y) if x = w j (x) and y = w (x) for some j = , 0 otherwise, (19)
where n(x, y) is the number of edges oriented clockwise in the counterclockwise arc from x to y in the inner decoration.

Then K -1 is an inverse of the Kasteleyn operator K on G F .

When k = 0, it is the unique inverse with bounded coefficients.

Alternatively, the coefficients of the inverse of the Kasteleyn operator admit the expression

K -1 x,y = ik 8π Cx,y f x (u + 2K)f y (u) e (x,y) (u)H(u)du + C x,y , (20)
where the function H is defined in (66)-(67), C x,y is a trivial contour oriented counterclockwise on the torus, not crossing Γ x,y and containing in its interior all the poles of g (x,y) and the pole of H, and C x,y is defined in [START_REF] Dubédat | Exact bosonization of the Ising model[END_REF].

Before we go on with the proof of this theorem, let us make a few comments about the formula of the inverse Kasteleyn matrix:

• As soon as x and y are not in the same decoration, or one of them is of type v, then the constant C x,y is zero, and the formula for K -1

x,y as a contour integral has the same flavour as the Green function of the Z-invariant massive Laplacian introduced in [10, Theorem 12].

• The constant C x,x is here to ensure that K -1

x,x is 0 if x is of type w (the integral is 0 when x is of type v as we shall see later).

• As one can expect, the full formula is skew-symmetric in x and y.

• To obtain the alternative expression (20) from [START_REF] De Tilière | Critical Ising model and spanning trees partition functions[END_REF], one can make use of a meromorphic multivalued function with a horizontal period of 1, like the function H defined in (66)-(67), originally introduced in [10] for k 2 ∈ (0, 1). Following this way, one may rewrite the integral as an integral over a contour bounding a disk, allowing one to perform explicit computation with Cauchy's residue theorem. One can add to H any elliptic function on T(k) without changing the result of the integral, given that C x,y encloses all the poles of the new integrand.

• Adding to the columns of K -1 functions in the kernel of K yield other inverses, with different behaviour at infinity. Such a function in the kernel is obtained by integrating g (x,y) (•) along a horizontal contour in T(k), see Remark 10. As a consequence, if we replace in [START_REF] De Tilière | Critical Ising model and spanning trees partition functions[END_REF] the contour Γ x,y by a contour winding a times vertically and b times horizontally, with a and b coprimes, and divide the integral by a, then we get a new inverse for the Kasteleyn operator K, which has an alternative expression as a trivial contour integral involving integer linear combinations of functions H and V , as defined in Appendix A.2.

• When k = 0, the "torus" T(k) is in fact a cylinder, with two points at infinity. The contour Γ x,y has infinite length. The function g x,y (u) decays sufficiently fast at infinity to ensure convergence of the integral. By performing the change of variable λ = -e iu in the integrals [START_REF] De Tilière | Critical Ising model and spanning trees partition functions[END_REF] or [START_REF] Duffin | Potential theory on a rhombic lattice[END_REF], one gets the adaptation to this variant of the Fisher graph of the formula for the inverse Kasteleyn operator in [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF], as an integral along a ray from 0 to ∞, or as an integral over a closed contour with a log.

We now turn to the proof of Theorem 11. We show that the operator K -1 with those coefficients satisfy KK -1 = Id and K -1 K = Id. These identities, understood as products of infinite matrices, make sense since K has a finite number of non zero coefficients on each row and column. Moreover, by skew-symmetry, it is enough to check the first one. When k = 0, it turns out that these coefficients for K -1 go to zero exponentially fast, see Theorem 13. This property together with K -1 K = Id imply injectivity of K on the space of bounded functions on vertices of G F , which in turn implies uniquess of an inverse with bounded coefficients.

The general idea for proving KK -1 = Id follows [START_REF] Kenyon | The Laplacian and Dirac operators on critical planar graphs[END_REF], but it is complicated by the fact that the Fisher graph G F itself is not isoradial. In this respect, the proof follows more closely that of Theorem 5 of [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF] with two main differences: we work with a different Fisher graph G F and more importantly we handle the elliptic case, making it a non-trivial extension. Section 3.3.1 corresponds to Sections 6.3.1 and 6.3.2 of [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF]. It consists in the delicate issue of encoding the poles of the integrand f x (u + 2K)f y (u) e (x,y) (u); for this question there are no additional difficulties so that we have made it as short as possible and refer to the paper [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF] for more details and figures. Section 3.3.2 consists in obtaining a sector s x,y on the horizontal axis of the torus T(k) from the encoding of the poles; this is then used to define the contour of integration Γ x,y . It corresponds to Section 6.3.3 of [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF] but requires adaptations to pass to the elliptic case. Section 3.3.3 is a non-trivial adaptation of Section 6.4 of [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF], handling a different Fisher graph G F and more importantly handling the elliptic case.

Preliminaries: encoding the poles of the integrand

Let G be an infinite isoradial graph and let G be the corresponding diamond graph. In order to encode poles of the integrand of K -1

x,y , we need the notion of minimal path which relies on the notion of train-tracks, see Section 2.2.1 for definition. A train-track is said to separate two vertices x, y of G if every path connecting x and y crosses this train-track. A path from x to y in G is said to be minimal if all its edges cross train-tracks that separate x from y, and each such train-track is crossed exactly once. A minimal path from x to y is in fact a geodesic for the graph metric on G . Since G is connected, it always exists. In general, there are several minimal paths between two vertices, but they all consist of the same steps taken in a different order.

For every pair of vertices x, y of G F , we now define an edge-path γ x,y of G encoding the poles of the integrand of K -1

x,y . Consider a minimal path from x to y and let e iα be one of the steps of the path, then the corresponding pole of the exponential function is α + 2K. Since e iα +2K = e iα +iπ = -e iα , this pole is encoded in the reverse step. As a consequence, poles of the exponential function are encoded in the steps {-e iα } of a minimal path from y to x.

We now have to add the poles of the functions f x (u + 2K) and f y (u). The difficulty lies in the fact that some of them might be canceled by factors in the numerator of the exponential function. By definition, the function f x (u + 2K) has either one or two poles {α j }, encoded in the edge(s) {e iα j } of the diamond graph G ; let T x = {T j

x } be the corresponding traintrack(s). Similarly, the pole(s) of f y (u) are at {α j + 2K} and are encoded in the edge(s)

{-e iα j } of G , and T y = {T j y } are the corresponding train-track(s). Let us start from a minimal path γ x,y from y to x. For every j, do the following procedure: if T j y separates y from x, then the pole α j + 2K is canceled by the exponential and we leave γ x,y unchanged. If not, this pole remains, and we extend γ x,y by adding the step -e iα j at the beginning of γ x,y . The path γ x,y obtained is still a path of G , denote by y the new starting point, at distance at most 2 from y.

When dealing with a pole of f x (u + 2K), one needs to be careful since, even when the corresponding train-track separates y from x, the exponential function might not cancel the pole, if it has already canceled the same pole of f y (u); this happens when T x and T y have a common train-track. The procedure to extend γ x,y runs as follows: for each j, if T j

x separates y and x and is not a train track of T y , then the pole α j is canceled by the exponential function, and we leave γ x,y unchanged. If not, this pole remains, and we extend γ x,y by attaching the step e iα j at the end of γ x,y . The path obtained in this way is still a path of G , starting from y. Denote by x its ending point, which is at distance at most 2 from x.

Obtaining a sector s x,y from γ x,y

Let x, y be two vertices of G F and let γ x,y be the path encoding poles of the integrand of K -1

x,y constructed above. Denote by {e iτ j } the steps of the path, seen as vectors in the unit disk; the corresponding poles of the integrand are {τ j }. Using these poles, we now define an interval/sector s x,y in the horizontal axis R/4KZ of the torus T(k). Given the sector s x,y , the contour of integration Γ x,y of K -1

x,y is then defined to be a simple closed curve winding once around the torus vertically, i.e., in the direction i, along which the second coordinate globally increases, and which intersects the horizontal axis in s x,y , see Figure 8 below.

General case. This case contains all but the three mentioned below. We know by Lemmas 17 and 18 of [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF] that there exists a sector in the unit circle, of size greater than or equal to π, containing none of the steps {e iτ j }. Equivalently, there exists a sector in the horizontal axis R/4KZ of the torus T(k), of size larger than or equal to 2K, containing none of the poles {τ j }. We let s x,y be this sector, it is represented in Figure 8. 2iK Here are the three cases which do not fit in the general situation.

Γ x,y C x,y s x,y Γ x,y
Case 1. The path γ x,y consists of two steps that are opposite. Then, the two poles separate the real axis of T(k) in two sectors of size exactly 2K, leaving an ambiguity. This can only happen when x = y = w j (x) and the two poles are {α j , α j + 2K}. In this case, the standard convention is that2 s w j ,w j = (α j , α j + 2K), and by Lemma 44 we have

ik 8π Γw j ,w j f w j (u + 2K)f w j (u)du = 1 4 . (21)
We choose the value of the constant C w j ,w j to compensate exactly the value of this integral so as to have K -1 w j ,w j = 0, that is C w j ,w j = -1 4 , and we recover the first line of the definition of C x,y of Equation [START_REF] Dubédat | Exact bosonization of the Ising model[END_REF].

Remark 12. It will be useful for the proof to consider also the non-standard convention with the complementary sector, defining a contour Γ w j ,w j . Returning to the definition of the function f w j (u), we have that the integral over Γ w j ,w j is equal to minus the one on the contour Γ w j ,w j , so that in order to have K -1 w j ,w j = 0, we set C w j ,w j = -C w j ,w j .

The two other cases correspond to situations when |T x T y | = 2. The corresponding path γ x,y does not enter the framework of Lemmas 17 and 18 of [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF]. They occur when x and y are equal or are neighbors in G F , and both of type 'v'. In other words, one has x = y = v j (x), or (x, y) = (v j (x), v (y)), with x ∼ y in G, j and being such that v j (x) ∼ v (y) in G F .

α j α j + 2K
Figure 9: Standard convention for the definition of the sector s w j ,w j .

Case 2. Suppose first that x = y = v j (x). Then the poles are {α j , α j + 2K, α j+1 , α j+1 + 2K} (the exponential function is equal to 1 and cancels no pole). If we take s v j ,v j = (α j , α j+1 ), then the integral is zero by symmetry. Indeed, the change of variable u → α j+1 + α j -u leaves the contour invariant (up to homotopy) and f v j (u + 2K) is changed into its opposite, whereas f v j (u) is invariant. Note that taking s v j ,v j = (α j + 2K, α j+1 + 2K) also gives a zero integral, because it is related to the previous one by the change of variable u → u + 2K. These two choices of sectors will be useful in the proof of Theorem 11, see Figure 10 (center).

Case 3. Suppose now that (x, y) = (v j (x), v (y)). Then f v j (x) and f v (y) induce twice the same poles {α j , α j+1 }. The exponential adds the poles {α j + 2K, α j+1 + 2K} and the numerator cancels one pair of {α j , α j+1 }, implying that there remains the poles {α j , α j +2K, α j+1 , α j+1 + 2K}. We set the convention given in Figure 10 (right).

α j + 2K α j+1 + 2K α j+1 α j α j + 2K α j+1 + 2K α j+1 α j α j + 2K α j+1 + 2K α j+1 α j v j e iα j+1 e iα j -e iα j+1 -e iα j v x y sv j ,v j sv j ,v
Figure 10: Definition of the sectors s v j ,v j (center) and s v j ,v (right).

Proof of the local formula for K -1 of Theorem 11

We need to prove that ∀ x, y ∈ V F , (KK -1 ) x,y = δ x,y .

We use the following notation. The vertex y is y = w j (y) or v j (y), for some vertex y of G and some j ∈ {1, . . . , d(y)}. If y = w j (y), it has four neighbors w j-1 (y), v j-1 (y), v j (y), w j+1 (y); if y = v j (y), it has three neighbors w j (y), w j+1 (y), v (y ), see Figure 11. We denote by x i the neighbors of x, with i ranging from 1 to 3 or 4.

y y w j v j v j-1 w j-1 w j+1 v w w +1
Figure 11: Notation for the cases where the general argument for proving (KK -1 ) x,y = δ x,y does not work and which have to be treated separately.

As long as the computation of (KK -1 ) x,y = i K x,x i K -1 x i ,y only involves terms K -1 x i ,y for which the constant C x i ,y is 0, and the sector s x i ,y defining the contour Γ x i ,y does not use any special convention, that is when • x is not in the same decoration as y, if y = w j (y),

• x / ∈ {w j (y), w j+1 (y), v (y )} {w (y ), w +1 (y ), v j (y)}, if y = v j (y), then by the argument of [START_REF] Kenyon | The Laplacian and Dirac operators on critical planar graphs[END_REF][START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF], all contours of integration Γ x i ,y can be deformed into a common contour Γ, crossing the horizontal axis in the nonempty intersection of the sectors i s x i ,y , so that by Proposition 9 (see also Remark 10) we have:

i K x,x i Γx i ,y f x i (u + 2K)f y (u) e (x,y) (u)du = Γ i K x,x i f x i (u + 2K)f y (u) e (x,y) (u)du = 0.
Let us check the remaining cases separately.

Suppose that y = w j (y). The degree of the vertex y is d(y) and indices below should be thought of as being modulo d. We have to handle all cases where the vertex x belongs to the decoration y, whether it is of type 'v' or 'w'.

• We first compute (KK -1 ) x,w j when x = v r (y) for some r ∈ {1, . . . , d(y)}. The vertex x has three neighbors w r (y), w r+1 (y) and a vertex of type 'v' in a neighboring decoration. We now omit the argument y from the notation.

When r ∈ {j -d + 1, . . . , j -2}, we are in the general case of the definition of Γ x i ,y ; when r = j -1, we choose the standard convention of Case 1, that is Γ w j ,w j and C w j ,w j ; when r = j, we choose the equivalent, non-standard convention of Case 1, that is Γ w j ,w j and C w j ,w j .

With these choices, the three sectors appearing in the expressions of K -1 x i ,w j have non-empty intersection, so that the contours Γ x i ,y in the three integrals can be continuously deformed into the same contour Γ, and thus the combination of the integral parts gives zero.

Since vertices of type 'v' have no constant contribution C x i ,y , we are left with proving that ∀ r ∈ {j -d + 1, . . . , j -1}, K vr,wr C wr,w j + K vr,w r+1 C w r+1 ,w j = 0, [START_REF] Duminil-Copin | Universality for the random-cluster model on isoradial graphs[END_REF] and that

K v j ,w j C w j ,w j + K v j ,w j-d+1 C w j-d+1 ,w j = 0. ( 23
)
Multiplying each of the equations of ( 22) by K vr,wr , and using that -K vr,wr K vr,w r+1 = K wr,w r+1 by the clockwise odd condition on triangles, we have that the first set of equations is equivalent to, for all r ∈ {j -d + 1, . . . , j -1}, C wr,w j = K wr,w r+1 C w r+1 ,w j , which in turn holds if and only if

C wr,w j = j-1 m=r K wm,w m+1 C w j ,w j = (-1) n(wr,w j ) C w j ,w j . (24)
Recalling that C w j ,w j = -1 4 and returning to the second line of the definition of C x,y , we see that this is indeed the case, whence ( 22) is proved.

We are left with proving that Equation ( 23) is satisfied. Doing the same steps as above, and using that C w j-d+1 ,w j = (-1) n(w j-d+1 ,w j ) C w j ,w j , this is equivalent to proving that C w j ,w j = C w j ,w j (-1) n(w j-d+1 ,w j ) K w j ,w j-d+1 .

Observing that (-1) n(w j-d+1 ,w j ) K w j ,w j-d+1 = -1 because of the clockwise odd condition on the inner circle of decorations, and recalling that C w j ,w j = -C w j ,w j (see Remark 12), we deduce that this equation is indeed true, thus ending the proof when x = v r (y).

• We now compute (KK -1 ) x,w j when x = w r (y) for some r ∈ {1, . . . , d(y)}. The vertex x = w r (y) has four neighbors: w r-1 (y), w r+1 (y), v r-1 (y) and v r (y), and we now omit the argument y.

Let us first handle the integral part. When r = j + 1, we are either in the general case of the definition of Γ x i ,y or in Case 1, and we choose the standard definition. When r = j + 1, we choose the non-standard definition of Case 1, that is Γ w j ,w j and C w j ,w j . With these choices, as long as r = j, the four sectors have non-empty intersection, so that the combination of the integral parts is equal to zero.

When r = j, then the four sectors enter the framework of the general case and are not compatible, see Figure 12.

A vertical contour Γ passing between α j and α j+1 is contained in the three sectors s w j-1 ,w j , s v j-1 ,w j , s w j+1 ,w j . If the fourth integral was taken along this contour, the combination of the

α j + 2K α j + 2K α j + 2K α j + 2K α j α j α j-1 α j-1 α j+1 α j+1 w j w j+1 w j-1 v j v j-1
sw j-1 ,w j sv j-1 ,w j sv j ,w j sw j+1 ,w j four would be zero. By adding and subtracting the integral for the pair (v j , w j ) along Γ , we have that the contribution of the integral part of (KK -1 ) w j ,w j is equal to

ik 8π K w j ,v j Γv j ,w j - Γ f v j (u + 2K)f w j (u)du. ( 25
)
The contour Γ v j ,w j -Γ is the (negatively oriented) boundary of a cylinder in the torus, which contains only one pole of the integrand, at u = α j . The function f w j (u) has no pole in the cylinder, and only the term involving f w j of the function

f v j (u + 2K) = K v j ,w j f w j (u + 2K) + K w j+1 ,v j f w j+1 (u+2K)
has a pole at u = α j . As a consequence, the contribution of the integral part is equal to

- ik 8π Γv j ,w j - Γ f w j (u+2K)f w j (u)du = - ik 8π Γ w j ,w j - Γw j ,w j f w j (u+2K)f w j (u)du = 1 2 ,
by continuously deforming the contours to those of Case 1, and using Equation ( 21) and Remark 12.

We now handle the constant part of (KK -1 ) wr,w j , keeping in mind that vertices of type 'v' have no constant contribution. As long as r + 1 / ∈ {j + 1, j + 2}, we have by Equation ( 24)

C wr,w j = K w r-1 ,wr K wr,w r+1 C w r+1,j , so that K wr,w r-1 C w r-1 ,w j + K wr,w r+1 C w r+1 ,w j = (K wr,w r-1 K w r-1 ,wr K wr,w r+1 + K wr,w r+1 )C w r+1 ,w j = 0.

When r = j + 1, recalling that we have chosen the non-standard definition from Case 1, factoring K w j ,w j+1 , using Equation ( 24) to write C w j+2 ,w j = C w j-d+2 ,w j , and finally remembering that C w j ,w j = -C w j ,w j , we have

K w j+1 ,w j C w j ,w j + K w j+1 ,w j+2 C w j+2 ,w j = K w j ,w j+1 1 + K w j ,w j+1 K w j+1 ,w j+2 (-1) n(w j-d+2 ,w j ) C w j ,w j ,
which is equal to 0 by the Kasteleyn orientation condition on inner cycles of decorations.

When r = j, using a similar argument, we obtain

K w j ,w j-1 C w j-1 ,w j + K w j ,w j+1 C w j+1 ,w j = -2C w j ,w j = 1 2 .
Wrapping up, we have proved that (KK -1 ) wr,w j is equal to 0 when r = j, and to 1 2 + 1 2 = 1 when r = j.

Suppose that y = v j (y). Note that since y is of type 'v', we always have C x i ,y = 0. We have to handle the cases where x ∈ {w j (y), w j+1 (y), v (y )} {w (y ), w +1 (y ), v j (y)}, and need to check whether the sectors defining the contours Γ x i ,y in the integral part of K x i ,y have non-empty intersections.

There are three values of x where one of the neighbors of x is y = v j (y): namely when x ∈ {w j (y), w j+1 (y), v (y )}. We now omit the arguments y, y from the notation. In these three cases, the sectors s w j ,v j and s w j+1 ,v j are compatible and intersect, either in the arc from α j to α j+1 (2 first cases), or from α j + 2K to α j+1 + 2K (last case). In all these situations, using the two possible definitions of Case 2 to write K -1 v j ,v j = 0 as the integral with a contour in that common sector, then by the general argument, we get that (KK -1 ) x,v j = 0.

We now need to check the remaining three cases where the combination uses K -1 v ,v j , corresponding to the situation where x ∈ {w (y ), w +1 (y ), v j (y)}.

In the two first situations, using the general Case and Case 3, we see that the sectors are compatible, and we can conclude with the general argument that (KK -1 ) x,v j = 0. Suppose now that x = y = v j . Its three neighbors are w j , w j+1 and v and the corresponding sectors are not compatible, see Figure 13.

α j + 2K α j + 2K α j α j α j+1 + 2K α j+1 + 2K α j+1 + 2K α j+1 α j+1 v j w j w j+1 v j-1 v sw j ,v j sw j+1 ,v j sv ,v j x x
Figure 13: Sectors s x i ,y when x = y = v j .

The two sectors for w j and w j+1 are compatible and intersect in the arc from α j to α j+1 , whereas according to the convention of Case 3, the one for v is the arc from α j + 2K to α j+1 +2K. A vertical contour Γ passing between α j and α j+1 is contained in the three sectors s w j ,v j , s w j+1 ,v j and s v ,v j . If the third integral was taken along this contour, the combination of the three would be 0. By adding and subtracting the integral of the pair (v , v j ) along Γ , and using Proposition 9 to write

K v j ,v f v (u + 2K) e y ,y (u) = -K v j ,w j f w j (u + 2K) + K v j ,w j f w j+1 (u + 2K) , we obtain (KK -1 ) v j ,v j = - ik 8π Γv ,v j - Γ K v j ,w j f w j (u + 2K) + K v j ,w j+1 f w j+1 (u + 2K) f v j (u)du.
By a change of variable u → u + 2K, the integral of the first term in the sum is

ik 8π Γ - Γv ,v j K v j ,w j f w j (u + 2K)f v j (u)du = - ik 8π K v j ,w j Γ +2K - Γ v ,v j +2K f v j (u + 2K)f w j (u)du,
which is exactly the same integral as the one computed in [START_REF] Grimmett | Universality for bond percolation in two dimensions[END_REF]. Indeed, Γ + 2K (resp. Γ v ,v j + 2K) is homologous to Γ v j ,w j (resp. to Γ ). Therefore it is equal to 1 2 . Using the same argument as for the computation of ( 25), we obtain that the integral of the second term in the sum

ik 8π Γ - Γv ,v j K v j ,w j+1 f w j+1 (u + 2K)f v j (u)du = - ik 8π Γ - Γv ,v j f w j+1 (u + 2K)f w j+1 (u)du = 1 2 .
Therefore (KK -1 ) v j ,v j = 1 2 + 1 2 = 1, which completes the proof. Note that the proof uses essentially the fact that the contour Γ x,y winds once vertically, but makes no use of the horizontal winding of the contour, which can be arbitrary. However, "verticality" of the contour plays a crucial role for the exponential decay of the coefficients of K -1 , as stated below in Theorem 13.

Asymptotics for the inverse Kasteleyn operator K

For any x, y ∈ G, define

χ(u) = 1 |x -y| log{e (x,y) (u + 2iK )}, (26)
with the exponential function introduced in Section 2.2.5. The main result of this section (Theorem 13) shows the exponential decay of the inverse Kasteleyn operator, with a rate that can be directly computed in terms of χ.

Since |x -y| will be typically large in this section concerned with asymptotic results, we are in the general case, according to Section 3.3.2. The poles of the exponential function are {τ j }; they belong to a sector of size strictly less than 2K, say τ j ∈ τ + (-K, K).

Theorem 13. Assume that k = 0. As |x -y| → ∞, one has

K -1 x,y = -f x (u 0 + 2K + 2iK )f y (u 0 + 2iK ) 4 2π|x -y|χ (u 0 ) e |x-y|χ(u 0 ) • (1 + o(1)),
where u 0 is the unique u ∈ τ + 2K + (-K, K) such that χ (u) = 0, and χ(u 0 ) < 0.

Proof. It consists in applying the saddle-point method to the contour integral [START_REF] De Tilière | Critical Ising model and spanning trees partition functions[END_REF]. It is very similar to the proof of Theorem 14 in [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF], which is devoted to the derivation of the asymptotics of the Green function [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]. Indeed, the integrands of ( 10) and ( 18) only differ by the prefactor function f x (• + 2K)f y (•) (as well as a constant multiplicative term). This prefactor function will affect the asymptotics by multiplying by its value at the saddle-point u 0 ± 2iK the asymptotics of [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]. Let us give some brief details.

• It follows from [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]Lemma 15] that the equation χ (u) = 0 has a unique solution u 0 in the interval τ + 2K + (-K, K), and moreover χ(u 0 ) < 0 (cf. [10, Lemma 16]). The point u 0 ± 2iK will be interpreted as the saddle-point.

• We then move the contour Γ x,y of ( 18) into a new one, denoted by Γ x,y , going through u 0 ± 2iK and satisfying some further properties. The validity of this change of contour is based on the fact that neither the exponential function nor the prefactor have poles in the sector s x,y , see Figure 8.

• We adjust the new contour Γ x,y so as to have, classically, a contribution exponentially negligible outside a neighborhood of u 0 ±2iK (this can be done by introducing suitable steepest descent paths).

• In the neighborhood of u 0 ± 2iK , we apply (a uniform version of) the saddle-point method, which eventually yields to the expansion written in Theorem 13.

Remark 14. Let us note that the constant -f x (u 0 + 2K + 2iK )f y (u 0 + 2iK ) in Theorem 13 is positive. Indeed, by [START_REF] Cimasoni | The critical temperature for the Ising model on planar doubly periodic graphs[END_REF], f(u) is the sum of one or two terms nc( u-α j

2 ). Due to the location of the poles described in Section 3.3.1, α j ∈ τ + (-K, K). Hence

u 0 ± 2iK -α j 2 ∈ K + iK + (-K, K),
and by Table 2, nc(

u 0 ±2iK -α j 2 ) = nc(K + iK + v 0 ) = ik -1 k cn(v 0 )
, with some v 0 ∈ (-K, K). Consequently, the constant -f x (u 0 + 2K + 2iK )f y (u 0 + 2iK ) equals minus the product of two (sums of) terms ik -1 k cn(v 0 ). The positivity follows from -i 2 = 1.

The case where G is Z 2 -periodic

In this section, we suppose moreover that the graph G is Z 2 -periodic, implying that the graph G F is also Z 2 -periodic. We consider the dimer model on the graph G F arising from a Z-invariant Ising model on G, and the Z-invariant rooted spanning forest model on G. Note that the weight function corresponding to each of the models is periodic. Consider the natural operators associated to the two models, that is the Kasteleyn operator K acting on C V F , arising from a periodic admissible orientation of the edges of G F3 ; and the massive Laplacian operator ∆ m acting on C V . Using Fourier techniques, see for example [START_REF] Cohn | A variational principle for domino tilings[END_REF], an important tool for understanding each of the models is the characteristic polynomial of the respective operators. In this section, we prove that the characteristic polynomials of the two models are equal up to an explicit constant. We state implications of this fact for the spectral curve of the dimer model on G F .

In the periodic case, we have two explicit expressions for the inverse operator K -1 . The one given by Theorem 11 and the one obtained using Fourier techniques. In Corollary 18, we prove that the two are equal.

These two facts are used in Section 3.6 for obtaining explicit expressions for the dimer model on G F and relating its free energy to that of the rooted spanning forest model.

Quasi-periodic functions

A natural toroidal exhaustion of G is {G n = G/nZ 2 } n 1 ; in a similar way {G F n = G F /nZ 2 } n 1 is a toroidal exhaustion of G F . The smallest graphs G 1 and G F
1 of the exhaustions are known as the fundamental domains.

We note with an addition sign the action of Z 2 on vertices, edges, faces of G and G F . Let γ x , γ y be two simple paths in G * connecting a given face f 0 , with f 0 + (1, 0) and f 0 + (0, 1) respectively.

To simplify notation, we also write γ x , γ y for the images of these paths when quotienting by the action of Z 2 , which are now non trivial cycles of G * 1 generating the first homology group of the torus on which G 1 is drawn. For (z, w) ∈ C 2 , denote by C V (z,w) the space of (z, w)-quasi-periodic functions on vertices of G:

∀ x ∈ V, ∀ (z, w) ∈ C 2 , f (x + (m, n)) = z -m w -n f (x).
For every vertex x of G 1 , define δ x (z, w) to be the (z, w)-quasi-periodic function equal to 0 on vertices which are not translates of x and to 1 at x. Then the collection {δ x (z, w)} x∈V 1 is a natural basis for C V (z,w) . Note that γ x , γ y can be deformed into directed cycles of the dual graph (G F 1 ) * (or of the diamond graph G 1 ). In a way similar to C V (z,w) we define C V F (z,w) , the space of (z, w)-quasiperiodic functions on vertices of G F , with basis {δ x (z, w)} x∈V F 1 .

Characteristic polynomials and spectral curves

Let K(z, w) be the matrix of the restriction of K to the space C V F (z,w) in the basis {δ x (z, w)} x∈V F 1 . The matrix K(z, w) is obtained from the Kasteleyn matrix of the fundamental domain G F 1 as follows: multiply coefficients of edges crossing γ x by z (resp. z -1 ) if the edge goes from the left of γ x to the right (resp. from the right of γ x to the left); coefficients of edges crossing γ y are multiplied by w or w -1 .

In a similar way, ∆ m (z, w) is the matrix of the restriction of ∆ m to the space C V (z,w) in the basis {δ x (z, w)} x∈G 1 .

The dimer characteristic polynomial of G F is the determinant of the matrix K(z, w),

P K (z, w) = det K(z, w).
The massive Laplacian characteristic polynomial of G is the determinant of ∆ m (z, w),

P ∆ m (z, w) = det ∆ m (z, w).
Consider a Laurent polynomial P (z, w) in two complex variables z, w. Then, the spectral curve of the polynomial, denoted by C(P ), is defined to be its zero locus:

C(P ) = {(z, w) ∈ (C * ) 2 : P (z, w) = 0}.
The following proves that the two characteristic polynomials are equal up to an explicit multiplicative constant. The constant is determined in Section 3.6.3. Proposition 15. There exists a nonzero constant c such that

P K (z, w) = cP ∆ m (z, w).
Proof. Let us first prove that P ∆ m (z, w) divides P K (z, w), using an argument similar to that of [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: the periodic case[END_REF]Lemma 11]. By [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]Proposition 21], any point of the spectral curve of the Laplacian is of the form (z(u), w(u)), with u ∈ T(k), and

z(u) = e iα ∈γx i √ k sc( u-α 2 ) , w(u) = e iα ∈γy i √ k sc( u-α 2 ) . (27)
The function g (•,y) (u) of Equation ( 14) is (z(u), w(u))-quasi-periodic, as it involves the function f • (u), which is invariant by translations, and the discrete massive exponential function, which is (z(u), w(u)) quasi-periodic. By Proposition 9 it is in the kernel of the Kasteleyn operator K. Therefore, P K (z(u), w(u)) = 0. As a consequence, P ∆ m (z, w) divides P K (z, w).

Now by [START_REF] Dubédat | Exact bosonization of the Ising model[END_REF] we know that, up to a constant, P K (z, w) is equal to the dimer characteristic polynomial on G Q . The graph G Q being bipartite, the corresponding spectral curve is a Harnack curve [START_REF] Kenyon | Dimers and amoebae[END_REF]. Hence, the characteristic polynomial on G Q , and thus P K (z, w), are irreducible.

The fact that P ∆ m (z, w) divides P K (z, w) and that P K (z, w) is irreducible implies that the two polynomials are equal up to a constant, and concludes the proof.

By Proposition 15, properties of the spectral curve of P ∆ m obtained in [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF], are automatically transferred to the spectral curve of the dimer characteristic polynomial P K .

Corollary 16. Let k 2 = 0.

• The spectral curve of the dimer model on G F is a Harnack curve of genus 1, with (z, w) ↔ (z -1 , w -1 ) symmetry.

• Every Harnack curve of genus 1 with (z, w) ↔ (z -1 , w -1 ) symmetry arises from such a dimer model.

• The characteristic polynomial P K (z, w) has no zero on the unit torus

{(z, w) ∈ C 2 : |z| = 1, |w| = 1}.
Proof. For k 2 > 0, this follows from our results of the paper [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]. For the last point we use the fact that (0, 0) does not belong to the amoeba of the spectral curve. For k 2 < 0, we prove later in Section 4 that the spectral curve is the same as for an elliptic parameter (k * ) 2 > 0, such that (k * ) k = 1.

Remark 17. For k = 0, the spectral curve of the dimer model on G F is still the spectral curve of the Laplacian with conductances tan θ. It is a Harnack curve of genus 0 with the same symmetry, and a double point at z = w = 1, see [START_REF] Kenyon | Planar dimers and harnack curves[END_REF][START_REF] Boutillier | The critical Z-invariant Ising model via dimers: the periodic case[END_REF][START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF].

Inverse Kasteleyn operator

Because of the periodicity of the graph, it is possible to define an inverse for the Kasteleyn matrix, by inverting in Fourier space the multiplication operator K(z, w): define the operator K -1 by its coefficients

K -1 x+(m,n),x +(m ,n ) = 1 (2πi) 2 {(z,w)∈C 2 :|z|=1,|w|=1} [Cof K(z, w)] x ,x P K (z, w) z m -m w n -n dz z dw w , (28)
for all x, x ∈ G F 1 , and m, m , n, n ∈ Z. By Corollary 16, we know that the characteristic polynomial P K (z, w) has no zero on the unit torus. This means that Proposition 5 of [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: the periodic case[END_REF] holds, and we have: Proposition 18. The inverse Kasteleyn operators K -1 and K -1 given by Equations ( 18) and (28), respectively, are equal.

Proof. For k = 0, there is no root of P K (z, w) on the unit torus, by Corollary 16. As a consequence, the quantities K -1 x+(m,n),x are the Fourier coefficients of an analytic periodic function, and as such, decay exponentially fast when |m| + |n| goes to ∞. In particular, these coefficients are bounded. By the uniqueness statement in Theorem 11, K -1 and K -1 are equal.

For k = 0, adaptating the computation of the asymptotics of the integral formula for K -1 from [8, Corollary 7] proves that these coefficients go to 0. The result then follows from [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: the periodic case[END_REF]Proposition 5], stating uniqueness of the inverse of the Kasteleyn operator with coefficients tending to 0 at infinity on a periodic Fisher graph.

Note that similarly to what has been done for the Green function of the Z-invariant Laplacian [10, Section 5.5.1], it is also possible to directly understand the transformation from the double integral expression [START_REF] Kasteleyn | The statistics of dimers on a lattice: I. the number of dimer arrangements on a quadratic lattice[END_REF] to the contour integral [START_REF] De Tilière | Critical Ising model and spanning trees partition functions[END_REF], by computing one integral (e.g., with respect to w) by residues, and then by perfoming the change of variable from z on the spectral curve to u ∈ T(k).

Dimer model on the graph G F

Consequences of the results of Sections 3.3-3.5 on the dimer model on G F are investigated: in Section 3.6.1, we prove an explicit local expression for a Gibbs measure, in the case where the graph G is periodic or not. In Section 3.6.2, we prove an explicit local expression for the free energy of the dimer model. Combining this with the high temperature expansion, we deduce an explicit local formula for the free energy of the Ising model, as a sum of contributions for each edge of the fundamental domain, similar to the one given by Baxter, see [5, (7.12.7)] and [START_REF] Baxter | Solvable eight-vertex model on an arbitrary planar lattice[END_REF].

Gibbs measure

When the graph is infinite, the notion of Boltzmann measure is replaced by that of Gibbs measure. A Gibbs measure on the set of dimer configurations of G F is a probability measure satisfying the DLR conditions: if one fixes a perfect matching in an annular region of G F , then perfect matchings inside and outside of this annulus are independent. Moreover, the probability of an interior perfect matching is proportional to the product of its edge-weights. Let F be the σ-field generated by cylinders of G F , a cylinder being the set of dimer configurations containing a fixed, finite subset of edges of G F . Following the argument of [START_REF] Cohn | A variational principle for domino tilings[END_REF][START_REF] De Tilière | Quadri-tilings of the plane[END_REF], we obtain the following.

Theorem 19. There is a unique probability measure P dimer on (M(G F ), F), such that the probability of occurrence of a subset of edges

E = {e 1 = x 1 y 1 , . . . , e n = x n y n } of G F in a dimer configuration of G F is: P dimer (e 1 , . . . , e n ) = n j=1 K x j ,y j Pf[ t (K -1 ) E ], ( 29
)
where K -1 is the inverse Kasteleyn operator whose coefficients are given by [START_REF] Duffin | Potential theory on a rhombic lattice[END_REF], and (K -1 ) E is the sub-matrix of K -1 whose rows and columns are indexed by vertices of E. The measure P dimer is a Gibbs measure.

Moreover, when the graph G F is Z 2 -periodic, the measure P dimer is obtained as weak limit of the dimer Boltzmann measures on the toroidal exhaustion G F n , and coefficients of K -1 are also given by K -1 of Equation [START_REF] Kasteleyn | The statistics of dimers on a lattice: I. the number of dimer arrangements on a quadratic lattice[END_REF].

Proof. Convergence of the Boltzmann measures in the periodic case follows the argument of [START_REF] Cohn | A variational principle for domino tilings[END_REF]. The delicate issue in the convergence comes from the possible zeros of the dimer characteristic polynomial on the torus {(z, w) : |z| = 1, |w| = 1}. By Corollary 16, we know that whenever k = 0, it has no zero, and the argument goes trough. When k = 0, it has a double zero at (1, 1), the argument is more delicate and has been done in [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: the periodic case[END_REF].

The argument in the non-periodic case follows that of [START_REF] De Tilière | Quadri-tilings of the plane[END_REF]. The key requirements are the convergence of the Boltzmann measure in the periodic case, uniqueness of the inverse Kasteleyn operator decreasing at infinity, and the locality property of the formula given by Theorem 11.

Example. Consider an edge e = v j (x)v (y) = v j v of G F corresponding to an edge e = xy of G with rhombus half-angle θ e . Then, the probability P dimer (e) that this edge occurs in a dimer configuration of G F , or equivalently the probability that this edge occurs in a high temperature polygon configuration of G, is equal to

P dimer (e) = K v j v K -1 v ,v j = 1 2 - 1 -2H(2θ e ) 2 cn θ e . ( 30
)
This is computed in Lemma 45 of Appendix B.

To compute probabilities of edges occurring in low temperature polygon configurations one uses the duality relation of Section 4.2.

Free energy

Suppose that the isoradial graph G is Z 2 -periodic. The free energy F of a model is defined to be minus the exponential growth rate of its partition function, that is

F k dimer = -lim n→∞ 1 n 2 log Z dimer (G F n , ν), F k Ising = -lim n→∞ 1 n 2 log Z Ising (G n , J), F k forest = -lim n→∞ 1 n 2 log Z forest (G n , ρ, m),
where ν, J, m 2 and ρ are given by ( 5), ( 2), ( 8) and [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: the periodic case[END_REF], respectively.

In Theorem 20 we prove an explicit local formula for the free energy F k dimer of the dimer model on the graph G F arising from a Z-invariant Ising model on G. From this we deduce an explicit formula for the free energy Proof. By Corollary 16 the dimer characteristic polynomial P K (z, w) has no zero on the unit torus. This implies [START_REF] Cohn | A variational principle for domino tilings[END_REF] that the free energy F dimer is equal to

F k dimer = -(|E 1 | + |V 1 |) log 2 2 + e∈E 1 - 1
F dimer = - 1 2 |z|=|w|=1 log P K (z, w) dz 2iπz dw 2iπw . ( 31
)
Following an idea of Kenyon [START_REF] Kenyon | The Laplacian and Dirac operators on critical planar graphs[END_REF], the next step consists in studying its variation as the embedding of the graph is modified by tilting the train-tracks. Note that the idea of tilting the train-tracks can already be found in the section free energy of [START_REF] Baxter | Solvable eight-vertex model on an arbitrary planar lattice[END_REF].

Let us consider a smooth deformation of the isoradial graph G, i.e., a continuous family of isoradial graphs (G(t)) t∈[0,1] obtained by varying the directions (α T (t)) of the train-tracks smoothly with t, in such a way that G(1) = G and G(0) = G flat , where G flat is an isoradial graph whose edges have half-angles θ flat equal to 0 or π 2 . Let G F (t) be the Fisher graph corresponding to the isoradial graph G(t), and let F dimer (t) be the corresponding free energy. We thus have F dimer = F dimer [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF].

As the angles of the train-tracks vary smoothly with t, recalling that P K (z, w) = det K(z, w), one has:

dF dimer (t) dt = - 1 2 |z|=|w|=1 d dt log det K(z, w) dz 2πiz dw 2πiw = - 1 2 |z|=|w|=1 x,y∈V F 1 ∂ log det K(z, w) ∂K(z, w) x,y dK(z, w) x,y dt dz 2πiz dw 2πiw = - 1 2 x,y∈V F 1 |z|=|w|=1 (K(z, w) -1 ) y,x dK(z, w) x,y dt dz 2πiz dw 2πiw = - e=xy∈E F 1 K -1 y,x dK x,y dt = - e=xy∈E F 1 P dimer (e) d log ν e dt .
In the penultimate line we used that for an invertible matrix M = (M i,j ), one has

∂ log det M ∂M i,j = (M -1 ) j,i .
In the last line we used: the explicit expression for K -1 given by Equation ( 28), Corollary 18 and Theorem 19 implying that K x,y K -1 y,x = K x,y K -1 y,x = P dimer (e), and the fact that ν e = |K x,y |.

The weight of an edge e of a decoration of G F 1 is equal to 1, i.e., is independent of the rhombus angle, so that it does not contribute to dF dimer (t) dt .

If e corresponds to an edge of G 1 having rhombus angle θ e , we have by Equations ( 30) and ( 3):

P dimer (e) =
F dimer (1) = F dimer (0) + 1 0 e∈E 1 dF dimer (t) dt dt = F dimer (0) - e∈E 1 θe θ flat e P(θ) d log ν(θ) dθ dθ. ( 33
)
We first consider the integral part of the above equation, and then the term F dimer (0). Replacing P(θ) using ( 32), we have:

-P(θ) d log ν(θ) dθ dθ = - 1 2 - 1 -2H(2θ) 2 cn θ d log ν(θ) dθ dθ = - 1 2 log ν(θ) + 1 -2H(2θ) 2 cn θ d log ν(θ) dθ dθ.
Using the explicit expression of ν(θ) (see again [START_REF] Kenyon | An introduction to the dimer model[END_REF]) and the formulas sn = cn dn, cn = -sn dn, cn 2 + sn 2 = 1, gives

dν(θ) dθ = (cn θ + cn 2 θ + sn 2 θ) dn θ (1 + cn θ) 2 = dn θ 1 + cn θ , which readily implies that d log ν(θ) dθ = dn θ sn θ = ds θ.
As a consequence,

1 -2H(2θ) 2 cn θ d log ν(θ) dθ = 1 -2H(2θ) 2 cn θ ds θ = 1 -2H(2θ) 2 d log sc θ dθ ,
using that sc = dn cn 2 ⇒ (log sc) = dn sn cn . Integrating by parts, we obtain

-P(θ) d log ν(θ) dθ dθ = - 1 2 log sn θ 1 + cn θ + 1 -2H(2θ) 2 log sc θ + 2H (2θ) log sc θ dθ.
We now need to compute -1 2 log sn θ 1+cn θ + ( 1-2H(2θ)

2

) log sc θ in the limit θ → θ flat for possible values of θ flat , i.e., for θ flat ∈ {0, K} since θ flat ∈ {0, π 2 }.

When θ → K. We have sn K = 1, cn K = 0 and H(2K) = 1 2 , see (71) for k 2 > 0 and the duality relation (67) for k 2 < 0; so that

- 1 2 log sn K 1 + cn K + 1 -2H(2K) 2 log sn K = 0. Moreover as θ → K, 1-2H(2K) 2 = O(θ -K) implying that lim θ→K 1-2H(2K) 2 log cn θ = 0. Wrapping up, we have lim θ→K - 1 2 log sn θ 1 + cn θ + 1 -2H(2θ) 2 log sc θ = 0.
When θ → 0. We have sn 0 = 0, cn 0 = 1,

H(0) = K π Z(0) = 0; implying that - 1 2 log 1 1 + cn 0 - 1 -2H(0) 2 log cn 0 = 1 2 log 2.
We are left with handling, (-

1 2 + 1 2 -H(2θ)) log sn θ = -H(2θ) log sn θ. As θ → 0, H(2θ) = O(θ), thus lim θ→0 H(2θ) log sn θ = 0. Wrapping up, we have lim θ→0 - 1 2 log sn θ 1 + cn θ + 1 -2H(2θ) 2 log sc θ = 1 2 log 2.
Plugging this into Equation ( 33), we obtain

F dimer (1) =F dimer (0) - log 2 2 |{e ∈ E 1 : θ flat e = 0}|+ + e∈E 1 - 1 2 log sn θ 1 + cn θ + 1 -2H(2θ) 2 log sc θ + θe θ flat e 2H (
2θ) log sc θ dθ .

Let us now compute F dimer (0), that is the free energy of the dimer model on the flat graph G F flat . By Fisher's correspondence, for every n 1, the number of dimer configurations of the toroidal graph G F n,flat = G F flat /nZ 2 is equal to 2 |V 1 |n 2 times the weighted number of polygon configurations of G n,flat with edge-weights sn θ 1+cn θ arising from the high temperature expansion. Let us describe G flat , see also [START_REF] Kenyon | The Laplacian and Dirac operators on critical planar graphs[END_REF][START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF]. Since the sum of the rhombus-angles around vertices of G flat is 2π, and since rhombus half-angles are equal to 0 or π 2 , there is around every vertex exactly two rhombi with rhombus half-angle π 2 , with corresponding high temperature expansion weight equal to 1. The other rhombi have rhombus half-angle 0; with corresponding weight equal to 0. The graph G n,flat thus consists of p disjoint cycles covering all vertices, for some p = O(n). A polygon configuration has even degree at every vertex so that for each such cycle there is exactly two polygon configurations. As a consequence,

Z dimer (G n,flat ) = 2 |V 1 |n 2 2 p , implying that F dimer (0) = -|V 1 | log 2.
From the geometric description of the graph G flat , we also know that

|{e ∈ E 1 : θ flat e = 0}| = |E 1 | -|{e ∈ E 1 : θ flat e = π/2}| = |E 1 | -|V 1 |.
As a consequence, the dimer free energy F dimer = F dimer (1) is equal to

-(|E 1 | + |V 1 |) log 2 2 + e∈E 1 - 1 2 log sn θ e 1 + cn θ e + 1 -2H(2θ e ) 2 log sc θ e + θe θ flat e 2H (2θ) log sc θ dθ ,
and the proof is concluded by recalling that sn θe 1+cn θe = sc θe 2 dn θe 2 .

We obtain the free energy of the Z-invariant Ising model, similar to the one given by Baxter, see [5, (7.12.7)] and [START_REF] Baxter | Solvable eight-vertex model on an arbitrary planar lattice[END_REF].

Corollary 21. The free energy of the Z-invariant Ising model is equal to

F k Ising = -|V 1 | log 2 2 -|V 1 | K 0 2H (2θ) log sc θ dθ + e∈E 1 -H(2θ e ) log sc θ e + θe 0 2H (2θ) log sc θ dθ .
Proof. From the high temperature expansion, see Equation (1), we have for every n 1,

Z Ising (G n , J) =   e∈En cosh(J e )   Z dimer (G F n , ν).
Let us compute cosh(J e ) for the Z-invariant coupling constants. By Equation (4), we have cosh(2J e ) = nc θ e , so that cosh(J e ) = 1 + cosh(2J e ) 2 = 1 + cn θ e 2 cn θ e = 1 + cn θ e 2 sn θ e sc θ e .

As a consequence,

F Ising = F dimer - e∈E 1 log cosh(J e ) (35)
=

F dimer + |E 1 | log 2 2 + e∈E 1 1 2 log sn θ e 1 + cn θ e - 1 2 log sc θ e = -|V 1 | log 2 2 +
F k Ising = -|V 1 | log 2 2 + F k forest 2 . ( 36
)
Proof. This is obtained by comparing the expressions for the free energies proved in Corollary 21 and [10, Theorem 2].

Remark 23. In the case k = 0, Corollary 21 is obtained in [8, Theorem 3]; the free energy becomes (the key point being that as k → 0, H(u) → u 2π , see Lemma 42)

F 0 Ising = -|V 1 | log 2 2 - e∈E 1 θ e π log tan θ e + 1 π L(θ e ) + L π 2 -θ e ,
where L is the Lobachevsky function: L(θ e ) = -θe 0 log(2 sin(t))dt. Moreover, for k = 0, Corollay 22 is derived in [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF], below Theorem 3.

Remark 24. It is possible to recover Baxter's "local" formula for the free energy when k = 0 without deformation argument, but rather by computing more directly the double integral [START_REF] Kenyon | The Laplacian and Dirac operators on critical planar graphs[END_REF], by first fixing w and evaluating the integral in z with residues: since P K is reciprocal, and has no root on the unit torus, then up to a multiplicative constant, one can rewrite P K (z, w) for |z| = 1 as

d j=1 z j (w) + 1 - z z + j (w) 1 - z - j (w) z ,
where z ± j (w), j = 1, . . . , d are the roots of P K (•, w), |z + j (w)| > 1 and z - j (w) = (z + j (w)) -1 . The log of this product can be expanded in series in z, whose term of degree 0 is the sum of logarithms of the roots z + j (w): this is the contribution we get when dividing by 2iπz and integrating over the unit circle. Ending here the computation for the square lattice yields the expression of the free energy of the Ising model on Z 2 given in [START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF]Equation (7.9.16)]. We can go even further, for general periodic isoradial graphs, if we perform the change of variable from w to u, as in [10, Section 5.5] to pass from the Fourier expression to the local expression of the massive Green function on periodic isoradial graphs. The roots (z, w) of P K where |z| > 1 and |w| = 1 form a closed curve on the spectral curve, which is mapped under log | • | to a horizontal segment joining the two connected components of the boundary of the amoeba, and lifts on T(k) as a contour Γ winding once vertically. We can then write

|w|=1 d j=1 log z + j (w) dw 2iπw = Γ log z(u) w (u) w(u) du 2iπ ,
where z(u) and w(u) are given by [START_REF] Kadanoff | Some critical properties of the eight-vertex model[END_REF]. As in [10, Section 5.3.3] when computing the area of the hole of the amoeba, the product of log z(u) and w (u)/w(u) can be rewritten as a sum over intersections of train-tracks on G 1 , i.e., over edges of the fundamental domain, thus yielding a local formula.

Computing the constant in

P K (z, w) = cP ∆ m (z, w)
With the computations of free energies, we can now compute the constant in Proposition 15.

Corollary 25. The dimer characteristic polynomial of the graph G F arising from the Zinvariant Ising model and the Z-invariant spanning forests characteristic polynomial are related by

P K (z, w) = 2 |V 1 |+|E 1 | e∈E 1 cn θ e 1 + cn θ e P ∆ m (z, w). (37)
Proof. In Proposition 15, we proved that P K (z, w) = cP ∆ m (z, w). Moreover, the dimer and spanning forests free energies can be computed using the characteristic polynomials:

F dimer = - 1 2 {|z|=|w|=1} log P K (z, w) dz 2iπz dw 2iπw , F forest = - {|z|=|w|=1} log P ∆ m (z, w) dz 2iπz dw 2iπw .
This implies that

log c = F forest -2F dimer = 2F Ising -2F dimer + |V 1 | log 2, by Corollary 22, = -2 e∈E 1 log cosh J e + |V 1 | log 2, by Equation (35), = |E 1 | log 2 + e∈E 1 log cn θ e 1 + cn θ e + |V 1 | log 2, by (34).
This is coherent with [8, Corollary 14], concerning the case k = 0.

Duality and phase transition in the Z-invariant Ising model

This section is about the behavior of the Z-invariant Ising model as the elliptic parameter k varies. Section 4.2 exhibits a duality relation in the sense of Kramers and Wannier, see also [START_REF] Chan | The Ising susceptibility scaling function[END_REF][START_REF] Mccoy | The anisotropic Ising correlations as elliptic integrals: duality and differential equations[END_REF]. In Section 4.3 we derive the phase diagram of the model and compare it to that of the Z-invariant spanning forests of [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]. In Section 4.4 we extend to all isoradial graphs a self duality relation proved by Baxter in the case of the triangular lattice. Finally in Section 4.5 we relate duality transformations to the modular group.

Dual elliptic modulus

Let k ∈ [0, 1) be a fixed elliptic modulus, and recall the notation k = √ 1 -k 2 for the complementary parameter. By definition, the dual parameter of k is:

k * = i k k , or equivalently k * = 1 k . ( 38
)
Notice that k 2 → k * 2 is an involutive bijection between [0, 1) and (-∞, 0]. As we shall see, the dual parameter k * parametrizes the dual temperature.

We need the following duality identities relating elliptic integrals and Jacobi elliptic functions with parameters k and k * :

√ k K(k) = √ k * K(k * ), [

Duality in the Z-invariant Ising model

Kramers and Wannier's duality [START_REF] Kramers | Statistics of the two-dimensional ferromagnet[END_REF][START_REF] Kramers | Statistics of the two-dimensional ferromagnet[END_REF] says the following. Consider an Ising model on a planar graph G with coupling constants J and an Ising model on G * with coupling constants J * . Perform the high temperature expansion of the first Ising model, and low temperature expansion of the second. Both expansions yield polygon configurations on G. The Ising models are said to be dual if both induce the same measure on polygon configurations. This is true if the coupling constants satisfy the following duality relation:

∀ e ∈ E, tanh(J e ) = e -2J * e * ⇐⇒ sinh(2J e ) sinh(2J * e * ) = 1, (41)
where e * is the dual edge of e. Duality maps a high temperature Ising model on a low temperature one and vice versa.

In this setting, the Z-invariant Ising model on G with parameter k and the one on G * with parameter k * are dual models, see also [START_REF] Chan | The Ising susceptibility scaling function[END_REF] for the case of the triangular/hexagonal lattices. Indeed, for the first and second model we respectively have, for every edge e of G and dual edge e * of G * , 38), ( 39) and ( 40), from which we deduce the duality relation [START_REF] Li | Critical temperature of periodic Ising models[END_REF] for Z-invariant Ising models:

     sinh(2J e ) = sinh(2J(θ e | k)) = sc θ e 2K(k) π k , sinh(2J * e * ) = sinh 2J π 2 -θ e k * = sc K(k * ) -θ e 2K(k * ) π k * . Moreover, sc K(k * ) -θ e 2K(k * ) π k * = 1 k * cs θ e 2K(k * ) π k * ,
sinh(2J(θ e |k)) sinh 2J π 2 -θ e k * = 1. ( 42
)
The elliptic parameters k and k * can be interpreted as parametrizing dual temperatures.

Note that this duality relation together with the computation of Equation ( 30) can be used to obtain the probability of an edge occurring in a polygon configuration arising from the low temperature expansion of the Ising model.

Range of the Z-invariant Ising model and phase transition

The following lemma shows that the Z-invariant coupling constants are analytic in k 2 and that the whole range of temperatures is covered as the parameter k varies.

together with the fact that sinh -1 , tan, am and K are all increasing functions. Moreover, as k → 1, am(u|k) → 2 arctan(e u ) -π 2 , see [1, 16.15.4] From the duality relation [START_REF] Lis | Phase transition free regions in the Ising model via the Kac-Ward operator[END_REF] and from the case k 2 ∈ [0, 1), we deduce that J(θ e |k) is increasing on k 2 ∈ (-∞, 0] and goes to 0 as k 2 → -∞.

The analyticity of J(θ e |k) is clear on (-∞, 0) and (0, 1). For the neighborhood of 0 we use the series expansion of sc in terms of the Nome q = exp(-πK (k)/K(k)), see [1, 16.23.9]:

sc θ e 2K(k) π |k = π 2 √ 1 -k 2 K(k) tan(θ e ) + 2π √ 1 -k 2 K(k) ∞ n=1
(-1) n q 2n 1 + q 2n sin(2nθ e ). ( 44)

The following proves a second order phase transition, which together with Lemma 26 allows to derive the phase diagram of the Z-invariant Ising model (see Figure 14). Comments on the result are given in Remark 28.

Corollary 27 (Phase transition). The free energy F k Ising of the Z-invariant Ising model on an isoradial graph G is analytic for k 2 ∈ R\{0}. The model has continuous (or second order) phase transition at k = 0. Namely, one has

F k Ising = F 0 Ising -|k| 2 log |k| -1 |V 1 | 2 + O(k 2 ).
Proof. We start from Corollary 22 relating the free energy of the Z-invariant spanning forests and of the Ising model, and from [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]Theorems 36 and 38]. The analyticity for k 2 > 0 and the phase transition when k 2 > 0 tends to 0 immediately follow from these results.

In the case k 2 < 0, the key point is the expression of the free energy of rooted spanning forests, derived in [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]Theorem 38] whenever k 2 > 0. It is an integral expression in terms of the function H. It turns out that a similar expression holds in the case k 2 < 0. Performing an asymptotic expansion (along the same lines as in the proof of [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]Theorem 38]) of the so-obtained expression of F k Ising when k tends to 0 then concludes the proof. In the case k 2 < 0, we could also use the forthcoming Corollary 30 relating F k Ising and F k * Ising . This allows us to express F k Ising in terms of the free energy associated with the elliptic modulus k * , whose square is positive. In this way we can again use [START_REF] Kramers | Statistics of the two-dimensional ferromagnet[END_REF] and [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]Theorem 38].

Remark 28.

• In the particular case of the square lattice, Corollary 27 is derived in [5, (7.12.7)],

proving criticality at k = 0 of the Z-invariant Ising model on Z 2 . Our result is thus a generalization of the latter to all isoradial graphs.

• Criticality for the Z-invariant Ising model has been proved in [START_REF] Li | Critical temperature of periodic Ising models[END_REF][START_REF] Cimasoni | The critical temperature for the Ising model on planar doubly periodic graphs[END_REF][START_REF] Lis | Phase transition free regions in the Ising model via the Kac-Ward operator[END_REF], with a different parametrization of the temperature and with different techniques: the authors multiply the Z-invariant weights at k = 0 by an inverse temperature parameter β, and prove that when β is equal to the particular value β c := 1, the model is critical. In their setting, the difference β -β c is to be compared to the first nonzero term in the expansion of J(θ|k)-J(θ|0) as k goes to zero. This expansion, which takes into account the fact that θ = θ 2K(k) π is itself a function of k, reads:

J(θ|k) = J(θ|0) + k 2 8 sin( θ) + O(k 4 ),
yielding that k 2 (β -β c ), which allows to compare the free energy as a function of β close to β c (as the one given by Onsager [START_REF] Onsager | Crystal statistics. I. A two-dimensional model with an order-disorder transition[END_REF]) and Corollary 27.

• What is remarkable and not present in the physics literature is that the phase transition of the Z-invariant Ising model is (up to a multiplicative factor 1 2 ) the same as the phase transition of the Z-invariant spanning forest model. As explained in the proof, this follows from fact that the free energies of the two models are related by a simple formula proved in Corollary 22.

We deduce the following phase diagram for the Z-invariant Ising model:

• k 2 = 0: critical Ising model,
• k 2 ∈ (0, 1): low-temperature Ising model,

• k 2 ∈ (-∞, 0): high temperature Ising model.

Note that the phase diagram is nicer when expressed with the complementary elliptic modulus (k ) 2 = 1 -k 2 ∈ (0, ∞), see Figure 14. In the rest of the paper we have nevertheless chosen to use the elliptic parameter k since it is the one classically used in the notation of elliptic functions. The domain (k ) 2 < 0 (or equivalently k 2 > 1) on Figure 14 corresponds to the reciprocal parameter, also named Jacobi's real transformation. Jacobi elliptic functions are still defined for k 2 > 1, see [1, 16.11]. However, one main difference is that the period K(k) is not real anymore. Accordingly, the angles θ e = θ e

2K

π have a non-zero imaginary part. We now examine the effect of this reciprocal transformation on the coupling constant J(θ e |k) defined in [START_REF] Au-Yang | Q-dependent susceptibilities in ferromagnetic quasiperiodic Z-invariant ising models[END_REF]. As it does not seem natural to use complex angles, we extend the formula [START_REF] Au-Yang | Q-dependent susceptibilities in ferromagnetic quasiperiodic Z-invariant ising models[END_REF] in the regime k 2 > 1 with only the real parts of the angles. Using [1, 16.11] in [START_REF] Au-Yang | Q-dependent susceptibilities in ferromagnetic quasiperiodic Z-invariant ising models[END_REF], one finds that for k 2 > 1

J(θ e |k) = 1 2 log 1 + 1 k sn(θ e 2K(1/k) π |1/k) dn(θ e 2K(1/k) π |1/k) . ( 45
)
Using similar arguments as in Lemma 26, we see that when k 2 goes from 1 to +∞, the coupling constant [START_REF] Martìnez | Multi-spin correlation functions for the Z-invariant Ising model[END_REF] goes decreasingly from +∞ to 0. So this is the same as for the classical regime k 2 ∈ (-∞, 1); this range of parameter does not allow to reach a new regime, as for instance the anti-ferromagnetic regime.

Self-duality for the Z-invariant Ising model

We now prove a self-duality relation for the Z-invariant massive Laplacian. From this and Corollary 22, we deduce a self-duality relation for the Z-invariant Ising model, see Remark 31.

Lemma 29. The Laplacian operators associated to k and k * satisfy the following self-duality relation:

√ k * ∆ m(k * ) = √ k ∆ m(k) ,
and hence the discrete massive harmonic functions are the same.

Proof. Due to the particular form of the Laplacian operator (6), Lemma 29 is equivalent to proving that √ k sc(θ e |k) = √ k sc( 2K(k) π θ e |k) and √ k n j=1 A(θ j |k) are self-dual. For the first quantity, this directly follows from [START_REF] Lawden | Elliptic functions and applications[END_REF] and [START_REF] Kuperberg | An exploration of the permanent-determinant method[END_REF]. For the second one, we write θ j = α j+1 -α j 2

and use n times the addition theorem (72) and finally the periodicity relation (73). We obtain

√ k n j=1 A(θ j |k) = - n j=1 { √ k sc(α j+1 |k)}{ √ k sc(α j |k)}{ √ k sc(α j+1 -α j |k)},
which is self-dual, for the same reasons as previously.

Corollary 30. The free energy of the Z-invariant Ising model on the graph G satisfies the following self-duality relation.

F k Ising + |V 1 | 2 log k = F k * Ising + |V 1 | 2 log k * .
Remark 31.

• In the case of the triangular lattice this is proved in [5, (6.5.1)], our result thus extends the latter to all isoradial graphs.

• There is no simple self-duality relation between the coupling constants J(θ e |k) and J(θ e |k * ) so that this result is not straightforward. Baxter's argument in the triangular case reads as follows: he transforms the Z-invariant Ising model with parameter k on the triangular lattice into the one on the honeycomb lattice with the same parameter k by using Y -∆ moves, from this he deduces that the partition functions differ by an explicit constant; then he uses Kramers and Wannier duality to map the partition function of the Z-invariant Ising model with parameter k on the honeycomb lattice into the one of the triangular lattice (dual graph) with parameter k * . Making the constants explicit allows to relate the free energies with parameters k and k * on the triangular lattice. It is not obvious that this argument should extend to general isoradial graphs, even though one can generically go from a periodic isoradial graph to its dual using Y -∆ moves4 : when working out the constants in Baxter's computation, there seems to be some cancellations that are specific to the triangular and honeycomb lattices.

• This self-duality relation and the assumption of uniqueness of the critical point is used in [5, (6.5.5)-(6.5.7)] to compute the critical temperature of the Ising model on the triangular and honeycomb lattice. Corollary 30 allows to extend this physics argument to all isoradial graphs.

Dualities of the Ising model and the modular group

Various changes of the elliptic modulus k are considered throughout this article. Besides the intrinsic complementary transformation k → k , we have seen the importance of the dual transformation k → k * , as many quantities are self-dual:

• √ k K(k), see (39)
;

• √ k sc( 2K(k) π u|k), see (40);
• the exponential function (9);

• the modified Laplacian operator, see Lemma 29;

• the modified free energy, see Corollary 30;

• the rescaled function H, as

H(K(k)u|k) = H(K(k * )u|k * ), see (67).
Moreover, in the proof of Theorem 13, we make use of the ascending Landen transformation

k → 2-k 2 -2 √ 1-k 2 k 2
(note, this does not appear explicitly in the proof, as we refer to the companion paper [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF] for the details).

Our aim in this paragraph is twofold: first we reformulate the self-duality as a parity property of expansions in terms of the Nome q, then we relate the various transformations of k to the modular group. Links between the Ising model and the modular group already exist in the physics literature, see in particular [START_REF] Krieger | Mathematically modeling the most everyday of physical phenomena[END_REF]Chapter 8] as well as [START_REF] Maillard | Modular invariance in lattice statistical mechanics[END_REF][START_REF] Bostan | The Ising model: from elliptic curves to modular forms and Calabi-Yau equations[END_REF].

Self-duality and expansions in terms of the Nome

Let us first mention that a function of the elliptic modulus k 2 is analytic at 0 if and only if it is analytic at 0 as a function of the Nome q = e -πK /K . Indeed, q is analytic in terms of k Accordingly, any generic quantity of our article admits an analytic expansion in terms of the Nome at 0. The following simple criterion translates the self-duality property as a matter of parity:

Lemma 32. Let f (k) be a function analytic in k 2 around 0. Then f is self-dual (i.e., f (k) = f (k * ))
if and only if its expansion in terms of the Nome is even.

Proof. When k is replaced by k * , the quarter-periods become K(k * ) = k K(k), see [START_REF] Kuperberg | An exploration of the permanent-determinant method[END_REF], and

K (k * ) = k (K (k) + iK(k)) = K (k) + iK(k) 1 -k * 2
, ([1, 17.4.17]).

It becomes obvious that q(k * ) = -q(k), Lemma 32 follows.

With Lemma 32 the question of finding expansions in terms of the Nome comes up. In fact, such expansions typically appear rather indirectly, when writing Fourier expansions of Jacobi functions or elliptic integrals; cf. [START_REF] Martìnez | Correlation functions for the Z-invariant Ising model[END_REF] for the Fourier expansion of the sc function, as well as [1, 16.23 and 16.38] and [START_REF] Lawden | Elliptic functions and applications[END_REF]Section 8.7] for a more systematic treatment.

On the other hand, the Ising weights [START_REF] Au-Yang | Q-dependent susceptibilities in ferromagnetic quasiperiodic Z-invariant ising models[END_REF], which are at the heart of our whole construction, are not self-dual. This default of duality is responsible for the non-analyticity (and in some sense of the phase transition, see Corollary 27) of the free energy F k Ising at 0.

Dualities and modular group

The modular group (see [START_REF] Lawden | Elliptic functions and applications[END_REF]Chapter 9] for an introduction) is the group generated by the transformations S(τ ) = -1/τ and T (τ ) = τ + 1, acting on the upper half-plane. This group is the set of all transformations τ → c + dτ a + bτ ,

with a, b, c, d ∈ Z such that ad -bc = 1.

Two pairs of complex vectors (1, τ ) and (1, τ ) generate exactly the same lattice Z + τ Z = Z + τ Z if and only if τ is obtained from τ by a modular transformation [START_REF] Mccoy | The anisotropic Ising correlations as elliptic integrals: duality and differential equations[END_REF].

The quantity τ should be interpreted as ratios of quarter-periods; for instance τ = iK K (and then the Nome is q = e iπτ ).

It is interesting to notice that both generators of the modular group correspond to a duality: S is the complementary duality and T the self-duality, see Table 1.

Elliptic modulus

Name of the change of modulus

Transformations of the modular group

k No change τ k Complementary S(τ ) = -1/τ k * (Self-)duality T (τ ) = τ + 1 2-k 2 -2 √ 1-k 2 k 2

Landen transformation 2τ

Table 1: Correspondance between changes of the elliptic modulus and transformations of the modular group.

The set of all transformations (47) with ad -bc 1 also forms a group, called the extended modular group. The quantity ad -bc is then named the order of the transformation [START_REF] Mccoy | The anisotropic Ising correlations as elliptic integrals: duality and differential equations[END_REF]. In our elliptic treatment of the Ising model we have also encountered higher order transformations: namely, the Landen ascending transformation (used in the proof of Theorem 13) has order 2, see again Table 1.

This short discussion suggests that combinatorial links could exist between any two Ising models associated with elliptic modulus whose τ 's are related by a transformation [START_REF] Mccoy | The anisotropic Ising correlations as elliptic integrals: duality and differential equations[END_REF].

5 The double Z-invariant Ising model via dimers on the graph G Q

In the whole of this section we consider the dimer model on the bipartite graph G Q arising from two independent Z-invariant Ising models defined on an infinite isoradial graph G. Edges of G Q are assigned the weight function ν of (5).

In Section 5.2 we introduce a one parameter family of functions in the kernel of the Kasteleyn operator K of this dimer model. This is the key object used in Section 5.3 to prove a local expression for an inverse of the operator K. In Section 5.4 we derive asymptotics of this operator, and in Section 5.5 we derive consequences for the dimer model Gibbs measure. We aslo give a few examples of computations.

Kasteleyn matrix/operator

Let us recall the construction of the bipartite graph G Q . Every edge of the graph G is replaced by a "rectangle" and the latter are glued together in a circular way using external edges. Each "rectangle" has two edges "parallel" to an edge of G and two edges "parallel" to the dual edge, see Figure 2. For instance, if G is the square lattice, G Q is the square-octagon lattice.

When the graph G is isoradial, so is the graph G Q with radii of circles being one half of those of G. The isoradial embedding of G Q is such that external edges have length 0 and "rectangles" are real rectangles; vertices of the rectangles are in the middle of the edges of the diamond graph G , and each rectangle is included in a rhombus of G , see Figure 15 or 16.

For the sequel it is useful to note that every vertex of G Q belongs to a unique rhombus of G .

The graph G Q being bipartite, its vertices can be split into white and black

V Q = W Q B Q .
In this case, the Kasteleyn matrix K has rows indexed by white vertices and columns by black ones. Following Kuperberg [START_REF] Kuperberg | An exploration of the permanent-determinant method[END_REF], instead of considering an admissible orientation of G Q as we have done for G F , one can assign phases (e iφ wb ) wb∈E Q to edges of G Q , in such a way that, for every face of G Q whose boundary vertices are w 1 , b 1 , . . . , w n , b n in counterclockwise order, we have

(-1) n-1 n j=1 e iφ w j b j e -iφ w j+1 b j = 1. (48)
In our case, we define the phasing of the edges to be:

e iφ wb =          1
if the edge wb is parallel to an edge e of G, i if the edge wb is parallel to the dual of an edge e of G, -ie -iθ if wb is an external edge and w belongs to a rhombus of G having half-angle θ.

The fact that Equation (48) holds is proved in [18, Lemma 4.1], see also [START_REF] Kenyon | The Laplacian and Dirac operators on critical planar graphs[END_REF].

Coefficients of the Kasteleyn matrix K are then given by, for every white vertex w and every

black vertex b of G Q , K w,b = e iφ wb ν wb ,
where ν is the dimer weight function [START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF], see also Figure 15. Before proving that this function is well defined, i.e., independent of the choice of path from b to w, we give some examples of computation that are useful for the sequel. 15. Let e iα and e iβ be the rhombus vectors of the rhombus of G Q containing the rectangle as in Figure 16, then the two rhombus vectors of (G Q ) of the edge:

• b 1 w are 1 2 e iα and 1 2 e iβ , implying that f (b 1 ,w) (u) = dc( u-α 2 ) dc( u-β 2 ),
• b 2 w are 1 2 e i(β-π) and 1 2 e iα , implying that

f (b 2 ,w) (u) = -ik nc( u-β+2K 2
) nc( u-α 2 ),

• b 3 w are 1 2 e iβ and 1 2 e i(β+π) , implying that

f (b 3 ,w) (u) = ie iθ dc( u-β 2 ) dc( u-β-2K 2
). We also compute f (b,w ) (u) for b ∈ {b 1 , b 2 }, where w is the white vertex facing w along the diagonal of the rectangle, see Figure 16. Then, the two rhombus vectors of (G Q ) of the edge:

• b 1 w are 1 2 e iβ and 1 2 e i(α+π) , implying that f (b 1 ,w ) (u) = -ik nc( u-β 2 ) nc( u-α-2K 2
),

• b 2 w are 1 2 e i(α+π) and 1 2 e i(β+π) , implying that

f (b 2 ,w ) (u) = dc( u-α-2K 2 ) dc( u-β-2K 2
).

Lemma 35. The function f is well defined, that is independent of the choice of path from b to w.

Proof. It suffices to check that when traveling around each face of the graph, the product of the contributions of the edges is 1. There are three types of faces to consider: rectangles which correspond to edges of the graph G (or G * ), faces corresponding to those of the graph G, and faces corresponding to those of the dual graph G * . As a consequence, for every k < (with cyclic notation for indices), we have

-1 j=k f (b j ,w j ) (u)f (w j ,b j+1 ) (u) = e -i 2 (α -α j ) sn( u-α 2 ) ns( u-α j 2 ). (53)
It is important to notice that the right-hand side of ( 53) is independent of the determination of the angles α and α j .

In particular, the product around the face is (with

α n+1 = α 1 + 4K) e -iπ sn( u-α n+1 2 ) ns( u-α 1 2 ) = -sn( u-α 1 2 -2K) ns( u-α 1 2 ) = 1.
Consider now a face of degree n of the graph G Q corresponding to a face of the graph G.

Using similar notation, the picture differs in that the vertex at the center of the face belongs to G, that black and white vertices are exchanged and that the angle 2θ j is at the center of the face, see Figure 17 (right). By (52), we have

f (b j ,w j ) (u)f (w j ,b j+1 ) (u) = e -iθ j sd( u-α j+1 2
) ds(

u-α j 2 ) = e -i 2 (α j+1 -α j ) sd( u-α j+1 2
) ds( u-α j

2 ), since we have θ j = 1 2 (α j+1 -α j ), see Figure 17 (right). As a consequence, for every k < (with cyclic notation for indices), we have

-1 j=k f (b j ,w j ) (u)f (w j ,b j+1 ) (u) = e -i 2 (α -α j ) sd( u-α 2 ) ds( u-α j 2 ). (54)
In particular, the product around the face is (with

α n+1 = α 1 + 4K) e -iπ sd( u-α n+1 2 ) ds( u-α 1 2 ) = -sd( u-α 1 2 -2K) ds( u-α 1 2 ) = 1.
Note that Equations ( 53) and ( 54) are used again in the proof of Lemma 40, which proves an alternative expression for the function f .

Next is the key proposition used in proving the local expression for an inverse of the Kasteleyn operator K. Multiplying these equations by the corresponding entries of the Kasteleyn matrix (namely, sn θ, i cn θ and -ie -iθ , see ( 49)) and summing them, one can group together terms having similar values of H. We are left with computing the sum of residues at 2iK . It turns out that this boils down to (57) for u = 2iK . Thus this sum equals 0. Therefore

Asymptotics of the inverse Kasteleyn operator

We first need to introduce some notation. For any b and w, there exists a path on the diamond graph G joining b and w, see In the same way as in Section 3.6, the inverse Kasteleyn operator K -1 can be used to obtain an explicit local expression for a Gibbs measure P Q dimer on dimer configurations of the infinite graph G Q arising from two independent Z-invariant Ising models. It can also be used to obtain an explicit local formula for the free energy of the model. By Dubédat [START_REF] Dubédat | Exact bosonization of the Ising model[END_REF] we know that this free energy is equal, up to an additive constant, to that of the dimer model on G F (since the characteristic polynomials differ by a multiplicative constant), so that we feel it presents no real interest to derive the formula, although it can be done using the approach of Theorem 20.

For the Gibbs measure, everything works out in exactly the same way so that we do not write out the details. We obtain that the probability of occurrence of a subset of edges E = {w 1 b 1 , . . . , w k b k } in a dimer configuration of G Q is:

P Q dimer (w 1 b 1 , . . . , w k b k ) =   k j=1 K w j ,b j   det[(K -1 ) E ],
where K -1 is the inverse Kasteleyn operator whose coefficients are given by ( 58) or ( 59) and (K -1 ) E is the sub-matrix of K -1 whose rows are indexed by b 1 , . . . , b k and columns by w 1 , . . . , w k .

A Useful identities involving elliptic functions

In this section we list required identities satisfied by elliptic functions.

A.1 Identities for Jacobi's elliptic functions Change of argument. Jacobi's elliptic functions satisfy various addition formulas by quarter-periods and half-periods, among which (cf. [1, Table 16.8]):

-u u ± K u + 2K u + iK u + 2iK u + K + iK sn -sn ± cd -sn Proof. First, it follows from the change of variable u → u + 2K and the above-mentioned property of nc that

1 2iπ Γ f(u + 2K)f(u)du = - 1 2iπ Γ-2K f(u + 2K)f(u)du = 1 2iπ Γ f(u + 2K)f(u)du,
where Γ is the contour Γ -2K crossed in the opposite direction of Γ. Further, using the 4iK -periodicity of the integrand, we deduce that

1 2iπ Γ f(u + 2K)f(u)du = 1 2 1 2iπ C f(u + 2K)f(u)du, ( 74
)
where C is the closed contour (Γ-2iK ) ( Γ-2iK ) S 1 S 2 , S 1 and S 2 being the horizontal segments joining (Γ -2iK ) and ( Γ -2iK ), see Figure 19. The main point is that the contour integral in the right-hand side of (74) can be computed with the residue theorem: the only pole (of order 1) is at α and has residue -2 k , see Table 2. Lemma 44 follows.

Γ C • • • • • 0 4K 4iK S 1

B.2 Inverse Kasteleyn operator at an edge

Here we compute the probability that a given edge e = xy of the isoradial graph G with rhombus half-angle θ e appears in the high temperature contour expansion of the Ising model. In terms of dimers, it corresponds to the probability that the corresponding edge e = v j (x)v (y) = v j v belongs to a dimer configuration of G F . By Theorem 19 this probability is given by P dimer (e) = K v j ,v K -1 v ,v j . Lemma 45. One has

K v j ,v K -1 v ,v j = 1 2 - 1 -2H(2θ e ) 2 cn θ e .
Proof. Instead of K v j ,v K -1 v ,v j we compute K v ,v j K -1 v j ,v . Both quantities are obviously equal, but the second one happens to be more convenient when applying the results of Section 3. in which we set, without loss of generality, α j = 0. This is equivalent to replace the function H(u) by H(u-α j ), which is possible because both functions satisfy the same jump conditions stated in Lemma 42. There are three residues, at 2θ e , 2K and 2i K . We thus have

(i K + K) nc(i K -θ e ), (78)
which may be slightly simplified, using that H(0) = 0. Thanks to (77), (78) and Table 2, we obtain that (76) equals The proof is complete.

K v ,v j K -1 v j ,v =
Using that lim k→0 H(u) = u 2π , see Lemma 42, we find that

lim k→0 K v j ,v K -1 v ,v j = 1 2 - π -2θ e 2π cos θ e .
This is in accordance with the computation of [8, Appendix A] for the case k = 0. Indeed, the dimer model considered in the latter paper corresponds to complementary polygon configurations meaning that the probability ( 30) is 1 minus the probability of the paper [START_REF] Boutillier | The critical Z-invariant Ising model via dimers: locality property[END_REF].

We now need to prove that when doing the inductive procedure around a cycle of G F , we recover the same angle mod 4π. There are two types of cycles to consider: inner cycles of decorations and cycles arising from boundary of faces of G.

Consider a cycle of a decoration corresponding to a vertex x of G of degree d, and let n be the number of edges of the inner cycle oriented clockwise. By definition of the angles, we have:

α d+1 -α 1 = d j=1 (α j+1 -α j ) = d j=1
2θ j + 2πn = 2π(n + 1).

Since the orientation of the edges is admissible, n is odd, thus proving that α d+1 -α 1 = 0 [4π].

Now consider a cycle C arising from the boundary x 1 , . . . , x m of a face of G, with vertices labeled in clockwise order, see also Figure 21. Up to a relabeling of vertices of the decorations, this cycle can be written as C = (w 2 (x 1 ), v 2 (x 1 ), v 1 (x 2 ), w 2 (x 2 ), . . . , v 1 (x 1 )). Using the definition of the angles within a decoration and in neighboring ones we deduce that, mod 4π, we have:

α 2 (x j ) -α 2 (x j+1 ) =
-π -2θ 1 (x j ) if co(w 2 (x j ), v 2 (x j ), v 1 (x j+1 ), w 2 (x j+1 )) is odd, π -2θ 1 (x j ) if co(w 2 (x j ), v 2 (x j ), v 1 (x j+1 ), w 2 (x j+1 )) is even.

Let n(C) denote the number of portions of the cycle C where co(w 2 (x j ), v 2 (x j ), v 1 (x j+1 ), w 2 (x j+1 )) is odd. Then, writing x m+1 = x 1 , we have: Since m j=1 (π -2θ 1 (x j )) is the sum of angles at the center of the cycle, it is equal to 2π. The orientation of the cycle being admissible, n(C) is odd, thus concluding the proof.

α 2 (x 1 ) -α 2 (x m+1 ) =

1 -

 1 2θ|k) log sc(θ|k) dθ + e∈E H(2θ e |k) log sc(θ e |k) + θe 0 2H (2θ|k) log(sc θ|k) dθ ,

Figure 1 :

 1 Figure 1: Left: a piece of a planar graph G and of a polygon configuration. Center and right: the corresponding Fisher graph G F and the two associated dimer configurations.

Figure 3 :

 3 Figure 3: Left: a piece of an infinite isoradial graph G (bold) with its circumcircles. Right: the diamond graph G .

Figure 4 :

 4 Figure 4: An edge e of G, the corresponding rhombus half-angle θ e and rhombus vectors e iαe , e iβ e .

Figure 6 :

 6 Figure 6: Notation for vertices of decorations, and rhombus vectors assigned to vertices. Since no confusion occurs, the argument x is omitted.

Figure 7 :Lemma 7 .Definition 3 . 1 .

 7731 Figure 7: Defining angles in neighboring decorations.

Figure 8 :

 8 Figure8: Left: the torus T(k) with the contours of integration Γ x,y and C x,y ; the poles of the integrand f x (u + 2K)f y (u) e (x,y) (u) are represented by black bullets and the pole 2iK of the function H by a black square. Right: in blue the sector s x,y used to define the contour of integration Γ x,y .

Figure 12 :

 12 Figure 12: Sectors s x i ,y when x = y = w j .

  F k Ising of the Z-invariant Ising model, see Corollary 21. Then in Corollary 22, we compare the latter to the free energy F k forest of Z-invariant spanning forests on G. Theorem 20. The free energy of the dimer model on the Fisher graph G F arising from the Z-invariant Ising model on G is equal to,

1 2 - 1 - 1 P

 1211 2H(2θ e ) 2 cn θ e := P(θ e ), ν e = sn θ e 1 + cn θ e := ν(θ e ). (θ e ) d log ν(θ e ) dθ e dθ e dt along the deformation, we have

e∈E 1 -Corollary 22 .

 122 H(2θ e ) log sc θ e + θe θ flat e 2H (2θ) log sc θ dθ . ) log sc θ dθ -|V 1 | K 0 2H (2θ) log sc θ dθ, since there are |V 1 | edges whose rhombus half-angle is π 2 in G flat . This concludes the proof. The free energy of the Z-invariant Ising model and spanning forests are related by

Lemma 26 .

 26 Let G be an isoradial graph and consider an edge e of G. Then the coupling constant J(θ e |k) = 1 2 log 1 + sn(θ e |k) cn(θ e |k) defined in (2), seen as a function of k 2 , is analytic on (-∞, 1) and increases from 0 to ∞ as k 2 goes from -∞ to 1. Proof. Jacobi's amplitude am(•|k), defined by am(u|k) = u 0 dn(v|k)dv, relates Jacobi and classical trigonometric functions. For the function sc we have sc(u|k) = tan(am(u|k)). The fact that J(θ e |k) is an increasing function of k 2 on [0, 1) comes from the relation sinh(2J(θ e |k)) = sc θ e 2K(k) π |k = tan(am(θ e 2K(k) π |k))

  , and thus am(θ e 2K(k) π |k) → π 2 . Using (43) then leads to lim k→1 J(θ e |k) = lim k→1 sinh(2J(θ e |k)) = ∞.

Figure 14 :

 14 Figure 14: Phase diagram of the Z-invariant Ising model on an isoradial graph G as a function of the complementary elliptic modulus (k ) 2 .

Figure 15 :

 15 Figure 15: Coefficients of the Kasteleyn matrix K around a rectangle face of G Q .

Example 5 . 1 .

 51 We compute f (b,w) (u) for b ∈ {b 1 , b 2 , b 3 }, where b 1 , b 2 , b 3 are the three black vertices incident to a white vertex w of G Q , see Figure

Figure 16 :

 16 Figure 16: Computation of f (b,w) (u) for b ∈ {b 1 , b 2 , b 3 }, and of f (b,w ) (u) for b ∈ {b 1 , b 2 }. To simplify the picture, the factor 1 2 is omitted in the notation of the rhombus vectors in the bottom part of the picture.

Figure 17 :

 17 Figure 17: Faces around dual (left) and primal (right) vertices.

Proposition 36 .

 36 Fixing a white base vertex w 0 of W Q , for every u ∈ T(k), the function f (•,w 0 ) (u), seen as a function on B Q , is in the kernel of the Kasteleyn operator K of the bipartite graph G Q .Proof. As we shall see, Proposition 36 follows from the identitysn(u + v) cn u -cn(u + v) dn v sn u -dn u sn v = 0, (55)which can be found in (iii) of Exercise 32 in [40,Chapter 2].By Equation (49), we need to prove that, for every white vertex w with neighbors b 1 , b 2 , b 3 as in Figure15, and every white base vertex w 0 , we have:Similarly, for the computation of K -1 b 2 ,w we have:2iπ{Res u=2K+α + Res u=β + Res u=2iK } 1 4iπ H(u)f (b 2 ,w) (u) = H(2K + α)(-ik ) -k cn(2K -θ) + H(β)(-ik ) -k cn θ + K π (-ik ) cn(iK -α/2) cn(K + iK -β/2).Finally, we obtain for K -1 b 3 ,w :2iπ{Res u=2K+β + Res u=4K+β + Res u=2iK } 1 4iπ H(u)f (b 3 ,w) (u) = -ie iθ H(2K + β) cd 0 + -ie iθ H(4K + β) cd(2K) + K π ie iθ cd(iK -β/2) cd(-K + iK -β/2).

  Terms with a H(2K + α) give:H(2K + α) sn θ -cd(K -θ) + cn θ cn(2K -θ) = H(2K + α) sn θ -sn θ + cn θ cn θ = 0.Similarly, those with a H(2K + β) give:H(2K + β) sn θ -cd(K + θ) + -ie -iθ (-ie iθ ) cd(0) = H(2K + β) -sn θ -sn θ -1 = 0.We group the terms in H(β) and H(4K + β), and use the fact that H(4K + β) = H(β) + 1, stated in Lemma 42:H(β) -ik i cn θ -k cn θ + H(4K + β) -ie -iθ (-ie iθ ) cd(2K) = H(β) (-1 + 1) + 1 = 1.

Figure 18 . 2 e 1

 1821 The first and the last edges of the path are half-edges of G , the other ones are plain edges. We call b = b 1 , and b n the black vertex adjacent to w. We further define b 2 , . . . , b n-1 as the successive black vertices in the middle of the edges of G joining b to w, see again Figure18. The n edges are equal to 1 iα 1 , e iα 2 , . . . , e iα n-1 , 1 2 e iαn . The α i are not well defined (in the sense that any multiple of 2π could be added to α i ), but the e α i are. The edges are orientated in such a way that b = b n .(60)We also define the points a j (j = 1, . . . , n -1) as the vertices of the diamond graph lying between b j and b j+1 . The notation of this paragraph is illustrated on Figure18.

Figure 18 :

 18 Figure 18: Notation for a path on the quad-graph G joining b and w, see (60).

Figure 19 :

 19 Figure 19: The contour C in (74) used in the proof of Lemma 44.

  e ) -H(2K) + H(2θ e -2K) 2 cn θ e + K k 2 2π sd θ e .The addition formula (69) results inH(2θ e -2K)-H(2K) = H(2θ e -4K) -K k 2 π sn K sn(θ e -K) sn(θ e -2K) if k 2 > 0 (-k 2 ) sd K sd(θ e -K) sd(θ e -2K) if k 2 < 0 = H(2θ e ) -1 -K k 2 πcn θ e sd θ e .

w 2 (x 1 ) v 2 (x 1 )v 1 (x 2 )w 2 (x 2 ) x 1 x 2 x 3 eCFigure 21 :

 212112212321 Figure 21: Notation for a cycle C arising from the boundary of a face of G.

m j=1 α 2

 2 (x j ) -α 2 (x j+1 ) = m j=1 (π -2θ 1 (x j )) -2πn(C).

  Fisher graph G F

	Let us compute the dimer weight function ν on G F corresponding to the Z-invariant Ising
	model on G with coupling constants J given by (2). For every edge e of G, we have
	tanh(J e ) =	e 2Je -1 e 2Je + 1	=	1+sn θe cn θe -1 1+sn θe cn θe + 1	=	1 + sn θ e -cn θ e 1 + sn θ e + cn θ e	=	sn θ e 1 + cn θ e	= sc	θ e 2	dn	θ e 2	,
	see [40, (2.4.4)-(2.4.5)] for the last identity.
	As a consequence of Section 2.1.3, the dimer weight function ν on the Fisher graph G F is
	ν e =	 1   1+cn θe = sc θe sn θe 2 dn θe 2	if e belongs to a decoration, if e corresponds to an edge e of G,
			  0				otherwise.

Table 2 :

 2 Addition formulas by quarter-periods and half-periods, taken from[1, 16.8].

	1 k ns	sn	1 k dc

In Equation (7.8.4), Baxter actually uses the complementary parameter k = √ 1 -k

and the parametrization, sinh(2Je) = -i sn(iθe|k ). The latter is equal to sc(θe|k) by[40, (2.6.12)].

When indicating sectors on the circle, the convention we adopt is that (α, β) represents the sector where the horizontal coordinate increases from α to β.

Such an orientation always exists by Theorem 3.1 of[START_REF] Cimasoni | Dimers on surface graphs and spin structures[END_REF], using that the number of vertices of the fundamental domain of G F is even.

If on the torus, there is a unique train-track with a given homology class, making it move across the torus through every pairwise intersection once by Y -∆ moves yields the same rhombic graph, but with the role of primal and dual vertices exchanged. If there are several of them, then moving them all to put the first one instead of the second one, the second instead of the third, etc., yields the result. Note though that this is not possible in the case when there are just two different homology classes, like in the case of the square lattice.

4iπ Γ b,w f (b,w) (u)du,(58)

i=1 K w,b i K -1 b i ,w = 1,thereby completing the proof of Theorem 37.

Note that K can also be seen as an operator mapping

More precisely, since every white vertex w has degree 3, denoting by b 1 , b 2 , b 3 its neighbors as in Figure 15, this relation can be rewritten as,

5.2 Functions in the kernel of the Kasteleyn operator K

We now define the function f in the kernel of the Kasteleyn operator K. Note that it generalizes to the elliptic case the function f introduced by Kenyon in [START_REF] Kenyon | The Laplacian and Dirac operators on critical planar graphs[END_REF] when the bipartite graph is G Q .

Rhombus vectors. In order to define the function f , we need to assign rhombus vectors to edges of the graph G Q . Since the graph G Q is isoradial, it also has a diamond graph (G Q ) ; note that rhombi of (G Q ) are obtained by cutting those of G in four identical rhombi, see Figure 16.

Consider an edge bw of G Q . Then we let 1 2 e iα and 1 2 e iβ be the two rhombus vectors of (G Q ) of the edge bw, where 1 2 e iα is on the right of the oriented edge (b, w). Some examples are given in Figure 16. Remark 33. The angles α and β above are defined so that β -α ∈ (0, 2K). Definition 5.1. For every edge bw of G Q and every u ∈ C, define

if bw is parallel to an edge e of G, -ik nc( u-α

2 ) nc( u-β 2 ) if bw is parallel to the dual of an edge e of G,

if bw is an external edge and w belongs to a rhombus of G having half-angle θ,

The function f :

. . , b n , w n = w be a path from b to w, then:

Remark 34. As the function g of Definition 3.2, the function f is meromorphic and biperiodic:

This comes from [START_REF] Nienhuis | Critical behavior of two-dimensional spin models and charge asymmetry in the Coulomb gas[END_REF] and from the addition formulas of cn and cd by 2K, see Table 2. We therefore also restrict the domain of definition to T(k) = C/(4KZ + 4iK Z).

Let us first check that this is true for rectangles, using the notation and computations of Example 5.1. Recalling that

We deduce that the function f is well defined around rectangles. We now turn to faces which are not rectangles, and do some preliminary computations. Thanks to Example 5.1 again, we have (see Table 2 for the various simplifications involving Jacobi elliptic functions in ( 51) and ( 52))

We also have

We have expressed the product Then, by [START_REF] Onsager | Crystal statistics. I. A two-dimensional model with an order-disorder transition[END_REF] we have

) ns(

see Figure 17 (left), implying that

) ns(

where Γ b,w is a vertical contour directed upwards on T(k), crossing the real axis outside of the sector of size 2K containing all the poles of f (b,w) .

Then K -1 is an inverse operator of K. For k = 0, it is the only inverse with bounded coefficients.

The quantity K -1 b,w in (58) can alternatively be expressed as

where H is related to Jacobi's zeta function and is defined in (66)-( 67), C b,w is a trivial contour on the torus, not crossing Γ b,w and containing in its interior all the poles of f (b,w) and the pole of H.

Proof. To show that the two expressions ( 58) and ( 59) indeed coincide, we use the same argument as in the proof of Theorem 11.

The structure of the proof of Theorem 37 is analogous to that of [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]Theorem 1]. Instead of using the form (58) as in the proof of Theorem 11, we use the alternative expression (59) of K -1 b,w . Indeed, from the computations done below to prove that (KK -1 ) w,w = 1, one can extract as a by-product the explicit probability of a given edge to be present in the random dimer configuration in the corresponding dimer model. b 1 ,w , K -1 b 2 ,w and K -1 b 3 ,w can be deformed into a common contour C without crossing any pole. Therefore, the entry (KK -1 ) w,w can be written as:

We now need to compute the entry (KK -1 ) w,w . This is done explicitly, via the residue theorem. In addition to the simple pole at u = 2iK coming from the function H (with residue 2K /π, see Lemma 42), there are other (simple) poles located at the zeros of the functions in the denominator of f , i.e., when the argument of the functions cd and cn is equal to K.

We shall successively compute

b 1 ,w is obtained from (the minus signs in the numerators in the right-hand side below come from the expansion of cd around K, see [1, Table 16.7])

Theorem 38. Let G be a quasicrystalline isoradial graph. When the distance |b -w| → ∞, we have

where θ is the rhombus-angle of the rhombus to which w belongs, u 0 is the unique u ∈ (-K, K) such that χ (u) = 0, and χ(u) < 0.

Theorem 38 should be compared to its genus 0 counterpart, i.e., Theorem 4.3 of [START_REF] Kenyon | The Laplacian and Dirac operators on critical planar graphs[END_REF] where polynomial decrease of coefficients of the inverse Kasteleyn operator is proved.

Remark 39. Contrary to Theorem 13, where we have shown that the constant in front of the exponential function is always positive (see Remark 14), we have less control on the constant in Theorem 38. First, it can have a phase, due to the terms e iθ and e -i 2 (αn-α 1 ) . The main point is that the quantity h(u 0 ± 2iK ), which is real, can be positive, negative and even 0. This follows from (65).

Note that the proof of Theorem 38 is similar to that of Theorem 13; therefore, we omit the details. We only need to give an expression for the function h appearing in the statement of Theorem 38. First of all, we prove that the function f (b 1 ,bn) looks very much like the exponential function e (a 1 ,a n-1 ) (we use the previous notation).

Lemma 40. The following formula holds:

where

It is important to note that taken independently, the factors e -i 2 (αn-α 1 ) and g(u) are not well defined: if α n (or α 1 ) is replaced by α n + 4K, these terms should be replaced by their opposite. However, the product e -i 2 (αn-α 1 ) g(u) is well defined, which suffices for our purpose.

Proof. With the previous notation we write f (b,bn

There are two cases for the computation of f (b j ,b j+1 ) , according to whether a j is a primal or a dual vertex: the identities ( 53) and (54) yield

)ds(

) if a j is primal, sn(

The term ±2K in (63) comes from the fact that the orientation of the rhombus vectors in ( 53)- [START_REF] Wannier | The statistical problem in cooperative phenomena[END_REF] and in (60) are reversed. The quantity (63) does not depend on the value of this sign (this is a consequence of the addition formulas by ±K for the sn and sd functions, see Table 2 in the appendix).

Due to the fact that the value of f (b j ,b j+1 ) depends on the type (primal or dual) of the vertex a j , there are four cases for the computation of f (b,bn) , according to the types of a 1 and a n-1 .

We write down the computations in the particular case where both a 1 and a n-1 are dual, the other cases would follow in a very similar manner. We have (with all signs ± = +):

2 ) e (a 1 ,a n-1 ) .

The third equality above uses addition formulas of Table 2. In the last line, we have applied the definition (9) of the exponential function. We also implicitly used the fact that n is even (because both a 1 and a n-1 are dual vertices). The first case of Equation ( 62), and thus of Lemma 40, is proved.

We now give a formula for f (b,w) for general vertices b and w.

Lemma 41. The following formula holds:

where θ is the rhombus-angle of the rhombus to which w belongs, and

Proof. We have f (b,w) = f (b 1 ,bn) f (bn,w) . The definition (50) of the function f provides

(The choice of the sign ± above comes from Remark 33.) Using that cd( u-αn∓2K

2

2 ) together with Lemma 40 ends the proof of Lemma 41.

A.2 Jacobi's epsilon, zeta and related functions

The explicit expressions of the inverse operators of Theorems 11 and 37 involve the function H defined as follows. In the remarks following Theorem 11, we also use the function V mentioned below.

where E is Jacobi's epsilon function: [1, 16.26.3]. For k 2 < 0, let

Properties of these functions are presented below.

Lemma 42. The functions H and V admit jumps in the horizontal and vertical directions, respectively:

In the fundamental rectangle [0, 4K] + [0, 4i K ], the function H (resp. V ) has a simple pole, at 2i K , with residue 2 K π (resp. 2iK π ). Moreover,

The following addition formulas hold:

Proof. We first prove the lemma in the case 0 < k 2 < 1. All properties concerning H are proved in [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]Lemmas 44 and 45]. The statements for V follow similarly.

A slightly different proof would consist in using a reformulation of H and V in terms of Jacobi's zeta function Z:

(Equation ( 71) is a consequence of (66) and [1, 17.4.28 and 17.3.13].) Then we could use the numerous properties satisfied by Z (see in particular the addition formula [1, 17.4.35]) to derive Lemma 42.

In the case k 2 < 0, we use the definition (67) of H and V , allowing to transfer all properties from the case k 2 > 0. In doing so we have to: consider K in (68) because in the case k 2 < 0 the quarter-period K is non-real, see [START_REF] Mccoy | The two-dimensional Ising model[END_REF], and use the transformation of the sn function into the sd one under the dual transformation, see [1, 16.10].

The Laplacian operator (6) uses the function A, which satisfies the following properties.

Lemma 43 (Lemma 44 in [START_REF] Boutillier | The Z-invariant massive Laplacian on isoradial graphs[END_REF]). The function A(•|k) is odd and satisfies the following identities:

B Some explicit integral computations

We gather here computations of some contour integrals appearing in the expression of the Kasteleyn operator, in the Fisher case of Section 3.

B.1 An important contour integral

The following result has been used when proving Theorem 11.

Lemma 44. One has

where Γ is a vertical contour on T(k) winding once vertically on the torus and crossing the horizontal axis in the interval [α, α + 2K] = {x : α x α + 2K} and f(u) = nc( u-α 2 ).

On the rectangle, the contour Γ is supposed to cross the horizontal axis inside of the interval [α, α + 2K]. If the vertical contour crosses the horizontal axis in the other interval [α + 2K, α(+4K)], the integral is equal to + 1 k , as it corresponds to changing α into α + 2K and nc( u-(α+4K) 2 ) = -nc( u-α 2 ). Note further that Lemma 44 is independent of the choice of the angle α mod 4K, and the integrand f(u + 2K)f(u) is.

We start from the expression of K -1 given in [START_REF] Duffin | Potential theory on a rhombic lattice[END_REF] of Theorem 11 (note that by [START_REF] Dubédat | Exact bosonization of the Ising model[END_REF], C v j ,v = 0):

Using the harmonicity property (16) of g (z,•) (u) enables us to rewrite

By definition of the function f v j , see ( 13), we thus have

Recalling that the orientation of the triangle (w j , v j , w j+1 ) is admissible, we moreover have K v j ,w j K v j ,w j+1 = -K w j ,w j+1 , and since

With (75) the quantity K v ,v j K -1 v j ,v is a sum of four terms. The first two ones are computed thanks to Lemma 44: with the choice of contour C v j ,v , we have

To compute the last two terms in (75), namely

recall that by ( 11)

), so that

) nc( u-α j 2 ).

We therefore focus on the term

C Proof of Lemma 7

This section consists in the proof of Lemma 7 stating that the angles (α j (x)) x∈V, j∈{1,...,d(x)} defined in Equations ( 11) and ( 12) are indeed well defined mod 4π.

Proof. We first need to check that angles around a rhombus corresponding to two neighboring decorations are well defined, see Figure 20. We want to prove that the following is equal to 0 mod 4π:

(α +1 -α ) -(α j+1 -α j ) + (α -α j ) -(α +1 -α j+1 ).

By definition of the angles within a decoration we have, mod 4π, (α +1 -α ) -(α j+1 -α j ) = 0 if co(w , w +1 ) + co(w j , w j+1 ) is even, 2π if co(w , w +1 ) + co(w j , w j+1 ) is odd.

By definition of the angles in neighboring decorations we have, mod 4π, (α -α j ) -(α +1 -α j+1 ) = 0 if co(w j+1 , v j , w j ) + co(w +1 , v , w ) is even, 2π if co(w j+1 , v j , w j ) + co(w +1 , v , w ) is odd. This implies that, mod 4π, we have:

if co(w j+1 , v j , w j , w j+1 ) + co(w +1 , v , w , w +1 ) is even, 2π if co(w j+1 , v j , w j , w j+1 ) + co(w +1 , v , w , w +1 ) is odd.

But since the orientation of the graph is admissible, we have that co(w j+1 , v j , w j , w j+1 ) and co(w +1 , v , w , w +1 ) are odd, implying that the sum is even, thus concluding the proof for angles around a rhombus.