On the Shuffling Algorithm for Domino Tilings

# On the Shuffling Algorithm for Domino Tilings

Eric Nordenstam
###### Abstract.

We study the dynamics of a certain discrete model of interacting particles that comes from the so called shuffling algorithm for sampling a random tiling of an Aztec diamond. It turns out that the transition probabilities have a particularly convenient determinantal form. An analogous formula in a continuous setting has recently been obtained by Jon Warren studying certain model of interlacing Brownian motions which can be used to construct Dyson’s non-intersecting Brownian motion.

We conjecture that Warren’s model can be recovered as a scaling limit of our discrete model and prove some partial results in this direction. As an application to one of these results we use it to rederive the known result that random tilings of an Aztec diamond, suitably rescaled near a turning point, converge to the GUE minor process.

## 1. Introduction

There has been a lot of work in recent years connecting tilings of various planar regions with random matrices. One particular model that has been intensely studied is domino tilings of a so called Aztec diamond. One way to analysing that model, [johansson:arctic_circle, johansson:discrete_orthogonal, johansson:gue_minors], is to define a particle process corresponding to the tilings so that uniform measure on all tilings induces some measure on this particle process.

In this article we will study the so called shuffling algorithm, described in [elkies:alternating_sign_matrices_II, propp:generalized_domino], which in various variants can be used either to count or to enumerate all tilings of the Aztec diamond or to sample a random such tiling.

The sampling of a random tiling by this method is an iterative process. Starting with a tiling of an order Aztec diamond, a certain procedure is performed, producing a random tiling of order . This procedure is usually described in terms of the dominoes which should be moved and created according to a certain procedure. We will instead look at this algorithm as a certain dynamics on the particle process mentioned above.

The detailed dynamics of the particle process will be presented in section 2 and how it is obtained from the traditional formulation of the shuffling algorithm is presented in section 4. For now, consider a process for , , , …, where . The quantity represents the position of the :th particle on line after steps of the shuffling algorithm have been performed. (The reason for the is technical convenience.) We will show that

###### Theorem 1.1.

For fixed , consider only the component from rescaled according to

 (1) ~Xni(t)=Xni(Nt)−12Nt12√N

and defined by linear interpolation for non-integer values of . The process converges to a Dyson Brownian motion with all particles started at the origin as , in the sense of convergence of finite dimensional distributions.

The full process has remarkable similarities to, and is we believe a discretization of, a process studied recently by Warren, [warren:dyson_brownian_motions]. It consists of many interlaced Dyson Brownian motions and is here briefly described in section 3. We will denote that process . There is reason to believe the following.

###### Conjecture.

Consider the process rescaled according to

 (2) ~Xni(t)=Xni(Nt)−12Nt12√N

and defined by linear interpolation for non-integer values of . The process converges to Warren’s process as , in the sense of convergence of finite dimensional distributions.

The key to our asymptotic analysis of the shuffling algorithm is that the transition probabilities of can be written down in a convenient determinantal form, see proposition 3.2. These formulas mirror beautifully formulas obtained by Warren.

As an application of our results we will use it to rederive an asymptotic result about random tilings near the point where the arctic circle touches the edge of the diamond. This result was first stated in [johansson:arctic_circle] and proved in [johansson:gue_minors].

Recall that the Gaussian Unitary Ensemble, or GUE for short, is a probability measure on Hermitian matrices with density where is a normalisation constant that depends on the dimension of the matrix. Let a GUE matrix and denote its principal minors by . Let be the eigenvalues of . Then is the so called GUE minor process.

###### Theorem 1.2 (Theorem 1.5 in [johansson:gue_minors].).

Let the -valued process be a rescaled version of with

 (3) ~Xji(t)=Xji(t)−t212√t.

Then as in the sense of weak convergence of probability measures.

To put this in perspective, let us note that a similar result for lozenge tilings is known from Okounkov and Reshetikhin [okounkov:birth]. They discuss the fact that, for quite general regions, that close to a so called turning point the GUE minor process can be obtained in a limit. A turning point is, just as in our situation, where the disordered region is tangent to the domain boundary.

## 2. The Aztec Diamond Particle Process

We will here content ourselves with stating the rules of the particle dynamics that we will study. The reader will in section 4 find a description the traditional formulation of the shuffling algorithm and how that relates to the formulas below.

