The domino shuffling algorithm and Anisotropic KPZ stochastic growth

TL;DR

This paper demonstrates that the domino shuffling algorithm models a stochastic growth process in the anisotropic KPZ class, revealing properties of growth speed, fluctuations, and discontinuities at specific slopes, using geometric and probabilistic methods.

Contribution

It establishes the anisotropic KPZ classification for the domino shuffling growth model and analyzes the discontinuities and fluctuation behavior without explicit formulas.

Findings

Growth speed depends on slope and weights, belonging to the anisotropic KPZ class.

Height fluctuations grow at most logarithmically in time.

Discontinuities in growth speed occur at finitely many smooth slopes.

Abstract

The domino-shuffling algorithm can be seen as a stochastic process describing the irreversible growth of a $(2+1)$-dimensional discrete interface. Its stationary speed of growth $v_{\mathtt w}(\rho)$ depends on the average interface slope $\rho$, as well as on the edge weights $\mathtt w$, that are assumed to be periodic in space. We show that this growth model belongs to the Anisotropic KPZ class: one has $\det [D^2 v_{\mathtt w}(\rho)]<0$ and the height fluctuations grow at most logarithmically in time. Moreover, we prove that $D v_{\mathtt w}(\rho)$ is discontinuous at each of the (finitely many) smooth (or "gaseous") slopes $\rho$; at these slopes, fluctuations do not diverge as time grows. For a special case of spatially $2-$periodic weights, analogous results have been recently proven in Chhita-Toninelli (2018) via an explicit computation of $v_{\mathtt w}(\rho)$. In the general case, such a computation is out of reach; instead, our proof goes through a relation between the speed of growth and the limit shape of domino tilings of the Aztec diamond.

Equations

$${\mathtt{w}}_k=\{(w_{i,j;k}^a,w_{i,j;k}^b,w_{i,j;k}^c,w_{i,j;k}^d):(i,j)\in {\mathbb{Z}}^2,(i+j)=k\mathrm{mod}2\}$$

$${\Delta }_{i,j;k}=w_{i,j;k}^a{w}_{i,j;k}^c+w_{i,j;k}^b{w}_{i,j;k}^d.$$

$$\begin{gathered} (w_{i,j;k+1}^a,w_{i,j;k+1}^b,w_{i,j;k+1}^c,w_{i,j;k+1}^d) \\ :=\left(\frac{w_{i,j+1;k}^a}{{\Delta }_{i,j+1;k}},\frac{w_{i+1,j;k}^b}{{\Delta }_{i+1,j;k}},\frac{w_{i,j-1;k}^c}{{\Delta }_{i,j-1;k}},\frac{w_{i-1,j;k}^d}{{\Delta }_{i-1,j;k}}\right) \end{gathered}$$

$$\frac{w_{i,j;k}^b{w}_{i,j;k}^d}{{\Delta }_{i,j;k}}$$

$$\frac{w_{i,j;k}^a{w}_{i,j;k}^c}{{\Delta }_{i,j;k}}$$

$${\eta }_k:=T_k\circ \cdots \circ T_1({\eta }_0),$$

$$h_{\eta }(f^{\prime })-h_{\eta }(f)=\sum _{e\sim C_{f\to f^{\prime }}}{\sigma }_e(\text{mathds}{1}_{e\in \eta }-1/4)$$

$$h_{T_{k+1}(\eta )}(f)-h_{T_{k+1}(\eta )}(f^{\prime })=h_{\eta }(f)-h_{\eta }(f^{\prime }).$$

$$h_{T_{k+1}(\eta )}(f)=h_{\eta }(f).$$

$$h_{T_{k+1}(\eta )}(f)-h_{\eta }(f)=\frac{H[\eta ]+H[T_{k+1}(\eta )]-V[\eta ]-V[T_{k+1}(\eta )]}{4},$$

$$H[\eta ]=\text{mathds}{1}_{e_1}(\eta )+\text{mathds}{1}_{e_3}(\eta )$$

$$V[\eta ]=\text{mathds}{1}_{e_2}(\eta )+\text{mathds}{1}_{e_4}(\eta );$$

