A (2+1)-dimensional Anisotropic KPZ growth model with a smooth phase

Stochastic growth processes in dimension $(2+1)$ were conjectured by D. Wolf, on the basis of renormalization-group arguments, to fall into two distinct universality classes, according to whether the Hessian $H_\rho$ of the speed of growth $v(\rho)$ as a function of the average slope $\rho$ satisfies $\det H_\rho>0$ ("isotropic KPZ class") or $\det H_\rho\le 0$ ("anisotropic KPZ (AKPZ)"class). The former is characterized by strictly positive growth and roughness exponents, while in the AKPZ class fluctuations are logarithmic in time and space. It is natural to ask (a) if one can exhibit interesting growth models with"smooth"stationary states, i.e., with $O(1)$ fluctuations (instead of logarithmically or power-like growing, as in Wolf's picture) and (b) what new phenomena arise when $v(\cdot)$ is not smooth, so that $H_\rho$ is not defined. The two questions are actually related and here we provide an answer to both, in a specific framework. We define a $(2+1)$-dimensional interface growth process, based on the so-called shuffling algorithm for domino tilings. The stationary, non-reversible measures are translation-invariant Gibbs measures on perfect matchings of $\mathbb Z^2$, with $2$-periodic weights. If $\rho\ne0$, fluctuations are known to grow logarithmically in space and to behave like a two-dimensional GFF. We prove that fluctuations grow at most logarithmically in time and that $\det H_\rho<0$: the model belongs to the AKPZ class. When $\rho=0$, instead, the stationary state is"smooth", with correlations uniformly bounded in space and time; correspondingly, $v(\cdot)$ is not differentiable at $\rho=0$ and we extract the singularity of the eigenvalues of $H_\rho$ for $\rho\sim 0$.


Introduction
A statistical physicist's view of stochastic interface growth models is that they describe the diverse phenomena of interface growth and crystal deposition [1]. These models are related mathematically to both interacting particle systems [24,30] and stochastic partial differential equations, in particular, the celebrated Kardar-Parisi-Zhang (KPZ) equation [21].
A (2 + 1)-dimensional stochastic growth process describes the time evolution, in three-dimensional space, of a two-dimensional interface, modeled mathematically as the graph of a function h from Z 2 (or some other twodimensional lattice) to R. The dynamics form an irreversible Markov chain and the transition rates, by which the interface height increases or decreases, are asymmetric and induce a non-trivial drift, that depends on the average local slope. Among the most interesting and challenging questions are those of describing the stationary states, of obtaining the deterministic PDE that describes on large scales the typical evolution of the height function (hydrodynamic limit) and of understanding the large-scale space-time structure of the height fluctuations.
In many natural cases, given a slope ρ ∈ R 2 , there exists a unique stationary state µ ρ with average slope ρ, i.e. µ ρ (h(x) − h(y)) = ρ · (x − y). What is stationary is actually only the law of the height gradients {h(x)−h(y)} x,y∈Z d , while the average height itself grows linearly in time, with a certain speed of growth v(ρ). The PDE describing the hydrodynamic limit is then expected to be of the Hamilton-Jacobi type with ψ = ψ(x, t), t ≥ 0, x ∈ R d representing the suitably rescaled height profile, ∇ denoting gradient w.r.t. the space variable and the solution of the PDE being interpreted in the vanishing viscosity sense.
As far as fluctuations are concerned, one usually introduces a roughness exponent α and a growth exponent β, that measure how height fluctuations grow in space and time in the stationary process, namely and vanishing exponents usually meaning logarithmic growth. (The exponents can also be negative, in which case the constant 1 dominates asymptotically; this is the case for the stochastic heat equation with additive noise in dimension d ≥ 3, see below.) Heuristically, one expects height fluctuations in the stationary process to be qualitatively described, on large space-time scales, by a stochastic PDE (KPZ equation) of the type ∂ t φ(x, t) = ∆φ(x, t) + ∇φ(x, t), H ρ ∇φ(x, t) + ξ(x, t), (1.4) where the diffusive Laplacian term tends to locally smooth out fluctuations, the 2 × 2 symmetric matrix H ρ in the non-linear term is the Hessian of the function v(·) computed at ρ and ξ(x, t) is a space-time noise that contains the randomness of the process. Equation (1.4) is singular if ξ is space-time white noise (Hairer's theory of regularity structures gives a meaning to it in space dimension d = 1 [18], where H ρ reduces to the scalar quantity v (ρ)). On the other hand we are interested in properties on large space-time scales and lattice models have a natural "ultraviolet" space cut-off (lattice spacing), so we can as well assume that the noise is smooth in space and correlated over distances of order 1.
Most of the known rigorous results on stochastic growth models have been proven in dimension (1 + 1), often in cases where the stationary measures of the interface gradients are of i.i.d. type. It is then expected (and mathematically proved in many examples, cf. e.g. [9,13,29] for recent reviews on the mathematical aspects of the KPZ equation and universality class in one dimension) that when v (ρ) = 0 the exponents α, β are given by α d=1 = 1/2, β d=1 = 1/3, as already predicted in the original work [21]. These exponents differ from those (α SHE = (2 − d)/2, β SHE = (2 − d)/4) of the linear stochastic heat equation with additive noise, obtained by dropping the non-linear term in (1.4). On the other hand, for (d + 1)-dimensional models, d ≥ 3, renormalization-group computations [1,21] applied to the stochastic PDE (1.4) predict that, if the non-linear term is sufficiently small, then it is irrelevant, meaning that the large-scale fluctuation properties and exponents α, β are the same as those of the stochastic heat equation (cf. [16,27] for recent mathematical progress in the λ 1 regime; let us also add that [21] predicts a transition at a critical non-linearity λ c = 0 but nothing is known rigorously).
In the (2 + 1)-dimensional case we consider here, the picture is different. On the basis of a renormalization-group analysis of (1.4) by D. Wolf [34], of numerical simulations [19,20,31] and on a few mathematically treatable models (see references in point (ii) below), the following conjectural picture has emerged (see also the introduction of [32] for a more detailed discussion): (i) If det H ρ is strictly positive, i.e. if v(·) is convex or concave with non-zero Gaussian curvature, the growth model is said to belong to the "KPZ (or Isotropic KPZ) class". In this case, the exponents α, β are strictly positive and, in particular, different from those of the two-dimensional stochastic heat equation, which are both zero.
The actual values of the exponents are known numerically to a high degree of precision [19,31], but we are not aware of any rigorous or even convincing arguments to predict them. (ii) If instead det H ρ ≤ 0 ("Anisotropic KPZ" or AKPZ class) one expects that α = β = 0 and that moreover the growth of (1.2), (1.3) is logarithmic, exactly like for the two-dimensional stochastic heat equation. This belief is supported by the mathematical analysis of various (2 + 1)-dimensional growth models [2,5,7,33] that share the following features: stationary states can be found explicitly and their height fluctuations behave on large space scales as a massless Gaussian Field (GFF), with α = 0 and logarithmic correlations; the speed of growth v(·) can be computed and det H ρ turns out to be negative; the growth exponent β is zero and the height variance grows at most logarithmically with time. This state of affairs naturally leads to two questions. The first is to understand what happens when the function v(·) is not differentiable, so that the Hessian matrix H ρ in (1.4) is ill-defined. Can singularities of v(·) lead to qualitatively different dynamic phenomena? Secondly, Wolf's picture predicts rough stationary states, i.e. such that the l.h.s. of (1.2) diverges at larges distances, either logarithmically (AKPZ class) or as a power law (KPZ class). On the other hand, it is well known that certain equilibrium, two-dimensional random interface models exhibit smooth phases, i.e., for some choices of the slope the l.h.s. of (1.2) is bounded uniformly in x, y. Classical example are Solid-on-Solid interfaces and +/− interfaces of the three-dimensional Ising model with Dobrushin boundary conditions, for slope ρ = 0 and sufficiently low temperature [3,4,10]. Is it possible to have non-trivial AKPZ growth models with a smooth phase and, if yes, how does this fit in Wolf's conjectured picture? (See also [25,Sec. 7] for a related phenomenon in a one-dimensional growth model and for the discussion of the corresponding "faceting transition").
In this work, through the analysis of an explicit model, we show that the two questions are closely related. The growth process we study is based on the so-called shuffling algorithm [12,28], that was devised to perfectly sample domino tilings of the Aztec diamond. Instead of working on the Aztec diamond, the dynamics are defined on the (toroidal) periodized lattice (Z/2LZ) 2 and eventually on the infinite lattice Z 2 , a framework that is more relevant for growth processes. The content of Theorem 3.4 is that the Gibbs measures [23] for the dimer model on Z 2 , with 2-periodic weights are translation invariant, stationary, non-reversible measures of the dynamics. The reason behind this choice of weights is that this is the simplest and best-studied [6,11] situation that admits a smooth phase: the Gibbs measure has O(1) fluctuations if the slope ρ is zero. (As a side remark, what we call here "rough" and "smooth" phases correspond to "liquid" and "gas" phases in the language of [23]). As far as growth of fluctuations in time, Theorem 3.8 yields that the l.h.s. of (1.3) grows at most logarithmically in t in the rough phase ρ = 0, as is typical for models in the AKPZ class, and is O(1) uniformly in time in the smooth phase ρ = 0. Most of the technical work is devoted to actually computing the speed of growth v(ρ). The "implicit" formula for the speed in terms of dimer occupation probabilities is very simple, see Eq. (3.10). On the other hand, in order to analyze its convexity properties and its singularity for ρ → 0, we need a more workable expression. This is the content of Theorem 3.10. Note that, while the speed vanishes for ρ = 0, the stationary measure is non-reversible also in this case.
Even with an explicit expression for v(·) like Eq. (3.17) at hand, formulas for its second derivatives are too complicated to check directly whether the sign of det H ρ is negative, as suggested by Wolf's picture. To bypass these difficulties, we found an equivalent but more complex-analytic type expression for v(·), see Theorem 3.12, that allows to prove det H ρ < 0 for every ρ = 0. Finally, for ρ → 0 we show that one eigenvalue of H ρ tends to −∞ and the other tends to 0, while their product tends to a negative constant.
It is natural to ask whether the latter behavior is common to other AKPZ growth models, when the slope approaches that of a smooth phase. The shuffling algorithm on Z 2 , with weights of period larger than 2, may provide examples where such conjecture could be tested, provided that the face weights evolve periodically with the shuffling algorithm (this is not true in general, see discussion in the next paragraph). Let us remark that if the space periodicity of the graph is high enough, there may exist several smooth stationary phases for different integer values of the slope ρ ∈ Z 2 [23]. It is presently unclear how to extend to this more general case the results of Theorems 3.10-3.13. The main obstacle is that part of our approach relies on an explicit computation of the speed of growth, a task which seems, at least to us, to increase in complexity as the periodicity increases.
To conclude, let us mention an intriguing connection between the growth process we study and a discrete dynamical system introduced by Goncharov and Kenyon [14]. Put simply, one of the results of [14] is that the shuffling algorithm on the torus, seen as a map on weighted periodic graphs, behaves similar to a classical Hamiltonian integrable dynamical system. More precisely, some quantities are conserved, such as the spectral curve of the associated dimer model, while (introducing a suitable algebraic structure) others evolve quasi-periodically in time, such as the face weights. The toroidal square grid with two-periodic weights, that we consider here, belongs to a special class of graphs where the evolution of the face weights is not just quasi-periodic but actually periodic in time.
The rest of the paper is organized as follows: in Section 2 we give the necessary notation for the paper while in Section 3, we describe the shuffling dynamics on the torus and give our main results. Proof that the dynamics are stationary is given in Section 4 and the proof of the formula for the speed (in terms of edge probabilities) and of the bound on fluctuation growth is given in Section 5. The explicit formula for the speed of growth for the twoperiodic weights is found in Section 6 and in Section 7 it is shown that the model is in the AKPZ class when ρ = 0.

Notations
We let T L be the square lattice Z 2 periodized, both horizontally and vertically, with period 2L. The graph T L is bipartite and the (2L) 2 vertices are alternately colored black and white. Faces of T L with a white top-right vertex are called "even", the others are called "odd". Denote Ω L to be the set of dimer coverings (perfect matchings) of T L . Let C 1 (resp. C 2 ) be a horizontal (resp. vertical) nearest-neighbor, oriented closed path of length 2L on faces, directed to the right (resp. upward). Given η ∈ Ω L , we let for i = 1, 2 where the sum is over edges crossed by the path, σ e is +1 if the edge e is crossed with the white vertex on the right and −1 if it is crossed with the white vertex on the left and 1 e∈η is the indicator that e is occupied by a dimer. Correspondingly, given ∆ = (∆ 1 , ∆ 2 ), let The "slope" (∆ 1 (η)/L, ∆ 2 (η)/L) belongs [22] to P, the closed square with vertices (−1, 0), (0, −1), (1, 0), (0, 1).  Edges of T L are assigned a positive weight taking values either 1 or a > 0 (without loss of generality, we let 0 < a ≤ 1). First, assign to each even face a label "1" or "a" in an alternate way, as in Fig. 1. Then, we establish that the weight of an edge is given by the label of the even face it belongs to. Note that weights are 2-periodic, i.e. they are left invariant by a horizontal or vertical translation of even length.
For ∆ := ∆ (L) = (∆ 1 , ∆ 2 ) such that ∆/L in • P (the interior of P), let π (L) ∆ be the Boltzmann-Gibbs measure on Ω L (∆ 1 , ∆ 2 ) with weight proportional to a Na(η) , with N a (η) the number of dimers on edges of type a. It is known [23] that, if then the measure π (L) ∆ (L) has a limit π ρ as L → ∞. The "infinite-volume Gibbs measure" π ρ is a probability measure on Ω Z 2 , the set of perfect matchings of the infinite lattice Z 2 , invariant and ergodic with respect to translations by vectors v ∈ (2Z) 2 , and convergence holds for all local bounded functions f . As discussed in the next section, the average "height function" under π ρ is affine with slope ρ. The limit measure π ρ has determinantal correlation functions, with an explicit kernel (that depends on ρ and a), see Section 6. In the same section we will recall also the following fact: with the nomenclature of [23], the Gibbs measure π ρ corresponds to a "liquid phase" (with powerlaw decaying correlations) if and to a gaseous phase (with exponentially decaying correlations) if a < 1, ρ = 0. In the following, we use the terms "rough" and "smooth" instead of "liquid" and "gaseous".

Dynamics
3.1. Spider moves and shuffling algorithm. The dynamics we study is a version of the shuffling algorithm for the Aztec diamond [8,12,28] but reworked for the torus. Its definition is based on two well-known transformations on weighted graphs, that we recall here in a general context; see also [14]. We will come back to the periodized square lattice with 2-periodic weights in Section 3.2.
Here are the two transformations of a graph with weighted edges into a new graph with new edge weights: (1) (Spider move) Suppose that a graph contains a square face f with positive edge weights a, b, c and d where the labeling is clockwise around the face starting from the topmost horizontal edge. We replace the face f by a smaller square with edge weights A, B, C and D (with the same convention as for f ), and add an edge, with edge weight 1, between each vertex of the smaller square and its original vertex. Call these added edges legs and see Fig.2 for a diagram of this move. Then, set A = c/(ac + bd), B = d/(ac + bd), C = a/(ac + bd) and D = b/(ac + bd). This transformation is called the spider move at face f . (2) (Edge contraction) Given a weighted graph, define a new weighted graph as follows: for any two-valent vertex with the incident edges having weight 1, contract the two incident edges. See Fig. 3. This procedure is called edge contraction.
To introduce the dynamics, we first define a sequence of transformations of the edge weights of the graph T L , built via spider moves and edge contractions. We label the faces using Cartesian coordinates (i, j) (both modulo 2L) so that even faces correspond to i + j mod 2 = 0 while odd faces correspond to i + j mod 2 = 1. For the face (i, j) with i, j ∈ {0, 1, . . . , 2L − 1}, the denotes the edge weight around the face (i, j) and at a time k, where we use the same convention a, b, c and d as given in Fig. 2. Note that w k is completely determined by the edge weights around even (resp. odd) faces; we keep track of the weights around both the even and odd faces for convenience.
The relation between w k and w k+1 is, by definition for k ≥ 0, i + j mod 2 = (k + 1) mod 2 and where we have used [r] 2L = r mod 2L for compactness of notation. Since the weights on T L are fully determined by the edge weights on either the odd faces or the even faces, then it follows that all weights {w k } k≥1 are determined by w 0 . Note also that the change of weights from w k to w k+1 is equivalent to first applying the spider move to every face of T L with the same parity (even or odd) as k and then applying the edge contraction transformation to the graph thus obtained. For k ≥ 0 an integer of even/odd parity σ, we define a random map Ω L η → F k (η) ∈ Ω L through the following three steps, cf.  edges, add two parallel vertical dimers to the face with probability w b i,j;k w d i,j;k /∆ i,j;k or two parallel horizontal dimers with probability w a i,j;k w c i,j;k /∆ i,j;k (the operations are performed independently for each (i, j) and k).
Let π (L) ∆;k denote the probability measure (on Ω L (∆)) of the dimer model on T L with weights w k and height change ∆ = (∆ 1 , ∆ 2 ). The following property is crucial for the study of the growth model defined in next section: The above proposition is implicitly known but we could not find an explicit statement in the literature.

3.2.
The growth process. Here we turn back to the case where the graph T L has 2-periodic edge weights with values 1 and 0 < a ≤ 1 as described in Section 2. Coherently with the conventions of Section 3.1, we call such weighting w 0 .
Let e 1 (resp. e 2 ) be the vector from the face (0, 0) to the face (1, 0) (resp. to the face (0, 1)) and τ n (η) be the translation of the dimer configuration η by n 1 e 1 + n 2 e 2 , n = (n 1 , n 2 ) ∈ Z 2 . Note that e 1 , e 2 are just the usual lengthone vectors pointing to the right and up, respectively. With some abuse of notation, given a face f or edge e we denote τ n (f ), τ n (e) the face/edge translated by n 1 e 1 + n 2 e 2 . Note that if n 1 , n 2 are both even, then the vector n connects two faces of the same parity and two edges of the same parity (i.e. both with white edge on top or both with black edge on top).
The random map T that defines our dynamics is given as follows: In other words, we apply first F 0 , then F 1 and then we shift the configuration one step down and one step to the left on T L . Given an initial condition η(0) ∈ Ω L , we let η(k), k ∈ N denote the configuration at time k, i.e. the result of the application of k i.i.d. copies of the map T to η(0).
We begin with the basic properties of the dynamics: ∆ . It is easy to deduce that stationarity holds also for the infinite dynamics. Indeed, note that the definitions of F 0 , F 1 and T make perfect sense as maps from Ω Z 2 to itself (with Ω Z 2 the set of perfect matchings of the infinite lattice Z 2 ). Then: Theorem 3.4. For every local function g and for any slope ρ in A dimer covering of Z 2 is naturally associated to a height function [22]; we briefly recall this definition and then discuss how the height function evolves under the dynamics. This will provide an interpretation of the Markov chain as a two-dimensional stochastic growth model. Given η ∈ Ω Z 2 one associates a height function h η : f ∈ (Z 2 ) * → h η (f ) (defined on the faces f of the graph) by fixing it as h η (f 0 ) = 0 on a given reference face f 0 , say the even face with coordinates (0, 0), and by establishing that its gradients are given as 1 where C f →f is any nearest-neighbor path from f to f (it is well-known that the r.h.s. of (3.4) is path-independent). By construction of the finite-volume measures π From Theorem 3.4, we know that if η(0) ∼ π ρ , the law of the height gradients of η(k) is time-independent. However, we have not yet determined how the additive constant h η(k) (f 0 ) changes with time from the initial value h η(0) (f 0 ) = 0. Assume from now on that the reference face f 0 is even. To motivate our convention given in Definition 3.6 below, we start from the following observation.
Proposition 3.5. If f 1 , f 2 are two even faces and η ∈ Ω Z 2 , then We make therefore the following choice: This way, the whole height functions of F 0 (η) and of T (η) are completely defined, including the overall additive constant, in terms of that of η.
Remark 3.7. From Proposition 3.5, it follows immediately that (3.8)-(3.9) hold for any even face f if they hold for f 0 ; in other words, Definition 3.6 is actually independent of the choice of the even face f 0 . One can check that (3.8)-(3.9) do not, in general, hold for odd faces. If they did, the height function of η(k) would be a deterministic function of the height function of η(k − 1), which is not the case.

Speed of growth and fluctuations.
Theorem 3.8 (Speed and fluctuations). Let f be any face of Z 2 , k ∈ N and let for lightness of notation h k (·) := h η(k) (·). Then, for ρ ∈ where e 1 is a horizontal edge with an "a" face above it and e 2 is a vertical edge with an "a" face to its left. One has for any face f Remark 3.9. As discussed at the beginning of Section 6, the function v(·) satisfies a number of simple reflection symmetries, namely: As a consequence, v(·) is fully determined by its restriction to the subset The function v(·) can be explicitly computed (this computation is rather involved and takes the whole Section 6): The speed of growth is differentiable, except if the weights are genuinely 2-periodic (i.e. a < 1) and the slope corresponds to that of the smooth phase: Theorem 3.11 (Asymptotic behavior of the speed close to the smooth phase). The function v(·) is C 2 on L (recall definition (2.3)). For a < 1 it is not differentiable at ρ = 0 and one has the following asymptotic expansion for ρ ∼ 0. Assume ρ + > 0, ρ − > 0 and let r := ρ − /ρ + ∈ (0, 1) (the other cases can be obtained by symmetry, see Remark 3.9). Then, as ρ + → 0 (under the assumption a < 1) When a < 1, one eigenvalue of H ρ diverges as ρ → 0, the other tends to zero and det(H ρ ) tends to a finite strictly negative limit.
The explicit expressions of f 1 , f 2 can be found in Section 6.3. From (3.19), one sees that as ρ + → 0, which is non-trivial and nowhere vanishing for c < 1/2. Hence ∂ ρ + v(ρ) depends on r as ρ → 0 and so ∇v(ρ) is not continuous at ρ = 0. As discussed in the introduction, of particular relevance for the large-scale fluctuation properties of the growth model is the sign of the determinant of the Hessian matrix H ρ , a negative or a vanishing determinant indicating that the model belongs to the AKPZ universality class. For a < 1, Theorem 3.11 gives that the sign is negative in a neighborhood of ρ = 0. For a = 1, the formulas for the derivatives of v(ρ) are relatively simple and an explicit computation of the determinant shows that for every ρ ∈ L (cf. Eq. (7.7) below). For a < 1, on the other hand, we see no direct way of proving (3.20), starting from the explicit expression (3.17) (we did try to simplify the resulting expressions using Mathematica).
To circumvent this problem, we first give an equivalent, complex-analytic characterization of the speed of growth. Assume that 0 < ρ − , ρ + < π and define with the branch cut as specified in (6.5). The mapping from Here, we adopt the convention that the argument arg(·) ranges from −π/2 to +3π/2.
From the characterization of v(·) provided by Theorem 3.12, we can prove that the growth model does indeed belong to the AKPZ universality class for ρ ∈ L, coherently with the logarithmic upper bound on growth of fluctuations provided by Eq. (3.12): Theorem 3.13 (AKPZ signature of the speed of growth). For every a ≤ 1 and ρ ∈ L, one has det(H ρ ) < 0. (3.23)

Stationarity of Gibbs measures
In this section, we prove Proposition 3.1 as well as Theorems 3.3 and 3.4.
Proof of Proposition 3.1. We consider the case k even, as the odd case follows from the same argument but interchanging even and odd below. Notice that the whole dimer configuration η is determined by dimers covering edges on the boundary of the even faces and that there is a height change crossing an even face horizontally (resp. vertically) if and only if there is exactly one vertical (resp. horizontal) dimer covering an edge on the boundary of that even face. It then follows immediately that the only step in the definition of F k which has an effect on the height function is the sliding step: since after sliding the single dimer on the boundary of the even face moves to an edge of opposite parity, the height change of F k (η) is the opposite as that of η. Take the torus T L with weights w k and apply the spider move to all the even faces, followed by edge contraction of all the resulting two-valent vertices; see Fig. 5 for a schematic of the transformation for the underlying graph. By construction, one obtains the graph T L with weights w k+1 . Explicitly, the weights around the even face with coordinates (i, j) are given by using the same labeling conventions as above.
−∆;k+1 , we proceed as follows. Call P k the transition matrix from Ω L (∆) to Ω L (−∆) of the random map F k and let w k e be the weight of the edge e as determined by w k . If we prove that, for every η ∈ Ω L (−∆), for some constant N independent of η , then the claim of the Proposition follows and actually N equals the ratio of partition functions of π ∆;k and π −∆;k+1 .
Call A(η ) the set of even faces such that η has no dimer along the boundary, D(η ) the set of even faces such that two parallel boundary edges are covered by dimers of η and R(η ) the set of even faces that are neither in A(η ) nor in D(η ). Note first of all that the set S(η ) := {η such that P k (η, η ) = 0} is the set of configurations such that: (a) for every face in A(η ), η has two parallel dimers (either vertical or horizontal) along the boundary of the face; (b) for every face in D(η ), η has no dimer along the face and (c) for every even face in R(η ), η has a single dimer along its boundary, on the edge that is opposed to the one on which η has a dimer. Secondly, observe that by the definition of the creation step, for all η ∈ S(η ), P k (η, η ) equals the product over (i, j) ∈ D(η ) of where (r i,j (η ), r i,j (η )) equals (a, c) or (b, d) according to whether η has two parallel horizontal or two parallel vertical dimers at even face (i, j). Next, note that for all η ∈ S(η ), one has where (r i,j (η), r i,j (η)) are as above, while p i,j (η) is a, b, c, d according to the position of the unique edge of η along the boundary of the even face (i, j). Multiplying by P k (η, η ) and then summing over η ∈ S(η ) (i.e. over the two possible orientations of dimers of η at faces in A(η )), we see that the l.h.s. of (4.2) equals e∈η w k+1 e ∆ i,j;k (4.4) where the product in the r.h.s. runs over all even faces of T L . We used the fact that, after moving a dimer in the sliding step and changing weights from w k to w k+1 , the weight assigned to the dimer is divided by ∆ i,j;k . In conclusion, (4.2) is proven with N = ∆ i,j;k .
Proof of Theorem 3.3. For the two-periodic weighting, w 0 is determined by We have just to prove that if we translate by −e 1 − e 2 the weights w 2 we obtain the original two-periodic weights w 0 , up possibly to an overall positive prefactor that multiplies all weights and is inessential in the definition of the measure. Remark 4.1. Note also that, modulo an overall multiplicative constant, the weights w 1 correspond to interchanging the positions of "a" and "1" faces in the original 2-periodic weighting w 0 .
As far as ergodicity is concerned, assume that there are n > 1 ergodicity classes within Ω L (∆), i.e. subsets C 1 , . . . C n of Ω L (∆) that are invariant for the dynamics. On the other hand, it is known that any two configurations in Ω L (∆) can be connected by a chain of elementary rotations (i.e. the flip of two parallel adjacent dimers from vertical to horizontal or vice-versa) . It follows that there exist two configurations, η ∈ C i , η ∈ C j with i = j that differ by a single elementary rotation at a face f . Assume first that f is an even face at single faces. Then, apply the transformation T to both η and η . Since in the "delete step" of F 0 the discrepancy at the face f disappears, the two updates can be coupled so that T (η) = T (η ). This contradicts the assumption that η, η belong to two different ergodicity classes. Assume that f is an odd face instead and call Ω η ⊂ C i , Ω η ⊂ C j the subset of configurations from which one can reach η, η via a single application of T . If i = j, then Ω η ∩ Ω η = ∅. Let σ ∈ Ω η and η = T (σ) = τ (−1,−1) [F 1 • F 0 ](σ); note that the two parallel dimers that η has around the face f have been added in the "addition" step of F 1 . On the other hand, the two parallel dimers at f could have been given (in the same addition step) the opposite orientation, with positive probability. In this case, the resulting configuration would be η instead, so that σ ∈ Ω η . This contradicts the assumption that Ω η ∩Ω η = ∅ and also contradicts that there is more than one ergodicity class within Ω L (∆).
Proof of Theorem 3.4. We give only a sketchy proof since the argument is standard. Call U f the support of f . The configuration T (η) restricted to U f depends only on the restriction of η toÛ f , the set of edges within distance d from U f , where d is an integer independent of f . Then, since the measure π (L) Lρ converges locally to π ρ [23], we can couple two random configurations sampled from the two measures in such a way that with high probability (as L → ∞) they coincide inÛ f . Then, (3.3) follows immediately from the stationarity of π (L) ∆ for the finite-volume dynamics.  In [33], one of us proved that the Gibbs measures π ρ of the dimer model on the infinite hexagonal graph are stationary for an irreversible Markov chain where particles can perform unbounded jumps. In that case, deducing stationarity for the infinite system from stationarity on the torus required non-trivial arguments since the dynamics are non local (i.e. in that case it is not true that the configuration at time 1 in a domain U depends only on the initial condition in a configuration-independent neighborhood of U ).

Proof of Theorem 3.8.
Proof of (3.10). Since the process is stationary, it is clear that the second expression in (3.10) is independent of k and of f . We will therefore take k = 1 and choose f to be an even face of type a. From Definition 3.6 and Remark 3.7 we see that the height increase in a step at f is Now we take expectation over η ∼ π ρ and obtain On the other hand, in fact, from Proposition 3.1 and Remark 4.1 we have that, if η ∼ π ρ , then τ (−1,0) [F 0 (η)] has the law π ρ with "a" and "1" faces interchanged and the dimer configuration shifted by −e 1 or, equivalently, has the law of a configuration sampled from π ρ and shifted by +e 2 . Altogether, (3.10) follows.
Remark 5.1. If we assumed that f is an even face of type "1" instead, we would get instead of (3.10) where c 3 (ρ) is the density of horizontal dimers with the 1 face above and c 4 (ρ) the density of vertical dimers with the 1 face on the left. Reassuringly, one can prove (via a change of variables in the integral kernels (6.2) giving the probabilities c i (ρ)) that Proof of (3.12). By stationarity of π ρ , the law of h k (f )−h 0 (f ) is independent of f , so we assume without loss of generality that f is an even face. Let Λ be the collection of the O( 2 ) even faces within distance from f and let We have clearly ≤ 2 V (k) Var Pπ ρ (K Λ (k)) + Var Pπ ρ (K Λ (k)). (5.11) We will prove in a moment the following Lemma 5.2. There exists C such that, for every k ≥ 0, ≥ 1, where σ(·) is as in (3.13).
Given this claim, it is easy to conclude the proof of (3.12). Indeed, we see that (with C i denoting positive constants) A simple inductive argument allows to deduce, for k = , (5.14) (Just check by induction that V (k) ≤ C 2 k 2 2 σ 2 ( ) for k = 0, 1, . . . , , if C 2 ≥ max(1, C 2 1 /4).) Now we are ready to prove (3.12). Write where o(1) tends to zero as u → ∞ uniformly in k, by Tchebyshev's inequality, thanks to (5.14). On the other hand In the second line, π ρ (h k (x) − h k (f )) is actually time-independent. The law of A (k) is time-independent by stationarity of π ρ and moreover: Therefore, by Tchebyshev's inequality, with probability 1+o(1) as u → ∞. Finally, we note from (5.18) that if event (5.20) holds and at the same time |Q Λ k (k) − E πρ Q Λ k (k)| ≤ √ uk 2 σ(k), since |Λ k | grows proportionally to k 2 one cannot have |Q f (k) − E πρ Q f (k)| ≥ uσ(k) for u large. Eq. (3.12) is then proven.
Proof of Lemma 5.2. Given that the sum in the definition of K Λ (k) runs over even faces, we apply (5.4) to write with e x (resp. e x ) the vertical edge connecting the top-right and the bottomright (resp. top-right and top-left) vertices of face x. Since the laws of η(k) and F 0 (η(k)) are stationary, the variance of K (j) Λ (k), j = 1, 2 are independent of k. One has with π(A; B) denoting the covariance of A, B. It is known [23] that dimerdimer correlations decay like the inverse distance square if π ρ is a rough phase (i.e. ρ ∈ L), and exponentially fast if π ρ is a smooth phase. Then, inequality (5.12) for K Proof of Lemma 5.3. This is an immediate consequence of the known fact that, if π ρ is a rough phase (as is the case for ρ = 0 or ρ = 0, a = 1) then the variance Var πρ (h(x) − h(y)) grows proportionally to log |x − y| when |x − y| → ∞, while if π ρ is a smooth phase then the variance is uniformly bounded in x, y [23].
In the rest of the section, we will therefore assume that |ρ 1 | = |ρ 2 |.
6.1. Kasteleyn matrix. The speed of growth v(ρ) is given by the difference between the probabilities of two events. We begin by recalling how to express probabilities of local events via the Kasteleyn matrix and how to rewrite its matrix elements as single integrals in the complex plane, as was done in [6]. We adopt a similar coordinate system to the one used in [6] but we have chosen to interchange the white and black vertices, so that we can keep the same height change conventions from a previous paper [5]. Namely, with reference to Fig. 8, where the graph has been rotated 45 degrees clockwise, we assign to each vertex in Z 2 coordinates x = (x 1 , x 2 ) with x 1 + x 2 ∈ 2Z + 1. Then, the set of black vertices is B = {(x 1 , x 2 ) ∈ Z 2 : x 1 mod 2 = 0, x 2 mod 2 = 1} and the set of white vertices is W = {(y 1 , y 2 ) ∈ Z 2 : y 1 mod 2 = 1, y 2 mod 2 = 0}. For i ∈ {0, 1}, we also let and W i = {(y 1 , y 2 ) ∈ W : y 1 + y 2 mod 4 = 2i + 1}. The fundamental domain of the two-periodic weighting is given by a 2 by 2 block consisting of a vertex from each W 0 , W 1 , B 0 and B 1 . We suppose that if the vertex b ∈ B 0 is in the fundamental domain, then so are the vertices b+ê 2 ∈ W 0 , b+ê 1 ∈ W 1 , and b+ê 1 +ê 2 ∈ B 1 , withê 1 = (1, 1) andê 2 = (−1, 1). We impose that the edges inside the fundamental domain have weight a and the edges crossing to another fundamental domain have weight 1; see Fig. 8. The Kasteleyn matrix K is a weighted version of the adjacency matrix of Z 2 . The matrix element K(y, x) is zero if x ∈ B and y ∈ W are not neighbors. If they are neighbors, then K(y, x) equals the weight (a or 1) of the edge from x to y, times i (the imaginary unit) if the edge is parallel toê 2 . The non-zero values of the matrix elements are therefore a, ia, 1, i.
Given positive real numbers r 1 , r 2 , we introduce the "inverse Kasteleyn matrix" K −1 r 1 ,r 2 as follows. Define first the 2 × 2 matrix and the characteristic polynomial Let x ∈ B ε 1 and y ∈ W ε 2 and let (u, v) ∈ Z 2 be such that the unique fundamental domain containing x (resp. y) has black vertex in B 0 with coordinates X (resp. Y ) satisfying Y − X = 2(uê 1 + vê 2 ). Then, we let where the integrals are taken in the complex plane and Γ r is a contour of radius r around the origin, oriented anti-clockwise.
The link between the Gibbs measure π ρ and the Kasteleyn matrix is as follows [23]. Whenever the slope ρ corresponds to a rough phase, i.e. ρ ∈ L (recall (2.3)) there exists a unique choice r = (r 1 , r 2 ) = r(ρ) (magnetic coordinates) such that, given any integer n and edges e i , i ≤ n with black/white vertices of coordinates x i and y i respectively, one has π ρ (e 1 , . . . , e n ∈ η) 3) The image of the curve P (z, w) = 0 in C 2 under the map (z, w) → (log |z|, log |w|) is called the amoeba of P , often denoted by A(P ). The set B = {r : r = r(ρ) for some ρ ∈ L} is the interior component of the amoeba and the correspondence r ∈ B ↔ ρ(r) ∈ L is a bijection [23]. In the rest of this section, we assume that r ∈ B. Note that this is not a restriction since ρ = 0 has already been taken into account. Both the speed v(ρ) (through (3.10)) and the slope ρ (through (3.4)) can be expressed as linear combinations of probabilities as in the l.h.s. of (6.3), with n = 1. Before working out the explicit expressions, let us recall a few known facts [6] (u 1 , u 2 ) = 2(1 + a 2 ) + a(u 1 + u −1 1 )(u 2 + u −1 2 ). Note thatc(u 1 , u 2 ) = P (u 1 /u 2 , u 1 u 2 ). We next state the following lemma without proof from [6]: Lemma 6.1 (Lemma 4.3 from [6]). For x = (x 1 , x 2 ) ∈ B ε 1 and y = (y 1 , y 2 ) ∈ W ε 2 with ε 1 , ε 2 ∈ {0, 1}, The proof of the above lemma, see [6], follows from a change of variables, which also explains the shift in contours of integration. The following are from [6] and are useful for later computations. Define for w ∈ where the logarithm takes arguments in (−π/2, 3π/2). Write 1/w 2 + 2c for the same function evaluated at 1/w. Set Note that, since c ≤ 1/2, the branch cut i[− √ 2c, √ 2c] is included in [−i, i]. The choice in branch cut also gives that √ w 2 + 2c = − (−w) 2 + 2c. (6.6) Let us note the following, for later use: . (6.7) If z → i √ 2c + 0 + (resp. z → −i √ 2c + 0 + ), then any limit point of The transformation u = G(ω) is a map between |u| < 1 and the whole plane (without the cut) in ω. For u = Re iθ with R < 1, the inverse of this transformation is given by Moreover, when R < 1 and θ runs from 0 to 2π, ω in (6.8) runs (clockwise) over an ellipse, that does not cross the cut. We denote by γ R the ellipse, oriented anti-clockwise. For x = (x 1 , x 2 ) ∈ B ε 1 and y = (y 1 , y 2 ) ∈ W ε 2 define for i ∈ {1, 2} (6.9) Define for two contours γ and γ (that avoid the branch cut) and k, l ∈ Z and also define We state the following lemma without proof from [6]. ). Let R 1 = r 1 /r 2 and R 2 = 1/ √ r 1 r 2 with R 1 , R 2 < 1. Then, we have The proof of this lemma involves using the change of variables u = G(ω) for each variable u 1 and u 2 . Note that the contours γ R 1 , γ R 2 encircle the branch cut i[− √ 2c, √ 2c] of G(ω j ), j = 1, 2.
We have the following lemma which is a culmination and extension of Lemmas 3.3 and 4.5 from [6]. Lemma 6.5. Let r ∈ B, let R 1 , R 2 be defined as in Lemma 6.3 as functions of r and assume R 1 , R 2 < 1. Then, there exists a unique ω c ∈ R >0 × iR >0 such that where the integration path from ω c to ω c (resp. from −ω c to −ω c ) stays in the half-plane with positive (resp. negative) real part, otherwise, (6.14) and R > 1.
Proof of Lemma 6.5. The function in the r.h.s. of (6.10) has a pole whenever ω 1 = 1/ω 2 . We show first of all: Claim 6.6. There exist four distinct values of ω 1 ∈ γ R 1 such that ω 2 := 1/ω 1 ∈ γ R 2 and there is exactly one of them for each of the four quadrants of C \ {R ∪ iR}. We call by convention ω c the one that is in R >0 × iR >0 , the others being −ω c , ω c , −ω c .
Proof. Since each ω i , i = 1, 2 runs along the clockwise-oriented ellipse γ R i defined by (6.8), we have