Consider the process for , , , …, where . It satisfies the initial condition

 (4) Xk(0)=¯xk

where for . At each time the process fulfils the interlacing condition

 (5) Xki(t)≤Xk−1i(t)

and evolves in time according to

 (6) X11(t) =X11(t−1)+β11(t) Xj1(t) =Xj1(t−1)+βj1(t) −1{Xj1(t−1)+βj1(t)=Xj−11(t)+1} for j≥2 Xjj(t) =Xjj(t−1)+βjj(t) +1{Xj1(t−1)+βj1(t)=Xj−1j−1(t)} for j≥2 Xji(t) =Xji(t−1)+βji(t) −1{Xji(t−1)+βji(t)=Xj−1i(t)+1} +1{Xji(t−1)+βji(t)=Xj−1i−1(t)} for j≥3 and 1

for , , …where all the are i.i.d. unbiased coin tosses, satisfying .

One way to think about this is that at each time , this is a set of particles on lines. The :th line has particles on it at positions , …, . At each time step each of these particles either stays or jumps one unit step forward independent of all others except that the particles on line can force particles on line to jump or to stay to enforce the the interlacing condition (5). Also note that the interlacing implies that at each time , i.e. two particles cannot occupy the same space at the same time.

As mentioned we can write down transition probabilities for this process on a particularly convenient determinantal form. Define such that if and otherwise. Let us first introduce some notation.

 (ϕ∗ψ)(x) =∑s+t=xϕ(s)ψ(t) (Convolution product) ϕ(0) =δ0 ϕ(n) =ϕ(n−1)∗ϕ for n=1,2,… Δϕ =(δ0−δ1)∗ϕ (Backward difference) Δ−1ϕ(x) =x∑y=−∞ϕ(y) ¯Δϕ =(−δ0+δ−1)∗ϕ (Forward difference) ¯Δ−1ϕ(x) =x−1∑y=−∞ϕ(y)

Let . For and , define

 (7) \mathpzcqnt((x,y),(x′,y′))=det[At(x,x′)Bt(x,y′)Ct(y,x′)Dt(y,y′)]

where

• is an -matrix where element is ,

• is an -matrix where element is ,

• is an -matrix where element is and

• is an -matrix where element is .

Let and for let

 (8) hn(x)=∏i

Finally, after all this notation, we can state a result.

###### Theorem 2.1.

The transition probabilities of from the process above are

 (9) \mathpzcqk,+t((x,y),(x′,y′)):=hk(y′)hk(y)\mathpzcqkt((x,y),(x′,y′))

that is

 (10) P[(Xk+1(s+t),Xk(s+t))=(x′,y′);(Xk+1(s),Xk(s))=(x,y)]=\mathpzcqk,+t((x,y),(x′,y′)).

A proof is given in section 6 and the reason I defined as opposed to defining directly will become obvious in the next section.

Given the exact expressions above it is a very straightforward computation to integrate out the component in expression (9). We find that the transition probabilities of from the process above is

 (11) \mathpzcpk,+t(y,y′):=hk(y′)hk(y)\mathpzcpkt(y,y′)

where given above. We recognise this as the transition probability for random walks conditioned never to intersect, a fact that is so important we state it properly.

###### Corollary 2.2.

The component of is a discrete Dyson Brownian motion of particles started at .

This fits nicely with theorem 1.1. The component from simply simple symmetric random walks conditioned never to intersect, their limit is Brownian motions conditioned never to intersect, which is exactly what from Warren’s process is.

## 3. Interlacing Brownian motions

We will now digress a bit and summarise Warren’s work in [warren:dyson_brownian_motions], so as to see the similarities between his continuous process and our discrete process. The reader is referred to that reference for more details of the construction. Consider an -valued stochastic process satisfying an interlacing condition

 (12) X1(t)≤Y1(t)≤X2(t)≤⋯≤Yn(t)≤Xn+1(t),

and equations

 (13) Yi(t) =yi+βi(t∧τ), (14) Xi(t) =yi+γi(t∧τ)+L−i(t∧τ)−L+i(t∧τ)

where

and are independent Brownian motions,

,

