On the Shuffling Algorithm for Domino Tilings

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.


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, [4][5][6], 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 [3,11], 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 n − 1 Aztec diamond, a certain procedure is performed, producing a random tiling of order n. 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 X (t) = (X 1 (t), . . . , X m (t)) for t = 0, 1, 2, . . . , where X k (t) = (X k 1 (t), . . . , X k k (t)) ∈ Z k . The quantity X j i (t) represents the position of the i:th particle on line j after t − j steps of the shuffling algorithm have been performed. (The reason for the t − j is technical convenience.) We will show that Theorem 1.1. For fixed k, consider only the component X k (t) from X (t) rescaled according to and defined by linear interpolation for non-integer values of N t. The processX n (t) converges to a Dyson Brownian motion with all particles started at the origin as N → ∞, in the sense of convergence of finite dimensional distributions.
The full process (X (t)) t=0,1,... has remarkable similarities to, and is we believe a discretization of, a process studied recently by Warren, [12]. It consists of many interlaced Dyson Brownian motions and is here briefly described in section 3. We will denote that process (X(t)) t≥0 . There is reason to believe the following.
Conjecture. Consider the process (X (t)) t=0,1,... rescaled according to (2)X n i (t) = X n i (N t) − 1 2 N t 1 2 √ N and defined by linear interpolation for non-integer values of N t. The processX (t) converges to Warren's process X(t) as N → ∞, 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 (X k , X k+1 ) 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 [6] and proved in [4].
Recall that the Gaussian Unitary Ensemble, or GUE for short, is a probability measure on Hermitian matrices with density Z −1 n e − Tr H 2 /2 where Z n is a normalisation constant that depends on the dimension n of the matrix. Let H = (h rs ) 1≤r,s≤n a GUE matrix and denote its principal minors by H j = (h rs ) 1≤r,s≤j . Let λ j = (λ j 1 , . . . , λ j j ) be the eigenvalues of H j . Then Λ = (λ 1 , . . . , λ n ) ∈ R n(n+1)/2 is the so called GUE minor process. Theorem 1.2 (Theorem 1.5 in [4].). Let the R n(n+1)/2 -valued processX (t) = (X 1 (t), . . . ,X n (t)) be a rescaled version of X (t) with ThenX (t) → Λ as t → ∞ 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 [10]. 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.

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 (X (t)) = (X 1 (t), . . . , X n (t)) for t = 0, 1, 2, . . . , where X k (t) = (X k 1 (t), . . . , X k k (t)) ∈ Z k . It satisfies the initial condition At each time t the process fulfils the interlacing condition and evolves in time according to for t = 1, 2, . . . where all the (β j i (t)) i,j,t are i.i.d. unbiased coin tosses, satisfying P[β 1 One way to think about this is that at each time t, this is a set of particles on n lines. The k:th line has k particles on it at positions X k 1 , . . . , X k k . 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 k can force particles on line k + 1 to jump or to stay to enforce the the interlacing condition (5). Also note that the interlacing implies that X k i < X k i+1 at each time t, 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 δ i : Z → Z such that δ i (x) = 1 if i = x and δ i (x) = 0 otherwise. Let us first introduce some notation.
Theorem 2.1. The transition probabilities of (X k , X k+1 ) from the process X above are that is A proof is given in section 6 and the reason I defined q k t as opposed to defining q k,+ t directly will become obvious in the next section.
Given the exact expressions above it is a very straightforward computation to integrate out the x component in expression (9). We find that the transition probabilities of (X k ) from the process X above is where p k t (y, y ′ ) := D t (y, y ′ ) 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.
This fits nicely with theorem 1.1. The component X k from X simply k simple symmetric random walks conditioned never to intersect, their limit is k Brownian motions conditioned never to intersect, which is exactly what X k from Warren's process X is.

Interlacing Brownian motions
We will now digress a bit and summarise Warren's work in [12], 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 R n+1 ×R n -valued stochastic process (Q(t)) t≥0 = (X(t), Y (t)) t≥0 satisfying an interlacing condition and equations where are twice the semimartingale local times at zero of X i − Y i and X i − Y i−1 respectively. This process can be constructed by first constructing the Brownian motions β i and γ i and then using Skorokhod's construction to push X i up from Y i−1 and down from Y i . The process is killed when τ is reached, i.e. when two of the Y i 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.
and t > 0 to be the determinant of the matrix [12]). The process (X, Y ) killed at time τ has transition densities q n t , that is Warren goes on to condition the Y i not to intersect via so called the Doob h-transform. The transition densities for the transformed process are given in terms of the those for the killed process by He also shows that you can start all the X i and Y i of the transformed process at the origin by giving a so called entrance law, that is, showing (lemma 4 of [12]) that this expression satisfies It is possible to integrate out the X components in that transition density and entrance law. The result is transition density and entrance law Now comes the interesting part. Let K be the cone of points where the (γ k i ) i,k are independent Brownian motions and L k,+ i and L k,+ i are continuous, increasing processes growing only when X k i (t) = X k−1 i (t) and X k i (t) = X k i−1 (t) respectively and the special cases L k,+ k and L k,− 1 are identically zero for all k. Think of this as essentially n(n + 1)/2 particles performing independent Brownian motions except that the k particles in X k can push the particles in X k+1 up or down to enforce the interlacing condition that the whole process should stay in K.
This full process process can be constructed inductively as follows.
(1) The process (X k ) has transition densities p k,+ t and entrance law µ k t . (2) The process (X k , X k+1 ) has transition densities q k,+ t and entrance law ν k t . (3) For k = 2, . . . , n − 1 the process (X k+1 ) is conditionally independent of (X 1 , . . . , X k−1 ) given (X k ). (4) This implies (by some explicit calculations) that (X k+1 ) has transition densities p k+1,+ t and entrance law µ k+1 t . This argument shows that the following. Proposition 3.2 (Warren). There exists such a process X(t) started at the origin and it satisfies that for k = 1, . . . , n− 1, the process (X k , X k+1 ) has entrance law ν n t and transition probabilities q k,+ .
It is this process X that is the continuous analog of our discrete process X .

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 n, denoted A n , is an area in the plane that is the union of those lattice squares [a, a + 1] × [b, b + 1] ⊂ R that are entirely contained in {|x| + |y| ≤ n + 1}. A n can be tiled in 2 n(n+1)/2 ways by dominoes of size 2 × 1. 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 [3], and very nicely explained and generalised in [11]. It is an iterative procedure that given a random tiling of A n and some number of coin-tosses, produces a random tiling of a diamond of A n+1 . You start with the empty tiling on A 0 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 A n . Destruction: All 2 × 2 blocks consisting of an S-domino directly above an N-domino are removed. Likewise all 2 × 2 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 A n+1 . The empty parts can be covered in a unique way by 2×2 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 x lines and exactly one of the y lines that passes through its interior. Indeed we can uniquely specify the location of a tile by giving its coordinates (x, y) and type (N, S, E or W). You can see that along the line y = k there are exactly k shaded tiles, for y = 1 . . . 8 where 8 is the order of the diamond. The obvious generalisation of that statement is true for tilings of A n for any n. 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 x j i be the x-coordinate of the i:th particle along the line y = j. It is clear from the definitions that these satisfy an interlacing criterion, 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 x 1 1 . As can easily be seen in figure 1, on the y = 1 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 x-coordinate will increase by one in the next step. Thus the first particle performs the simple random walk were γ i j (t) are independent coin tosses, i.e. P [γ j i (t) = 1] = P [γ j i (t) = 0] = 1 2 , for t, j = 1, . . . and 0 ≤ i ≤ j.
Consider now the particles on row y = 2. For x 2 1 , while x 2 1 (t) < x 1 1 (t) it performs a random walk independently of x 1 1 , at each time either staying or adding one with equal probability. However, when there is equality, x 2 1 (t) = x 1 1 (t), 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 Symmetry completes our analysis of this row with the relation For the third row, our previous analysis applies to the first and last particle.
On y = 3 between x 2 1 and x 2 2 there must be first a sequence of zero or more E dominoes, then x 3 2 , then a sequence of zero or more N dominoes. While x 3 2 is in the interior of this area it performs the customary random walk. It must interact with x 2 1 and x 2 2 in the same way as we have seen other particles interacting above. So The same pattern repeats itself evermore.

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 [12].
Consider the W n+1,n -valued process process (Q n (t)) = (X(t), Y (t)) with components X 1 (t), . . . , X n+1 (t) and Y 1 (t), . . . , Y n (t), satisfying the equations (36) where α i (t) and β i (t) are i.i.d. coin tosses, s.t. P[α i (t) = 0] = P[α i (t) = 1] = 1 2 . They evolve until the stopping time τ = min {t : Y i (t) = Y i+1 (t) for some i ∈ {1, . . . , n − 1}}. At the time τ the process is killed and remains constant for all time after that. This is a very simple dynamics, each Y i either stays or increases one independently of all others. The X i do the same but are sometimes pushed up or down by a Y i−1 or Y i respectively so as to stay in the cone W n,n+1 . This is the discrete analog of the process Q defined in section 3 of this paper.  ((x, y)).
Proof. Let m = 2n + 1 and z 1 = x 1 , z 2 = y 1 , . . . , z m−1 = y n , z m = x n+1 . Equation (42) in [12] states that to both sides of that equality turns the left hand side into q n 0 ((x, y), (x ′ , y ′ )) and the right hand side into 1{z 1 = z ′ 1 , We want of course to prove that F and G 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 The master equation satisfied by G is 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 y i -particles coincide, this corresponds to the event t = τ , which does not contribute to the expectation in (41). Thus G(t+1, ·) is uniquely determined from G(t, ·) using the recursion equation and boundary values above. It follows that G is uniquely defined by the recursion equation (43) and the boundary conditions (42,44,45,46).
Observe that, all functions g : Z → R, satisfy Using this identity many times on (43) shows that (48) which can be rewritten as The boundary conditions can be rewritten in this notation as well, equations (44,45,46) can be rewritten to Now let us look at F . The observation (47) gives that φ * ψ = (1 + 1 2 ∆)ψ. In particular, which shows that F satisfies the same recursion (49) as G. Now let us take a look at the boundary values.
q n t ((x, y), (x ′ , y ′ )) is zero when y i = y i+1 because two of its rows are then equal. When y i = x i for some i then∆ xi q n t ((x, y), (x ′ , y ′ )) = 0 because two rows will be equal when you take the difference operator into the determinant. The same argument shows that∆ xi+1 q n t ((x, y), (x ′ , y ′ )) = 0 when y i = x i+1 . Applying this knowledge to the sum F , shows that F (t, (x, y)) = 0 when y i = y i+1 (53)∆ xi F (t, (x, y)) = 0 when Since F and G satisfy the same recursion equation with the same boundary values, they must be equal. Again, following the example of Warren, we observe that it is possible to condition the processes never to leave W n,n+1 via a so called Doob h-transform. See for example [8] for details about h-transforms for discrete processes. Let The h-transform of the process above has transition probabilities Just for the sake of notation, call the transformed process (Q n,+ (t)).
The idea now is to stitch together the process X from processes Q k,+ for k = 1, . . . , n − 1, just like Warren does in the continuous case. For this we need to establish some auxiliary results about Q n and Q n,+ .
One observation to make is that it is possible to integrate out the x variables from q n,+ t .
(58) p n,+ t (y, y ′ ) := where dx ′ is counting measure on W n+1 = {x ∈ Z n+1 : x 1 < · · · < x n+1 }. The reader might recognise p n,+ t , as a h-transformed version of the transition probabilities from the Lindström-Gessel-Viennot theorem. Thus this can be seen as the transition probabilities for a process on W n , where all n particles perform independent random walks but are conditioned never to intersect, i.e. never to leave W n . We state this as a proposition.
Now for a technical lemma. For x ∈ W n , let W n (x) = {y ∈ R n : x 1 ≤ y 1 < · · · ≤ y n < x n+1 } ⊂ Z n and for y ∈ W n (x) let It is not a difficult calculation to show that λ n (x, ·) is a probability measure on W n (x). Just rewrite h n (y) as a Vandermonde matrix and perform the summation over all y.
Proof. An elementary calculation given the explicit formula for q n,+ t .
Our proof of this is very similar to the proof of proposition 5 in [12].
Proof. Recall the the transition probabilities for (Q n,+ (t)) are q n,+ t .

Transition probabilities for the Aztec Diamond Process
Let us now return to the process (X (t)) that came from the shuffling algorithm. Observe in the recursions (6), the formulas that define X k+1 contain X k but not X j for j < k. Thus X k is conditionally independent of (X 1 , . . . , X k−1 ) given X k . Also, the dependence of X k+1 on X k is the same as the dependence of X on Y in Q k and in Q k,+ , see (36). This, together with theorem 5.5 lends itself to an inductive procedure for constructing the process X .
(1) The process (X k (t), t = 0, 1, . . . ) is started at X k (0) =x k and has transition probabilities governed by p k,+ for k = 1, 2, . . . , k, (2) By proposition 5.3, X k can be considered as the Y component of the process Q k,+ . By the observation above, the pair of processes (X k (t), X k+1 (t)) has the same distribution as Q k,+ started at (x k ,x k+1 ) and are thus governed by transition probabilities q k,+ , (3) The process X k+1 is conditionally independent of (X 1 , . . . , X k−1 ) given X k . (4) By theorem 5.5, the process X k+1 is governed by transition probabilities p k+1,+ and started at X k+1 (0) =x k+1 . This proves theorem 2.1.

Asymptotics
We shall now see some results that lend support to our conjecture that Warren's process X can be recovered as a scaling limit of X , the process from the Aztec diamond. Let us rescale time byt = N t and space byx i = 1 2 N t + 1 2 √ Nx i for i = 0, . . . , n in the above processes. First let us show how to recover the entrance law for Warrens process. Recall that the Aztec diamond process the discrete process starts at X n (0) =x n and X n+1 (0) =x n+1 .
as N → ∞ where ν n t (x, y) is given by (19). Proof. Observe that we can write (62) where dx n is counting measure on the space W (x n+1 ) which has only one element. Then we apply lemma 5.4.
which can be written explicitly as which can be evaluated by a formula of Krattenthaler (Theorem 26 of [9]). Applying Stirling's approximation to the result shows our theorem with Z n = (2π) n/2 j<n j!. Likewise, Warren's expression for q n t can be recovered as a scaling limit from our expression for q n t . By Stirling's approximation, uniformly on compact sets as N → ∞ where uniformly on compact sets as N → ∞.
Proof. Just insert the limit relations given above for φ and ϕ in the explicit expression for q n t . Finally, one of the main results of this article is the following. Theorem 7.3. The process (X k (t), X k+1 (t)) from X , extended by interpolation to non-integer timest, converges in the sense of finite dimensional distributions to the process (X k (t), X k+1 (t)) from X.
Proof. For times t 1 , . . . , t m , and compact sets A 1 , . . . , A m ∈ W n,n+1 , we need to study (69) lim ,ȳ), (x 1 ,ỹ 1 ))q ñ t2−t1 ((x 1 ,ỹ 1 ), (x 2 ,ỹ 2 )) · · · q ñ tm−tm−1 ((x m−1 ,ỹ m−1 ), (x m ,ỹ m )) where of course dx i dỹ i is point measure on W n,n+1 . This is a case of a Riemann-sum converging to an integral. By the uniform convergence of q ñ t in lemmas 7.1 and 7.2, we can interchange the order of integration and taking the limit. This gives y 1 ), (x 2 , y 2 )) · · · q n t k −t k−1 ((x k−1 , y k−1 ), (x k , y k )) which proves our theorem. Theorem 1.1 follows from theorem 7.3 by just restricting to (X k ). We feel that given this theorem together with the fact that X k+1 is conditionally independent of X 1 , . . . , X k−1 given X k , lends a lot of credibility to the conjecture in the introduction.
We can also say something about the limit at a fixed time. Let K be the cone of points For each x n ∈ W n we will denote by K (x n ) the set of all (x 1 , . . . , x n−1 ) such that (x 1 , . . . , x n−1 , x n ) ∈ K . The number of points in K (x n ) is It follows from the characterisation in section 6 that at a fixed time the distribution of X n−1 (t) given X n (t) is λ n−1 (X n (t), ·). Together with the conditional independence noted in that section this implies that the distribution of (X 1 (t), . . . , X n−1 (t)) given X n (t) is uniform in K (X n (t)). So the probability distribution of X (t) is (73) m n t (x) = p n t (x n , x n ) χ(x 1 , x 2 ) . . . χ(x n−1 , x n ) card(K (x n )) where χ(x k , x k+1 ) is one iff x k+1 i ≤ x k i < x k+1 i+1 for all i = 1, . . . , k − 1 and zero otherwise. For x n ∈ W n , define K(x n ) as the set of (x 1 , . . . , x n−1 ) where x k ∈ R k satisfying x k+1 i+1 . The n(n − 1)/2-dimensional volume of K(x n ) is (74) vol(K(x n )) = h n (x n ) k<n k! . Theorem 7.4. Consider the process (X (t)) t∈Z + under the rescaling As N → ∞,X → Λ weakly where Λ has distribution (76) µ n 1 (x n ) χ(x 1 , x 2 ) . . . χ(x n−1 , x n ) vol(K(x n ) .
The expression for µ n 1 is given in equation (22). Proof. Put in the correct rescaling in (73) and perform a computation that is practically the same as that in the proof of 7.1.
This distribution Λ also happens to be the distribution of X(1), which is consistent with our conjecture. It has been studied [1] and is the distribution of the GUE minor process mentioned in the introduction. This proves theorem 1.2.

Closing Remarks
Looking at the expression of transition probabilities q n,+ t it is natural to ask the question, what happens if we plug in a different φ than 1 2 (δ 0 + δ 1 ) into that determinantal formula? It turns out that for many other φ this gives a valid transition probability, although we do not fully understand why the Doob h-conditioning still works in that case. It would be interesting to see sufficient and necessary conditions on φ for this construction to work.
Another reference worthy of attention is [7], by Johansson. He considers only geometric jumps and studies only the top particles (X 1 1 (t), X 2 2 (t), . . . , X n n (t)) from our model, with a slight change of variables that is of no real importance. He not only writes down transition probabilities, but also recovers the top particles Warren's process X as the limit of his process properly rescaled.
All the results proved in this article can be generalised to φ = pδ 0 + qδ 1 where p + q = 1. It is also not very difficult given my calculations to write down transition probabilities for the top particles and to rescale that to obtain the top particles in Warren's continuous process, analogous to Johansson [7] but with Bernoulli jumps. We have not written included that calculation here since we don't think it ads much to our knowledge of these processes, but it is a fact that adds to the plausability of the conjecture of this article.