(6.15)
Set u 1 = ie iθ 1 and u 2 = ie iθ 2 . Then, Set w = 1/(u 1 u 2 ) and z = u 1 /u 2 and note |z| = |w| = 1. This gives Comparing with (6.1) we conclude that 1 − ω 1 ω 2 = 0 is equivalent to having P (zr 1 , wr 2 ) = 0. Recall from Remark 6.4 that z 0 , w 0 are not both real. We also have that they cannot be both purely imaginary, since it is clear from (6.1) that P (z 0 r 1 , w 0 r 2 ) would not be zero in that case. Therefore, we can assume that at least one among z 0 /w 0 and z 0 w 0 is not real, and assume that the former is the case (in the other case, the following argument works the same exchanging the roles of u 1 and u 2 ). Recalling the definitions of z, w in terms of u 1 , u 2 we see that the condition ω 1 ω 2 = 1 implies that u 1 can have only one of the four possible values ± z 0 /w 0 , ± z 0 /w 0 . The four values are distinct and there is one of them in each of the four quadrants of C \ {R ∪ iR}. As a consequence, Claim 6.6 immediately follows.
We resume the proof of Lemma 6.5. Via the change of variables ω 1 → ω −1 1 , whereγ R 1 represents the contour from γ R 1 under the image ω → 1/ω, oriented anti-clockwise. Claim 6.6 indicates thatγ R 1 intersects γ R 2 at four points, ω c , −ω c , ω c , and −ω c . Since R 1 , R 2 < 1, the contourγ R 1 is outside of γ R 2 on the real axis and inside γ R 2 on the imaginary axis. Note also that the branch cut i(−∞, Fixing ω 2 we integrate over ω 1 , by deforming the contourγ R 1 to a circle Γ 1/R with small 1/R < R 2 < 1. The deformation crosses a pole ω 1 = ω 2 if and only if ω 2 is in the portion of γ R 2 from −ω c to ω c or from ω c to −ω c . Fig. 9 shows these deformations.

The integral becomes theñ
where the first line is the contribution from the pole. For the double integral in the above equation, we can deform γ R 2 to Γ R because 1/R < 1 < R, so that the branch cut of G(·) remains inside the contour. We can then apply the change of variables ω 1 → 1/ω 1 and so the double integral becomes , (6.20) where the orientation of Γ R is as usual anti-clockwise. To show that the above integral is equal toS(k, l), notice that G(ω) behaves like √ c/( √ 2ω) as |ω| → ∞ and when k ≥ 0 or l ≥ 0, we can push one or the other contour through infinity to get that the integral is zero. Also when l = k = −1, we can compute the residue at infinity which gives 1/a and hence the double integral in the above equation is equal toS(k, l). Remark 6.7. Note that, as a consequence of arg ω c ∈ (0, π/2), one has arg G(ω c ) ∈ (π/2, π). An easy way to see this is to observe that log G(ω) is analytic in the open quadrant R >0 × iR >0 , so that its imaginary part arg(G(ω)) is harmonic. On the other hand, along the positive real axis arg G(ω) = π, while along the positive imaginary axis, arg G(ω) ∈ [π/2, π], as follows from Remark 6.2. Given that at infinity G(ω) ∼ −c/ω whose argument is also in [π/2, π] and since ω c is in the interior of the quadrant, we deduce the claim. A similar argument gives arg G(ω −1 c ) ∈ (π, 3π/2). The following result reduces by symmetry the case of R 1 , R 2 not both smaller than 1 to the case R 1 , R 2 < 1: Lemma 6.8. Let r ∈ B and assume x = (x 1 , x 2 ) ∈ B ε 1 , y = (y 1 , y 2 ) ∈ W ε 2 with ε 1 , ε 2 ∈ {0, 1}. If R 1 = r 1 /r 2 = 1 and R 2 = 1/ √ r 1 r 2 = 1. Then, where for i ∈ {0, 1}γ (6.22) and Proof. This is immediate from making the change of variables u 1 → u −1 1 or u 2 → u −1 2 depending on whether R 1 > 1 or R 2 > 1 (or both) in the formula for K −1 r 1 ,r 2 given in Lemma 6.1.
Remark 6.9. Expanding the formulas for the slopes given in (6.31), (6.32) in terms of their integral formulas given in Lemma 6.1, shows that R 1 = 1 or R 2 = 1 is equivalent to |ρ 1 | = |ρ 2 |; see Lemma 6.12 below. For the remaining computations for the speed, it is sufficient to consider R 1 < 1 and R 2 < 1 (recall that we already computed the speed for |ρ 1 | = |ρ 2 | at the beginning of this section) and use symmetry.

6.2.
Computation of the speed of growth. We first express the speed as a function of ω c , δ 1 and δ 2 , which are themselves functions of R 1 , R 2 and therefore of r 1 , r 2 , see Lemma 6.10. Then, we find a formula (Lemma 6.11) for the slope ρ also in terms of ω c , δ 1 and δ 2 . Finally, in Lemma 6.13 we put together the two results to obtain the speed as function of the slope and formula (3.17).
Lemma 6.10. For r ∈ B, let v = v(r) be defined as v := v(r) = π ρ(r) (e 1 ∈ η) − π ρ(r) (e 2 ∈ η) (6.24) (compare with Eqs. (3.10) and (3.11)). For R 1 , R 2 , δ 1 , and δ 2 as given in Lemma 6.8 with R 1 , Here we have defined v to be a function of r whereas in (3.10) we have defined it as a function of ρ. It turns out that this is equivalent because the correspondence r ∈ B ↔ ρ(r) ∈ L is a bijection [23], as already mentioned.
We also have using (6.46) and (6.48) Taking the real part of the above equation yields The equations (6.41), (6.42) for arg G(ω c ) and arg G(ω −1 c ), together with angle addition formulas followed by factorizing and simplifying, give the result (6.45).

AKPZ signature of the growth model
In this section we prove Theorem 3.12 and 3.13.
V (z) V (1/z) = α (7.5) for some α > 0. Now we use the explicit form of G(z), which implies that The l.h.s. of (7.5) is real on the positive real axis. One can check from Remark 6.2 that its argument is in [0, π/2] when the imaginary axis is approached from the first quadrant. To be precise, if z tends to iy, y > 0, the argument tends to π/2 if y > 1/ √ 2c or y < √ 2c and to 0 if √ 2c < y < 1/ √ 2c; if y = √ 2c or 1/ √ 2c the limit need not exist but the limit points are in [0, π/2]. Also, for |z| → ∞ or |z| → 0 (with z ∈ Q + ) the argument of the l.h.s. of (7.5) is in (0, π/2) because the ratio is z √ 2c + o(z) and z/ √ 2c + o(z) respectively. Since arg(V (z)/V (1/z)) is a harmonic function in the open set Q + with boundary values in [0, π/2] and not identically 0, it vanishes nowhere in Q + , so (7.5) cannot hold for any α > 0. By the way, the same harmonicity argument gives that the real parts of V (z) or V (1/z) vanish nowhere in Q + , as mentioned above.
Once we know that the map z → (X, Y ) is a local diffeomorphism, we can conclude that it is a global diffeomorphism, via a theorem by Hadamard [15,17], if we prove that it is proper, i.e. for every sequence z n ∈ Q + such that z n → z ∈ ∂Q + , (X(z n ), Y (z n )) tends to a point on the boundary of (π/2, π) × (π, 3π/2) 2 . Here, when we write z n → z ∈ ∂Q + , we mean that either z is purely real, or it is purely imaginary, or |z n | tends to infinity. When z n tends to a real limit, it is clear that arg G(z) tends to π while when |z n | → ∞, arg(G(1/z)) tends to π. When instead z n tends to iy, then Remark 6.2 shows that arg G(z n ) tends to π/2 if y > √ 2c while arg(G(1/z n )) tends to 3π/2 if y < 1/ √ 2c. Since √ 2c ≤ 1/ √ 2c, we proved our claim that in all cases (X(z n ), Y (z n )) approaches the boundary of (π/2, π) × (π, 3π/2).