and

 (15) L+i(t) =∫t01(Xi(s)=Yi(s))dL+i(s) L−i(t) =∫t01(Xi(s)=Yi−1(s))dL−i(s)

are twice the semimartingale local times at zero of and respectively.

This process can be constructed by first constructing the Brownian motions and and then using Skorokhod’s construction to push up from and down from . The process is killed when is reached, i.e. when two of the meet.

Warren then goes on to show that the transition densities of this process have a determinantal form similar to what we have seen in the previous section. Let and . Let .

Define for , and to be the determinant of the matrix

 (16) [At(x,x′)Bt(x,y′)Ct(y,x′)Dt(y,y′)]

where

is an -matrix where element is ,

is an -matrix where element is ,

is an -matrix where element is and

is an -matrix where element is .

###### Proposition 3.1 (Prop 2 in [warren:dyson_brownian_motions]).

The process killed at time has transition densities , that is

 (17) qnt((x,y),(x,y))dx′dy′=Px,y[Xt∈dx′,Yt∈dx′;t<τ]

Warren goes on to condition the not to intersect via so called the Doob -transform. The transition densities for the transformed process are given in terms of the those for the killed process by

 (18) qn,+t((x,y),(x′,y′))=hn(y′)hn(y)qnt((x,y),(x′,y′)).