$$(K(z,w))_{xy}=\left\{\begin{array}[]{lll}\mathrm{wt}(x,y)z^{a}&\text{if}&\mbox{$(x,y)$ is a horizontal edge,}\\ \mathrm{i}\,\mathrm{wt}(x,y)w^{b}&\text{if}&\mbox{$(x,y)$ is a vertical edge,}\\ 0&\text{if}&\mbox{$(x,y)$ is not an edge}\end{array}\right.$$

$$a=\left\{\begin{array}[]{ll}1&\mbox{if $x=(2n-1,k)$ and $y=(0,k)$}\\ -1&\mbox{if $x=(0,k)$ and $y=(2n-1,k)$}\\ 0&\mbox{otherwise}\end{array}\right.$$

$$b=\left\{\begin{array}[]{ll}1&\mbox{if $x=(l,2n-1)$ and $y=(l,0)$}\\ -1&\mbox{if $x=(l,0)$ and $y=(l,2n-1)$}\\ 0&\mbox{otherwise}\end{array}\right.$$

$$N(P)=\text{mbox}convexhull\{(j,k)\in {\mathbb{Z}}^2|z^j{w}^k\text{mbox}isamonomialinP(z,w)\}\subset {\mathbb{R}}^2.$$

$$\prod _{e\in \eta \cap \Lambda }{\mathtt{w}}_0(e)$$

$$R(B)=\frac{1}{(2\pi \mathrm{i}{)}^2}\int {\int }_{\begin{array}{c}|z|=e^{B_1}\\ |w|=e^{B_2}\end{array}}\log |P(z,w)|\frac{dz}{z}\frac{dw}{w}$$

$$\sigma (\rho )=\underset{B\in {\mathbb{R}}^2}{\sup }(-R(B)+\rho \cdot B).$$

$${\mathrm{Var}}_{{\pi }_{\rho ,{\mathtt{w}}_0}}(h_{\eta }(f)-h_{\eta }(f^{\prime }))\sim \frac{1}{{\pi }^2}\log |f-f^{\prime }|$$

$$\displaystyle\mathcal{S}\subset\Bigl{[}\stackrel{{\scriptstyle\circ}}{{N(P)}}\cap\;\mathbb{Z}^{2}\Bigr{]};$$

$$Q=\{(x_1,x_2):|x_1|+|x_2|\le 1/(2n)\}.$$

$${\hat{h}}_{\eta }(\hat{f}):=\frac{1}{N}{h}_{\eta }(f),$$

$${\hat{h}}_{\eta }(\hat{f})\stackrel{N\to \infty }{=}{\psi }_{\partial Q}(x):=\frac{n}{2}\left(|x_1|-|x_2|\right).$$

$$\displaystyle\lim_{N\to\infty}\pi_{\mathtt{w},N}\Bigl{(}\exists\hat{f}\in\hat{F}_{N}:|\hat{h}_{\eta}(\hat{f})-\psi_{\mathtt{w}}(\hat{f})|>\delta\Bigr{)}=0.$$

$${\int }_Q\sigma (\nabla \psi )dx$$

$$v_{\mathtt{w}}(\rho ):=\underset{k\to \infty }{\lim }\frac{1}{k}\sum _{j=1}^k\left({\mathbb{E}}_{{\pi }_{\rho ,\mathtt{w}}}(h_{{\eta }_j}(f))-{\mathbb{E}}_{{\pi }_{\rho ,\mathtt{w}}}(h_{{\eta }_{j-1}}(f))\right)$$

$$\underset{k\to \infty }{\lim }\frac{1}{k}{\mathbb{E}}_{{\pi }_{\rho ,\mathtt{w}}}(h_{{\eta }_k}(f)-h_{{\eta }_0}(f))=v_{\mathtt{w}}(\rho ).$$

$$v_{\mathtt{w}}(\rho )={\psi }_{\mathtt{w}}(x_{\mathtt{w}}(\rho ))-x_{\mathtt{w}}(\rho )\cdot \rho .$$

$$D{v}_{\mathtt{w}}(\rho )=-x_{\mathtt{w}}(\rho ).$$

Full text

1. Introduction

In the realm of stochastic interface growth [2], dimension $(2+1)$ (i.e., growth of a two-dimensional interface in three-dimensional physical space) plays a distinguished role. In $(1+1)$ dimensions, one finds a non-trivial KPZ growth exponent $\beta=1/3$ as soon as the growth process is genuinely non-linear, while in dimension $(d+1),d\geq 3$ a phase transition is expected [18] between a regime of small non-linearity, where the process behaves qualitatively like the stochastic heat equation (SHE) with additive noise, and a regime of large non-linearity, characterized by new growth and roughness critical exponents. See the recent [25, 12, 11] for mathematical progress on the small non-linearity regime of the KPZ equation for $d\geq 3$. On the other hand, dimension $(2+1)$ is the “critical” or “marginal” case: here, the critical exponents are expected to depend not so much on the intensity of the non-linearity, but rather on its structure. In fact, in this case, the existence of two different universality classes has been conjectured [35, 2] (see [33] for a recent mathematical review). The first, called Anisotropic KPZ (or AKPZ) class, is characterized by logarithmic growth of height fluctuations in space and time, like the two-dimensional SHE with additive noise. The second, called KPZ class tout court, has universal and non-trivial roughness and growth exponents, $\alpha_{KPZ}\simeq 0.39$ and $\beta_{KPZ}\simeq 0.24$ respectively (these values are known only numerically, cf. e.g. [32, 17]). Conjecturally, the universality class of a model is determined by the properties of the average speed of growth $v(\rho)=\lim_{t\to\infty}\frac{1}{t}\mathbb{E}[h(t,x)-h(0,x)]$ of the interface height function $h$, where $\rho$ is the macroscopic slope of the initial condition. Namely, a model is expected to belong to the AKPZ class if and only if $\det(D^{2}v(\rho))\leq 0$, where $D^{2}v(\rho)$ is the $2\times 2$ Hessian matrix. From the mathematical point of view, the understanding of the AKPZ universality class has remarkably progressed lately but it is still limited to a few special cases (see Section 1.1 for references). For the KPZ class, very interesting recent developments (in a somewhat different direction) concern the weak non-linearity (or weak-disorder) regime [6, 5]: if non-linearity is scaled to zero as $\hat{\beta}/\sqrt{|\log\epsilon|}$, with $\epsilon\to 0$ a noise regularization parameter and provided $\hat{\beta}$ is smaller than a precisely identified critical value $\hat{\beta}_{c}$ [5], then the KPZ equation scales to the SHE with additive noise. In this regime, the non-trivial exponents $\alpha_{KPZ},\beta_{KPZ}$ do not emerge.

In the present work, we focus on the so-called “domino shuffling algorithm”. This is a discrete-time Markov chain on perfect matchings (or “domino tilings”) of $\mathbb{Z}^{2}$, that was originally devised [13, 30] as a way to exactly sample and to count perfect matchings of certain special two-dimensional domains (Aztec diamonds). When this algorithm is run on the infinite square grid, it can be seen also as a $(2+1)$-dimensional growth model, and it is from this point of view that we consider it here. The shuffling algorithm is actually an infinite-dimensional family of growth processes, indexed by the edge weights $\mathtt{w}$, that we only assume to be positive and periodic in both lattice directions, with some period $2n\in 2\mathbb{N}$. Along the dynamics, the edge weights also evolve (deterministically) in time. In fact, the evolution $\{\mathtt{w}_{k}\}_{k\geq 0}$ of edge weights under the shuffling algorithm (or “spider moves”) has a a remarkable interest in itself, as a classical integrable dynamical system [16]. Its trajectories are in general not time-periodic.

For generic edge weights of period $2n$, there are $2n(n-1)+1$ special values for the slope (“smooth” or “gaseous” slopes), that correspond to “cusps” of the surface free energy $\sigma(\rho)$ of domino tilings with weights $\mathtt{w}$. The slopes at which $\sigma$ is smooth are instead referred to as “rough slopes” (the reason for the nomenclature smooth/rough is reminded in Section 2.2). We let $\mathcal{S}$ (resp. $\mathcal{R}$) denote the set of smooth (resp. rough) slopes.

Our main result is that the domino shuffling algorithm (with general weights $\mathtt{w}$) belongs to the AKPZ class, and that the speed of growth is singular at each of the smooth slopes (see Theorem 2.3 and Section 2.4.1 for more precise statements):

** Main Theorem**** (Informal version).**

For $\rho\in\mathcal{R}$, the speed of growth function $\rho\mapsto v_{\mathtt{w}}(\rho)$ is $C^{\infty}$ and $\det[D^{2}v_{\mathtt{w}}(\rho)]<0$. On the other hand, the gradient $Dv_{\mathtt{w}}(\rho)$ is discontinuous at each of the finitely many slopes $\rho\in\mathcal{S}$. For $\rho\in\mathcal{R}$, the height fluctuations grow logarithmically in space (they scale to a Gaussian Free Field) and at most logarithmically in time. For $\rho\in\mathcal{S}$, the variance of the height fluctuations is uniformly bounded in space and time.

In a special case of $2$-periodic weights ($n=1$) analogous results have been proven recently in [8]. In that case, there is a single smooth slope ($|\mathcal{S}|=1$) and the explicit computation of $v_{\mathtt{w}}(\rho)$ is doable, though rather involved, via Kasteleyn theory. In the general case we are considering here, computing $v_{\mathtt{w}}(\rho)$ directly using Kasteleyn theory seems very complicated, and we do not proceed that way. The first key point in the proof of the theorem is a simple relation (cf. (2.16)) between $v_{\mathtt{w}}(\cdot)$ and the limit shape $\psi_{\mathtt{w}}$ of the dimer model with edge weights $\mathtt{w}$ in the Aztec diamond. The limit shape is nothing but the solution of the Euler-Lagrange equation [20] associated to the dimer model’s surface tension, with weights $\mathtt{w}$ and boundary conditions determined by the geometry of the domain. This relation allows to translate analytic properties of $v_{\mathtt{w}}(\cdot)$ into analytic properties of the limit shapes, for which we use results from [1, 31]. In particular, singularities of $v_{\mathtt{w}}(\cdot)$ are in bijection with the facets (flat portions) of $\psi_{\mathtt{w}}$ that do not touch the boundary of the Aztec diamond or, equivalently, with the holes of the amoeba of the spectral curve [22]. In [8], the discontinuity of $Dv_{\mathtt{w}}(\rho)$ at the unique smooth phase was found via the explicit formula, but the connection with the facet of the limit shape was not realized. Another point we wish to emphasize is that, since edge weights change non-periodically with time as $\mathtt{w}=\{\mathtt{w}_{k}\}_{k\geq 0}$, it is a priori not obvious that an asymptotic speed of growth even exists (the connection with the limit shape shows that it does, because the limit shape $\psi_{\mathtt{w}_{k}}$ is actually independent of $k$).

Let us conclude this section by mentioning a recent article [36], that proves a hydrodynamic limit for the domino shuffling dynamics, in the form of the convergence of the rescaled height profile to the viscosity solution of the non-linear Hamilton-Jacobi PDE $\partial_{t}\phi=v_{\mathtt{w}}(\nabla\phi)$. The result of [36] is stated for the case of edge weights with space periodicity $1$, but the same proof presumably works for general periodic edge weights, as in the framework of the present article.

1.1. Related works on AKPZ growth models

Historically, the first rigorous result we are aware of, on a $(2+1)$-dimensional growth model in the AKPZ class, is [29], that computed the speed of growth of the Gates-Westcott model [15], verified that $\det(D^{2}v(\rho))<0$ and proved that stationary states are only logarithmically rough, in agreement with the above conjecture (growth of flucutations in time was not studied there). More recently, a growth model that is a $(2+1)$-dimensional, discrete, analog of Hammersley’s process has been introduced in [3]. Besides the computation of the speed of growth and the verification of $\det(D^{2}v(\rho))<0$, rigorous results on this model include the proof that height fluctuations grow at most logarithmically in space and time [3, 34], the study of stationary states [34], hydrodynamic limits for the height profile [3, 24], determinantal formulas for certain space-time correlations [3] and a CLT on scale $\sqrt{\log t}$ for height flucutations under special initial conditions [3]. Some of these results have been extended to an AKPZ growth process defined in terms of the dimer model on the square grid, see [7].

Apart from the above references, that deal with specific models, let us mention [4], that gives a sufficient condition for a $(2+1)$-dimensional growth model to belong to the AKPZ class. In simple terms, [4, Th. 2.1] states that if the hydrodynamic equation $\partial_{t}\phi=v(\nabla\phi)$ preserves solutions of the Euler-Lagrange equations associated to some strictly convex surface tension function $\sigma(\cdot)$, then $\det(D^{2}v(\rho))\leq 0$. This condition can be verified on several growth models, e.g. the one defined in [3], and it is related to the fact that these stochastic processes preserve a certain “local Gibbs property” ($\sigma$ is then the surface tension corresponding with such Gibbs potential).

The rest of the paper is organized as follows. In Section 2, we introduce the dimer model on $\mathbb{Z}^{2}$ for general weights, give some dimer model theory and give a precise version of our theorem. In Section 3, we prove the existence of the speed and its formula, while the main properties of the speed are proven in Section 4.

2. Model and results

In this section, we introduce the shuffling algorithm for the dimer model on $\mathbb{Z}^{2}$ for general weights, some of the basic dimer model theory and we precisely formulate our results.

2.1. Shuffling algorithm for the dimer model on $\mathbb{Z}^{2}$ (general weights)

The vertices of the graph $\mathbb{Z}^{2}$ are colored black and white in a bipartite way and they are assigned Cartesian coordinates, that is the neighbouring vertices which share a common edge with the vertex $(0,0)$ are $(1,0),(0,1)$, $(-1,0)$ and $(0,-1)$. We label a face $(i,j)\in\mathbb{Z}^{2}$ if its center has coordinates $(i+1/2,j+1/2)$; see Fig. 1 for an example on a $4\times 4$ torus graph.

The discrete time index of the Markov chain will be denoted $k=0,1,\dots$. We will say that a face $(i,j)$ is even if its bottom-left vertex is white, and odd otherwise. In the dynamics defined below, the colors of the vertices will interchange at each time step $k$ and we assume that initially the vertex $(0,0)$ is white. Therefore, a face with coordinates $(i,j)$ will be even at time $k$ if $i+j=k\mod 2$ and odd otherwise.

Given a weighting $\mathtt{w}$ of the edges, i.e. an assignment of a strictly positive weight to each edge, we first define a deterministic sequence $\{\mathtt{w}_{k}\}_{k\geq 0}$ of edge weightings with $\mathtt{w}_{0}:=\mathtt{w}$. To this purpose, note first that the weighting is uniquely defined if we specify the weights of edges on the boundary of every even face. We write then

[TABLE]

where the 4-tuple of positive numbers $(w_{i,j;k}^{a},w_{i,j;k}^{b},w_{i,j;k}^{c},w_{i,j;k}^{d})$ denotes the edge weights around the face $(i,j)$ at time $k$, where $a,b,c$ and $d$ are the edges labelled clockwise around the face, with $a$ being the topmost horizontal edge on that face. Also for $(i,j)\in\mathbb{Z}^{2}$ and $k\geq 0$, set

[TABLE]

The relation between $\mathtt{w}_{k}$ and $\mathtt{w}_{k+1}$ is, by definition,

[TABLE]

for $k\geq 0$ and $(i+j)=k+1\mod 2$; see Fig. 2.

We are now ready to define the shuffling algorithm. This is a discrete-time Markov chain on $\Omega$, the set of dimer coverings, or perfect matchings, of $\mathbb{Z}^{2}$. That is, each $\eta\in\Omega$ is a subset of edges of $\mathbb{Z}^{2}$, such that each vertex is contained in exactly one of them. Each edge contained in $\eta$ will be said to be “occupied by a dimer”. The chain is not time-homogeneous, since the transition rates depend on the time index $k$, via the edge weights $\mathtt{w}_{k}$. For $k\geq 0$, we define a random map $\Omega\ni\eta\mapsto{T_{k+1}}(\eta)\in\Omega$ through the following four steps, cf. Fig. 3 (only the third one is actually random):

1. (Deletion step)

All pairs of parallel dimers of $\eta$ covering two of the four boundary edges of any face that is even (at time $k$) are removed. 2. (Sliding step)

For every even face (at time $k$) with only one boundary edge covered by a dimer of $\eta$, slide this dimer across that face. 3. (Creation step)

For each face that is even at time $k$ (call $(i,j)$ its coordinates), if there are no dimers of $\eta$ covering any of its four boundary edges, add two parallel vertical dimers to the face with probability

[TABLE]

or two parallel horizontal dimers with probability

[TABLE]

(the operations are performed independently for each $(i,j)$ and $k$). 4. (Interchange step)

Interchange the white and black colors of vertices of the graph.

It is well known, and easy to check, that $T_{k}(\eta)\in\Omega$ if $\eta\in\Omega$. The swapping of colors at each step, that may seem to be pointless at this stage, will appear more natural in the discussion below of the evolution of the height function.

The maps $T_{k}$ are independent but not identically distributed, since the edge weights depend on $k$. Iteratively applying these maps and letting

[TABLE]

one obtains the desired Markov chain $\{\eta_{k}\}_{k\geq 0}$ on $\Omega$.

2.1.1. Height function and its evolution

Each dimer configuration $\eta\in\Omega$ is in one-to-one correspondence (up to a height offset) with a height function $h_{\eta}(\cdot)$ which is defined on the faces of $\mathbb{Z}^{2}$ [19]. That is, one fixes the height to be zero at some reference face $f_{0}$ and one defines the height gradients as

[TABLE]

where the sum runs over the edges crossed by a nearest-neighbor path $C_{f\to f^{\prime}}$ from $f$ to $f^{\prime}$, $\mathds{1}_{e\in\eta}$ is the indicator that $e$ is occupied by a dimer and $\sigma_{e}=+1$ or $-1$ according to whether $e$ is crossed with the white vertex on the right or left. The r.h.s. of (2.2) is well-known to be independent of the choice of $C_{f\to f^{\prime}}$.

In order for the shuffling algorithm to define a Markovian evolution of the height profile, we have to complement the definition of the maps $T_{k}$ with a prescription of how the height offset evolves as time $k$ increases. The convention that we adopt here is slightly different from that of [36, 8]. We start with the following observation, which is immediately verified from the definition of $T_{k}$ and of the height function (recall that vertex colors are swapped at each step). Let $f,f^{\prime}$ be any two faces that are odd at time $k$, i.e. they have coordinates $(i,j)$ and $(i^{\prime},j^{\prime})$ respectively, with $i+j=k+1\mod 2$ and $i^{\prime}+j^{\prime}=k+1\mod 2$; then,

[TABLE]

Therefore, we make the following choice:

Definition 2.1**.**

If $f$ is an odd face at time $k$, then

[TABLE]

This convention fixes unambiguously the whole height function of $\eta_{k+1}$ and in particular the value of $h_{T_{k+1}(\eta)}(f)$ for even faces $f$. Namely, let $f$ be any face and let $\eta|_{\partial f}$ (resp. $T_{k+1}(\eta)|_{\partial f}$) be the restriction of the dimer configuration $\eta$ (resp. $T_{k+1}(\eta)$) to the four boundary edges of $f$. Then, one may check by direct inspection starting from the definition of $T_{k+1}$ that, if $f$ is even at time $k$, then

[TABLE]

where (denoting $e_{1},\dots e_{4}$ the four boundary edges of $f$, labeled clockwise from the top one),

[TABLE]

and

[TABLE]

see Fig. 4.

2.2. Periodic weights

In this section, we introduce briefly some of the main aspects of the dimer model machinery needed for the formulation of the main result. Since we are interested in stochastic growth in a translationally invariant situation, here we specialize to the case where the edge weights are periodic in both directions of space. Let the fundamental domain $D_{0,0}$ of size $2n\times 2n,n\in\mathbb{N}$, consist of the vertices $\{(i,j):0\leq i,j\leq 2n-1\}$, half of which are black and half white. For $0\leq j\leq 2n-1$, the vertices $(2n,j)$ are identified with the vertices $(0,j)$ but are on the fundamental domain $D_{1,0}$ (obtained from $D_{0,0}$ via a horizontal translation by $2n$), while for $0\leq i\leq 2n-1$, the vertices $(i,2n)$ are identified with the vertices $(i,0)$ but on the fundamental domain $D_{0,1}$. The edge weights are chosen on all edges on $D_{0,0}$ and its boundary edges and then extented by periodicity to the whole graph. Call this weighting $\mathtt{w}_{0}$.

Underlying the dimer model theory is the characteristic polynomial $P$. To define $P$, consider $D_{0,0}$ embedded on a $2n\times 2n$ torus as above and let $\mathrm{wt}(x,y)$ denote the weight of the edge $(x,y)$ for two vertices $x$ and $y$ of $D_{0,0}$. Given $z,w\in\mathbb{C}$, define $K(z,w)$ to be the Kasteleyn matrix with rows indexed by white vertices and columns indexed by black vertices of $D_{0,0}$, with

[TABLE]

where $x$ is a white vertex and $y$ is black vertex in $D_{0,0}$, and

[TABLE]

and

[TABLE]

for $0\leq k,l\leq 2n-1$. The Laurent polynomial $P(z,w)=\det K(z,w)$ is called “characteristic polynomial” [21]. Of course, $P$ depends on $n$ and on the weights.

From [21], the Newton Polygon (depending on $n$) is defined to be

[TABLE]

One can check, for the $K(z,w)$ specified above, that $N(P)$ is the (closed) square with vertices $(\pm n,0),(0,\pm n)$.

A probability measure $\mu$ on $\Omega$ is said to be an ergodic Gibbs measure (corresponding to the edge weights $\mathtt{w}_{0}$) if:

• •

it is invariant and ergodic with respect to horizontal/vertical translations by multiples of $2n$;

• •

it satisfies the following Dobrushin-Lanford-Ruelle (DLR) property. Given any finite subset $\Lambda$ of edges and any dimer configuration $\bar{\eta}\in\Omega$, let $\Omega_{\Lambda,\bar{\eta}}$ be the (finite) set of dimer configurations $\eta\in\Omega$ that coincide with $\bar{\eta}$ outside $\Lambda$. Then, conditionally on $\eta=\bar{\eta}$ outside $\Lambda$, the $\mu$-probability of a configuration $\eta$ is proportional to the product

[TABLE]

of $\mathtt{w}_{0}$-weights of the edges in $\Lambda$ occupied by dimers.

Thanks to translation invariance, one may associate to each ergodic Gibbs measure $\mu$ an average slope $\rho=(\rho_{1},\rho_{2})$. Here, $\rho_{1}$ (resp. $\rho_{2}$) is the expected height difference between a face in $D_{0,0}$ and its translate in $D_{1,0}$ (resp. $D_{0,1}$). The slope $\rho$ is contained in the Newton polygon $N(P)$. Moreover, provided that $\rho$ belongs $\stackrel{{\scriptstyle\circ}}{{N(P)}}$, the interior of $N(P)$, there exists a unique Gibbs measure with slope $\rho$ [21] and we denote it by $\pi_{\rho,\mathtt{w}_{0}}$. This measure is known to be determinantal, in the sense that the probability that $r$ given edges $e_{1},\dots,e_{r}$ belong to $\eta$ is given by the determinant of an $r\times r$ matrix, whose entries are elements of the so-called inverse Kasteleyn matrix.

Define the Ronkin function associated to $P$ as

[TABLE]

for $B=(B_{1},B_{2})\in\mathbb{R}^{2}$. From [21], $R$ is the Legendre transform of the so-called surface tension $\sigma$ of the dimer model with the given periodic weights, i.e.

[TABLE]

We will recall later the relation between $\sigma$ and the “limit shapes” of the dimer model.

We write $\stackrel{{\scriptstyle\circ}}{{N(P)}}$, the interior of the Newton polygon, as the disjoint union of $\mathcal{R}$ (rough region) and $\mathcal{S}$ (smooth region), whose definition we recall now. (Rough (resp. smooth) phases are called “liquid” (resp. “gaseous”) phases in [21].) From [21], it is known that if $\rho\in\stackrel{{\scriptstyle\circ}}{{N(P)}}$, two cases can occur:

• •

either the measure $\pi_{\rho,\mathtt{w}_{0}}$ is rough, meaning that height fluctuations of $h_{\eta}(f)-h_{\eta}(f^{\prime})$ grow logarithmically w.r.t. the distance between the faces $f,f^{\prime}$. More precisely,

[TABLE]

as the distance $|f-f^{\prime}|$ between $f$ and $f^{\prime}$ diverges. Moreover, the scaling limit of the height profile is a Gaussian Free Field [19]. We call $\mathcal{R}$ the set of such “rough slopes”;

• •

or the measure $\pi_{\rho,\mathtt{w}_{0}}$ is smooth, meaning that height flucutations of $h_{\eta}(f)-h_{\eta}(f^{\prime})$ have uniformly bounded variance. In this case, $\sigma$ has a cone singularity (i.e. $\nabla\sigma$ is discontinuous) at this value of $\rho$. The set of “smooth slopes” is denoted $\mathcal{S}$.

In both cases, $\sigma(\cdot)$ is strictly convex.

From [21], it is further known that $\mathcal{S}$ is a finite set and moreover

[TABLE]

for generic edge weights $\mathtt{w}_{0}$, $\mathcal{S}$ actually coincides with the whole $\stackrel{{\scriptstyle\circ}}{{N(P)}}\cap\;\mathbb{Z}^{2}$, which contains $2n(n-1)+1$ points. However, this may fail for particular choices of weights: for instance, when all edge-weights are equal, then it is known that $\mathcal{S}=\emptyset$.

2.2.1. Shuffling algorithm with periodic weights

A remarkable feature of the shuffling algorithm is the following (see for instance [8, Proposition 3.1] and [36, Prop. 2.2]):

Proposition 2.2**.**

If the initial condition $\eta_{0}$ at time [math] is drawn from $\pi_{\rho,\mathtt{w}_{0}}$ (i.e., $\eta_{0}\sim\pi_{\rho,\mathtt{w}_{0}}$), then at time $k$ one has $\eta_{k}\sim\pi_{\rho,\mathtt{w}_{k}}$.

If we had not swapped vertex colors at each step, the slope $\rho$ would swap to $-\rho$ at each step. There are two observations that we will need going forward. The first one is that the characteristic polynomial only changes by a multiplicative constant factor when the weights $\mathtt{w}_{k}$ are replaced by $\mathtt{w}_{k+1}$ [16] In particular, in view of (2.7) and (2.8), this implies that the surface tension $\sigma(\cdot)$ for weights $\mathtt{w}_{k}$ equals that for weights $\mathtt{w}_{k+1}$, up to an additive constant. Another consequence is the following: since the rough or smooth nature of $\pi_{\rho,\mathtt{w}}$ depends on $\mathtt{w}$ only through the zeros of the characteristic polynomial $P(z,w)$ on the torus $\{z,w\in\mathbb{C}:|z|=|w|=1\}$ [21], we deduce that $\pi_{\rho,\mathtt{w}_{k+1}}$ is rough (resp. smooth) iff $\pi_{\rho,\mathtt{w}_{k}}$ is. In other words, the condition $\rho\in\mathcal{R}$ does not depend on $k$.

Another important observation is the following: even though weights $\mathtt{w}_{0}$ (and therefore $\mathtt{w}_{k}$) are periodic in space, the sequence $\{\mathtt{w}_{k}\}_{k\geq 0}$ is in general not periodic w.r.t. the time index $k$. Time-periodicity can, however, hold for special choices of $\mathtt{w}_{0}$ and indeed the cases studied in [8, 36] are time-periodic.

2.3. The Aztec diamond

The Aztec diamond $A_{N}$ of size $N$ is the subset of the graph $\mathbb{Z}^{2}$ whose vertices have Cartesian coordinates $(x_{1},x_{2})$ satisfying the condition $|x_{1}-1/2|+|x_{2}-1/2|\leq N$. We let $E^{+}_{N}$ denote the set of edges outgoing from $A_{N}$, $F^{+}_{N}$ the set of faces not in $A_{N}$ but neighboring $A_{N}$ and $F_{N}$ the set of internal faces of $A_{N}$.

Let $\pi_{\mathtt{w},N}$ be the probability measure on $\Omega_{N}$, the set of perfect matchings of $A_{N}$, where the weight of a configuration is proportional to the product of the $\mathtt{w}$-weights of edges occupied by dimers. Since all vertices of $A_{N}$ are matched among themselves, all edges in $E^{+}_{N}$ are empty and therefore the height difference between two faces in $F^{+}_{N}$ is independent of the choice of $\eta\in\Omega_{N}$. We assume that the coloring of the vertices is such that the vertex of coordinates $(-N+1,1)$ is white. We fix the height offset as in Fig. 5, by setting the height to $+N/4$ on the leftmost face of $F^{+}_{N}$; then, the boundary height ranges from $-N/4$ to $+N/4$.

The height function in $A_{N}$ satisfies a limit shape phenomenon (or law of large numbers) as $N\to\infty$. Namely, rescale the lattice mesh by $1/(2nN)$ and call $\hat{A}_{N}$ the rescaled Aztec diamond (and correspondingly denote $\hat{E}_{N}^{+},\hat{F}_{N}^{+},\hat{F}_{N}$ the analog of $E_{N}^{+},F_{N}^{+},F_{N}$). The union of the faces of $\hat{A}_{N}$ tends to the square

[TABLE]

Define the rescaled height function $\hat{h}_{\eta}:\hat{F}_{N}\mapsto\mathbb{R}$ as

[TABLE]

with $f\in F_{N}$ the face of $A_{N}$ that corresponds to $\hat{f}$ before rescaling. Thanks to the factor $2n$ in the rescaling, $\hat{h}_{\eta}$ is a Lipschitz function whose gradient is contained in the Newton polygon $N(P)$. Note that, if $\hat{f}=\hat{f}_{N}$ is a face in $\hat{F}^{+}_{N}$ whose center tends to $x=(x_{1},x_{2})\in\partial Q$ as $N\to\infty$, then

[TABLE]

The limit shape theorem (cf. [9] for the model with uniform weights and [23] for the general periodic case) states that there exists a Lipschitz function $\psi_{\mathtt{w}}:Q\mapsto\mathbb{R}$ that coincides with $\psi_{\partial Q}$ on $\partial Q$, such that for every $\delta>0$,

[TABLE]

Here, with some abuse of notation, $\psi_{\mathtt{w}}(\hat{f})$ means $\psi_{\mathtt{w}}$ computed at the center of the face $\hat{f}$. The “limit shape” $\psi_{\mathtt{w}}$ is characterized by being the unique minimizer of the surface tension functional

[TABLE]

among Lipschitz functions that equal $\psi_{\partial Q}$ on the boundary. While the boundary condition does not depend on $\mathtt{w}$, the limit shape does (through the surface tension), but $\psi_{\mathtt{w}_{k+1}}=\psi_{\mathtt{w}_{k}}$ because, as we already mentioned, $\sigma$ changes only by an additive constant when $\mathtt{w}_{k}$ is changed into $\mathtt{w}_{k+1}$.

2.4. Statement of Main theorem

Our main result concerns the average speed of growth for the Markov process in the infinite graph, started from $\pi_{\rho,\mathtt{w}}$. By definition, this is given by the limit (provided it exists)

[TABLE]

with $f$ any face of $\mathbb{Z}^{2}$ and $\mathbb{E}_{\nu}$ the law of the process started from the probability measure $\nu$. Note that every second term in the sum is zero because every second time the face $f$ is odd.

Since $\eta_{j}\sim\pi_{\rho,\mathtt{w}_{j}}$ with $\mathtt{w}_{0}\equiv\mathtt{w}$, each non-zero term in the sum could in principle be computed via (2.4) and Kasteleyn theory, using the determinantal structure of the measure $\pi_{\rho,\mathtt{w}_{k}}$. Following this route, however, it is not clear how to get any manageable expression or to prove that the limit $k\to\infty$ in (2.14) even exists. One reason is that, for generic periodic weights, it is hard to invert the infinite-volume Kasteleyn matrix explicitly. Fortunately, an alternative way exists, that leads to:

Theorem 2.3**.**

For every $\rho\in\stackrel{{\scriptstyle\circ}}{{N(P)}}$ and positive periodic weighting $\mathtt{w}$, there exists $v=v_{\mathtt{w}}(\rho)$ such that, for any face $f\in\mathbb{Z}^{2}$,

[TABLE]

The speed $v_{\mathtt{w}}(\cdot)$ is determined as follows: let $\psi_{\mathtt{w}}(\cdot)$ be the limit shape for the dimer model in the Aztec diamond with weights $\mathtt{w}$ and let $x_{\mathtt{w}}(\rho)\in\stackrel{{\scriptstyle\circ}}{{Q}}$ (the interior of the unit square in (2.11)) be a point such that $\nabla\psi_{\mathtt{w}}(x_{\mathtt{w}}(\rho))=\rho$. Then

[TABLE]

On the rough region $\mathcal{R}$, $v_{\mathtt{w}}(\cdot)$ is $C^{\infty}$and $\det(D^{2}v_{\mathtt{w}})<0$. On the other hand, $Dv_{\mathtt{w}}$ is discontinuous at every $\rho\in\mathcal{S}$.

A few comments are in order:

• •

the existence of $x_{\mathtt{w}}(\rho)$ is part of the statement. Uniqueness in general fails (the limit shape $\psi_{\mathtt{w}}$ may have “facets”, i.e. open regions where it is affine) but for $\rho\in\mathcal{R}$, the point $x_{\mathtt{w}}(\rho)$ is unique (see Section 4).

• •

Using smoothness of $v_{\mathtt{w}}(\cdot)$ on $\mathcal{R}$ and (2.16), one sees that

[TABLE]

Note that the r.h.s. of (2.16) looks like (minus) the Legendre transform of $\psi_{\mathtt{w}}$, except that there is no infimum over $x$ and in fact neither $v_{\mathtt{w}}$ nor $\psi_{\mathtt{w}}$ have any definite convexity.

• •

It was observed in [20] that the Euler-Lagrange equation satisfied by the limit shape $\psi_{\mathtt{w}}$ of dimer models can be written (in the “rough region” where the limit shape is $C^{2}$) in terms of a first-order PDE (“complex Burgers equation”) for a complex pair $(z,w)$ related by the relations $P(z,w)=0$ and $\pi\nabla\psi_{\mathtt{w}}=(-{\rm arg}(w),{\rm arg}(z))$. Locally, these relations give a bijection between $z$ and $\rho=\nabla\psi_{\mathtt{w}}$. Then, using [4, Sec. 3], the above Theorem 2.3 can be complemented by the following statement:

[TABLE]

• •

For a special case of two-periodic weights ($n=1$), it was found via explicit computation in [8, Th. 3.11] that the behavior of $v_{\mathtt{w}}$ near the unique gas slope $\rho=0$ is of the type

[TABLE]

The absence of the $|\rho|^{2}$ term can be given an interesting interpretation. In fact, this is a simple consequence of formula (2.16) plus the fact that, if $x$ approaches a point $x_{0}$ on the boundary of the “facet” where $\nabla\psi_{\mathtt{w}}\equiv 0$, then generically $\psi_{\mathtt{w}}(x)-\psi_{\mathtt{w}}(x_{0})$ vanishes as $|x-x_{0}|^{3/2}$ [20] (this behavior is referred to as “Pokrovsky-Talapov law” [28]).

2.4.1. Fluctuations

One can further prove that height fluctuations grow slowly (at most logarithmically) in time, as is typical for growth models in the AKPZ universality class. In fact, one has uniformly in $k\geq 1$

[TABLE]

for some constant $c$, where $g(k)=\sqrt{\log(k+1)}$ if $\rho\in\mathcal{R}$ and $g(k)\equiv 1$ if $\rho\in\mathcal{S}$. The proof of this fact works the same as in [8] so we will not add details (the speed of convergence $O(u^{-2})$ was not explicitly stated in [8], but it can be immediately extracted from the proof). Note in particular that

[TABLE]

and since the r.h.s. is summable in $k$, one can upgrade (2.15) to the almost-sure convergence, with respect to the joint law of the initial condition and of the process,

[TABLE]

3. Identification of the speed of growth

In this section, we prove existence of the speed and formula (2.16).

3.1. General properties of the dynamics

We need two general facts: the dynamics is monotone (it preserves stochastic ordering among height profiles) and it is local (information travels at most ballistically through the system).

Let us start with monotonicity. Given two dimer configurations $\eta,\eta^{\prime}$, we say that $h_{\eta}\preceq h_{\eta^{\prime}}$ if $h_{\eta}(f)\leq h_{\eta^{\prime}}(f)$ for every face $f$. Given two initial configurations $\eta_{0},\eta^{\prime}_{0}$, we can couple the two Markov chains $\{\eta_{k}\}_{k\geq 0},\{\eta^{\prime}_{k}\}_{k\geq 0}$ in the following way (global monotone coupling): for any face $f$, if in both configurations $\eta_{k-1}|_{\partial f}=\eta^{\prime}_{k-1}|_{\partial f}=\emptyset$, then in the “creation step” of the shuffling map $T_{k}$ we choose the same randomness to decide whether we add two vertical or two horizontal dimers around $f$. Then, the following statement holds (it implies the preservation of stochastic order mentioned above): if $h_{\eta_{0}}\preceq h_{\eta^{\prime}_{0}}$, then the same holds at all later times $k$ [36, Lemma 2.4].

As far as locality is concerned the point is that, by the definition of the shuffling algorithm, the value of $h_{\eta_{k}}(f)$ is completely determined by the height at time $k-1$ at the face $f$ and at its four neighbors (this determines the dimer configuration $\eta_{k-1}$ on $\partial f$), plus the randomness used to create parallel dimers at $f$, if the face is even and $\eta_{k-1}|_{\partial f}=\emptyset$. From this, it is immediate to deduce the following:

Proposition 3.1**.**

Let $\eta_{0},\eta^{\prime}_{0}$ be two dimer configurations whose height coincides on all faces at $\ell_{1}$-distance up to $N+1$ from a given face $f$. Couple the Markov chains started from $\eta_{0},\eta^{\prime}_{0}$ via the global monotone coupling. Then, $h_{\eta_{k}}(f)=h_{\eta^{\prime}_{k}}(f)$ for every $k\leq N$.

Let us also describe in some more detail how the shuffling algorithm works on the Aztec diamond (this is the framework where the algorithm was originally introduced [13, 30]). In a step of the algorithm, a dimer configuration $\eta$ on $A_{N}$ is mapped to a configuration $\eta^{\prime}$ on the larger domain $A_{N+1}$. Suppose that we have $\eta_{N}\in\Omega_{N}$, i.e. a dimer configuration on the diamond of size $N$. We can also view $\eta_{N}$ as a subset of edges of $A_{N+1}$ (but not a perfect matching, since the boundary vertices are necessarily unmatched). To construct $\eta_{N+1}$, apply the map $T_{N+1}$ in $A_{N}$ to $\eta_{N}$ (with weights $\mathtt{w}_{N}$ as above). Note that the faces in $A_{N+1}$ that are closest to the boundary, i.e. the faces in $F^{+}_{N}$, are even. It is well known that the resulting dimer configuration $\eta_{N+1}$ is a perfect matching of $A_{N+1}$. Due to the swapping of colors, at the next step the faces in $F^{+}_{N+1}$ are again even and the procedure goes on.

The analog of Proposition 2.2 in the Aztec diamond is the well known fact that, if we start at time zero with a configuration $\eta_{0}$ on $A_{N}$ such that $\eta_{0}\sim\pi_{\mathtt{w}_{0},N}$ for certain periodic weights $\mathtt{w}_{0}$, then at time $k$ one has $\eta_{k}\sim\pi_{\mathtt{w}_{k},N+k}$.

There is an important point to be discussed: when we introduced the shuffling algorithm on the infinite lattice, we fixed the evolution of the height offset via Definition 2.1. On the other hand, on the Aztec diamond the height offset is fixed by the requirement that the left-most face in $F^{+}_{k}$ has height $k/4$. These two conventions must be compatible, i.e., if we adopt the convention (2.4) for the evolution of the height function, then the height on the left-most face of $F^{+}_{k}$ must be $k/4$ deterministically. This is easily seen inductively in $k$, as explained in the caption of Fig. 6.

3.2. The speed of growth

Here we prove the following:

Proposition 3.2**.**

Let $\rho\in\stackrel{{\scriptstyle\circ}}{{N(P)}}$ and assume that there exists $x_{\mathtt{w}}(\rho)$ in the interior of $Q$, such that $\psi_{\mathtt{w}}(\cdot)$ is $C^{1}$ in a neighborhood of $x_{\mathtt{w}}(\rho)$ and $\nabla\psi_{\mathtt{w}}(x_{\mathtt{w}}(\rho))=\rho$. Then, the limit in (2.15) exists and (2.16) holds.

The existence of $x_{\mathtt{w}}(\rho)$ for every $\rho$ in the interior of the Newton polygon will be proved in the next section.

Proof.

For an integer $N$, let $\bar{f}_{N}$ be a face of $A_{N}$ whose center is at minimal distance from $(2nN)x_{\mathtt{w}}(\rho)$. One should think of $N$ as being a large multiple of $k$, the time in (2.15), with $\epsilon:=k/N$ that will be sent to zero at the end. For later convenience, we let $\Lambda_{N,\epsilon}$ be the square box of side $2k+1$ centered at the face $\bar{f}_{N}$. Recall that $A_{N}$ denotes the $N\times N$ Aztec diamond and take the edge weights to be given by $\mathtt{w}$. We run the shuffling dynamics in the Aztec diamond, starting at time zero with the domain $A_{N}$ and with an initial condition $\tilde{\eta}_{0}$ sampled from $\pi_{\mathtt{w},N}$. We denote $\tilde{\eta}_{k}$ the configuration at time $k$, where the tilda is used just to distinguish this from the evolution in the infinite graph. The height function of $\tilde{\eta}_{0}$ is concentrated at the limit shape $\psi_{\mathtt{w}}(\cdot)$. In particular, from (2.13) with $\delta=\epsilon^{2}$ we have

[TABLE]

As before, we identify with some abuse of notation a face $f$ with the point at its center. As observed in Section 3.1, at time $k$, the configuration $\tilde{\eta}_{k}$ has law $\pi_{\mathtt{w}_{k},N+k}$ and we still have (3.1) with $N$ replaced by $N+k$. Altogether, we see that

[TABLE]

where we used that $\nabla\psi_{\mathtt{w}}(x_{\mathtt{w}}(\rho))=\rho$ and the error term $o_{\epsilon}(1)$ vanishes as $\epsilon\to 0$, since the limit shape is $C^{1}$ around $x_{\mathtt{w}}(\rho)$. We also used the fact that $|h_{\tilde{\eta}_{0}}|/N,|h_{\tilde{\eta}_{k}}|/N$ are uniformly bounded for $k\leq\epsilon N$, to deduce from (3.1) a statement about their average.

Our goal now is to prove a statement analogous to (LABEL:eq:perA) for the dynamics $\{\eta_{k}\}_{k\geq 0}$ on the infinite graph. By Proposition 3.1, the evolution of $h_{\eta_{j}}(\bar{f}_{N}),j\leq k$ is not influenced by the height function of $\eta_{0}$ outside $\Lambda_{N,\epsilon}$. Recall that $\eta_{0}$ is sampled from the infinite-volume measure $\pi_{\rho,\mathtt{w}}$. Under this probability measure, the height function is essentially linear, with slope $\rho$ and sub-linear fluctuations. More precisely,

[TABLE]

where once more we have identified a face with its center and the factor $1/(2n)$ is there because $\rho$ is the average height change per fundamental domain. To get (3.3), observe first that

[TABLE]

uniformly in $f\in\Lambda_{\epsilon,N}$ (the error term is there because $f$ is not necessarily an exact translation of $\bar{f}_{N}$ in a different fundamental domain). Also, recall that the fourth centered moment of $h_{\eta_{0}}(\bar{f}_{N})-h_{\eta_{0}}(f)$ under $\pi_{\rho,\mathtt{w}}$ grows at most like $(\log|\bar{f}_{N}-f|)^{2}$ for $|\bar{f}_{N}-f|$ large (see [21, Section 4] for a $O(\log|\bar{f}_{N}-f|)$ bound on the variance; higher moments are treated analogously). Then, a union bound over $f\in\Lambda_{N,\epsilon}$ and an application of Chebyshev’s inequality leads to (3.3).

Note that we have not yet specified the height offset $h_{\eta_{0}}(\bar{f}_{N})$ at time zero. We will fix it in such a way that, with high probability (w.h.p.) as $N\to\infty$,

[TABLE]

For this, note that (3.3) implies that w.h.p.

[TABLE]

while (3.1) and $C^{1}$ continuity of the limit shape implies that w.h.p.

[TABLE]

with $R_{N,\epsilon}/N\epsilon=o_{\epsilon}(1)$. Then, (3.4) holds provided we choose

[TABLE]

By monotonicity of the dynamics and Proposition 3.1 we see that $h_{\eta_{k}}(\bar{f}_{N})\leq h_{\tilde{\eta}_{k}}(\bar{f}_{N})$ and therefore, w.h.p.,

[TABLE]

where we used (3.1) in the last step and $k=\epsilon N$. Note that $[h_{\eta_{k}}(f)-h_{\eta_{0}}(f)]/k$ is deterministically bounded by $1$, so we can turn the statement w.h.p. into a statement in average and obtain that

[TABLE]

where we used also (LABEL:eq:perA). Note that the face $\bar{f}_{N}$ depends on the time $k=N\epsilon$. However, since the measure $\pi_{\rho,\mathtt{w}}$ is invariant by translations of multiples of $2n$ and the height function has bounded Lipschitz constant, we have (3.9) also for any fixed face $f$. Finally, we let $\epsilon\to 0$.

A lower bound is proven in the very same way and altogether the statements (2.15) and (2.16) follow. ∎

With similar arguments, we also obtain the following result, that will be useful later:

Proposition 3.3**.**

If there exists $x$ in the interior of $Q$ such that $\psi_{\mathtt{w}}$ is $C^{1}$ in a neighborhood of $x$ and $\nabla\psi_{\mathtt{w}}(x)=\rho$ with $\rho=(\rho_{1},\rho_{2})$ at one of the four corners of the Newton polygon (i.e., $\rho=(\pm n,0)$ or $\rho=(0,\pm n)$) then

[TABLE]

Proof.

Assume to fix ideas that $\rho=(n,0)$. As above, let $\bar{f}_{N}$ be the face of $A_{N}$ closest to $(2nN)x$, let $k=N\epsilon$ and $\Lambda_{N,\epsilon}$ be the square of side $2k+1$ centered around $\bar{f}_{N}$. One has, in analogy with (LABEL:eq:perA) and with the same argument,

[TABLE]

On the other hand, let $F$ be the collection of faces in $\Lambda_{N,\epsilon}$ (there are approximately $4k^{2}$ of them). Write $F=F^{(+,j)}\cup F^{(-,j)}$, where $F^{(+,j)}$ contains the faces that are even at time $j$ and $F^{(-,j)}$ all the others. Because of (2.13) and the fact that $\nabla\psi_{\mathtt{w}}=(n,0)+o_{\epsilon}(1)$ in an $\epsilon$-neighborhood of $x$, from the definition of height function we see that, with probability $1-o_{\epsilon}(1)$, a proportion $1-o_{\epsilon}(1)$ of the dimers of $\tilde{\eta}_{0}$ in $\Lambda_{N,\epsilon}$ occupy a vertical edge with bottom white vertex. The same holds for $\tilde{\eta}_{j},j\leq k$, because $\tilde{\eta}_{j}$ has the same limit shape as $\tilde{\eta}_{0}$. Therefore, a proportion $1-o_{\epsilon}(1)$ of the faces in $F^{(+,j)}$ have a single vertical dimer of $\tilde{\eta}_{j}$ along their boundary. From (2.4) we see that each such even face contributes $-1/2$ to the height change from time $j$ to $j+1$. Since $|F^{(+,j)}|/|F|=1/2+o_{\epsilon}(1)$, the l.h.s. of (3.11) equals also $-1/4+o_{\epsilon}(1)$ and (3.10) follows. ∎

3.3. The limit shape

Here we give some analytic properties of the limit shape $\psi_{\mathtt{w}}(\cdot)$ and prove the existence of $x_{\mathtt{w}}(\rho)$:

Theorem 3.4**.**

There exists a non-empty, open subset $\mathcal{F}$ of the rescaled Aztec diamond $Q$ (cf. (2.11)) where $\psi_{\mathtt{w}}(\cdot)$ is $C^{1}$ and the gradient $\nabla\psi_{\mathtt{w}}(\cdot)\in\stackrel{{\scriptstyle\circ}}{{N(P)}}$. For every $\rho\in\stackrel{{\scriptstyle\circ}}{{N(P)}}$, there exists $x_{\mathtt{w}}(\rho)\in\mathcal{F}$ such that

[TABLE]

Proof of Theorem 3.4.

For $x=(x_{1},x_{2})\in Q$, let

[TABLE]

and note that $\psi^{-}$ (resp. $\psi^{+}$) is the mimimal (resp. maximal) Lipschitz function with gradient in $N(P)$ that equals $\psi_{\partial Q}$ on $\partial Q.$ We let

[TABLE]

It is easy to see the following (the proof is given below):

Lemma 3.5**.**

The set $\mathcal{F}_{0}$ is non-empty.

We need some regularity properties of the limit shape $\psi_{\mathtt{w}}$, and for this we appeal to [31, 1]. Let us compactify the Newton polygon by introducing a continuous map $H:N(P)\mapsto S^{2}$ (the two-dimensional sphere) in such a way that $\partial N(P)$ is mapped to a point of $S^{2}$ while $H$ is a homeomorphism between $\stackrel{{\scriptstyle\circ}}{{N(P)}}$ and $H(\stackrel{{\scriptstyle\circ}}{{N(P)}})$. Then one has:

Proposition 3.6**.**

[31*, Th. 4.1 and Th. 1.3]**

The map $x\mapsto H(\nabla\psi_{\mathtt{w}}(x))$ is continuous in the interior of $Q$. Moreover, $\psi_{\mathtt{w}}$ is $C^{1}$ in $\mathcal{F}_{0}$.*

Define further the open set

[TABLE]

that is the one appearing in the statement of Theorem 3.4. Decompose $\mathcal{F}$ as the union of the open set

[TABLE]

and the closed set

[TABLE]

Proposition 3.7**.**

The set $\mathcal{F}_{R}$ is non-empty.

Let us assume for the moment Proposition 3.7 (the proof is given below) and let us proceed with the proof of Theorem 3.4. In general, $\mathcal{F}_{S}$ consists of a collection of disjoint, simply connected sets (these were called “bubbles” in [20]); on each bubble, the gradient $\nabla\psi_{\mathtt{w}}$ is constant and belongs to one of the finitely many slopes in $\mathcal{S}$. It is also known [26] that, on $\mathcal{F}_{R}$, the limit shape $\psi_{\mathtt{w}}$ is not just $C^{1}$ but actually $C^{\infty}$, since the surface tension $\sigma(\rho)$ is $C^{\infty}$ for $\rho\in\mathcal{R}$. Therefore, in particular, the map $D:x\mapsto\nabla\psi_{\mathtt{w}}(x)$ is a $C^{1}$ map from $\mathcal{F}_{R}$ to $\mathcal{R}$. The next step requires the following:

Theorem 3.8**.**

[1*]**

The map $D:x\mapsto\nabla\psi_{\mathtt{w}}(x)$ is a proper map111 Properness does not hold for general domains $Q$ and boundary values $\psi_{\partial Q}$. Theorem 3.8 holds for the Aztec diamond because in this case $Q$ is a convex polygonal domain with sides perpendicular to the sides of the Newton polygon, and $\psi_{\partial Q}$ in (2.12) is a “natural boundary value” for $Q$. The notion of “natural boundary value” is defined in [1] and it requires in particular that, if the side $\ell$ of $Q$ is perpendicular to the side $[p_{i},p_{i+1}]$ of the Newton polygon with $p_{i},p_{i+1}$ two of its adjacent corners, then the derivative of $\psi_{\partial Q}$ along $\ell$ equals $\langle t_{\ell},p_{i}\rangle$ with $t_{\ell}$ the tangent vector to $\partial Q$ along $\ell$. from $\mathcal{F}_{R}$ to $\mathcal{R}$ (i.e. the pre-image of every compact subset of $\mathcal{R}$ is compact). *

Let us prove that the Jacobian $\det(J(x))$ of the map $D$ is everywhere non-positive on $\mathcal{F}_{R}$ and not identically zero. The Jacobian matrix equals

[TABLE]

On the other hand, on $\mathcal{F}_{R}$, $\psi_{\mathtt{w}}(\cdot)$ satisfies the Euler-Lagrange equation

[TABLE]

with $\sigma_{ab}$ the derivative of $\sigma(\rho)$ w.r.t. the arguments $\rho_{a},\rho_{b}$, computed at $\rho:=\nabla\psi_{\mathtt{w}}(x)\in\mathcal{R}$. For $\rho\in\mathcal{R}$, the matrix $\{\sigma_{ab}\}_{a,b=1,2}$ is strictly positive definite, in particular $|\sigma_{12}|<\sqrt{\sigma_{11}\sigma_{22}}$. From this, we deduce that

[TABLE]

In fact, assume first that $\partial^{2}_{x_{1}x_{2}}\psi_{\mathtt{w}}(x)=0$. Then, $\partial^{2}_{x_{1}}\psi_{\mathtt{w}}(x)\partial^{2}_{x_{2}}\psi_{\mathtt{w}}(x)\geq 0$ (because $\det J(x)\geq 0$) but on the other hand (3.18) reduces to

[TABLE]

Since both $\sigma_{11},\sigma_{22}$ are strictly positive, the only possibility is that $\partial^{2}_{x_{1}}\psi_{\mathtt{w}}(x)=\partial^{2}_{x_{2}}\psi_{\mathtt{w}}(x)=0$. On the other hand, assume (by contradiction) that $\partial^{2}_{x_{1}x_{2}}\psi_{\mathtt{w}}(x)\neq 0$, so that $\partial^{2}_{x_{1}}\psi_{\mathtt{w}}(x)\partial^{2}_{x_{2}}\psi_{\mathtt{w}}(x)>0$. Then

[TABLE]

which is a contradiction because the second inequality is strict. Altogether, (3.19) follows. From this, we see that $\det(J(\cdot))$ can vanish identically on $\mathcal{F}_{R}$ only if $\psi_{\mathtt{w}}(\cdot)$ is affine, which is clearly not possible in view of Proposition 3.8.

We have that the map $D$ is proper and its Jacobian is non-negative and not identically vanishing. Then, by [27, Th. 1], we deduce that the map $D$ is onto: for every $\rho\in\mathcal{R}$, there exists $x_{\mathtt{w}}(\rho)\in\mathcal{F}_{R}$ with $\nabla\psi_{\mathtt{w}}(x_{\mathtt{w}}(\rho))=\rho$.

It remains to show the existence of $x_{\mathtt{w}}(\rho)$ for every $\rho\in\mathcal{S}$. Let $\{\rho_{i}\}$ be a sequence of slopes in $\mathcal{R}$ that converges to $\rho$. Any limit point $\bar{x}$ of $x_{\mathtt{w}}(\rho_{i})$ is in $\mathcal{F}_{0}$ (because of Proposition 3.8). Due to Proposition 3.6, the slope of $\psi_{\mathtt{w}}$ at $\bar{x}$ is $\rho$, so we can set $x_{\mathtt{w}}(\rho):=\bar{x}$. ∎

We conclude this section by proving the two technical results, Lemma 3.5 and Proposition 3.7 that were stated above.

Proof of Lemma 3.5.

Since $\psi^{-}(x)<\psi^{+}(x)$ for every $x$ in the interior of $Q$ and $\psi_{\mathtt{w}}$ is continuous, we have just to exclude that $\psi_{\mathtt{w}}\equiv\psi^{-}$ or $\psi_{\mathtt{w}}\equiv\psi^{+}$. Assume for instance that $\psi_{\mathtt{w}}\equiv\psi^{+}$; we are going to exhibit a function $\psi$, with the right boundary value, such that

[TABLE]

For this purpose, let for $\epsilon>0$ small

[TABLE]

It is immediate to see that $\psi(x)=4\epsilon n^{2}x_{1}^{2}+(1/4-\epsilon)$ in

[TABLE]

and $\psi(x)=\psi^{+}(x)$ in $Q\setminus S_{\epsilon}$, so in particular $\psi$ equals $\psi_{\partial Q}$ on $\partial Q$. The difference between the r.h.s. and the l.h.s. of (3.23) is then

[TABLE]

Since $\sigma(\cdot)$ is strictly convex, one has

[TABLE]

Therefore, using continuity of $\sigma(\cdot)$, for $\epsilon$ small enough the difference (3.24) is strictly positive and, as a consequence, the minimizer $\psi_{\mathtt{w}}$ of the surface tension functional cannot coincide with $\psi^{+}$. ∎

Proof of Propoposition 3.7.

We begin by making an observation on the shape of $\mathcal{F}_{0}$. Recall from (3.12) the definition of $\phi_{a},a\in\{N,E,S,W\}$ and define the (possibly empty) regions

[TABLE]

$Q_{N}$ belongs to the triangle $\{x\in Q:x_{2}\geq 0\}$, otherwise $\psi_{\mathtt{w}}$ would exceed the maximal function $\psi^{+}$; similar statements hold for $Q_{S},Q_{E},Q_{W}$. See Fig. 7. Also, it follows from [31, Th. 4.2] that $\partial Q_{N}\cap\stackrel{{\scriptstyle\circ}}{{Q}}$ is the graph of a concave function; analogously, $\partial Q_{a}\cap\stackrel{{\scriptstyle\circ}}{{Q}}$ for $a\in\{E,S,W\}$ is the graph of a concave function in a reference frame rotated clockwise by $\pi/4,\pi/2$ and $3\pi/4$ respectively.

Because of the definition of $\psi^{\pm}$, we see that

[TABLE]

Note that $\mathcal{F}_{0}$ is convex.

Before proving that $\mathcal{F}_{R}$ is non-empty, let us show that $\mathcal{F}$ is non-empty. Let $\rho^{(a)},a\in\{N,E,S,W\}$ be the the gradient of $\phi_{a}(\cdot)$ (these are also the four corners of $N(P)$) and $\ell^{(i)},i\in\{NE,SE,SW,NW\}$ the open segment connecting $\rho^{(N)}$ to $\rho^{(E)}$ etc. Remark that if $x\in\mathcal{F}_{0}$, then $\nabla\psi_{\mathtt{w}}(x)$ cannot coincide with any of the slopes $\rho^{(a)},a\in\{N,E,S,W\}$. In fact, thanks to Proposition 3.3, in this case one would have

[TABLE]

which contradicts the fact that $x\in\mathcal{F}_{0}$. Therefore, we have that

[TABLE]

with

[TABLE]

In general, the region $\mathcal{F}$ is a proper subset of $\mathcal{F}_{0}$, see Fig. 8.

Using also the second statement in Proposition 3.6, we conclude that if (by contradiction) $\mathcal{F}$ is empty, then necessarily $\mathcal{F}_{0}$ must coincide with one of the four sets $\mathcal{F}^{(i)}$. To fix ideas, say that $\mathcal{F}_{0}=\mathcal{F}^{(NW)}$, i.e. everywhere in $\mathcal{F}_{0}$, $\nabla\psi_{\mathtt{w}}$ is a non-trivial convex combination of $\rho^{(W)}=(-n,0)$ and $\rho^{(N)}=(0,-n)$. Let $\gamma$ be the curve along $\partial\mathcal{F}_{0}$ from point $A$ to point $B$, as in Fig. 9, and let $t_{p}$ be the tangent vector at a point $p\in\gamma$.

From the definition of $Q_{S},Q_{W}$ one has that the directional derivative of $\psi_{\mathtt{w}}$ in direction $t_{p}$ equals $t_{p}\cdot g_{p}$, with $g_{p}\in\ell^{(SE)}$. On the other hand, if $\gamma^{\prime}$ is a curve from $A$ to $B$ that runs slightly inside $\mathcal{F}_{0}$ at distance $\delta$ from $\gamma$, we have that the directional derivative along $\gamma^{\prime}$ at a point $p^{\prime}$ equals $t^{\prime}_{p^{\prime}}\cdot\nabla\psi_{\mathtt{w}}(p^{\prime})$, with $\nabla\psi_{\mathtt{w}}(p^{\prime})\in\ell^{(NW)}$, because $\mathcal{F}_{0}=\mathcal{F}^{(NW)}$ by assumption. Taking $\delta\to 0$, one easily sees that these two facts are not compatible with $\psi_{\mathtt{w}}$ being continuous along $\gamma$. This proves that $\mathcal{F}$ is not empty.

Finally, the fact that $\mathcal{F}_{R}\neq\emptyset$ follows easily from $\mathcal{F}\neq\emptyset$. In fact, if $\mathcal{F}_{R}$ were empty, then $\nabla\psi_{\mathtt{w}}(x)$ would belong to $\mathcal{S}$ for every $x\in\mathcal{F}$ and (because of Proposition 3.6) it would actually take a constant value $\bar{\rho}$ on $\mathcal{F}$. If $\mathcal{F}=Q$, this is a contradiction since the affine function with slope $\bar{\rho}$ cannot match the boundary datum $\psi_{\partial Q}$. If on the other hand $Q\setminus\mathcal{F}\neq\emptyset$, then take a sequence of points $x_{i}\in\mathcal{F}$ and a sequence $y_{i}\in Q\setminus\mathcal{F}$ that have the same limit in the interior of $Q$. One has $\nabla\psi_{\mathtt{w}}(x_{i})=\bar{\rho}$ while $\nabla\psi(y_{i})\in\partial N(P)$, which contradicts Proposition 3.6. ∎

4. Properties of $v_{\mathtt{w}}(\rho)$

We start with the following statement, whose proof is given below:

Proposition 4.1**.**

The function $\rho\mapsto v_{\mathtt{w}}(\rho)$ is $C^{\infty}$ on $\mathcal{R}$.

Remark 4.2**.**

We know that the determinant of the Hessian matrix $J(x)$ of $\psi_{\mathtt{w}}$ is negative or zero on the rough region $\mathcal{F}_{R}$; if we knew that the inequality is everywhere strict, $C^{\infty}$ continuity of $v_{\mathtt{w}}(\cdot)$ would easily follow from formula (4.2) below and from further derivation w.r.t. $\rho$. On the other hand, non-vanishing of $J(x)$ in the rough region is not a general property of macroscopic shapes of dimer models. For instance, for the dimer model on the honeycomb graph with uniform weights, one can verify from the explicit solution [10] that the macroscopic shape $\psi$ in a hexagonal domain has a Hessian with strictly negative determinant in the whole rough region, except at a single point (the center of the domain), where all entries of the Hessian matrix are zero. To overcome this problem, for the proof of Proposition 4.1 we will not rely directly on analytic properties of the limit shapes, bur rather on the definition (2.14) of the speed and on the properties of the dimer measure $\pi_{\rho,\mathtt{w}_{j}}$ under the dynamics $\{\mathtt{w}_{j}\}_{j\geq 0}$ of the edge weights (“spider move dynamics”).

From Proposition 4.1 and the formula (2.16) for the speed, we deduce

[TABLE]

By the way, this shows that $x_{\mathtt{w}}(\rho)$ is unique for $\rho$ in the rough region. This formula also allows to prove that the speed is not $C^{1}$ at smooth slopes. Indeed, we know from Theorem 3.4 that for every $\bar{\rho}\in\mathcal{S}$, there exists $x_{\mathtt{w}}(\bar{\rho})$ in the interior of $Q$, where the slope of $\psi_{\mathtt{w}}$ is $\bar{\rho}$. Moreover, it is known [1] that, since the boundary condition $\psi|_{\partial Q}$ is “natural” (cf. footnote 1), the set $B_{\bar{\rho}}:=\{x\in Q:\nabla\psi_{\mathtt{w}}(x)=\bar{\rho}\}$ is a closed set with non-empty interior. Letting $x\in\mathcal{F}_{R}$ approach different points of $B_{\bar{\rho}}$ (so that $\nabla\psi_{\mathtt{w}}(x)$ approaches $\bar{\rho}$, by continuity of $x\mapsto\nabla\psi_{\mathtt{w}}(x)$), we see from (4.1) that $Dv_{\mathtt{w}}(\rho)$ does not have a unique limit as $\rho\to\bar{\rho}$.

From (4.1) we see also that, for $\rho\in\mathcal{R}$,

[TABLE]

where the $2\times 2$ Jacobian matrix $J(\cdot)$ is as in (3.17). We already know that $\det(J(x))\leq 0$, and the fact that the speed is $C^{2}$ means that the inequality is strict. In particular,

[TABLE]

as wished.

Proof of Proposition 4.1.

Let $f$ be an even face. From (2.14) and (2.4) one has, with $\mathtt{w}\equiv\mathtt{w}_{0}$,

[TABLE]

On the other hand, recall from (2.5) and (2.6) that $H(\eta),V(\eta)$ are sums of dimer indicator functions. From the determinantal structure of the measures $\pi_{\rho,\mathtt{w}}$, one has an explicit expression for the probability that an edge $e$ is occupied. Assume that the white endpoint of $e$ is in the fundamental domain $D_{m_{1},m_{2}}$ (that is the translation of $D_{0,0}$ by $2m_{1}n$ in the horizontal direction and by $2m_{2}n$ in the vertical one) and that, modulo this translation, it is equivalent to the white vertex $x$ of the fundamental domain $D_{0,0}$. Similarly, assume that the black endpoint is in $D_{\ell_{1},\ell_{2}}$ and that it is equivalent to the black vertex $y$ in $D_{0,0}$. Then,

[TABLE]

where $\mathbb{K}_{\mathtt{w}}(e)$ equals the $\mathtt{w}$-weight of $e$, times the complex unit $\mathrm{i}$ if the edge is vertical, while

[TABLE]

We recall that $K(z,w)$ is the $2n^{2}\times 2n^{2}$ Kasteleyn matrix of the fundamental domain $D_{0,0}$ (recall Section 2.2) and $B=B(\rho)=(B_{1}(\rho),B_{2}(\rho))$ is the value that realizes the supremum in (2.8). For $\rho=(\rho_{1},\rho_{2})\in\mathcal{R}$ the maximizer is unique and the relation between $\rho$ and $B(\rho)$, through

[TABLE]

is a $C^{\infty}$ diffeomorphism between $\mathcal{R}$ and $A(P)\subset\mathbb{R}^{2}$ (the amoeba of $P$, $A(P)$, defined as the image of the curve $P(z,w)=0$ in $\mathbb{C}^{2}$ under the map $(z,w)\mapsto(\log|z|,\log|w|)$) [22]. We will prove:

Lemma 4.3**.**

The r.h.s. of (4.6) is a $C^{\infty}$ function of $B$.

As a consequence, (4.5) and therefore the sum in (4.4), for every fixed $k$, are $C^{\infty}$ functions of $\rho$. To conclude the proof of the proposition, we will prove:

Lemma 4.4**.**

Let $\mathtt{w}=\mathtt{w}_{j}$. The derivatives (of any order) of (4.5) w.r.t. $B$ can be bounded uniformly w.r.t. the index $j$.

The smoothness claim for $v_{\mathtt{w}}$ then easily follows from (4.4). ∎

Proof of Lemma 4.3.

Assume without loss of generality (by translation invariance) that $\ell_{1}=\ell_{2}=0$. Write

[TABLE]

with $P(z,w)=\det K(z,w)$ the characteristic polynomial and $Q(z,w)$ (that is also a Laurent polynomial in $z,w$) the cofactor $(x,y)$ of $K(z,w)$, so that (4.6) reduces to

[TABLE]

The prefactor of the integral is smooth and will be dropped; also, we write $\tilde{Q}$ for $Q\times e^{i\theta m_{1}+i\phi m_{2}}$. If $B=B(\rho)$ as in (4.7) with $\rho\in\mathcal{R}$, it is known that $(\theta,\phi)\mapsto P(e^{B_{1}+i\theta},e^{B_{2}+i\phi})$ has two distinct simple zeros [21], call them $(\theta^{\omega},\phi^{\omega}),\omega=\pm$. Write

[TABLE]

where $P^{\omega}_{1}$ is the first-order Taylor expansion around $(\theta^{\omega},\phi^{\omega})$. The zeros $(\theta^{\omega},\phi^{\omega})$ and also $a^{\omega},b^{\omega}$ are real analytic functions of $B_{1},B_{2}$, and the ratio $a^{\omega}/b^{\omega}$ is not real. Write

[TABLE]

where

[TABLE]

and $\chi:\mathbb{R}\mapsto[0,1]$ is a $C^{\infty}$ function that equals $1$ (resp. [math]) when its argument is smaller than $\epsilon$ (resp. larger than $2\epsilon$), with $\epsilon$ sufficiently small so that the supports of $f^{\pm}$ are disjoint. The integral of

[TABLE]

is $C^{\infty}$ w.r.t. $B$. Now look at the integral of $f^{\omega}\tilde{Q}/P$. Suppose we want to prove it is $C^{k}$ w.r.t $B$. Write

[TABLE]

with $\tilde{Q}^{\omega}:=\tilde{Q}(\theta^{\omega},\phi^{\omega})$ and $\hat{Q}^{\omega}:=\tilde{Q}-\tilde{Q}^{\omega}$. Write

[TABLE]

Since the ratio $a^{\omega}/b^{\omega}$ is not real, the Jacobian of the change of variables $(\theta,\phi)\leftrightarrow(X,Y)$ is non-singular. One has then

[TABLE]

which is zero by symmetry. Next look at

[TABLE]

where, with some abuse of notation, we write

[TABLE]

The constant prefactor has a $C^{\infty}$ (in fact, real analytic) dependence on $B$. Also, $\hat{Q}$ is a polynomial with real analytic coefficients and it vanishes at least linearly when $(X,Y)$ tends to zero. Then, it is easy to deduce that (4.18) is a $C^{\infty}$ function of $B$. Finally, we look at

[TABLE]

with the same convention as in (4.19). Since $R^{\omega}$ is at least quadratic for $X,Y$ close to zero, the derivatives of order $k$ (w.r.t. the components of $B$) of the integrand are upper bounded by

[TABLE]

uniformly for $B$ in compact sets of the amoeba $A(P)$. The function (4.21) is integrable and the claim of the Lemma easily follows. ∎

Proof of Lemma 4.4.

We have seen that for each choice of $\mathtt{w}$, the derivatives of (4.5) w.r.t. $B$ are bounded. Now we let $\mathtt{w}=\mathtt{w}_{j}$ and we need to show uniformity of the bounds w.r.t. $j$. It is immediate to see that uniformity follows if all edge weights stay bounded away from [math] and $\infty$, uniformly in $j$.

Let us recall that the probability measure $\pi_{\rho,\mathtt{w}}$ depends on the edge weights only modulo gauge transformations [19, Sec. 3.2]. That is, if edge weights $\mathtt{w}$ are changed as $\mathtt{w}(e)\mapsto\mathtt{w}(e)f(b)g(w)$, with $e$ the edge with black/white endpoints $b/w$ and $f/g$ two non-vanishing functions defined on black/white vertices, then the measure is unchanged. In the ($2n\times 2n$) periodic setting with fundamental domain $D_{0,0}$ as in the present work, the knowledge of the edge weights modulo gauge is equivalent to the knowledge of:

1. (1)

the “face weights”: for each of the $4n^{2}$ faces $f$ of the fundamental domain $D_{0,0}$, one lets $\mathtt{w}(f)$ be the alternate product

[TABLE]

with $e_{1},\dots,e_{4}$ the four boundary edges of $f$ labeled cyclically clocwise, with $e_{1}$ chosen such that it is clockwise oriented from white to black endpoint. Actually, the product of face weights over all faces gives $1$, so we need to know only $4n^{2}-1$ of them. 2. (2)

the “magnetic coordinates”, i.e. the alternate product $\mathtt{W}_{1}$ (resp. $\mathtt{W}_{2}$) of the weights of the edges belonging to a cycle on $D_{0,0}$ with winding number $(1,0)$ (resp. $(0,1)$).

If the face weights, as well as $\mathtt{W}_{1},\mathtt{W}_{2}$, are all bounded away from [math] and $+\infty$, then there exists a suitable gauge such that edge weights are also all bounded away from [math] and $+\infty$.

When the weights $\mathtt{w}$ evolve along the sequence $\{\mathtt{w}_{j}\}_{j\geq 0}$ associated to the shuffling algorithm, the magnetic coordinates $\mathtt{W}_{1},\mathtt{W}_{2}$ stay constant [16]. This is related to the fact that the measure $\pi_{\rho,\mathtt{w}_{j}}$ is mapped to $\pi_{\rho,\mathtt{w}_{j+1}}$ and the slope $\rho$ is unchanged, recall Proposition 2.2. On the other hand, the face weights do change with $j$: in general they are not periodic in time but only quasi-periodic, they stay in a compact set (that depends on the initial weights $\mathtt{w}_{0}$) and they approach neither zero nor infinity. This can be extracted from the classical integrability of the dynamics of the face weights under the spider moves [16] (cf. also [14] and [22, Sec. 3]). More explicitly, the spider move preserves the spectral curve and for positive-real-valued edge weights, the common level set of the Hamiltonians is homeomorphic to a finite cover of the product of the compact ovals of the spectral curve. ∎

Acknowledgements We would like to thank Sanjay Ramassamy for discussions on the Goncharov-Kenyon dynamical system, and Erik Duse for sharing the results of [1] before publication and for discussions on the variational principle. We would also like to thank the referees for their careful reading and constructive comments. F.T. was partially supported by the CNRS PICS grant 151933, by ANR-15-CE40-0020-03 Grant LSD and ANR-18-CE40-0033 Grant DIMERS. S.C acknowledges the support of EPSRC grant EP /T004290/1.