Non-integrable dimers: Universal fluctuations of tilted height profiles

We study a class of close-packed dimer models on the square lattice, in the presence of small but extensive perturbations that make them non-determinantal. Examples include the 6-vertex model close to the free-fermion point, and the dimer model with plaquette interaction previously analyzed in \cite{A,AL,GMT17a,GMT17b}. By tuning the edge weights, we can impose a non-zero average tilt for the height function, so that the considered models are in general not symmetric under discrete rotations and reflections. In the determinantal case, height fluctuations in the massless (or `liquid') phase scale to a Gaussian log-correlated field and their amplitude is a universal constant, independent of the tilt. When the perturbation strength $\lambda$ is sufficiently small we prove, by fermionic constructive Renormalization Group methods, that log-correlations survive, with amplitude $A$ that, generically, depends non-trivially and non-universally on $\lambda$ and on the tilt. On the other hand, $A$ satisfies a universal scaling relation (`Haldane' or `Kadanoff' relation), saying that it equals the anomalous exponent of the dimer-dimer correlation.


Introduction
The question of universality, that is the independence of the critical properties of macroscopic systems from the microscopic details of the underlying model Hamiltonian, is a central issue in statistical physics, whose mathematical understanding is largely incomplete. A convenient framework where it can be studied is that of planar dimer models, which exhibit a rich critical behavior: algebraic decay of correlations, conformal invariance, ... The dimer model on a bipartite planar lattice is integrable and, more precisely, determinantal (also said 'free fermionic'): its correlation functions are given by suitable minors of the so-called inverse Kasteleyn matrix [29]. The model is parametrized by edge weights t and has a non-trivial phase diagram. By varying t, one can impose an average non-zero tilt ρ for the height field. A central object of the dimer model is the so-called characteristic polynomial P (z, w), where z, w are complex variables. For instance, the infinite-volume free energy is given by an integral of log |P (z, w)| over the torus T = {|z| = |w| = 1}. Also, the large-distance decay of correlations is dictated by the so-called spectral curve, i.e. the algebraic curve C(P ) = {(z, w) ∈ C 2 : P (z, w) = 0}. When the edge weights are such that the spectral curve intersects T transversally one is in the "liquid" or "massless" phase, where the two-point dimer-dimer correlation of the model decays like the inverse distance squared. Correspondingly the height field scales to a Gaussian Free Field (GFF) and the variance grows like the logarithm of the distance times 1/π 2 . Remarkably, this pre-factor is independent of the weights t and of the specific choice of the bipartite periodic planar lattice. This is related [31] to the fact that C(P ) is a so-called Harnack curve. Summarizing, in the massless phase the scaling limit of height fluctuations of the dimer model is universal, in a very strong sense: the limit is always Gaussian, with logarithmic growth of the variance; moreover, the pre-factor in front of the logarithm in the variance is independent of the details of the underlying microscopic structure (edge weights and lattice).
The previous results heavily rely on the determinantal structure of the model, but universality is believed to hold much more generally. Motivated by this, we consider weak, translation-invariant, perturbations of the dimer model (for simplicity, we restrict to the square lattice). Generically, as soon as we switch on the perturbation, the determinantal structure provided by Kasteleyn's theory breaks down. Two particular examples of perturbed, non-determinantal, models that we consider are: the 6-vertex model with general weights a 1 , . . . , a 6 , in the disordered phase, close to, but not exactly at, the free-fermion point; and the dimer model with plaquette interaction, originally introduced in [26] and recently reconsidered in [2,3,35] in the context of quantum dimer models. There is a basic difference between two such cases: the 6-vertex model, even if non-determinantal, is still solvable via Bethe Ansatz (BA), see [4] and reference therein (note that the BA solution is not as explicit as the Kasteleyn solution of standard dimers: only a few thermodynamic functions can be explicitly computed). On the other hand, dimers with plaquette interaction are believed not to be solvable, i.e., not even the basic thermodynamic functions admit an explicit representation. From the exact solution, one finds that some of the critical exponents of the 6-vertex model depend continuously on the vertex weights 1 ; they differ, in general, from those of the standard dimer model. On the other hand, the existence of non-trivial critical exponents in the dimer model with plaquette interaction, as well as in other planar models in the same 'universality class' (such as coupled Ising models, Ashkin-Teller and 8-vertex models) can be proved by constructive Renormalization Group (RG) methods [6,8,22,32], which allow one to express them as convergent power series in the interaction strength.
In this setting, it is natural to ask whether the height fluctuations are still described by a GFF at large scales and, in case, whether the pre-factor in front of the logarithm still displays some universal features. The very fact that the critical exponents depend non-trivially on the interaction strength suggest that universality cannot then be true in the naive, strong, sense that 'large-scale properties are independent of the microscopic details of the model': in fact, in this context, a weaker form of universality is expected, in the form of a number of scaling relations, originally proposed by Kadanoff [28], which allow one to determine all the critical exponents of the critical 1 More precisely, the limit of the critical exponents of the 8-vertex model as the additional vertex weights a7 = a8 tend to zero have a non-trivial continuous dependence on the remaining vertex weights a1, . . . , a6, see [4,Eqs.10.12.23 and 10.12.27]. theory in terms of just one of them; this form of universality is often referred to as 'weak universality', see e.g. [4,Section 10.12]. Support for the Kadanoff scaling relations comes from the so-called bosonization picture, see e.g. [24] for a basic introduction. Only some of these universality relation have been rigorously proven [6,8]; an example is the identity X c X e = 1 [28,Eq.(13b)], relating the "crossover exponent" X c and "energy exponent" X e , see [6,Eq.(1.10)]. The proof in [6] covers both solvable and non-solvable models, but only works for scaling relations involving the critical exponents of the "local observables", i.e., those that admit a representation in terms of a local fermionic operator. Other scaling relations, involving the critical exponents of non-local observables (e.g. monomer-monomer correlations in dimer models, or spin-spin correlations in the Ashkin-Teller model) remained elusive for many years. In particular, the relation X p = X e /4 [28,Eq.(13a)], relating the energy exponent X e to the "polarization exponent" X p in the AT model, remained unproven at a rigorous level.
In this paper, we prove the stability of the Gaussian nature of the height fluctuations for non-integrable perturbations of the dimer model, with logarithmic growth of the variance in the whole liquid region. The pre-factor A in front of the logarithm depends, in general, non-trivially on the strength of the perturbation (see Remark 4 below) and on the dimer weights, so it is not universal in a naive, strong, sense. The non-trivial dependence of A on the interface tilt has been also verified numerically for the 6-vertex model [25]. Nevertheless, A satisfies a scaling relation, that connects it with the critical exponent of the dimer-dimer correlations.
Main Theorem. In a weakly perturbed dimer model with perturbation of strength λ, the variance of the height difference between two faraway points grows like the logarithm of the distance, with a pre-factor A/π 2 , where A = 1 + O(λ) is an analytic function of λ and of the dimer weights. Moreover, the pre-factor satisfies the scaling relation where 2ν is the anomalous decay exponent of the dimer-dimer correlation. Higher cumulants of the height difference between two points are bounded uniformly in their distance, that is, the fluctuations of the height difference are asymptotically Gaussian.
For a more precise statement, see Theorem 2 and the remarks and comments that follow it. Note that in the un-perturbed case, λ = 0, the dimerdimer correlation decays at large distances like (dist.) −2 in the whole liquid phase, i.e., its decay exponent is equal to 2 (so that ν = 1), irrespective of the specific choice of the dimer weights. In this case, of course, our result reduces to the one of [31], A = 1. Note also that our result covers both integrable models, such as 6-vertex, and non-integrable ones, in the spirit of the universality picture.
Scaling relations involving exponents and amplitudes were conjectured by Haldane [27] and proved by Benfatto and Mastropietro [11,12] in the context of quantum one-dimensional models. Even if formulated in different notations, the scaling relation (1.1) is strictly related to one of those proposed by Kadanoff, in particular to the above-mentioned, elusive, identity X p = X e /4 [28,Eq.(13a)]. In fact, there is a duality (called 'discrete bosonization' in [17]) between the 6-vertex model, which is part of the class of perturbed dimer models considered in this paper, and the AT model; the duality implies non-trivial identities between the correlations of 6-vertex model and those of AT, see [17,Section 2.6]. In particular, the two-point correlation of the polarization operator in AT equals the 'electric correlator' e iπ(hx−hy) 6V of the 6-vertex model, see [17,Section 2.6] 2 , while the energy critical exponent of AT equals the anomalous decay exponent of the arrow-arrow correlations of 6-vertex 3 . Given these identities, (1.1) implies that X p = X E /4 [28,Eq.(13a)], provided that e iπ(hx−hy) 6V ∼ e − π 2 2 (hx−hy) 2 6V ∼ e − A 2 log |x−y| (1.2) at large distances, as suggested by the asymptotic Gaussian behavior of the height difference 4 .
To prove our results, we start by periodizing the non-integrable dimer model on the toroidal graph of size L. Then we map it into a system of interacting two-dimensional lattice fermions, by rewriting its moment generating function as an integral over Grassmann variables, with non-quadratic action. At this point, we apply tools from the so-called constructive fermionic RG to control the L → ∞ limit of the correlation functions. In particular, we need a very sharp asymptotic description of the large-distance behavior of the dimer-dimer correlation function (cf. Theorem 1). The large-scale logarithmic behavior of height correlations, as well as the validity of the 'Haldane' scaling relation (1.1), rely on non-trivial identities (cf. (2.43)) between the coefficients appearing in the large-distance asymptotics of the dimer-dimer correlation function. In turn, (2.43) is the result of so-called Ward identities, i.e. exact relations between the correlation functions of the interacting lattice fermionic model, which the dimer model maps into.
The analogs of Theorems 1 and 2 have been proven in our previous works [23,24] for the specific case of plaquette interaction and uniform edge weights t ≡ 1. In this case, the average tilt of the height field is just ρ = 0 and the model has all the discrete symmetries of the lattice Z 2 . The extension to the general case, achieved here, is non-trivial: the loss of discrete rotation and reflection symmetries results, in the RG language, in the emergence of four new running coupling constants (two "Fermi velocities" and two "Fermi points"), whose flow, along the multi-scale integration procedure, has to be controlled via the choice of suitable counter-terms. Another consequence of the loss of rotation and reflection symmetry is that the cancellation at the basis of the logarithmic growth of the variance does not follow simply from 2 Here hx is the height function of the 6-vertex model at face x and · 6V is the corresponding statistical average; the factor π at the exponent depends on our definition of height function, which differs by a multiplicative factor 2π from that of [17]. 3 In the dimer formulation of 6-vertex, the arrow-arrow correlations translate into the dimer-dimer correlations. 4 As discussed in [23,Remark 2], our method allows us to compute the average of exp{iπ(hx − hy)} only after coarse-graining the height difference at exponent against a smooth test function.
the basic symmetries of the model, as it was the case in [23,24]: the proof of the key identity, (2.43), now requires the use of a lattice Ward Identity for the dimer model, in combination with an emergent Ward Identity for an effective continuum model, which plays the role of 'infrared fixed point' of the RG flow. Quite surprisingly, the loss of rotation and reflection symmetry plays a role also in the technical control of the thermodynamic limit of correlations: in [23,24], in order to simplify the analysis of the finite-size corrections to the critical correlation functions, we first studied a modified, slightly massive model of mass m > 0 (the modification consisted in adding a modulation of size m on the horizontal dimer weights; in the tilt-less case, this was enough to guarantee that the modified correlations decayed exponentially with rate m), and then we took the massless limit m → 0 after the thermodynamic limit. However, this strategy fails for general dimer weights: in this case, neither a modulation of the dimer weights nor other simple modifications of the model produce a mass; therefore, in the present paper, we directly derive quantitative estimates on the corrections to the thermodynamic limit of the massless correlations, by a careful control of the finite-size effects in the multi-scale procedure.
1.1. Related works. Let us conclude this introduction by mentioning some recent related works. While most literature on dimer models focuses on the determinantal case, there have been recently various attempts to go beyond the exactly solvable situation [34]. As far as "limit shape phenomena" (i.e. laws of large numbers for the height profile) for non-solvable random interface models are concerned, let us mention for instance [16,33,14]. Closer in spirit to our results is [15], that provides a central limit theorem for height fluctuations of ∇φ-interface models with continuous heights and strictly convex potential. This work uses the Helffer-Sjöstrand formula, that is not available for discrete-height model like the dimer model. Let us mentional also [1], that obtains convergence to the GFF (with interactiondependent amplitude) for a dimer model with a special non-local interaction that makes it integrable, although not determinantal. Finally, a very interesting recent development is [13]: while in this work the convergence to the GFF is proven only for the non-interacting dimer model, the method of proof, that goes through Temperley's bijection and Wilson's algorithm rather than via Kasteleyn's theory, might prove robust enough to allow for extensions to some non-determinantal situations.

1.2.
Organization of the article. The rest of this work is organized as follows. The dimer model is defined in Section 2. There, we recall the large-scale behavior of the integrable model and we state our results for the non-integrable one. In Section 3 we give the Grassmann representation of the interacting dimer model and its lattice Ward identities. In Section 4 we recall the continuum reference model that plays the role of infrared fixed point of interacting dimers. Theorems 1-2 are proven in Section 5, conditionally on technical results, based on the multi-scale expansion, whose proofs are postponed to Section 6.