He also shows that you can start all the and of the transformed process at the origin by giving a so called entrance law,

 (19) νnt(x,y):=n!Zn+1t−(n+1)2/2exp{−∑ix2i/(2t)}{∏i

that is, showing (lemma 4 of [warren:dyson_brownian_motions]) that this expression satisfies

 (20) νnt+s(x′,y′)=∫Wn,n+1νns(x,y)qn,+t((x,y),(x′,y′))dxdy.

It is possible to integrate out the components in that transition density and entrance law. The result is transition density

 (21) pn,+t(y,y′):=h(y′)h(y)detDt(y,y′)

and entrance law

 (22) μnt(y):=1Znt−n2/2exp{−∑iy2i/(2t)}{∏i

Now comes the interesting part. Let be the cone of points where . Warren defines a process taking values in such that

 (23) Xki(t)=xki+γki(t)+Lk,−i(t)−Lk,+i(t)

where the are independent Brownian motions and and are continuous, increasing processes growing only when and respectively and the special cases and are identically zero for all .

Think of this as essentially particles performing independent Brownian motions except that the particles in can push the particles in up or down to enforce the interlacing condition that the whole process should stay in .

This full process process can be constructed inductively as follows.

1. The process has transition densities and entrance law .

2. The process has transition densities and entrance law .

3. For the process is conditionally independent of given .

4. This implies (by some explicit calculations) that has transition densities and entrance law .

This argument shows that the following.

###### Proposition 3.2 (Warren).

There exists such a process started at the origin and it satisfies that for , …, , the process has entrance law and transition probabilities .

It is this process that is the continuous analog of our discrete process .

## 4. Shuffling algorithm

We will now show how relate some well known facts about sampling random tilings of an Aztec diamond before showing how to get the particle dynamics in section 2.

The Aztec diamond of order , denoted , is an area in the plane that is the union of those lattice squares that are entirely contained in . can be tiled in ways by dominoes of size . We will be interested picking a random tiling. By random tiling in this article we will always mean that all possible tilings given the same probability.

A key ingredient of almost all results concerning tilings of this shape is the realization that one can distinguish four kinds of dominoes present in a typical tiling. The obvious distinction to the casual observer is the difference between horizontal and vertical dominoes. These can be subdivided further. Colour the underlying lattice squares black and white according to a checkerboard fashion in such a way that the left square on the top line is black. Let a horizontal domino be of type N or north if its leftmost square is black, and of type S or south otherwise. Likewise let a vertical domino be of type W or west if its topmost square is black and type E or east otherwise. In figures 1 and 2 the S and E type dominoes have been shaded for convenience.

One way of sampling from this measure is the so called shuffling algorithm, first described in [elkies:alternating_sign_matrices_II], and very nicely explained and generalised in [propp:generalized_domino]. It is an iterative procedure that given a random tiling of and some number of coin-tosses, produces a random tiling of a diamond of . You start with the empty tiling on and you repeat this process until you have a tiling of the desired size. It is a theorem that this procedure gives all tilings equal probability, provided that the coin-tosses we have made along the way are fair.

The algorithm works in three stages. Start with a tiling .

Destruction:

All blocks consisting of an S-domino directly above an N-domino are removed. Likewise all blocks of consisting of an E-domino directly to left of a W-domino are removed.

Shuffling:

All N, S, E and W-dominoes respectively move one unit length up, down, right and left respectively.

Creation:

The result is a tiling of a subset of . The empty parts can be covered in a unique way by 22 squares. Toss a coin to fill these with two horizontal or two vertical dominoes with equal probability.

Figure 1 illustrates the process. In the leftmost column there are tilings of successively larger diamonds. From column one to column two, the destruction step is carried out. From there to the third column, shuffling is performed. These figures contain several dots which will concern us later in this exposition. The creation step of the algorithm applied to a diamond in the last column gives (with positive probability) the diamond in the first column on the next row.

To study more detailed properties of random tilings it is useful to introduce a coordinate system suited to the setting and a particle process such that the possible tilings correspond to particle configurations.

In the left picture in figure 2, the S and E type dominoes are shaded and a coordinate system is imposed on the tiling. For each tile there is exactly one of the lines and exactly one of the lines that passes through its interior. Indeed we can uniquely specify the location of a tile by giving its coordinates and type (N, S, E or W). You can see that along the line there are exactly shaded tiles, for where is the order of the diamond. The obvious generalisation of that statement is true for tilings of for any . We shall call the occurrence of a shaded tile a particle. The right picture in figure 2 is the same tiling but with dots marking the particles.

Just to fix some notation, let be the -coordinate of the :th particle along the line . It is clear from the definitions that these satisfy an interlacing criterion,

 (24) xji≤xj−1i≤xji+1.

We will now see how the shuffling algorithm described above acts on these particles.

It turns out that the positions of the particles is uniquely determined before the creation stage of the last iteration of the shuffling algorithm, and we have marked these with dots in the last column in figure 1. As can be seen in that figure, running the shuffling algorithm to produce tilings of successively larger Aztec diamonds imposes certain dynamics on these particles. That is the central object of study in this article.

Let us first consider the trajectory of . As can easily be seen in figure 1, on the line there are always a number of W-dominoes, then the particle, then a number of N-dominoes. Depending on whether the creation stage of the algorithm fills the empty space in between these with a pair of horizontal or vertical dominoes, either the particle stays or its -coordinate will increase by one in the next step. Thus the first particle performs the simple random walk

 (25) x11(t)=x11(t−1)+γ11(t).

were are independent coin tosses, i.e. , for and .

Consider now the particles on row . For , while it performs a random walk independently of , at each time either staying or adding one with equal probability. However, when there is equality, , then the particle must be represented by a vertical (S) tile. Thus it does not contribute to growth of the west polar region, thus the particle will remain fixed. In order to represent this as a formula, we subtract one if the particle attempts to jump past .

 (26) x21(t)=x21(t−1)+γ21(t)−1{x21(t−1)+γ21(t)=x11(t−1)+1}

Symmetry completes our analysis of this row with the relation

 (27) x22(t)=x22(t−1)+γ22(t)+1{x22(t−1)+γ22(t)=x11(t−1)}.

For the third row, our previous analysis applies to the first and last particle.

 (28) x31(t) =x31(t−1)+γ31(t)−1{x31(t−1)+γ31(t)=x21(t−1)+1} (29) x33(t) =x33(t−1)+γ33(t)+1{x33(t−1)+γ33(t)=x22(t−1)}

On between and there must be first a sequence of zero or more E dominoes, then , then a sequence of zero or more N dominoes. While is in the interior of this area it performs the customary random walk. It must interact with and in the same way as we have seen other particles interacting above.

So

 (30) x32(t)=x32(t−1)+γ32(t) −1{x32(t−1)+γ32(t)=x22(t−1)+1} +1{x32(t−1)+γ32(t)=x21(t−1)}.

The same pattern repeats itself evermore.

 (31) xj1(t) =xj1(t−1)+γj1(t)−1{xj1(t−1)+γj1(t)=xj−11(t−1)+1} (32) xjj(t) =xjj(t−1)+γjj(t)+1{xjj(t−1)+γjj(t)=xj−1j−1(t−1)} (33) xji(t) =xji(t−1)+γji(t)−1{xji(t−1)+γji(t)=xj−1j(t−1)+1} (34) +1{xji(t−1)+γji(t)=xj−1j−1(t−1)}.

with initial conditions for and .

In order to analyse this situation it is suitable to perform a change of variables,

 (35) Xji(t)=xji(t−j),

which gives the equations given in section 2.

## 5. Transition probabilities on two lines

In order to analyse the dynamics just described we follow Warren’s example and first consider just two lines at a time. What we do in this section is very similar to section 2 of [warren:dyson_brownian_motions].

Consider the -valued process process with components and , satisfying the equations

 (36) Yi(t+1) =Yi(t)+βi(t) X1(t+1) =X1(t)+α1(t)−1{X1(t)+α1(t)=Y1(t+1)+1} Xi(t+1) =Xi(t)+αi(t)+1{Xi(t)+αi(t)=Yi−1(t+1)} −1{Xi(t)+αi(t)=Yi(t+1)+1} Xn+1(t+1) =Xi(t)+αn+1(t)+1{Xn+1(t)+αn+1(t)=Yn(t+1)}

where and are i.i.d. coin tosses, s.t. . They evolve until the stopping time . At the time the process is killed and remains constant for all time after that. This is a very simple dynamics, each either stays or increases one independently of all others. The do the same but are sometimes pushed up or down by a or respectively so as to stay in the cone . This is the discrete analog of the process defined in section 3 of this paper.

###### Lemma 5.1.

For any ,

 (37) ∑(x′,y′)∈\mathpzcWn+1,n\mathpzcqn0((x,y),(x′,y′))f((x′,y′))=f((x,y)).
###### Proof.

Let and , , …, , . Equation (42) in [warren:dyson_brownian_motions] states that

 (38) det{1{zi≤z′j}i≥j−1{zi≤z′j}i

for . Applying the operator to both sides of that equality turns the left hand side into and the right hand side into , , …, . ∎

###### Proposition 5.2.

, for , are the transition probabilities for the process , i.e. for , ,

 (39) \mathpzcqnt((x,y),(x′,y′))=P(x,y)[X(t)=x′,Y(t)=y′;t<τ]
###### Proof.

Take some test function . Let

 (40) F(t,(x,y)):=∑(x′,y′)∈\mathpzcWn+1,n\mathpzcqnt((x,y),(x′,y′))f(x′,y′)

and

 (41) G(t,(x,y)):=E(x,y)[f(Xt,Yt);t<τ]

We want of course to prove that and are equal and we will do this by showing that they satisfy the same recursion equation with the same boundary values. By lemma 5.1 we already know that

 (42) F(0,⋅)≡G(0,⋅)≡f(⋅).

The master equation satisfied by is

 (43) G(t+1,(x,y))=122n+1∑ai,bi∈{0,1}G(t,x1+a1,y1+b1,x2+a2,…,yn+bn,xn+1+an+1).

This formula simply encodes the dynamics that each particle either stays or jumps forward one step. This needs to be supplemented with some boundary conditions that have to do with the interactions between particles.

When two of the -particles coincide, this corresponds to the event , which does not contribute to the expectation in (41). Thus

 (44) G(\mathpzct,…,yi−1=z,xi,yi=z,…):=0.

Also, the particle cannot jump past ,

 (45) G(\mathpzct,…,xi=z+1,yi=z,…):=G(\mathpzct,…,xi=z,yi=z,…)

and must not drop below ,

 (46) G(\mathpzct,…,yi=z,xi+1=z,…):=G(\mathpzct,…,yi=z,xi+1=z+1,…).

is uniquely determined from using the recursion equation and boundary values above. It follows that is uniquely defined by the recursion equation (43) and the boundary conditions (42,44,45,46).

Observe that, all functions , satisfy

 (47) 12(g(x)+g(x+1))=g(x)+12¯Δg(x).

Using this identity many times on (43) shows that

 (48) G(t+1,(x,y))=(1+12¯Δx1)(1+12¯Δy1)(1+12¯Δx2)⋯(1+12¯Δyn)(1+12¯Δxn+1)G(t,x1,y1,…,yn,xn+1)

which can be rewritten as

 (49) ¯Δ