Model and main results
2.1. Dimers and height function. A dimer covering, or perfect matching, of a graph Γ is a subset of edges that covers every vertex exactly once. The set of dimer coverings of Γ is denoted Ω Γ . We color the vertices of the bipartite graph Z 2 black and white so that neighboring vertices have different colors. A white vertex is assigned the same coordinates x = (x 1 , x 2 ) as the black vertex just at its left. The choice of coordinates is such that the vector e 1 is the one of length √ 2 and angle −π/4 w.r.t the horizontal axis, while e 2 is the one of length √ 2 and angle +π/4. The finite graph T L denotes Z 2 periodized (with period L) in both directions e 1 , e 2 . See Fig.  1. For simplicity we assume that L is even. Black/white sites are therefore indexed by coordinates x ∈ Λ = {(x 1 , x 2 ), 1 ≤ x i ≤ L}. An edge e = (b, w) of T L is said to be of type r ∈ {1, 2, 3, 4} if its white endpoint w is to the right, above, to the left or below the black endpoint b. If e = (b, w) is an edge of type r and x(b) is the coordinate of b then If Γ is planar and bipartite, the height function allows us to interpret a dimer covering as a two-dimensional discrete surface. Let us recall the standard definition of height function for the infinite lattice Z 2 . Given M ∈ Ω Z 2 , the height function h(·) := h M (·) is defined on the dual lattice (Z * ) 2 , i.e. on the faces η of Z 2 . We set h(η 0 ) := 0 at a given reference face η 0 , and we let its gradients be given by where η, η are any two faces, 1 e denotes the dimer occupancy, i.e., the indicator function that e is occupied by a dimer in M , while C η→η is any nearest-neighbor path on the dual lattice (Z * ) 2 from η to η (the right side of (2.2) is independent of the choice of C η→η ). The sum runs over the edges crossed by the path and σ e = +1/ − 1 depending on whether the oriented path C η→η crosses e with the white site on the right/left.

2.2.
Definition of the model. We define here both the non-interacting dimer model [29] and the interacting one. Both are probability measures on Ω L := Ω T L , denoted P L,t and P L,λ,t respectively, where λ ∈ R is the interaction strength and t are the edge weights. For lightness of notation, the index t will be dropped.
2.2.1. The non-interacting dimer model. We assign a positive weight to each edge. More precisely, an edge of type r ∈ {1, 2, 3, 4} is given a weight t r > 0. Then, the weight of a configuration M ∈ Ω L is with N i (M ) the number of dimers on edges of type i in configuration M . Since the total number of dimers is constant, we can rescale all weights by a common factor and we will set t 4 ≡ 1 from now on. It is known that the free energy per site has a limit as L → ∞ (the infinite volume free energy): Note that The "characteristic polynomial" mentioned in the introduction is P (z, w) := µ(−i log z, −i log w). Also, the measure P L itself has a limit P as L → ∞, in the sense that the probability of any local event converges. The non-interacting model is integrable, and both the measure P L and its limit P admit a determinantal representation, recalled in Section 3.1.
In the special case where t 1 = t 3 =: t and t 2 = 1, i.e. assigning weight t to horizontal edges and 1 to vertical ones, one recovers the model originally solved by Kasteleyn [29]. For general weights t 1 , t 2 , t 3 , the model is equivalent to Kasteleyn's model with different weights for horizontal and vertical edges, and a non-zero average slope ρ = ρ(t 1 , t 2 , t 3 ) ∈ R 2 for the height function, i.e., where E denotes the average with respect to P. In fact, the weights t i are chemical potentials by which one can fix the densities of the four types of edges. Then, the slope ρ is obtained as a function of the four densities using the definition (2.2) of height function.
Another special case is obtained letting e.g. t 3 → 0: then, the model reduces to the closed-packed dimer model on the hexagonal graph with weights 1, t 1 , t 2 for the three types of edges.
Note that the condition µ(k) = 0 gives that determines the intersections of two circles in the complex plane. We will make the following important assumption: The parameters t are such that µ(·) has two distinct simple zeros, that we call p + and p − , on [−π, π] 2 (i.e. the two circles intersect transversally). In view of (2.7), one has p + + p − = (π, π).
Under Assumption 1, it is known [31] that the infinite-volume measure has power-law decaying correlations (in the language of [31], the dimer model is said to be in a "liquid phase"). With the nomenclature of condensed matter theory, the zeros p ± are called "Fermi points".
2.3. The interacting dimer model, and relation to the 6-vertex model. In order to study the effect of the breaking of integrability we introduce interacting dimer measures of the following form: (2.11) and the interaction potential W L is given as where f is some fixed local function of the dimer configuration and τ x M denotes the configuration M translated by x 1 e 1 + x 2 e 2 . We do not require f (·) to be symmetric under reflections or rotation by π/4. Let us mention two interesting particular examples of interaction W L (M ). The first one is the plaquette interaction that was considered in our works [23,24] and previously in the theoretical physics literature [2] in the context of quantum dimer models. Namely, In this case, the interaction W L (M ) in (2.13) is modified in that the sum runs only over one of the two sub-lattices of T * L (the subset of faces with black top-right vertex). Then, it is known that this interacting dimer model is equivalent to the 6-vertex model [5,18,19]. Recall that configurations of the 6-vertex model are assignments of orientations (arrows) to the edges of Z 2 such that at each vertex there are two incoming and two outgoing arrows. There are 6 possible arrow configurations at any vertex, each being assigned a positive weight a 1 , . . . , a 6 (see Fig. 3) and the weight of a configuration is the product of the weights over all vertices. By multiplying all weights by a common factor, one can reduce e.g. to a 3 = 1. Moreover, on the torus, the number of vertices of type 5 equals the number of vertices of type 6, so one can set without loss of generality a 5 = 1. One is left with four positive weights a 1 , a 2 , a 4 , a 6 and the model can be mapped to the interacting dimer model with weights t 1 , t 2 , t 3 , interaction (2.15) and interaction parameter λ such that (2.16) More precisely, as in Fig. 4, the dimer model lives on a square grid rotated by 45 degrees w.r.t. the lattice of the 6-vertex model. The mapping is obtained by associating to the arrow configuration at a vertex x of G 6v a dimer configuration at the even face of G d containing x, as in Fig. 5. The map is one-to-many because arrow configurations of type 6 are mapped to with dotted edges, while the dimer model lives on the square grid G d with full edges. Faces of G d containing a vertex of G 6v are called "even faces" and the others "odd faces".
or Figure 5. The local arrow-to-dimer mapping two possible dimer configurations. However, it is easily checked that the partition functions of the two models are equal provided the parameters are identified as in (2.16). Moreover, the height function of the dimer model, restricted to odd faces of G d , equals (up to a global prefactor) the canonical height function of the 6-vertex model [37]. The 6-vertex model is known to be free-fermionic (i.e. determinantal) if and only if It is immediately checked that this condition is equivalent to λ = 0 for the interacting dimer model.

2.4.
Non-interacting model: dimer-dimer correlations and logarithmic height fluctuations. It is known [31] that, under the infinitevolume measure P of the non-interacting model, dimer-dimer correlations decay like the inverse distance squared and the height field behaves on large scales like a massless Gaussian field. We briefly recall the basic facts here, since they serve to motivate our main result for the interacting dimer model.
For ω = ±, we let where p ± are the two zeros of µ(·), as in Assumption 1. (The complex numbers α ω , β ω are called "Fermi velocities" in the jargon of condensed matter.) Define also (2.20) Remark 2. Under Assumption 1 on the weights t, the complex numbers α ω and β ω are not colinear, as elements of the complex plane [31], i.e. α ω /β ω is not real. Therefore, φ ω is a bijection from R 2 to the complex plane. More precisely, one has that In fact, parametrize the weights t 1 , t 2 , t 3 as For B 1 , B 2 = 0 it is immediately checked that p + = (0, 0), p − = (π, π) and that (2.21) holds. On the other hand, once t is fixed, it is known [31] that the set B t of values of B = (B 1 , B 2 ) for which Assumption 1 holds is a connected subset of R 2 on which Im(β + /α + ) vanishes nowhere, and it is therefore everywhere positive.
Remark 3. Note that the prefactor 1/π 2 is independent of the weights t.
In [31], such universality is related to the fact that the spectral curve, i.e. the algebraic curve defined by the zeros on C 2 of the polynomial P (z, w) := µ(−i log z, −i log w), is a so-called Harnack curve.
It is useful to recall the key points of the proof of (2.22) in order to understand the main new features posed by the presence of the interaction. From the definition of height function, This correlation function has an exact expression involving the inverse Kasteleyn matrix of the infinite lattice; at large distances, it can be expressed as where: • the edge e (resp. e ) is of type r = r(e) (resp. r = r(e )) and the coordinate of its black endpoint is x = x(e) (resp. x = x(e )); • K ω,r = K r e −ip ω ·vr (see (2.1) for the definition of v i ) with Note that, since p ± are distinct by assumption, the complex exponential in the definition of B r,r is genuinely oscillating. For simplicity, assume that the paths C η 1 →η 2 , C η 3 →η 4 are a concatenation of elementary steps in direction ± e 1 and ± e 2 , connecting faces of the same parity: e.g., assume that an elementary step s(x, 1) in direction + e 1 'centered at x' consists in crossing the two bonds (•, •) = (x, x + v 3 ) and (•, •) = (x, x + v 4 ) with the white vertex on the right, while an elementary step s(x, 2) in direction + e 2 centered at x consists in crossing the two bonds (•, •) = (x, x) and (•, •) = (x − v 4 , x) with the white vertex on the right. A simple but crucial observation is that where ∆ j φ ω denotes the discrete gradient in direction e j of the affine function φ ω defined in (2.20). By inserting (2.27) in (2.26) one can see that the contribution from R r,r is subdominant, and the same for B r,r due to the oscillating complex exponential. As for the dominant contribution to (2.26), coming from the term A r,r , one sees using (2.29) that it approximately equals the integral in the complex plane whose explicit evaluation gives the main term in the r.h.s. of (2.22).
2.5. The interacting case: main results. In the presence of the interaction, λ = 0, Kasteleyn theory is not valid anymore, so that one cannot rely on an explicit computation of the dimer correlations to check the validity of the asymptotic Gaussian behavior of the height function. However, dimer correlations can be written as a renormalized expansion based on multiscale analysis. From now on, we will assume that the interaction is small: and all claims above hold if ε is small enough (uniformly in L). Our first result is: Theorem 1. Given a local function g of the dimer configuration, the limit exists. The infinite-volume dimer-dimer correlations are given by where: • r = r(e) is the type of the edge e, x = x(e) is the coordinate of the black site of e, and similarly for r , x ; (2.38) these are all analytic functions of λ and satisfy the symmetries Finally,R r,r (x, x ) = O(|x − x | −5/2 ) (the exponent 5/2 could be replaced by any δ < 3 provided λ is small enough).
[A warning on notation: given a quantity (such as α ω , φ ω ) referring to the non-interacting model, the corresponding λ-dependent quantity for the interacting model will be distinguished by a bar, such asᾱ ω , etc. On the other hand, we denote by z * the complex conjugate of a number z.] Note that the interaction modifies the decay rate of the correlation, producing a non-trivial ('anomalous') critical exponent ν. The analytic functions appearing in (2.37) are expressed as convergent power series but, due to the complexity of the expansion, the coefficients can be explicitly evaluated only at the lowest orders. This makes impossible to verify directly the validity of relations like (2.29), which were essential for the proof of large-scale Gaussian behavior of the height field in the non-interacting case. However, we can prove non-perturbatively that the parameters appearing in (2.34) are not independent, but related by exact relations, which are the central result of the present work:

43)
where ν = ν(λ) is the same as the critical exponent in Theorem 1. Here, s(x, j) is the elementary step in direction + e j centered at x, thought of as a collection of two bonds, as defined before (2.29). As a consequence, (the exponent 1/2 could be replaced by any δ < 1 provided λ is small enough; as in (2.22), when η i = η j ,φ + (η i ) −φ + (η j ) has to be read as 1).
Note that the result contains two non-trivial pieces of information: first, the sum of σ eKω,r(e) along a step in direction e i is proportional to the discrete gradient ofφ ω in the same direction; second, the coefficient of proportionality is related in an elementary way to the critical exponent ν that appears in (2.36). The latter relation immediately implies the identity (cf. (2.44)) between height fluctuation amplitude and critical exponent ν and is a form of universality.

Remark 4.
Recall that for the non-interacting model ν = 1, in particular it is independent of the weights t i . This is not true anymore for the interacting model. Indeed, an explicit calculation of ν at first order in λ for the model with plaquette interaction shows a non-trivial dependence both on λ and on the weights.
Theorem 2 follows from a combination of exact relations among correlation functions of the interacting dimer model ("lattice Ward identities") together with chiral gauge symmetry emerging in the continuum scaling limit; it is remarkable that such a symmetry, valid only in the continuum limit, implies nevertheless exact relations for the coefficients of the lattice theory.
Remark 5. The analog of Theorem 1 has been proven in [23], [24] in the special case t 1 = t 2 = t 3 = 1 and with plaquette interaction as in (2.14), which has the same discrete symmetries as the lattice. In that case, for symmetry reasons one obtains automatically that the ratiosK ω,r Kω,r are independent of r, ω and thatᾱ ω αω =β ω βω , the ratios being again ω-independent. Then, the analog of Theorem 2 is trivial in that case.
Let us add also that, in the works [23,24], the existence of the L → ∞ limit of the measure P L,λ itself was not proven: instead, we modified the measure P L,λ by an infra-red cut-off m > 0 (mass) and then we took the limit where first L → ∞ and then m → 0. We explain in Section 6 how the need of the cut-off m can be bypassed.
To upgrade Theorem 2 into a statement of convergence of the height field to a Gaussian Free Field with covariance − ν 2π 2 log(φ + (x) −φ + (y)), one needs to complement (2.44) with the statement that higher cumulants are negligible, i.e. that, for n > 2 and some θ > 0, In turn, this requires an analog of (2.34) for multi-dimer correlation functions. This can be done following the ideas of Sections 5 and 6 below but, in order to keep this work within reasonable length, we decided not to develop this point. The interested reader may look at [23,Theorem 3 and Sec. 7], where the precise statements on multi-dimer correlations and on the convergence to the GFF are given in detail for the model with edge weights t ≡ 1 and interaction (2.13).

Grassmann integral representation
3.1. Kasteleyn theory. For the statements of this section and more details on Kasteleyn theory, we refer the reader for instance to [30,31].
The partition function and the correlations of the non-interacting model (2.3) can be explicitly computed in determinantal form, via the so-called Kasteleyn matrix K. This is a square matrix of size L 2 ×L 2 with rows/columns indexed by black/white vertices b/w of T L , as follows. If b, w are not neighbors, then K(b, w) = 0. Otherwise, if (b, w) is an edge of type r one sets K(b, w) = K r , cf. (2.28). We actually need four Kasteleyn matrices K θ , θ = (θ 1 , θ 2 ) ∈ {0, 1} 2 , where the two indices label periodic/anti-periodic boundary conditions (depending on whether the index is 0/1) in the directions e i . To obtain K θ from K, one multiplies by (−1) θ 1 (resp. by (−1) θ 2 ) the matrix elements corresponding to edges (b, w) where w has first coordinate equal L and b has first coordinate equal 1 (resp. w has second coordinate equal L and b has second coordinate equal 1). See Fig. 1. Of course, K 00 = K. We have then [29,30] that where c θ ∈ {−1, +1} and, moreover, three of the c θ have the same sign and the fourth one has the opposite sign. More precisely, for the square grid, with our choice of Kasteleyn matrix, one finds (recall that we are assuming that L is even). The matrices K θ are diagonalized in the Fourier basis and where µ(·) is as in (2.5) and The matrices K θ are not necessarily invertible (e.g., if t i ≡ 1 then K 00 is not because µ(0) = 0) and this question will play a role in Section 6. When the four matrices K θ are invertible, the correlation functions of the non-interacting measure can be written as where the edge e j has black/white vertex b j /w j . The inverse of the matrix K θ can be computed explicitly as where w x (resp. b y ) is the white (resp. black) site with coordinate x (resp. y). Provided that it is easy to see that Condition (3.7) can fail for some values of L and of θ. For this reason, in Section 6 the values k ± θ ∈ P(θ) that are closest to the zeros of µ will be treated separately, see in particular Sections 6.1 and 6.5.
Due to the two zeros of µ, the matrix element g(x, y) decays only as the inverse distance between w x and b y . More precisely where r(x, y) = O(1/|x − y| 2 ) and φ ω was defined in (2.20).

3.2.
Grassmann representation of the generating functions. We refer for instance to [20] for an introduction to Grassmann variables and Grassmann integration; here we just recall a few basic facts. It is well known that determinants can be represented as Gaussian Grassmann integrals. For our purposes, we associate a Grassmann variable ψ + x (resp. ψ − x ) with the black (resp. white) site indexed x. We denote by Dψf (ψ) the Grassmann integral of a function f and since the variables ψ ± x anti-commute among themselves and there is a finite number of them, we need to define the integral only for polynomials f . The Grassmann integration is a linear operation that is fully defined by the following conventions: the sign of the integral changes whenever the positions of two variables are interchanged (in particular, the integral of a monomial where a variable appears twice is zero) and the integral is zero if any of the 2|Λ| variables is missing. We also consider Grassmann intergrals of functions of the type f (ψ) = exp(Q(ψ)), with Q a sum of monomials of even degree. By this, we simply mean that one replaces the exponential by its finite Taylor series containing only the terms where no Grassmann variable is repeated. It is well known that the definition of Grassmann integration allows one to write the determinant of a matrix as the integral of the exponential of the associated Grassmann quadratic form (such integral will be called a "Gaussian Grassmann integral", for the obvious formal analog with usual Gaussian integrals). In particular, where 3.12) and the index (θ) below the integral means that one has to identify ψ ± (L+1,x 2 ) := (−1) θ 1 ψ ± (1,x 2 ) and similarly ψ ± (x 1 ,L+1) := (−1) θ 2 ψ ± (x 1 ,1) . More compactly we write where the sum runs over edges of T L and, if e is an edge (b, w), Our goal here is to express, via a Grassmann integral, the partition function of the interacting dimer model, and more generally the generating function W Λ (A) defined by (3.14) where the product runs over the edges of T L and A e ∈ R. Note that e W L (0) is the partition function and that any multi-dimer truncated correlation function of the type E L,λ (1 e 1 ; . . . ; 1 e k ) can be obtained by differentiating W Λ (A) with respect to A e 1 , . . . , A e k and setting A ≡ 0.
Recall that the perturbed probability weight p L,λ depends on the local 'energy function' f via (2.11)-(2.12). Without loss of generality, we can assume that (2.12) holds with where c s are real constants, n is an integer, P s are finite collections of edges such that no space translation of P s coincides with a P s , s = s and 1 Ps = e∈Ps 1 e is the indicator that all edges in P s belong to M . Again without loss of generality we assume that each P s contains at least 2 edges (if P s consists in just one edge, its effect is just to modify the weights t). Under these assumptions, the following representation holds. Proposition 1. Let λ be small enough. Then, one has The first sum runs over all edges of T L and E e is as in (3.13). In the second sum, γ are finite subsets of disjoint edges of T L such that |γ| ≥ 2, and c(γ) is a real constant satisfying translation invariance (c(γ) = c(τ x γ)) and the bound for some constants a, b > 0, independent of L, and δ(γ) the tree distance of γ, that is, the length of the shortest tree graph on Λ containing γ (the precise definition of c(γ) is given below).
Remark 6. Both S(ψ) and V (ψ, A) are invariant under the following symmetry transformation of the Grassmann fields: where c → c * indicates that all the constants appearing in S(ψ) and V (ψ, A) are mapped to their complex conjugates. Also, we used the notation (−1) x := (−1) x 1 +x 2 . It is straightforward to check that, under this transformation, E e → E e , for all the edges e, which clearly shows that the considered transformation is in fact a symmetry of the Grassmann action. This symmetry will play a role in Section 6, in reducing the number of independent running coupling constants arising in the multiscale computation of the Grassmann generating function.
Proof of Proposition 1. The proposition has been proven in [23] in the case of constant weights t i ≡ 1 and plaquette interaction as in (2.14); the extension to the present situation is rather straightforward, so we will be concise. Let S = {τ x P s , s = 1, . . . , n, x ∈ Λ} and remark that by assumption all elements of S are distinct and contain at least two edges. If B ∈ S is a space translation of P s , set We start by writing and P . By manipulating the sum in the r.h.s. of (3.21), one can rewrite it as where the term n = 0 has to be interpreted as equal to 1 and the sum * is over non-empty, mutually disjoint subsets γ i of edges of T L . The constant c(γ) is given as follows. Let Σ γ be the set of all collections of the type The partition function Z corresponds to a non-interacting dimer model, with edge-dependent weights t e e Ae . Then, as in (3.1) and (3.11) we have Using expression (3.25) in (3.24) one readily concludes, as in [23], that (3.17) holds with If λ is small enough, it is easy to see that the bound (3.18) holds.
For the 6-vertex model with interaction (2.15), the potential V is exactly quartic in the fields ψ: indeed, c(γ) = 0 only if γ is the pair of edges γ = {e 1 , e 2 } or γ = {e 3 , e 4 } as in Fig. 2 or a translation thereof. For the plaquette model with interaction (2.14), instead, c(γ) is non-zero only if γ is a collection of |γ| ≥ 2 adjacent parallel edges, in which case In the following (in the comparison between the discrete lattice model and the continuum reference model) we will also need the generating function for mixed dimer and fermionic correlations. Namely, let {φ + x , φ − x } x∈Λ be Grassmann variables that anti-commute among themselves and with the ψ ± variables. Then, we let Dψ e S(ψ)+V (ψ,A)+(ψ,φ) (3.27) and We define g L (e 1 , . . . , e k ; x 1 , . . . , x n ; y 1 , . . . , y n ) as the truncated correlations associated with the generating function 5 W L (A, φ): g L (e 1 , . . . , e k ; x 1 , . . . , x n ; y 1 , . . . , y n ) Two cases that will play a central role in the following are k = 0, n = 1 (the interacting propagator), and k = n = 1 (the interacting vertex function), which deserve a distinguished notation.
Interacting propagator: is the Grassmann counterpart of the dimer observable at e, and e is an edge of type r with black site labelled z, then where the semicolon indicates truncated expectation.
In the following we will also need a distinguished notation for the twopoint dimer-dimer correlation: if e 1 , e 2 are two edges of type r, r , and black sites labelled x, y, respectively, we let Note that all the correlations g L (e 1 , . . . , e k ; x 1 , . . . , x n ; y 1 , . . . , y n ) are well defined for any finite L, despite the fact that the Kasteleyn matrix K θ may not be invertible for some choices of θ, L. The multipoint correlations, g L (e 1 , . . . , e k ; x 1 , . . . , x n ; y 1 , . . . , y n ), admit a thermodynamic limit as L → ∞, as shown in Section 6; the limit can be expressed as a convergent multiscale fermionic expansion and will be denoted g(e 1 , . . . , e k ; x 1 , . . . , x n ; y 1 , . . . , y n ).
In particular, the thermodynamic limit of the two-point dimer-dimer correlation will be denoted by G (0,2) r,r (x, y), while the L → ∞ limit of the interacting propagator and vertex function will be denoted G (2) (x, y) and G (2,1) (z, x, y).

3.3.
Lattice Ward Identity. The generating function W L (A, φ) has a gauge symmetry property that implies certain identities (lattice Ward identities) involving its derivatives. These identities were derived in [24] for the model with t i ≡ 1 and they hold (with the same proof) also for the general model studied here. We recall here, without giving the proof, the Ward Identity for the 'vertex function', but similar relations can be easily derived for higher point correlations: for any finite L, with δ x,y the Krokecker delta, see [24, Eq.(4.9)-(4.10)]. By taking the difference between these two equations, we get (see [24,Eq.(4.17) is the (un-normalized) discrete derivative acting on the x variable. By taking the limit L → ∞, we get the infinite volume version of (3.33)-(3.35).
In Fourier space, we definê Then, the infinite-volume limit of (3.33)-(3.35) can be rewritten as In the following the asymptotic behavior at large distances of the interacting propagator and vertex function will be computed in terms of a reference continuum model, see next section, which plays the role of the 'infrared fixed point' of our lattice dimer model in its Grassmann formulation.

The infrared fixed point theory
In order to introduce the "infra-red fixed point" of our theory (referred to in the following as "the continuum model" or "the reference model"), we need a couple of preliminary definitions. First, we let M be the 2 × 2 matrix with unit determinant whereᾱ j ,β j ∈ R, j = 1, 2 and ∆ =ᾱ 1β2 −ᾱ 2β1 > 0 (for the moment, these are free parameters; eventually, they will be the real and imaginary parts of the functionsᾱ ω ,β ω that appear in Theorem 1). Also, given L > 0 (the system size), an integer N (ultra-violet cut-off) and Z > 0, we introduce a Grassmann Gaussian integration 6 P [≤N ] Z (dψ) on the family of Grassmann 6 We recall (cf. e.g. [20,Sec. 4]) that, given a family {ψ − x , ψ + x }x∈I of Grassmann variables and a |I| × |I| matrix g, the "Grassmann Gaussian integration with propagator g", denoted sometimes Pg(dψ) . . . in the following, is the linear map acting on polynomials of the Grassmann variables, such that Pg(dψ)ψ − defined by the propagator where: that is equal to 1 if its argument is smaller than 1 and equal to 0 if its argument is larger than 2; Observe that, since we are assuming ∆ > 0, we have that While K is an infinite set, we effectively have only a finite number of nonzero Grassmann variablesψ ± k,ω , because χ N (k) is non-zero only for a finite number of values of k in K.
Note that, setting q = M −1 k, the r.h.s. of (4.3) equals In the language of Quantum Field Theory, in the limit lim L→∞,N →∞ , (4.6) is just the propagator of chiral massless relativistic fermions. It is convenient to define, for x ∈ R 2 , the Grassmann variables Note that ψ ± x,ω has anti-periodic boundary conditions on and that The generating functional W L,N (J, φ) of the continuum model is non-singular, one can write more explicitly x,ω } j=1,2 ω=±, x∈Λ are external "sources" (real-valued test functions) and φ = {φ σ x,ω } σ,ω=± x∈Λ are "external Grassmann sources", i.e. φ σ x,ω is a Grassmann variable. Also, we used the notation and Finally, the interaction V in (4.9) is and v 0 (·) is a smooth rotationally invariant potential, exponentially decaying to zero at large distances, normalized as We emphasize that, while this expression seems to depend on an uncountable set of Grassmann variables {ψ ± x,ω , φ ± x,ω } x∈Λ , writing everything in Fourier space there is only a finite number of non-zero Grassmann variables.
Remark 7. In order to recognize the equivalence of the model (4.9) with α ω = −i,β ω = ω and the one in, e.g., [11, Section 3] (or, analogously, the one in [6, Section 3]), one needs to set to zero some of the external fields, rotate the coordinate system and rescale some constants. More precisely, if W L,N (J (1) , φ) denotes the generating functional used in [11] with J (1) where the constant is independent of J (1) , φ (so that it does not influence the correlation functions; it depends upon ∆ and is due to the rescaling of the Grassmann fields), while and we denoted explicitly the dependence of the generating function on λ ∞ . This immediately implies obvious relations between the correlation functions G R,ω (x, y) and S (j,j) R,ω,ω (x, y), defined below, and the analogous ones of [11].
The peculiarity of the continuum model is that its correlations can be computed exactly. This is because, as compared to its lattice counterpart, the continuum model is "chiral gauge invariant", which means that the correlation functions satisfy two hierarchies of Ward Identities, distinguished by the choice of the 'chirality index' ω, see (4.23) below. These additional symmetries, together with other identities among correlation functions (the so-called Schwinger-Dyson equations), allow one to get closed equations for correlations functions. In this sense, the infrared fixed point theory can be regarded as "integrable".
We shall use the following definitions: if x, y, z are distinct points of Λ, From the construction of the correlation functions of the model, see e.g. [6, Section 3 and 4], one obtains in particular the existence of the following limits where cut-offs are removed: Away from x = 0, the so-called "density-density" correlation S (1,1) is given by [24, Eq. (5.12)] On the other hand, the "mass-mass correlation" S (2,2) satisfies (see [24, Eq.(6.14)]) whereB is an analytic function of λ ∞ , Z,ᾱ ω ,β ω , which is equal to 1 at see [24,Eq.(6.15)] and [24, Appendix C], and R 2 is a correction term such that |R 2 (x)| ≤ C|x| −2−θ , for some θ > 0 that, e.g., can be chosen θ = 1/2. We will not need the explicit form of G R,ω (x, 0) and G (2,1) R,ω ,ω (x, 0, z); let us just mention that they diverge as x, z tend to zero but they are locally integrable functions (see, e.g., the expression of the interacting propagator in [6, eq.(4.18)]) and therefore admit Fourier transforms in the sense of distributions 7 , For later use, let us mention that the small-momenta behavior ofĜ R,ω are (see, e.g., [9, Theorem 2]) if p, k are both of order c → 0. A very useful consequence of the exact solution of the continuum model is that the "propagator" and the "vertex function" satisfy the following Ward Identity (see [24, Eq.(5.9)]): Note that this identity resembles formally the lattice Ward identity (3.40) of the dimer model, with the crucial difference that (4.23) are actually two identities (one for each choice of ω).

5.
Comparison between lattice and continuum model, and proof of Theorems 1-2 The reason why the continuum model plays the role of the "infrared fixed point theory" for our interacting dimer model is that the large distance behavior of the dimer correlation functions can be expressed in terms of linear combinations of the correlations of the continuum model, for a suitable choice of the parameters Z, λ ∞ ,ᾱ ω ,β ω . Let us spell out the explicit relation between correlation functions of the two models, in the special cases of the dimer interacting propagator, the vertex function and the dimer-dimer correlation. The result is a consequence of the multi-scale analysis described in Section 6 (see in particular Section 6.6) and can be stated as follows.
By replacing (5.4) into (5.7) and recalling that ν We claim that 4 r=1K ω,r = 0 (we shall prove this fact in a moment): therefore, the first equation can be rewritten asK ω,1 +K ω,4 = i √ νᾱ ω . In terms of the 'elementary steps' s(x, j) in direction e j centered at x, introduced before (2.29), the two equations in (5.10) become which are the desired identities. In order to complete the proof of (5.11)-(5.12), we need to prove that 4 r=1K ω,r = 0, as claimed above. For this purpose, we consider (3.39), and combine it with By using (5.8) and the fact thatv(0) = 1, this becomeŝ Now, using (4.21), (4.22) and taking the limit c → 0, one obtains that for any sequence p j along which the ratioD ω (p j )/D −ω (p j ) admits a limit. Note that, in general, the limit depends upon the chosen subsequence. For On the other hand, these two values cannot be equal since we know that the ratio α ω /β ω is not real (cf. (4.5)). In conclusion, we find that 4 r=1K ω,r = 0 that, in light of (5.4), is equivalent to 4 r=1K ω,r = 0, as desired. With the identities (5.11)-(5.12) at hand, we can easily prove (2.44), by repeating the analogue of the discussion leading, in the non-interacting case, to (2.31). We will be very sketchy since the analogous argument has been given in detail in [23] in the case of the model with weights t ≡ 1. We start from the very definition of the covariance of the height difference: (5.16) where C η 1 →η 2 and C η 3 →η 4 are two lattice paths connecting η 1 with η 2 , and η 3 with η 4 , respectively. For simplicity, we assume that η 1 and η 2 have the same parity, and similarly for η 3 and η 4 : in this way, it is possible to choose the two paths C η 1 →η 2 and C η 3 →η 4 to be concatenations of 'elementary steps' s(x, j) in directions ± e j , see the discussion after (2.28) above. For simplicity, let us also assume that the mutual distances between the faces η 1 , . . . η 4 are all comparable, i.e.
In this case, we choose the two paths C η 1 →η 2 and C η 3 →η 4 to be of length at most C max i =j |η i − η j | and to be at mutual distance We now insert (2.34) into (5.16) and, by repeating the discussion of [23, Section 3.2], we find that the dominant contribution comes fromĀ r,r (the contribution fromB r,r is sub-dominant due to the oscillating pre-factors): where r(e) is the type of the edge e, x(e) is the coordinate of the black site of e. By using the explicit expression ofĀ r,r , (2.35), and by decomposing the two paths C η 1 →η 2 , C η 3 →η 4 , into a sequence of elementary steps, we obtain (denoting the generic elementary step in C η 1 →η 2 , resp. C η 3 →η 4 , by s(x, j), resp. s(x , j )) We now use (5.11)-(5.12), the symmetryφ ω =φ * −ω and realize that the dominant term in (5.21) is the Riemann sum approximation to the following integral: whose explicit evaluation gives the main term in the r.h.s. of (2.44). Putting together the error terms, we obtain the statement of Theorem 2, as desired.

Renormalization Group analysis
In this section we discuss the multiscale analysis of the dimer model and the comparison with the continuum model, which leads us to the results spelled out in the last two sections.
The goal is to obtain sharp estimates on W (θ) L (A, φ), see (3.27), as L → ∞, for all θ ∈ {0, 1} 2 . These will then be combined as in (3.28), to finally obtain the control of the large-scale behavior of the correlation functions of the interacting dimer model. From now on, C, C , . . . , and c, c , . . . , denote universal constants, whose specific values might change from line to line.
6.1. Preliminaries. As a preliminary step, we rewrite the quadratic part S of the action in (3.27) as a "dressed" term S 0 plus a "counter-term" N = S − S 0 , whose role is to fix the location of the interacting Fermi points and Fermi velocities. Namely, letting as usual we write: In this equation: (1)p ω =p ω (λ), with ω = ±, are points in [−π, π] 2 , such thatp + +p − = (π, π), and they will be fixed via the multiscale construction. A posteriori they can be interpreted as "dressed Fermi points"; they are the same functions appearing in Theorem 1. (2) a ω = a ω (λ) ∈ C and b ω = b ω (λ) ∈ C are such that a ω = −a * −ω and b ω = −b * −ω ; they will also be fixed via the multiscale construction. A posteriori, their choice fixes the "dressed Fermi velocities" via the following relations: whereᾱ ω ,β ω are the same functions appearing in Theorem 1. (3) the functionχ 0 is defined as:χ 0 (k ) =χ(|M −1 k |), where: (1) M is the same matrix as (4.1), withᾱ 1 andᾱ 2 (resp.β 1 andβ 2 ) the real and imaginary parts ofᾱ + (resp.β + ); (2)χ : R + → [0, 1] is a C ∞ cut-off function in the Gevrey class of order 2 (see [23,Appendix C]) that is equal to 1, if its argument is smaller than c 0 /2, and equal to 0, if its argument is larger than c 0 ; here c 0 is a small enough constant, such that in particular the support ofχ 0 (· −p + ) is disjoint from the support ofχ 0 (· −p − ). For later reference, we also let for h a negative integerχ From the properties just stated ofp ω , a ω , b ω andχ(·), we see that µ 0 ((π, π) − k) = µ * 0 (k). (6.7) In the integration over ψ in (3.27), the Fourier modes k that are the closest to the zeros of µ 0 (·) play a somewhat special role, so they have to be treated separately, at the very last step of the multi-scale procedure (cf. Section 6.5). Namely, given θ ∈ {0, 1} 2 , let k ± θ be the values of k ∈ P(θ) that are closest top ± and note that k + θ = (π, π) − k − θ [If there is more than one momentum at minimal distance fromp ± (there are at most four), any arbitrary choice will work]. Next, we decompose the quadratic action S 0 (ψ) as a sum of a term depending only on k ± θ plus a term depending only on the modes in P (θ) := P(θ) \ {k + θ , k − θ }, and we rewrite (3.27) as . We multiply and divide by (the product is non-zero since we singled out the possibly zero modes k ± θ ) and, lettingΨ we rewrite the generating function as with P g 0 the Grassmann Gaussian integration (cf. footnote 6) with propagator (6.14) From this point, we proceed as follows. First, we perform in a multi-scale way the integration over the Grassmann variables ψ, i.e. over the Fourier modes except k ± θ : the inductive integration procedure, including the definition of the running coupling constants (RCC), is described in Section 6.2; the outcome of the construction can be conveniently expressed in terms of a Gallavotti-Nicolò tree expansion, similar to the one described in [23, Section 6.2]. The main definitions (and the main differences compared to the case treated in [23]) are summarized in Section 6.3; in the same section, we also state the bounds satisfied by the kernels of the effective potential, see Proposition 2, under the assumption that the RCC are uniformly bounded in the infrared, see condition (6.64). The proof that the RCC remain in fact bounded under the iterations of the renormalization group map is given in Section 6.4; the flow of the RCC can be controlled only if their initial data are properly fixed: as shown there, the choice of the initial data fixes the dressed Fermi pointsp ω and the dressed Fermi velocitiesᾱ ω ,β ω , as anticipated after (6.3). In Section 6.5 we describe the integration of the last two modes and prove the existence of the thermodynamic limit for the correlation functions, with explicit bounds on the speed of convergence as L → ∞. Finally, in Section 6.6, we compute the fine asymptotics of the correlations functions, via a comparison of the tree expansion of the dimer model with that of the continuum model of Section 4; in particular, we show how to obtain (5.1), (5.2), (5.3), relating the dimer correlations with those of the reference model, thus concluding the proofs of Theorems 1-2.

6.2.
Multi-scale analysis. In this section we describe the multi-scale computation of W (θ) L (A, φ, Ψ), see (6.13). We consider explicitly only the case φ = 0; the general case can be treated analogously but we will not belabor the details in this paper.
The procedure is based on a systematic use of the 'addition principle' for Gaussian Grassmann integrals, namely the following property [20,Sec. 4]: if P g (dψ) is the Grassmann Gaussian integration with propagator g and g = g 1 + g 2 then P g (dψ)F (ψ) = P g 1 (dψ 1 )P g 2 (dψ 2 )F (ψ 1 + ψ 2 ). (6.15) We apply this formula to P g 0 , in connection with the following decomposition of the propagator g 0 (x, y): and, if P ω (θ) = {k : k +p ω ∈ P (θ)}, By using the decomposition (6.16) and (6.15), we rewrite (6.13) as where ψ (0) + ψ (≤−1) + Ψ is a shorthand notation for (6.20) P (0) is the Grassmann Gaussian measure with propagator g (0) (x, y), while P (≤−1) is the Grassmann Gaussian measure with propagator Since the cutoff functionχ −1 in (6.18) is a Gevrey function of order 2, the propagator g (0) has stretched-exponential decay at large distances: for suitable L-independent constants C, κ > 0, if |x − y| is the distance on the torus Λ. This is seen by writing g (0) via the Poisson summation formula as a sum of Fourier integrals, as in [23, App. A]; each integral decays in the desired way because it is the Fourier transform of a Gevrey function [36]. Next, we denote by V (0) (·, J) the combination N (·)+V (·, A), re-expressed in terms of the variables J = {J x,r } x∈Λ, 1≤r≤4 , instead of A: here, if b is the bond of type r and black site x, we let J x,r := e A b − 1. The result of the integration over ψ (0) is rewritten in exponential form: where [20,Sec. 4] with E T 0 the truncated expectation 8 w.r.t. the Grassmann Gaussian integration P (0) (dψ (0) ), and E (−1) , S (−1) (·) are fixed by the condition S (−1) (0) = 0, V (−1) (0, J) = 0. The series in the r.h.s. is absolutely summable, for λ sufficiently small (independently of L), see [20,Sec. 4.2]. The reason is that the propagator g (0) has a fast decay in space, uniformly in L, as in (6.21).
Here n,m;ω,r are absolutely convergent series of (ν 0,ω , a 0,ω , b 0,ω , λ 0 ), that each term in the expansion admits a limit as L → ∞ (as one can check by inspection) and that they satisfy uniform bounds as L → ∞, see (6.29), implies that their infinite volume limits exist and satisfy the same bounds. For later reference, the infinite volume limit of W (−1) n,m;ω,r will be denoted by W (−1),∞ n,m;ω,r , and similarly for its Fourier transform.
After this first integration step, we still need to integrate ψ (≤−1) out, see (6.19). Let us first informally explain how this is done, before giving the precise inductive procedure in Sections 6.2.1-6.2.3. The idea is to repeat the same procedure as above: we rewrite (via the addition principle) ψ ) is a Grassmann field with propagator supported, in momentum space, on momenta k ∈ P ω (θ) with |k | ∼ 2 −1 (resp. |k | 2 −2 ); we integrate ψ (−1) ω out; we exponentiate the result of the integration, thus defining the effective potential on scale −2, in analogy with (6.22)-(6.23); and so on. One after the other, we integrate the fields ψ (−2) , . . . , ψ (h+1) out, define the effective potential V (h) on scale h (which involves fields ψ (≤h) with momenta k ∈ P ω (θ) that belong to the support ofχ h (·) (cf. (6.6)), and continue until we reach the 'last scale', h L , fixed by the finite volume L, which induces a natural infrared cut-off. More precisely, h L is fixed as the smallest (in absolute value) negative integer h such that the support ofχ h (·) has empty intersection with P ω (θ). Note that, since all momenta in P ω (θ) are at distance at least π/L fromp ω , we have h L ∼ − log 2 L for L large. The result of the integration of the Grassmann fields ψ (≤h L ) gives the generating function W In order for the bounds on the generating function to be uniform in L, we need to improve the procedure roughly described here: at each step, before integrating the field on the next scale, we actually need to isolate and re-sum a certain selection of potentially dangerous contributions to the effective potential, the so-called marginal and relevant terms. We refer, e.g., to [23,Sec. 5], see in particular [23, Section 5.2.2] for a dimensional classification of the divergent terms arising in a 'naive' multiscale scheme. As discussed there, see [23, Eq. (5.8)] and following lines, the scaling dimension of the kernels with n external fields of type ψ and m external fields of type J is 2 − n/2 − m; in the renormalization group jargon, positive scaling dimension (that is, 2 − n/2 − m > 0 ⇔ (n, m) = (2, 0)) corresponds to relevant contributions, vanishing scaling dimension (that is, 2 − n/2 − m = 0 ⇔ (n, m) = (4, 0), (2, 1)) corresponds to marginal contributions, and negative scaling dimension corresponds to irrelevant ones. In order to cure the potential divergences associated with the terms with (n, m) = (2, 0), (4, 0), (2, 1), at each step of the multiscale construction we properly 'localize' and re-sum these terms, via an iterative procedure that we now describe. 6.2.1. The inductive statement. Let us inductively assume that the fields ψ (0) , ψ (−1) , . . . , ψ (h+1) , h ≥ h L , have been integrated out, and that after their integration the generating function has the following structure, analogous to the one at scales 0, −1: for suitable real constants E (h) , Z h , and suitable 'effective potentials' S (h) (J), V (h) (ϕ, J), to be defined inductively below, and fixed in such a way that V (h) (0, J) = S (h) (0) = 0. In the second line, whereD ω (k) =ᾱ ω k 1 +β ω k 2 and is a remainder of order O(k 2 ) for k small. Finally, P (≤h) (dψ) is the Grassmann Gaussian integration with propagator (diagonal in the index ω) We will also prove inductively that: (1) V (h) (ϕ, J) has the same structure as (6.25), with the upper index (−1) in the kernels replaced by (h); (2) the kernels of V (h) (ϕ, J) satisfy the following symmetry: Remark 9. It is important to emphasize right away that we will view the kernels W (iii) the irrelevant part of V (−1) , denoted by RV (−1) . The running coupling constants, as well as the irrelevant part of the effective potentials, will be defined along the iterative procedure.
6.2.3. The inductive step. We assume that (6.31) holds with V (h) satisfying the properties specified in the inductive statement, and we discuss here how to get the same representation at the next scale h − 1. First, we split V (h) into its local and irrelevant parts: n,m;ω,r the infinite volume limit ofŴ

Remark 10.
A few remarks about this definition are in order: (1) The existence of the limit ofŴ (h) n,m;ω,r as L → ∞ is a corollary of the inductive bounds on the kernels of V (h) , which are uniform in L, as it was the case for h = −1, cf. with Remark 8. More details on the inductive bounds on the kernels of V (h) are discussed below.
(2) The reason why, in the second line of (6.36), we only include terms where the Grassmann fields have the same index ω, is that the terms with opposite ω indices give zero contribution to the generating function, due to the support properties of the Grassmann fields. In fact, in (6.31) we need to compute V (h) at Grassmann fieldsψ (≤h)± k,ω that, in momentum space, have the same support asĝ (≤h) ω (k), i.e., |M −1 k| ≤ c 0 2 h (note that the support properties ofĝ (≤h) ω are the same as those ofχ h (cf. (6.6)), and these were discussed in the third item after (6.3)). If h ≤ −1 and c 0 is sufficiently small, quadratic terms of the formψ k+p ω −p −ω ,−ω would involve two fields that cannot both satisfy this support property.
(3) Due to the Grassmann anti-commmutation rules and the anti-symmetry of the kernels, the quartic term in (6.36) can be rewritten as Along the induction step, we will need a function W The function W (h),R 2,0;(ω,ω) will be shown to satisfy both the identity (6.34) and the extra symmetries (in Fourier space) (ω,ω) , h ≥ h has been already shown to satisfy (6.38) and below we explain how to prove the same at scale h − 1.
In order to define the running coupling constants on scale h, we decompose the term containing ∂ kŴ for some real constant z h . We now combine this term with the Grassmann Gaussian integration P (≤h) (dψ), and define: (6.43) and e L 2 t h is a constant that normalizesP (≤h) (dψ) to 1: By using (6.41), we rewrite the Grassmann integral in the right side of (6.31) as and the running coupling constants at scale h are defined as (h),s 2,0;(ω,ω) (0), (ω 1 ,ω 2 ),r (0,p ω 1 −p ω 2 ). Thanks to the symmetry (6.34) of the kernels (that by inductive hypothesis holds at step h) the running coupling constants satisfy the following:  For later reference, we rewrite the local part of V (h) (ϕ, J) as (for the definitions of F ν;ω (ϕ), F a;ω , F b;ω , etc., compare (6.49) with the first two lines of (6.46)). We now decompose the propagator (6.42) as with g (≤h−1) ω as in (6.33) and g (h) ω as in (6.35). To see that this decomposition holds, note thatZ h−1 (k) ≡ Z h−1 on the support ofχ h−1 (·).
To conclude the proof of the induction step, it remains to prove that the kernels of V (h−1) satisfy (6.34) and that (6.38) holds, at scale h − 1. The proof of the former statement is very similar (but not identical) to the argument used in Section 6.2.2 to prove (6.34) at scale h = −1 starting from the symmetries of V (0) . Namely, thanks to (6.34) at scale h, the potential x,−ω together with complex conjugation of the kernels. Then, the claim follows from the representation (6.52), together with the fact that the propagator g (h) (defined in (6.35)) satisfies the symmetry −ω (x, y). (6.53) As for (6.38) at scale h − 1, the proof uses the symmetries of the relativistic propagator (6.37), together with the fact that λ h is real. See Appendix A.
Remark 11. Note that, if the function z h in (6.43) is sufficiently small for all the scales h ≤ h ≤ −1, say |z h | ≤ uniformly in L, h , then e −c |h| ≤ Z h ≤ e c |h| . As a consequence, g ω satisfies a bound analogous to (6.21), namely In fact, note that on the support of f h (·) (which is concentrated on k : |k| ∼ 2 h ),Z h−1 (k)/Z h−1 = 1 + O( ) and recall that r ω (·) is quadratic for small values of its argument, so that R,ω satisfies the same estimate as (6.54), while the difference g R,ω satisfies an estimate that is better by a factor 2 h . 6.2.4. The Beta function. The iterative integration scheme described above allows us to express the kernels of V (h) and, in particular, the running coupling constants (RCC) at scale h, as functions of the sequence of RCC on higher scales, where B # h,· , h ≤ −1, is the so-called Beta function. One has to think of B # h,· as a function of the RCC on scales h with h ≤ h ≤ 0. Note that the first four equations makes sense also with h = 0, in which case they express the relation between (ν −1,ω , a −1,ω , b −1,ω , λ −1 ) and (ν 0,ω , a 0,ω , b 0,ω , λ 0 ), see (6.28). Note also that by construction the beta function B # h,· depends on Z h only via the combinations Z h /Z h −1 = (1 + z h ) −1 , with h < h < 0. For later reference, we rewrite the definition of z h , (6.40), in a form analogous to (6.55), z h−1 = B z h , h ≤ 0, (6.56) where the right side is thought of as a function of (λ h , z h ) h≤h ≤0 , with the convention that z 0 = z −1 = 0 (the latter is because W (−1),R 2,0;(ω,ω) ≡ 0). Remark 12. The components of the beta function for ν h,ω , a h,ω , b h,ω , λ h are independent of Y h ,r,ω , h > h. Therefore, we can first solve the flow equation for ν h,ω , a h,ω , b h,ω , λ h and then inject the solution into the flow equation for Y h,r,ω .
Before we proceed in describing the dimensional bounds satisfied by the kernels of the effective potential, let us comment on their structure. We have proven inductively that V (h) has, in momentum space, the same structure as (6.25). If one writes V (h) in real space, due to iterative action of the R operator in the inductive procedure explained above, the structure that emerges naturally is that of a polynomial with spatial derivatives acting on some of the Grassmann fields ϕ ± x,ω . For an explanation of why this is the case see [23, Section 6.1.4] and Appendix B below, where finite-size effects associated with the action of R are also discussed. Correspondingly, V (h) can be represented as V (h) (ϕ, J) = n,m≥0: n even, n≥2 x, y, ω, r, i,q W (h) n,m,i,q;ω,r (x, y) × (6.57) ×∂ q 1 i 1 ϕ (≤h)+ x 1 ,ω 1 · · ·∂ qn in ϕ (≤h)− xn,ωn J y 1 ,r 1 · · · J ym,rm . The main difference between this formula and (6.24), besides the scale label h replacing −1, is the presence of the indices i = (i 1 , . . . , i n ) ∈ {1, 2} n and q = (q 1 , . . . , q n ) ∈ {0, 1, 2} n and the operators∂ q i acting on the Grassmann fields: this is a differential operator, dimensionally equivalent to a derivative of order q in direction i. Let us stress that the representation in (6.57) is not unique: the claim is that there exists such a representation, with the kernels satisfying natural dimensional estimates, discussed below.
In order for the iterative construction to allow us to compute the thermodynamic and correlation functions, we need to prove that: (i) the RCC ν h,ω , a h,ω , b h,ω , λ h , z h are small, uniformly in the scale (say, smaller than a sufficiently small constant ε), provided the functionsp ω , a ω , b ω (see (6.3)) have been properly fixed; (ii) the kernels of the effective potential are all well defined (i.e. the sums (6.52) are convergent uniformly in L), quasilocal (i.e., fast decaying, with a stretched-exponential behavior) and satisfy natural scaling properties, i.e., for suitable constants C, c, κ, independent of L, h.
The boundedness of the flow of the RCC and the validity of the dimensional bounds for the kernels of the effective potential will, in fact, be the final outcome of our analysis. The logic of proof goes as follows: one first proves the validity of the dimensional bounds on the kernels, under the assumption that the RCC remain small. These bounds will, in particular, imply that the components of the beta function are well defined and satisfy bounds that are uniform in L and h. This part of the proof is pretty standard: it follows from a representation of the effective potential in terms of Gallavotti-Nicolò (GN) trees, see Section 6.3 below, and an iterative application of the Battle-Brydges-Federbush-Kennedy (BBFK) determinant formula, see, e.g., [23,Lemma 3].
Next, we prove that the RCC remain bounded, by studying the flow generated by the beta function. The key point is that, as long as the RCC on scales larger than h are small, then the beta function on scale h is well defined, and can be used to control the evolution of the RCC for another step. This opens the way to an inductive proof of the smallness of the RCC. Of course, the fact that RCC remain small at all scale requires a specific (modeldependent) structure of the beta function. In our case, we are lucky enough that the beta function has structure which maintains the RCC small at all scale, provided the initial data are small, and thatp ω , a ω , b ω are properly fixed, see Section 6.4 below. It is not just a matter of luck, of course: a key point in the analysis is played by the comparison of the λ-component of the beta function of our dimer model, with the corresponding quantity for the reference continuum model (the two functions are the same at dominant order). The exact solvability of the reference model implies the validity of a remarkable cancellation in the λ-component of the beta function for the reference model and, therefore, a posteriori, for our dimer model, as well. 6.3. Tree expansion for the effective potential. As anticipated above, the detailed structure of the kernels of V (h) , arising from the iterative construction described in the previous section, can be conveniently represented in terms of GN trees. The definition of GN trees, of their values, and the procedure leading to their introduction have been discussed at length in several previous papers and will not be repeated here, see e.g. [20]; in particular, we refer to [23, Section 5.2.1 and 6.2] for a description of the GN tree expansion in a context very similar to the present one, i.e., in the case of isotropic, 'tilt-less', interacting dimer models with weights t ≡ 1 and plaquette interaction. The present case differs from the one treated in [23] for the fact that here the model is anisotropic (and, correspondingly, the height has an average slope that is different from zero). Technically, this means that in the present case the expansion involves more running coupling constants than those considered in [23]: the RCC ν h,ω , a h,ω , b h,ω are identically zero in the tilt-less case. In particular, the trees involved in our construction are characterized by the following features, slightly different from those listed in [23, Section 6.2]: has root on scale h and can have endpoints (either 'normal' or 'special', which are those represented as black dots or white squares, respectively, in [23], see, e.g., [23, ), depending on its type (recall that the monomials F λ , F ν;ω , etc., were defined in (6.49)); in this case, the node v immediately preceding v on τ , of scale h v = h v − 1, is necessarily a branching point. If v is of type RV (−1) , then h v = 0, and v is associated with (one of the monomials contributing to) RV (−1) (ϕ (≤−1) , 0); in this case, the node immediately preceding v on τ , of scale h v − 1, is not necessarily a branching point.
(3) A special endpoint v on scale h v ≤ 0 can be either local, or non-local.
If v is local, then it is associated with for some r ∈ {1, 2, 3, 4}, ω = (ω 1 , ω 2 ) ∈ {±} 2 ; if ω 1 = ω 2 , we shall say that v is a 'density endpoint', while, if ω 1 = ω 2 , that v is a 'mass endpoint'. Note that the factors Z h v −1 in (6.60) simplify: the summand equals Y h v ,r,ω F Y ;r,ω (ϕ (≤h v ) , J); in (6.60), these factors are kept just for uniformity of notation with the cases in the previous item. In the case that v is local, the node v immediately preceding v on τ , of scale h v = h v − 1, is necessarily a branching point. If v is non-local, then h v = 0, and v is associated with (one of the monomials contributing to) −1) , 0); in this case, the node immediately preceding v on τ , of scale h v − 1, is not necessarily a branching point. In addition to the items above, let us recall that each vertex of the tree that is not an endpoint and that is not the special vertex v 0 (the leftmost vertex of the tree, immediately following the root on τ ) is associated with the action of an R operator.
In terms of the tree expansion, we can express the effective potential and the single-scale contributions to the free energy and generating function as Eq.(6.62) is the analogue of [23, (6.64)], and we refer the reader to that paper for the notation and a sketch of the proof (in this formula, the indices i, q replace the multi-indices β ∈ B T [23, (6.64)]). We recall that P ψ v 0 and P J v 0 are two sets of indices that label the Grassmann external fields and the external fields of type J, respectively; moreover, Clearly, the kernels in (6.57) are obtained by summing W τ,P,T,i,q ( Nn,Ns , under the constraint that the number of external fields of type ψ and J is equal to n and m, respectively, that the elements of i are the same as i, etc. The bound (6.58) is a corollary of the following fundamental bound on the weighted L 1 norm of W τ,P,T,i,q , which is the analogue of [23,Proposition 8] and of [10, (3.110)]; for the proof, we refer the reader to [10,23]. See also Appendix B below for some technical details.
where: |I ψ v 0 | = v e.p. |P ψ v | is the total number of Grassmann fields associated with the endpoints of the tree; the first product in the second line runs over the special endpoints, while the second over all the vertices of the tree that are not endpoints. Since the exponents 2 (6.65) are all strictly negative, one can sum (6.65) over τ ∈ T (h) Nn,Ns , over T ∈ T, and over P ∈ P τ , under the constraint that |P ψ v 0 | = n and |P J v 0 | = m, we get the bound (6.58); see also the discussion after [23,Proposition 8]. Similarly, if we sum (6.65) over τ ∈ T (h) Nn,Ns , T ∈ T, P ∈ P τ , with |P ψ v 0 | = n, |P J v 0 | = m, under the additional constraint that τ has at least one node on scale k > h, then we get a bound that is the same as (6.58) times an additional gain factor 2 θ (h−k) , where θ is a positive constant, smaller than 1 (estimates are not uniform as θ → 1 − ; from here on, we will choose θ = 3/4). This is the so-called short memory property, see Remark 16 after [23, Proposition 8].
6.4. The flow of the running coupling constants. As explained above, as long as the RCC ν h ,ω , a h ,ω , b h ,ω , λ h , z h stay small, for all h > h, in the sense of (6.64), the beta function controlling the flow of the same constants on all scales larger or equal to h, see (6.55)-(6.56), can be represented in terms of a convergent GN expansion, induced by the one of the kernels of the effective potential discussed above. The goal is then to fix the initial data ν 0,ω , a 0,ω , b 0,ω , in such a way that the resulting flow of ν h,ω , a h,ω , b h,ω , λ h , z h driven by the beta function stays uniformly small in the scale index. For this purpose, not only we have to make a careful choice of the 'counterterms' ν 0,ω , a 0,ω , b 0,ω , but we also need to exploit a number of remarkable cancellations, some of which follow from the exact solution of the reference model of Section 4. Let us emphasize that we have the right to fix the counter-terms, that up to now were chosen arbitrarily in (6.3) (recall (6.28)), but we cannot change λ 0 = λ, that enters the definition of the model.
Given a sequence λ := (λ h ) h≤−1 satisfying (6.68), we construct the solution of the beta function equation iteratively in h, starting from h = 0 [where z := (z h ) h≤−1 and, of course, the right side only depends on the components of λ, z of scale index larger or equal to h; note that by definition B z h does not depend on λ 0 ]. We denote this solution by z * (λ). By using the tree expansion of the beta function, we now show that z * (λ) is a Cauchy sequence, differentiable in λ; more precisely, we prove that, for λ 0 fixed, such that |λ 0 | ≤ ε , and λ satisfying (6.68), . Let us prove the first inequality in (6.71), inductively in h. Note that at the first step, h = 0, the inequality is trivially true, simply because z * 0 (λ) = z * −1 (λ) = 0. We assume that |z * h −1 (λ) − z * h (λ)| ≤ C 0 (ε ) 2 2 θh , for all scales h < h ≤ 0, and we want to prove that the same bound holds for h = h. Note that, for ε sufficiently small, the first inequality in (6.71) also implies that |z * h (λ)| ≤ ε ≤ ε, ∀h ≤ h ≤ −1, uniformly in λ. Recall that the definition of B z h is induced by (6.40). The kernel W (h),R 2,0;(ω,ω) can be written as a sum over GN trees, analogous to (6.62), and the contribution associated with each tree can be bounded as in Proposition 2. Therefore, B z h itself can be written as a sum over trees that, by definition, have only endpoints of type λ (and, more precisely, at least two such endpoints): where z * = z * (λ) (recall that B z h depends only on the components of z * with scale index ≥ h, which have already been inductively defined), and B z h (λ, z * ; τ, P, T ) can be bounded in a way analogous to (6.65): (6.74) Here we used the fact that the endpoints are all of type λ, and that their values are bounded as in (6.69). We now split B z h (λ, z * ) as follows: By definition, the difference in square brackets is expressed in terms of a sum over trees that have at least one endpoint on scale 0, while B z h (λ, z * ) λ −1 =0 is a sum over trees that have no endpoints on scale 0. By using the short memory property (see comments after the statement of Proposition 2), we find An important remark is that by rescaling h → h + 1, we can re-express where S is the shift operator, namely, (Sλ) h := λ h−1 , and similarly for Sz * .
In conclusion, . We want to bound the difference in the second line as The beta function B z is O(( ) 2 ) because N ≥ 2 in (6.73), so we have just to get the extra factor 2 θh . The left-hand side can be rewritten as with λ(t) := λ + t(Sλ − λ), and similarly for z * (t). B z h+1 can be written in terms of its tree expansion, see (6.73), so that, when the derivative w.r.t. t acts on it, it can act on the factors λ h (t) associated with the endpoints v of the tree, or on the factors z * h (t) associated with the propagators and with the branches of the tree. If it acts on an endpoint of type λ, whose value is λ h (t), its effect is to replace it by λ h − λ h −1 , which is bounded by ε 2 θh , see (6.68); if it acts on a factor z * h (t), its effect is to multiply the tree value by z * h (λ) − z * h −1 (λ), which is bounded by C 0 (ε ) 2 2 θh , thanks to the inductive hypothesis. Using these facts and the short memory property, (6.79) follows. Putting this together with (6.76), we get the desired bound on z * h−1 (λ) − z * h (λ). The proof of the second inequality (6.71) is completely analogous: it can be proved inductively in h (at the first step is trivially valid, again because z * −1 (λ) ≡ 0), by using the tree representation of the beta function, (6.73), and the short memory property. The details are left to the reader. Remark 13. The limiting value of z * h (λ) as h → −∞, which certainly exists, due to the first of (6.71), only depends on λ −∞ (λ) := λ 0 + h≤0 (λ h−1 − λ h ). This fact follows from the recursive equation for z * , (6.70), from the tree representation of the beta function, (6.73), and from the short memory property.
6.4.2. The solution of the flow equation as a fixed point. Given |λ 0 | ≤ ε and λ satisfying (6.68), we fix z = z * (λ) as described in the previous subsection, and plug it into the flow equations for ν h,ω , a h,ω , b h,ω , λ h : these are coupled equations, whose beta functions are thought of as functions of λ 0 and u, with u as in (6.72). In order to find the desired solution to these flow equations, we first note that the equations for ν h,ω , a h,ω , b h,ω in (6.55) imply that, for k < h ≤ 0, [Clearly, B · j,ω (λ 0 , u) actually depends only on the the components of u on scales larger than j.] If we send k → −∞ in (6.81) and impose the desired condition on the exponential decay of ν h,ω , a h,ω , b h,ω , see (6.67), we get Regarding λ h , we study its flow equation by extracting the first order contribution in (λ 0 , u) from the beta function. By inspection, one verifies that the first order contribution does not depend on u: therefore, we can write whereB λ h is at least of second order in (λ 0 , u) and c λ h can be computed in terms of first order perturbation theory. Note that the GN trees that contribute to it have only a normal endpoint at scale 0, of type RV (−1) . Then, due to the short memory property, for someC > 0. By iterating the beta function equation for λ h , we get: where C λ h = 1 + 0 j=h c λ j . In conclusion, given a sufficiently small λ 0 , we look for initial data ν 0,ω , a 0,ω , b 0,ω , depending on λ 0 , such that the corresponding flow satisfies, for all scales with u satisfying (6.67), (6.68). The system (6.85) will be viewed as a fixed point equation u = T (u) for a map T on the space of sequences 86) see (6.67), (6.68). In this equation, and from now on, we let ε be sufficiently small, and we fix θ = 1/2 and ε = ε/K, K = max 1, where C 1 , C 1 are the constants in (6.90) and (6.96) below, whose explicit values can be computed in terms of the first order contributions in λ to B ν h,ω , B a h,ω , B b h,ω . We now want to prove that T is a contraction on X ε , with respect to the metric d(u, u ) := u − u , where More precisely, we intend to prove that the image of X ε under the action of T is contained in X ε , and that T (u) − T (u ) ≤ (1/2) u − u for all u, u ∈ X ε . If this is the case, then T admits a unique fixed point in X ε , which corresponds to the desired initial data ν 0,ω , a 0,ω , b 0,ω , generating a flow satisfying conditions (1)-(2) discussed at the beginning of this section.

6.4.3.
Invariance of X ε under the action of T . In this subsection we show that T (X ε ) ⊆ X ε , i.e. T (u) ≤ ε under the condition that whereC is the same as in (6.83). Note that, in order to prove that T (X ε ) ⊆ X ε , it is enough to show that for some K-independent constants C 1 , C 2 (in order to see that (6.90)-(6.91) imply T (u) ≤ ε, it is enough to plug them in the right side of (6.85) and use (6.87) and (6.89)).
1. The bound on B ν h,ω (λ 0 , u). We start by proving the bound on B ν h,ω (λ 0 , u) in (6.90). Recall that the definition of B ν h,ω is induced by the first of (6.47), combined with the first of (6.55). As for the case of B z h discussed in Section 6.4.1, B ν h can be written as a sum over trees: and B ν h,ω (λ 0 , u; τ, P, T ) can be bounded in a way analogous to (6.65): |F v | × (6.93) × v not e.p.
if v is of type ν, a, b, or λ, see item (2) in the list of properties of trees in Section 6.3; if v is of type RV (−1) , then F v is a kernel of RV (−1) (the one associated with the given choice of P v ), and |F v | is its norm (6.30) (more precisely, it is its un-weighted counter-part, i.e., the case κ = 0), which is bounded as in (6.29). We now split B ν h,ω in a dominant plus a subdominant contribution, in the same spirit as the decomposition (6.39): includes the sum over the trees whose endpoints are all of type λ and all the single-scale propagators g (k) ω /Z k−1 have been replaced by g (k) R,ω /Z k−1 , see (6.37); B ν,s h,ω is the remainder, which includes the sum over trees that have at least one endpoint of type a, b, ν or RV (−1) (the scale k of the endpoints of type λ, a, b, ν satisfies h < k ≤ 0, while the scale of the endpoints of type RV (−1) is necessarily k = 0), or at least one 'remainder propagator' on some scale k between h and 0, (g The key observation is that B ν,R h,ω = 0: in fact the definition of B ν,R h,ω is induced by the first of (6.47), withŴ (h),∞ 2,0;(ω,ω) (0) replaced byŴ (h),R 2,0;(ω,ω) (0) which is zero, as follows immediately from (6.38).
The subdominant contribution, B ν,s h,ω , is not zero, but it is easy to bound. We further distinguish various contributions to it. (1) Let us start with the contributions from trees with at least two endpoints, one of which is of type ν, a, b, RV (−1) and is on scale k ∈ [h + 1, 0]: these are bounded proportionally to ε 2 2 θ (h−k) 2 θk , where θ = 3/4 > θ; the factor 2 θ (h−k) is due to the short memory property (see the comment after (6.66)), while the factor ε2 θk comes from the norm |F v | associated with the endpoint of type ν, a, b, RV (−1) on scale k, and the other ε from the second endpoint. Summing the bound over k in [h + 1, 0], we get const × ε 2 2 θh , with the constant being independent of K. (2) Next, we consider the contributions from trees with exactly one endpoint, of type RV (−1) (and, therefore, on scale 0). Recalling that the norm of the value of the endpoint, |F v |, is proportional to |λ 0 | ≤ ε/(2K), we find that the total contribution from these trees is O(εK −1 2 θh ), the factor 2 θh coming from the short memory property, the proportionality factor being independent of K. (3) Finally, we are left with the contributions from trees whose endpoints are all of type λ and at least one remainder propagator on some scale k between h and 0, (g R,ω )/Z k−1 . Recalling from Remark 11 that the dimensional bound of the remainder propagator is better by a factor 2 θk , as compared to the bound of g (k) R,ω /Z k−1 , we find that the these contributions are bounded by const × (ε/K) 0 k=h 2 θ (h−k) 2 θk ≤ const × (ε/K)2 θh (once again, the factor 2 θ (h−k) is due to the short memory property, and the constant is independent of K). Putting things together, we obtain the desired estimate on B ν h,ω . R,ω )/Z k−1 on some scale k between h and 0. By proceeding as in the previous item, in particular in the discussion of the bound on B ν,s h,ω , we find that |B a h,ω | ≤ const × (ε/K) 0 k=h 2 θ (h−k) 2 θk , which is the desired estimate, and similarly for |B b h,ω |.
is a sum of trees with two or more endpoints, given that we have extracted the first-order contribution c λ h λ 0 . The non-trivial issue is to show that the bound is proportional to 2 θh . For this purpose, we splitB λ h into a dominant and a subdominant part, following once again the same logic: we writẽ R,ω /Z k−1 , see (6.37); B λ,s h is the remainder, which includes the sum over trees that have at least one endpoint of type a, b, ν or RV (−1) , or at least one 'remainder propagator' on some scale k between h and 0, (g The key observation is that the dominant term, B λ,R h = B λ,R h (λ) is the same as the one of the reference model discussed in Section 4: by this, we mean that B λ,R h (λ) is the same that we would get in the reference model, by applying the same multi-scale integration procedure. On the other hand, for the continuum model it is known that, if we denote by λ * 1 the constant sequence (λ * 1) h ≡ λ * , then for λ * sufficiently small, see [7,Theorem 3.1]. Moreover, once the bound (6.94) is known for the beta function computed on the constant sequence λ * 1, by using the short memory property, we find that the same bound holds for more general sequences: more precisely, we find that B λ,R h (λ) ≤ (const.)ε 2 2 θh , for any Cauchy sequence λ satisfying λ θ ≤ ε, as desired.
We are left with the subdominant term, B λ,s h (λ), that, non surprisingly, can be bounded in a way similar to the subdominant contribution B ν,s h,ω ; the result is, once again, |B λ,s h (λ)| ≤ (const.)ε 2 2 θh (details left to the reader). This concludes the proof of (6.90)-(6.91) and, therefore, that T (X ε ) ⊂ X ε .

6.4.4.
T is a contraction on X ε . We now show that T (u) − T (u ) ≤ (1/2) u − u , for all pairs of sequences u, u ∈ X ε (here · is the norm in (6.88)). We consider the component at scale h of T (u) − T (u ), In order to prove that T is a contraction, it is enough to show that, if u, u ∈ X ε and λ 0 satisfies (6.89), then the analogues of (6.90)-(6.91) hold, namely: for some K-independent constants C 1 , C 2 . The proof of (6.96) is very similar to the one of (6.90)-(6.91): in order to illustrate the main ideas, let us focus on B ν h,ω (u) − B ν h,ω (u ), the other components being treatable in a similar manner. By using (6.92), we rewrite the difference under consideration as a sum over trees: We further rewrite the difference in parentheses in the right side in a way similar to (6.80): with u(t) := u + t(u − u ). When the derivative w.r.t. t acts on B ν h,ω (u(t); t, P, T ), it can act on the modified running coupling constants ν h ,ω (t), a h ,ω (t), b h ,ω (t), λ h (t) associated with the endpoints v of the tree, or on the modified constants z * h (λ(t)) associated with the propagators and with the branches of the tree. If, e.g., it acts on an endpoint v of type ν, which is associated with ν h ,ω (t), its effect is to replace it by ν h ,ω − ν h ,ω ; when bounding the norm of the tree value, the endpoint v comes with a factor |ν h ,ω − ν h ,ω |, which leads to a factor u − u ; this has to be compared with the 'standard' factor |ν h ,ω | appearing in the bound of the un-modified tree value, which led to a factor u ≤ ε in (one of the contributions to) the first of (6.90): therefore, the bound on d dt B ν h,ω (u(t); τ, P, T ) is qualitatively the same as the one on B ν h,ω (u; τ, P, T ), up to an additional factor u − u /ε. The terms in which the derivative w.r.t. t acts on other RCCs, or on z * h ((t)) are treated similarly. In light of these considerations, recalling the bound |B ν h,ω (u)| ≤ C 1 K −1 ε2 θh on the un-modified ν-component of the beta function, we obtain the bound in the first line of (6.96) with # = ν. The other components are treated similarly, but we do not belabor further details here.
6.4.6. The flow of Z h and its critical exponent η. Once the initial data are fixed as in (6.99) and the corresponding flow of RCC is bounded and exponentially convergent, we immediately find that , as easily follows from (6.71). Note that, while z −∞ depends only on λ −∞ , A −∞ depends on the whole sequence. For future reference, we let η = η(λ) = log 2 (1 + z −∞ (λ)) be the so-called critical exponent of the wave function renormalization. Then, (ω,ω) and, therefore, it is the same as we would get in a multiscale expansion of the reference model of Section 4: as a consequence, the critical exponent η(λ) is the same as the one of the reference model, η R (λ ∞ ), provided that the bare coupling λ ∞ of the reference model is fixed in such a way that the infrared limit λ −∞,R = λ −∞,R (λ ∞ ) of its coupling equals that of the dimer model, This equation is analytically invertible w.r.t. λ ∞ , as one can show by repeating the study of the flow of λ h for the reference model: in that case, λ h,R satisfies the analogue of the fourth equation in (6.85), which reads where B λ h,R is given by a convergent tree expansion, and satisfies |B λ h,R (λ ∞ , u R )| ≤ (const.)|λ ∞ | 2 2 θh . With respect to the dimer model (cf. (6.82)), note that there is no linear term in the beta function of λ: this is because the interaction potential of the continuum model is exactly quartic in the Grassmann fields. From this, one finds that λ −∞ = λ ∞ + f λ,R (λ ∞ ), where f λ,R is analytic in λ ∞ and of second order in λ ∞ ; in particular, λ −∞,R (λ ∞ ) is analytically invertible with respect to λ ∞ . In conclusion, (6.102) can be inverted into λ ∞ = λ −1 −∞,R λ −∞ (λ) ; this choice guarantees that the asymptotic couplings as h → −∞ of the dimer and reference models are the same. Finally, by inspection of second order perturbation theory, it turns out [9, Th. 2] that η R (λ ∞ ) = aλ 2 ∞ + O(λ 3 ∞ ), for a suitable a > 0. Therefore, by fixing λ ∞ as in (6.102) and recalling that 6.4.7. The flow of Y h,r,ω . On scale −1, one sees by direct inspection of the non-interacting dimer model that Y −1,r,(ω 1 ,ω 2 ) := −K r e −ip ω 2 ·vr +O(λ). Once the fixed point sequence u has been determined, we can plug it into the beta function equation for Y h,r,ω , where Y r = (Y h,r,ω ) h≤−1, ω∈{±} 2 . Note that, by definition, B Y h,r,ω is linear in Y r . This equation can be solved iteratively in h, via a procedure analogous to the one used to compute z given λ, see Subsection 6.4.1. In particular, B Y h,r,ω admits a tree expansion, by using which (6.103) can be rewritten as h,k;ω,ω is the relativistic contribution, i.e., it is expressed as a sum over trees whose endpoints are all of type λ and all the propagators have been replaced by relativistic ones (it is easy to check, by inspection, that B Y,R h,k;ω,ω is independent of r), and B Y,s h,k;r,ω,ω is the remainder. Thanks to the short memory property, and the known bounds on the components of the fixed point sequence u, we find that We now let y h,r,ω = Y h−1,r,ω /Y h,r,ω − 1, and iteratively compute y h,r,ω for h ≤ −1 from (6.104), starting from h = −1. Proceeding by induction, as in subsection 6.4.1, we find that y h,r,ω is a Cauchy sequence, whose limiting value as h → −∞, y −∞,ω (λ), only depends on λ −∞ (λ), see Remark 13. This limiting value defines new critical exponents, η ω (λ) := log 2 (1 + y −∞,ω (λ)). By using the same considerations in Remark 14, we conclude that η ω (λ) are the same as the corresponding exponents in the continuum model, provided that the bare coupling λ ∞ is fixed in such a way that λ −∞,R (λ ∞ ) = λ −∞ (λ). Thanks to the symmetries of the reference model, it is known that η (ω 1 ,ω 2 ) (λ) are real and only depend on the product ω 1 ω 2 ; we denote by η 1 (λ), resp. η 2 (λ), the critical exponent with ω 1 = −ω 2 , resp. ω 1 = ω 2 . Remarkably, it is known also that η 2 (λ) = η(λ), see [9,Theorem 3]. On the other hand, an explicit computation shows that Remark 14). In terms of these critical exponents, we can rewrite Y h,r,ω in a way analogous to (6.101) for suitable complex constants B h,r,ω , C h,r,ω , such that B h,r,−ω = B * h,r,ω and C h,r,−ω = C * h,r,ω . The critical exponent ν of Theorems 1 and 2 is given in terms of η(λ), η 1 (λ) by the simple relation ν(λ) = 1 + η(λ) − η 1 (λ).
(6.107) 6.5. Thermodynamic limit for the correlation functions. In the previous sections, we have obtained a convergent expansion for the effective potentials, valid for |λ| small enough and a suitable choice ofp ω =p ω (λ), α ω =ᾱ ω (λ),β ω =β ω (λ). In particular, after the integration of all the scales h > h L we obtain 9 from (6.31 where we defined are given by the convergent tree expansion discussed above. In order to obtain the Grassmann generating function with θ boundary conditions, W (θ) L (A, 0), we need to integrate out Ψ, see (6.11); finally, the dimer generating function is obtained by taking a linear combination of e W (θ) L (A,0) , see (3.28). Using (6.108) we write: In order to study the thermodynamic limit for correlations, it is important to characterize how E (h L ) , S (h L ) (J) and V (h L ) (Ψ, J) depend on the system size L and on the boundary conditions θ. For illustrative purposes, we start by considering the case A = J = 0, in which case the generating function reduces to the partition function. As shown in Appendix C, Z θ := e W (θ) can be rewritten as where: ∆ is analytic in λ, independent of L and of the boundary conditions; s θ (λ) depends on L, θ and is of order O(λ), uniformly in L, θ; , uniformly in L, θ. From now on, for lightness of notation, we drop the argument λ in u θ,ω (λ), v θ (λ), s θ (λ). The integration of Ψ is elementary, and gives (recall (6.12)) Recall that hL is the first scale at which the support ofχ h has empty intersection with P ω , so that (cf. (6.33)) at scale hL one can remove in (6.31) the integration w.r.t. P (≤h L ) and just replace ψ with 0.
or, equivalently, We now let θ 0 be the boundary conditions for which k ω θ is at the largest distance fromp ω ; if L is large enough, and for a suitable L-independent constant c − . Moreover, for a suitable L-independent constant c + , for all choices of θ, θ (see Appendix D.1). We now multiply and divide the termZ 0 θ L −2 σ θ in (6.113) by Z 0 θ 0 and rewrite it as By using (6.116)-(6.118), we immediately conclude that σ θ,θ 0 = O(λ), uniformly in L, θ. If we now take the appropriate linear combination of Z θ , we obtain the partition function of the interacting dimer model that, in light of the previous considerations, can be written as We now let Q 0 L = 1 2 θ c θ Z 0 θ ; we recall that the constants c θ are either 1 or −1, depending on θ and on the parity of L/2, see the definition after (3.1). A simple computation shows that If we use these inequalities in (6.120), we get Z L = e L 2 ∆(λ) Q 0 L (1 + r L (λ)), (6.123) where the error term r L (λ) is of order O(λ), uniformly in L.
Let us now adapt the previous discussion to e W (θ) L (A,0) , in the presence of the external field A. In this case, the analog of (6.112) (the proof being based on a similar reasoning, see also Appendices B and C) is for some C > 0. Moreover, in the second line of (6.124), V (h L ) (Ψ, J) admits the following explicit expression: or, equivalently, using thatĴ p,r = y∈Λ J y,r e ipy , for suitable translationally invariant functionsw θ,L m;r (y) (the summand with m = 0 should be interpreted as being equal to σ θ ), such that Finally, we take the appropriate linear combination of e W (θ) L (A,0) in order to obtain the generating function of the interacting dimer model: where (recall (6.123)) .
(6.132) By using the properties described above for S θ (J) andw θ,L m;r (y), it is easy to see thatS L (J) admits a representation analogous to (6.127), with w θ,L m;r replaced by a modified kernelw L m;r , satisfying the same estimate as (6.128). Thanks to these estimates, and to the explicit form of S L (J), we conclude, as desired, that the thermodynamic limit of the correlations of the interacting dimer model exist and are given by (we let e i be the edge of type r i with black vertex x i , and we assume the m edges e 1 , . . . e m to be all different from each other): A similar discussion can be repeated for mixed dimer/Grassmann field correlations, but we will not belabor further details here.
6.6. Asymptotic behavior of the dimer correlation functions. In order to complete the proof of our main theorems, we are left with proving that the large distance behaviour of the (thermodynamic limit of the) interacting propagator, vertex function and dimer-dimer correlation can be expressed in term of linear combinations of the appropriate correlations of the reference model, as stated in Section 5. We limit ourselves to the discussion of the asymptotic behaviour of the two-point dimer-dimer correlation, x 2 ), i.e., to the proof of (5.3), and we leave the analogous discussion for the propagator and vertex function (leading to (5.1), (5.2)) to the reader.
We use a strategy analogous to the one of [23, Section 7.1 and 7.2] and we refer the reader to those sections for further details. The starting point is (6.133) with m = 2, which we write as type ν, a, b, RV (−1) , if located at scale h ≤ k ≤ 0, is of the order O(2 θk λ); by the short memory property, we get a factor 2 θ (h−k) , θ < θ < 1 and a sum over k ≥ h finally produces the factor 2 θh in the right side of (6.137).
The contributions with one remainder propagator on scale k ≥ h are treated analogously.
We are left with the dominant parts, S (j),d r 1 ,r 2 ;ω 1 ,ω 2 , j = 1, 2, and to prove (5.3) we need to connect them to the correlation functions of the continuum model of Section 4. Let us fix the coupling constant λ ∞ of the continuum model in such a way that λ −∞;R (λ ∞ ) = λ −∞ (λ), so that the critical exponents η(λ), η 1 (λ) are the same as for the dimer model and one has, in analogy with (6.101), denoting by Z h;R the wave function renormalization of the reference model, withÃ −∞ an analytic function of λ ∞ (and therefore of λ) that equals 1 for λ = 0. The form of the special endpoints is different for the dimer and the continuum model, simply because the external fields J of the continuum model have the structure J (j) x,ω instead of J x,r . In fact, when the multiscale construction is applied to the continuum model, the value of a special endpoint of type j = 1 (density) or j = 2 (mass) is of the form (to be compared with (6.60)) Note that biliner terms such as ν h ,ω F ν;ω (ϕ), a h ,ω F a;ω (ϕ), b h ,ω F b;ω (ϕ) in (6.49) would break the above-mentioned symmetries, unless the coefficients ν h ,ω , etc. are zero. On the other hand, terms such as λ h F λ (ϕ) are invariant because we know by induction that λ h , h > h is real. Next, let us show that (6.38) implies (6.40). Write for lightness of notation ζ ω :=D ω (k), ζ * + = −ζ − . Since k · ∂ kŴ (h),R 2;(+,+) (0) is linear in k, we can write it as f (ζ + , ζ − ) = cζ + + c ζ − for some complex constants c, c . Note that the transformation k → A −1 σ 3 Ak (resp. k → A −1 σ 1 Ak) implies (ζ + , ζ − ) → (−iζ − , iζ + ) (resp. (ζ + , ζ − ) → (−ζ − , −ζ + )). By linearizing the second and third equation in (6.38) we get: which readily imply that c = 0, c ∈ R. This is the desired formula for ω = +. By using the first of (6.38), we get the desired formula for ω = −.

Appendix B. Finite size corrections and bounds on RV (h)
The bounds on the kernels of the effective potential arising in the multiscale procedure, such as Proposition 2, as well as the reason why the action of R, responsible for the factors 2 −z(Pv) in (6.65), makes the renormalized perturbation theory convergent, have been discussed several times in the literature in similar models, see e.g. [23, Section 6.1.4]. In particular, finite-size details have been discussed in [10], but the definition of the L, R operators given there is different from the one proposed in this paper: in [10] the action of L on the kernels of V (h) explicitly depends on the size L of the box, see [10, eq.(2.74)], while in the present case it only depends on the L → ∞ limit of the kernels, see (6.36). This new definition simplifies some technical aspects of the multi-scale construction: for instance, the flow of the running coupling constants is independent of L in the present work. The goal of this appendix is to discuss the modifications induced by the new definition of L, R on the proof of the bounds on the kernels of RV (h) . Familiarity with [23, Sec. 6] is assumed.
For illustrative purposes, we restrict our attention to the part of RV (h) , denoted RV (h) 4 , that is quartic in the Grassmann fields, has no derivative terms∂ϕ, and is independent of J. A similar discussion applies to the terms quadratic in the Grassmann fields (either independent of or linear in J), but we shall not belabor the details here. For the quartic term (using the same notation for the kernel as in (6.57)), we have from (6.36): where the * indicates the constraint that the field and coordinate labels associated with the external fields must match with the prescribed values of ω 1 , . . . , ω 4 , x 1 , . . . , x 4 . Among the various contributions to the right side of (B.2), there are those from trees such that |P v | > 4, for which the action of R on all its vertices v > v 0 (i.e. for all vertices of v that are descendants of v 0 along the tree), is trivial; we letT In this appendix, we limit our discussion to RV  We recall that, if τ ∈T  (ψ)), the action of R on a vertex v1, if non trivial, can interfere with the one on a vertex v < v1 preceding it. Such interference does not cause any conceptual extra difficulty, but it complicates the explicit form of the corresponding tree values must be expressed in an inductive form, rather than by a formula as explicit as (B.4). For a discussion of these issues, see e.g. [ which is the analogue of (6.65) (the factors z(P v ) are absent because R acts on none of the the vertices v > v 0 ; this bound on "non-renormalized trees" has been discussed in several previous papers, see e.g. [20, Section 6]). After summation over τ, P, T , this leads to the bound W (h) 4,0,0;ω κ,h ≤ Cε, uniformly in L: in particular, the bound applies to the kernel of RV On the other hand, by suitably taking into account cancellations between the two terms in the right side of (B.3), one can find an improved bound on RV (h) 4 , which we now discuss. 1. Let us first consider the terms in the right side of (B.3) such that either the argument ofW (h) 4,0,0;ω , (x 1 , x 2 , x 3 , x 4 ), or the argument ofW (h),∞ 4,0,0;ω , (x 1 , x 2 , x 3 , x 4 ), have tree-distance (i.e. length of the shortest tree on Λ including the four points) larger than L/4 (this is the first 'finite size correction' that we intend to discuss in this appendix). Recall that each of the trees contributing to these kernels comes with a a product of propagators 'along the spanning tree', see the factor v not e.p. ∈Tv g (hv) in the right side of (B.4). Therefore, by using the stretched exponential decay of the propagators in (6.54), we find that each of these contributions can be bounded by the right side of (B.5) times an additional, exponentially small, factor e −(κ/4) √ 2 h L = e −(const.)2 (h−h L )/2 , which is the desired dimensional gain.

2.
After having estimated the terms in the previous item, we are left with the terms with tree distance smaller than L/4, which can be rewritten as In the two products in the second line, ε i := (−1) i−1 ; notice that the Grassmann variables in the first product are computed at x i , while in the second product they are computed at x 1 . The term in the second line is the 'usual', infinite volume, renormalized term, which can be treated as discussed in, e.g., [23, Section 6.1.4]; we refer to that section for a discussion of why these terms have the 'usual' dimensional gains leading to the factors 2 −z(Pv) in (6.65). The term in the first line is, instead, the second 'finite size correction' that we intend to discuss in for some constants u 2 , u 4 depending on λ, L, θ. Using the dimensional estimates (see (6.58)), it is easy to deduce that |u n | ≤ C n |λ|2 h L (2−n/2) (C. 3) uniformly in θ, which implies the desired estimates on u θ,ω , v θ , because h L ∼ − log 2 L. Let us now prove (C.1). From the multiscale computation of the effective potential, it follows that where t h was defined in (6.44), andẼ h is the sum of the vacuum diagrams with smallest scale label equal to h, namelỹ which can be represented as a sum over trees, see (6.61)-(6.62). Let us start by discussing the contribution from t h ; using the definition (6.44), we rewrite The term with m = 0, which we denote by t 0,h , is L, θ independent and satisfies |t 0,h | ≤ C|λ|2 2h . (C. 8) To see this, observe that the area of the support ofχ h is O(2 2h ) and recall that r ω (k) = O(k 2 ), that |z h | ≤ C|λ| uniformly in h and that Z h = O(2 −ηh ) (see (6.101)), with η(λ) that tends to zero for λ → 0. The sum of the terms with m = 0, which we denote by t 1,h , is bounded from above as the stretched-exponential decay coming from the fact that the integrand is a function in the Gevrey class of order 2, by assumption onχ h . Finally, recalling thatχ h (k ω θ −p ω ) = 1 for all h > h L and that 1 + z h = Z h−1 /Z h , we find that, if h > h L , the sum in the second line of (C.6) can be rewritten as − 2L −2 log(Z h−1 /Z h ) + t 2,h , |t 2,h | ≤ CL −2 |λ|2 h(1−|η|) . (C.10) Putting things together we write: (C.11) The second term in the right side contributes to ∆(λ): it is L, θ independent and, thanks to (C.8), it is bounded by C|λ|. The term in brackets contributes to L −2 log(1+s θ (λ)): thanks to (C.8) and (C.10), it is bounded by CL −2 |λ|, as we wanted.
We are left with the sum over scales ofẼ h , see (C.4)-(C.5). As mentioned after (C.5),Ẽ h can be written as a sum over trees, where E(τ ), τ ∈ T (h) N,0 , is bounded as in (6.65), with |P ψ v 0 | = |P J v 0 | = |q| = 0, namely (C.13) We now rewriteẼ h as a sum of two terms: the first, which we denote bỹ E 0,h , is the sum over trees of the thermodynamic limit of the tree values (where sums over lattice points in Λ are replaced by sums on Z 2 and singlescale propagators g (h ) ω are replaced by their infinite-volume counterparts g (h ),∞ ). The second is the finite-size remainder, which we denote byẼ 1,h . By construction,Ẽ 0,h is L, θ independent, and it is bounded by the sum over trees of the right side of (C.13), which gives |Ẽ 0,h | ≤ C|λ|2 2h . (C.14) The finite size remainder admits an improved dimensional bound of the form (C. 16) The first term in the right side contributes to ∆(λ): it is L, θ independent and, thanks to (C.14), it is bounded by C|λ|. The term in brackets contributes to L −2 log(1 + s θ (λ)): thanks to (C.15), it is bounded by CL −2 |λ|, as desired. This concludes the proof of (C.1), with the desired bounds on ∆(λ), s θ (λ).