Shuffling algorithm for coupled tilings of the Aztec diamond

In this article we define a generalization of the domino shuffling algorithm for tilings of the Aztec diamond to the interacting $k$-tilings recently introduced by S. Corteel, A. Gitlin, and the first author. We describe the algorithm both in terms of dynamics on a system of colored particles and as operations on the dominos themselves.


Introduction
Domino tilings of the Aztec diamond were first introduced by Elkies, Kuperberg, Larsen, and Propp [11] in their study of alternating-sign matrices.See Figure 1 for an example of a domino tiling of rank 3.In their work, the authors introduced domino shuffling, an algorithm by which one can generate a tiling of the Aztec diamond of rank (N + 1) from a tiling of the Aztec diamond of rank N .Using this they were able to derive a recursive formula for the number of tilings.Solving this recursion they found the beautiful result: Theorem 1.1.The number of tilings of the Aztec diamond of rank N is given by 2 ( N +1 2 ) .
The shuffling algorithm has since proved very useful in the study of these tilings.An immediate benefit is that shuffling allows for efficient exact sampling with arbitrary weights [22].Furthermore, it has also been used as a tool for asymptotic analysis, in the following way: One central result in the study of tilings of the Aztec diamond is the arctic circle theorem [15].It states that for large N , a uniformly random tiling exhibits a brickwork pattern in four regions (called frozen regions or polar regions), one adjacent to each corner of the Aztec diamond, whose union is approximately the region outside of the largest circle (called the arctic circle) that can be inscribed in the Aztec diamond.The strategy used in the original proof of this fact was a careful analysis of the shuffling algorithm [15].
There are many ways to view domino shuffling.Originally it was described using sequences of moves that one must perform with the dominos on the tilings themselves, see [14,22].Alternatively, one can see it as an example of dynamics on a certain space of particle configurations, as first observed by Nordenstam [20].With a restricted class of weights, these particle dynamics can be re-derived using the algebraic structure of the Schur process [5,6,2,3].Furthermore, when viewed as a deterministic discrete time dynamical system on the weights, domino shuffling is an example of a cluster integrable system [13].This is a consequence of the fact that the shuffling algorithm consists of a collection of structure preserving local moves called spider moves.
In this article, we describe a generalization of the shuffling algorithm.Recently, in [7] the authors described a model of k interacting tilings of the Aztec diamond.The authors computed the partition function of the model by relating it the LLT polynomials of Lascoux, Leclerc, and Thibon [18].We generalize the shuffling algorithm to these interacting k-tilings.
More precisely, a k-tiling of the Aztec diamond is a collection of k domino tilings of the Aztec diamond.We consider the tilings to be indexed by colors, which are ordered.Thinking of the different tilings as being overlaid one on top of the other, we define an interaction between two of the tilings as an instance of a local configuration of the form , , , or where above blue is the smaller color in our ordering and red the larger.We assign a weight to the k-tilings given by t # interactions .Even restricting to just k = 2, this distribution on 2-tilings is an integrable one parameter deformation of the double dimer model, which is obtained by setting the interaction strength t = 1.In the large N limit, the model appears to exhibit many of the same phenomena as the dimer model, including arctic curves and limit shapes.On the other hand, the properties of these limit shapes and arctic curves appear to be very different from those observed in the dimer model.See Appendix B for a brief discussion and several simulations.
In addition to questions about limit shapes and arctic curves, there are many other natural questions one can ask about the coupling between the 2 (or k) interacting tilings in the scaling limit.For example, near the arctic curve we expect to observe a one parameter deformation of k independent Airy processes, in which each color's edge fluctuations are coupled in a nontrivial way.It would also be interesting to study the global height fluctuations, and in particular how the fluctuations of different colors are coupled together.We expect that the efficient sampling algorithm provided by our main theorem below could be an important tool for the investigation of these questions.
The following is a special case of the main result.
Theorem 1.2.The following algorithm generates a random k-tiling of the rank-N Aztec diamond with probability proportional to its weight.Algorithm: Start with a rank-0 Aztec diamond.To get from a rank-(T − 1) to rank-T k-tiling, 1. Slide and destroy as in the normal domino shuffle, independently for each color.
2. Fill in empty 2 × 2 squares according to the rule: (a) For the smallest color put two horizontal dominoes with probability t #1(1) 1 + t #1 (1)   where # 1 (l) = # colors m > l that locally have or or creation.Do all of these first.
This gives a k-tiling of the Aztec diamond whose rank has increased by one.
Repeat steps (2) and (3) until you get a rank-N Aztec diamond.
The main tools we use are the Cauchy and branching identities for the LLT polynomials, as these allow us to apply a general construction of Borodin and Ferrari [2].In fact, a standard bijection allows one to view a k-tiling as a tuple of interlaced particle arrays.The aforementioned construction (after some calculation) prescribes explicit transition probabilities for these particles, such that if the initial particle positions correspond to a random rank-N tiling, then the update generates a random k-tiling of rank-(N + 1).Using the bijection to interpret the dynamics as local moves on dominos, we obtain the shuffling algorithm in Theorem 4.1 in the text, which gives Theorem 1.2 by setting the weights to be uniform.This Markov chain on colored particle arrays generalizes the Markov chain on a single particle array described in [3, Section 2], which corresponds to the usual shuffling algorithm.
In addition to the proof of Theorem 1.2 using LLT polynomials, in Appendix A we give an alternative proof which employs a local resampling procedure which generalizes the resampling coming from the spider move.The resampling relies on a set of relationships between local partition functions, which are listed in Lemma A.1.Contrary to the one color case, these relations are not sufficient to produce a shuffling algorithm for k-tilings with arbitrary weights.However, they can still be used to construct a shuffling algorithm for certain choices of weights, including uniform weights, which is the setting of Theorem 1.2.
The paper is organized as follows: 1.In Section 2, we give a brief review of background material.We begin by reviewing the Aztec diamond and stating some fundamental results.We focus on highlighting the relationship with interlacing partitions, interlacing arrays of particles, and Schur polynomials which will be useful in later sections.We then define the k-tilings.We state some fundamental results from [7].We also state the necessary properties of the LLT polynomials that will be used in the subsequent sections.Of particular importance to what follows is the bijection between tilings and particle configurations.
2. In Section 3, we present the Markov chain on colored interlacing particle configurations which preserves a class of probability measures on colored particle arrays called 'LLT processes'.First, we review the one color case, which is the Schur case, and then we describe the generalization to multiple colors, which is powered by LLT polynomials.
3. In Section 4, we interpret the particle dynamics described in the previous Section as an operation on dominos.We review the shuffling algorithm for a single tiling of the Aztec diamond before describing the corresponding result for the k-tiling.The main result is an algorithm for generating random k-tilings with probability proportional to their weight.

4.
In Section 5, we summarize our results and give some possible avenues of future research.

5.
In Appendix A, we give an alternate description of our shuffling algorithm in terms of a generalization of the 'spider move' on the underlying double dimer model.
6. Finally, in Appendix B, we present some simulations of the k-tilings generated using our shuffling algorithm.As noted above, the coupled tilings appear to exhibit limit shapes and arctic curves.We give a discussion of the apparent features.
Acknowledgements.The authors would like to thank Alexei Borodin, Sylvie Corteel, and Ananth Sridhar for many useful discussions.

Background
2.1 Tilings of Aztec diamond

The Aztec diamond
Let A N +1 be the union of faces of Z 2 which are entirely contained in the region |x| + |y| ≤ N + 1.A tiling of the Aztec diamond of rank N is a tiling of the region A N +1 with 2 × 1 or 1 × 2 dominos.See Figure 1 for an example.Label the faces by (i, j) ∈ (Z + 1 2 ) 2 , and label the diagonals from 0 to 2N by declaring that the face (i, j) is on diagonal j − i + N .
We assign a (position dependent) weight to each domino in a tiling.Let C N = (c 1 , . . ., c N ), B N = (b 1 , . . ., b N ) be two tuples of real numbers.The domino weights we are interested in are given by: • Suppose a horizontal domino D is occupying the two squares (i, j), (i+1, j) with (i, j) on diagonal 2m − 1.Then the weight of this domino is c m .
• Suppose a horizontal domino D in a tiling is occupying the two squares (i, j), (i+1, j) with (i, j) on diagonal 2m.Then the weight of D is b N −m+1 .
• The weights of vertical dominos are 1.
The weight of a whole tiling tiling T is given by the product of the weights of each domino, wt(T ) = dominos D∈T wt(D) .
Define the rank-N Aztec diamond partition function with these weights as The probability of a random rank-N tiling is given by wt(T ) 6]).The partition function of the Aztec diamond of rank N with the above weights is given by One of the ways this theorem can be proved is via the machinery of Schur polynomials, the basics of which we briefly review in the next subsection.

Schur polynomials and interlacing partitions
A partition λ = (λ 1 , λ 2 , λ 3 , . ..) is a non-negative sequence of integers such that We draw our diagrams in French notation, in the first quadrant, as shown in the below example: We refer to the elements in D(λ) as cells.The cell labelled above has coordinates (1,3).The content of a cell u = (i, j) in the i-th row and j-th column of the Young diagram is c(u The size of a partition is the number of cells in its Young diagram and is denoted by |λ|.The above partition has size |λ| = 4 + 2 + 1 = 7. See Figure 2 for another example.
Given two partitions, λ and µ, we say that µ is contained in λ and write µ ⊂ λ if the Young diagram of µ is contained within λ.Given that µ ⊂ λ we can define the skew diagram λ/µ as the Young diagram of λ with the cells from the Young diagram of µ removed.
To a partition we can associate an infinite sequence of particles and holes by assigning a particle to every vertical edge on the boundary of its Young diagram, and a hole to every horizontal edge.This is known as the Maya diagram of the partition.See Figure 2. Note that the Maya diagram has a unique content line such that the number of particles to the right of this line is equal to the number of holes to the left.We call this the zero-content line and view it as the center of our Maya diagram.If we place the Maya diagram on Z + 1 2 centered at zero then the position x i of the i-th particles (counting from right to left) is given by For every partition λ we can associate a second partition λ known as the conjugate of λ.The conjugate λ is defined as the partition whose Young diagram is given by reflecting the Young diagram of λ across its zero-content line.For example, λ = (4, 3, 2, 2, 1), the partition in Fig. 2, has conjugate λ = (5, 4, 2, 1).
Given two partitions λ and µ say that λ and µ interlace if and write λ µ.Say that λ and µ co-interlace if their conjugate partitions interlace and write λ µ.Note that µ λ implies that µ ⊂ λ.Recall that given two partitions µ ⊂ λ and variables X n = (x 1 , . . ., x n ) the skew-Schur polynomial is defined as Here the sum is over all semi-standard Young tableaux of shape λ/µ, where a semi-standard Young tableaux is a filling of the cells of the diagram by the integers 1, . . ., n such that they are weakly increase along the rows and strictly increasing up the columns.

The Aztec diamond and Schur processes
There is a bijection between tilings of the Aztec diamond of rank N and sequences of interlacing partitions ∅ λ (1)  µ (2) . . .λ Given a tiling of the Aztec diamond of rank N assign particles and holes to the dominos according to the rules , , , .
Along each diagonal slice of the Aztec diamond view the resulting sequence of particles and holes as the Maya diagram of some partition by extending it infinitely to the South-West with particles and infinitely to the North-East by holes.Let us index the slices starting from 0, such that µ (i) is the partition along slice 2i−2 and λ (i) is the partition along slice 2i−1.Note that µ (1) = µ (N +1) = ∅ is forced.Figure 3 gives an example of our notation and this bijection.One can check [6,16] that the requirement that these partitions come from a valid tiling is exactly the interlacing condition given in Eqn.(1).
Given the bijection from tilings to sequences of interlacing partitions λ (1) , µ (2) , λ (2) , . . ., µ (N ) , λ (N )   described above, one can write the weight of the tiling in terms of Schur polynomials.The weight of a tiling is given by where here we used the notation (λ/µ) := λ /µ .Remark 2.3.A probability measure on sequences of partitions of the form with u 1 , . . ., u N , v 1 , . . ., v N ∈ R >0 is a particular case of a Schur Process.Schur processes are well-studied.Using them one can derive exact determinantal formulas for correlation functions and study various statistics of random tilings asymptotically as N → ∞, see [4] for a survey and see also [21].
By repeated applications of the identities in Proposition 2.2, the above simplifies to the product in Theorem 2.1.
It will often be more convenient to consider only the particle positions.From this point of view, the tiling becomes an array of interlacing particles.Let n } be the position of particles corresponding to λ (n) , and n−1 } those corresponding to µ (n) .A set of particle positions corresponds to a tiling if and only if they satisfy the interlacing conditions x and the bounds for each n = 1, . . ., N .See Figure 4 for several examples.

k-tilings of the Aztec Diamond
In this section we define the interacting k-tilings of the Aztec Diamond.See [7] for a more detailed discussion.
Consider an Aztec diamond of rank N .A k-tiling T = (T (1) , . . ., T (k) ) is a k-tuple of tilings of the Aztec diamond.We will often draw the different tilings in different colors (see Figure 4) and refer to tiling T (a) as being color a.We order the colors so that color a is smaller than color b if a < b. If , each of the tilings T (a) has its own weight wt(T (a) ) by giving the horizontal dominos weights c m , b N −m+1 on diagonals 2m − 1 and 2m, respectively, as described in subsection 2.1.1.Now we define an interaction between pairs of tilings.Consider two tilings T (a) and T (b) with a < b.Let blue be the smaller color and red the larger color.We define an interaction between the two tilings to be any instance of the local configuration , , , or when the two tilings are superimposed on top of one another.
Remark 2.5.A k-tiling with these interactions is called the "white-pink" model in [7].
Now we can define the weights we use for k-tilings.
Definition 2.6.The weight of a k-tiling is As usual we will study the probability measure on k-tilings where the probability of each k-tiling is proportional to its weight: .

The partition function Z (k)
AD (C N , B N ; t) can be written in a simple product form.

Theorem 2.7 ([7]
).The partition function Note that when t = 1 the tilings are independent and we have that is, the partition function is a product of k copies of the partition function for a single tiling of the Aztec diamond, as one would expect.More surprisingly, when t = 0 we have that is, the partition function of the k-tiling is equal to the partition function of a single tiling.See [7] for a bijective proof of this fact.

LLT polynomials
While Schur polynomials were valuable in studying the single tilings of the Aztec diamond, for the k-tilings it is the LLT polynomials that are useful.
LLT polynomials are certain symmetric polynomials introduced by Lascoux, Lecler, and Thibon [18] as the generating functions of semi-standard ribbon tableaux counting a spin statistic.Recently, a version of these polynomials called coinversion LLT polynomials were constructed as the partition function of a class of integrable lattice models [1,8,9,12].In [7], the authors used the lattice model formulation of the coinversion LLT polynomials to compute the partition function of the interacting k-tilings of the Aztec diamond that are the focus of this paper.
Here we collect the relevant definitions and properties of the coinversion LLT polynomials.

x ≤ y ≤ z.
While y must be in the Young diagram of λ (a) /µ (a) , we let x = 0 and z = ∞ if they are not in the Young diagram of λ (b) /µ (b) .
Define the coinversion LLT polynomial by where the sum is over k-tuples of semi-standard Young tableaux with shapes given by λ/µ and coinv(σ) is the number of coinversion triples in the filling.We will also need a 'dual' version of the polynomials.Define L by where the sum is over k-tuples of tableaux that are weakly increasing up columns and strictly increasing across rows, and inv(σ) counts the number of triples of the form x z y where 1. y is an entry of σ (a) and x, z are entries in σ (b) , with 1 ≤ a < b ≤ k.
2. y and z lie along the same content line.

x < y < z.
As for coinversions, y must be in the diagram but we let x = 0 and z = ∞ if they are not in the diagrams.Example 2.8 gives an example of both a coinversion and inversion triple.
Remark 2.9.Note that here our notation differs slightly from that of [7].Our L are the same as their L P .Also, our d(λ, µ) in the Cauchy identity, Prop.
The polynomials L and L satisfy the following properties: 1.They are symmetric in the x i .
2. When k = 1 we have 3. When t = 1 we have In addition, we have the following propositions.
Proposition 2.10 (Branching rule, [12]).The L and L satisfy the branching rules Proposition 2.11 (Cauchy identity, [12]).The L and L satisfy the Cauchy identity where d(λ, ν) has an explicit formula in terms of the parts of the partitions.If ν = µ = 0, then this simplifies to See [12] for proofs of these properties via integrable vertex models.
Remark 2.12.Note that when k = 1 this reduces to the dual Cauchy identity for Schur polynomials.
We will also need the following definitions: The size of a k-tuple of partitions is denoted |λ| and is given by the sum of the sizes of each partition in the tuple.We say that two k-tuples of partitions λ and µ (co-)interlace if for each i = 1, . . ., k we have that λ (i) and µ (i) (co-)interlace.We write λ µ or λ µ, for interlacing and co-interlacing, respectively.
See Figure 4 for an example.

Proposition 2.13 ([7]
).In terms of the corresponding Maya diagrams λ (1) , µ (2) , . . ., µ (N ) , λ (N ) , the weight of a k-tiling of the Aztec diamond of rank N can be written as Note that when have only a single variable, the LLT polynomials can be written In particular, they are monomials.It will be useful in the proof of Prop.3.9 to know precisely how each partition contributes to the powers of t in each of these monomials.
Lemma 2.14.Fix two colors b > a and consider the i-th row of color a in λ.
1.There is an interaction between this row and row j of color b whenever min(λ to the total power of t in L λ/µ (x; t).

2.
There is an interaction between this row and row j of color b whenever in which case its contributes a single power of t to the total power of t in t d(λ,µ) Lλ/µ (y; t).
Proof.This follows from the discussion in Section 4.4 of [7].
As with the single tiling of the Aztec diamond, we can associate to each k-tiling a particle configuration.In this case, we have particles of k colors, one for each tiling.For each color, the particles must satisfy the required interlacing conditions (2) and (3).Again, see Figure 4 for an example.

Markov chains on colored interlacing arrays
In this section we define a Markov chain on colored interlacing particle arrays which will be equivalent to the shuffling algorithm under the identification between particle arrays and k-tilings described in Section 2. The essential ingredients are the Cauchy identities and branching rule for LLT polynomials.Using these, we apply a construction whose original form was introduced in [10], and which was further developed in the case of random tilings in [2,3].We first elaborate on the (well-known) construction in the one color Schur case, and then describe the corresponding generalization to the LLT case.

Markov Chains on Schur Processes
Recall that under the bijection to interlacing partitions, the probability measure on domino tilings becomes As mentioned in Remark 2.3, this is an example of a Schur process.Now we introduce transition kernels which are building blocks for Markov chains that map Schur processes to other Schur processes.For generalizations of this construction and more details, see [4,3] and references therein.Suppose that b, c 1 , c 2 , . . .are positive real numbers, and denote the tuple c [1,k] = (c 1 , . . ., c k ).Define transition probabilities by where which can be computed by the Cauchy identity.We have where the second equality follows from the branching rule.Note that p ↑ λ→µ is nonzero only if λ ⊂ µ and we view this as the partition growing from λ to µ.Similarly, p ↓ λ→µ is nonzero only if µ ⊂ λ and we view this as the partition shrinking.
A key property of the transition kernels p ↑ , p ↓ is their commutation: This property also follows from the skew-Cauchy and branching identities.Using these, we will define transition probabilities out of which our Markov chain will be built.Now, let c 1 , . . ., c N , b 1 , . . ., b N be positive real numbers.Define c [i,j] := (c i , . . ., c j ).Given partitions λ and ν, define the transition probabilities to a new partition λ by .
2. For each k = 1, . . ., N + 1, update λ (k) → λ(k) with transition probabilities Now we note that we may write the Schur process corresponding to a rank-N domino tiling of the Aztec diamond as This form is useful in the following proposition.
Then after the Schur parallel update, the updated partitions are distributed according to the Schur process This proposition is a special case of well-known facts (again see e.g.[4,3]), but we include a proof for clarity of exposition, and in particular because this proof generalizes immediately to the multi-colored case.
Proof.Using Eqn. ( 14), we must show Indeed, first note that since μ(i) = λ (i−1) deterministically, we are not summing over the λ (i) .In particular, the numerators of the transition probability and thus can be pulled out of the sum.Moreover, this product is exactly the Schur process we desire.We are left to show that the denominator of the transition probabilities cancels with what remains in the numerator.Consider the sum over µ (2) .The relevant term is Using the commutation relation (13) we have This is precisely the denominator of P2 and we see that the terms cancel.
More generally, for k = 2, . . ., N + 1, one can see that the sum over µ (k) cancels with the denominator of Pk through the commutation relation Computing the transition probabilities in terms of particle positions respectively, one observes that we can define the Schur parallel update in terms of the particles as follows: 1.For each n = 2, . . ., N + 1, set ỹ(n) = x (n−1) .
These transitions are independent for 1 ≤ i ≤ n ≤ N + 1.
Remark 3.4.That the Markov chain has this form follows from the discussion of the Aztec diamond shuffling algorithm in [3].It will also follow from our generalization to multiple colors in the next section.

Markov Chains on LLT Processes
Now we generalize the Markov chain above to a Markov chain on k-tuples of interlacing partitions.While the definition in terms of interlacing partitions will follow directly from machinery of LLT polynomials, the interpretation as particle dynamics will require careful computation.Let c 1 , . . ., c N , b 1 , . . ., b N ∈ R >0 .Define the LLT process to be a probability measure on arrays of tuples of partitions given by P (λ (1) , µ (2) , λ (2) , . . ., Recall that this probability measure describes random k-tilings (Proposition 2.13).
The transition kernels from which we will build the Markov chain are, for c = (c 1 , . . ., c l ) and b, c ∈ R, These transition kernels satisfy a similar commutation relation to Equation (13): In the exact same way as in the Schur case, we can write the LLT process as P (λ (1) , µ (2) , λ (2) , . . ., We define the following update step that we will use for transitioning from rank N to rank (N + 1).Definition 3.7 (LLT parallel update).Suppose we are given a k-tuple of sequences of interlaced partitions 0 λ (1)  µ (2) . . .
Then the dynamics on the multi-colored interlacing particle arrays defined below is equivalent to the LLT parallel update defined in Definition 3.7 above.
2. Given x (n) and ỹ(n) , x(n) can be sampled according to the following rules: For each l = 1, . . ., k with probability 1.
x (n,l) i otherwise where These transitions are independent for each In (18), we allow j = 1, . . ., n, and we use the convention that ỹ(n,m) Proof.First let us recall that for the Markov chain on LLT processes, at level n the transition probabilities are proportional to the product of the transition kernels To simplify notation, in what follows we let c = c n , b = b N −n+2 , and λ = λ(n) , λ = λ (n) , μ = μ(n) .We see that the probability of a particular λ under Pn ( λ|λ, μ) is proportional to where we have only kept factors depending on λ.It is not hard to see from the definition of LLT polynomials that the quantity Lλ /λ (b; t) will be 0 unless λ corresponds to a configuration where all particles jump by at most 1.Furthermore, if λ, μ correspond to a particle configuration where there are particles forced to jump or stay, then L λ/ μ(c; t) will be 0 unless λ corresponds to a configuration where all of these particles do in fact jump or stay.Therefore, the possible k-tuples λ correspond exactly to the possible outcomes of particle jumps in the proposition.
Thus, to prove the proposition, it suffices to show that for these k-tuples λ, the ratios of their transition probabilities are equal to the ratios of the particle transition probabilities described in the proposition.It is enough to compare the ratio for pairs of x which differ by a single non-forced particle jump, as any ratio of the particle transition probabilities can be written as a product of such simple ratios.For concreteness, suppose the jump was made by particle i of color l.
Two particle configurations differing by a single particle jump are equivalent to two tuples of partitions differing by the corresponding single cell in one of their Young diagrams.Define δ(l, i) = (δ(l, i) (1) , . . ., δ(l, i) (k) ), where From the above discussion, we see that we must show that where # i (n, l) is defined in (18), and λ(l To show that the powers of c and b are correct in the equality (19), recall Eqn.(10): for LLT polynomials with a single variable we have Since the Young diagram of ( λ + δ(l, i))/λ has exactly one more cell than that of λ/λ, the numerator will have exactly one extra factor of b, and similarly one extra factor of c.
We are left to show the powers of t match on both sides of (19).The powers of t on the LHS of ( 19) come from two sources: interactions between color l and colors m > l, and interactions between color l and colors m < l.We show in detail that the contributions when the color m is larger give exactly the first term in (18).A similar analysis can be done to show that the contributions when the color m is smaller give the second term in (18).
Looking at the first term of ( 18), note the particle position inequalities ỹ(n,m) (n,m) j can be written in terms of the parts of the partitions as μ(m) We will show that in the LHS of (19) we get an extra power of t in the numerator that is not present in the denominator exactly when there exists an m and j where the above inequality holds.
To do so we do an exhaustive check over all possible relative positions of the corresponding particles and in each case use Lemma 2.14 to determine the powers of t on the LHS of ( 19).This casework is listed below: i − i there are three subcases for the contribution from Eqn. ( 11): − j then there is a power of t 0 in the numerator and no contribution to the denominator.
− j there is no contribution to either the numerator or the denominator.
For the contribution of Eqn.(12), note there is no contribution to the denominator since λ(m) In the numerator there is no contribution since λ (l) Overall, we get a net power of t in the ratio in case (a) and none in the other cases.

If λ(m)
i +i in the denominator and a t i +i in the numerator.Eqn.(12) does not contribute to either for the same reason as in case 1.
Overall, we have a net power of t in the ratio.

If λ(m)
in both the numerator and the denominator.
Eqn. (12) does not contribute to either, again, for the same reason as in case 1.
Overall, there is no net power of t in the ratio.

If λ(m)
i − i then for the contribution from Eqn. (11) we have two subcases: i +i+1 in both the numerator and the denominator.(b) Otherwise, there is no contribution to either.
Eqn. (12) still does not contribute to either for the same reasons as above.
Overall there is no net power of t in the ratio.Overall, there is no net power of t in the ratio if λ(m) , and a net a power of t otherwise.

If λ(m)
Then for the contribution from Eqn. (11) we get a t λ(m) Eqn. (12) does not contribute to the numerator for the same reason as in case 5.The contribution to the denominator from Eqn. ( 12) can be split into two subcases:  x (2)   ỹ(3) Figure 5: Shown is a possible configuration of particles of the Maya diagrams on diagonals 2, 3 and 4, after step one of the update.We use the convention blue = 1 < red = 2 < green = 3.The solid red arrow denotes a forced jump for the red particle to preserve the interlacing of the red particles.Supposing we have c i = b i ≡ 1, the dashed red arrow denotes a jump which will happen with probability t 2 1+t 2 , because ỹ(2,1) . Similarly, the dashed green arrow denotes a jump which has probability t 1+t . ( − j + 1, it contributes a single power of t to the denominator.(b) Otherwise, it contributes nothing to the denominator.
Overall, there is no net power of t in the ratio if λ(m) − j + 1, and a net a power of t otherwise.
One can check the ratio of the powers of t is given by exactly as we desired.Summing over all m > l and all j gives the first term in Eqn.(18).
As a simple example of the particle dynamics, one may set all parameters b i ,c i equal to 1 and consider x (1) , that is, the bottom level of particles.Along this row we have exactly one particle of each color, and the marginal evolution of these particles is itself a Markov chain.At each step each particle independently either stays in place or jumps by 1 to the right.The probability for the particle of color l, at position x    In other words, if at time T particles are ordered 0, 1, . . ., n − 1 from right to left, breaking ties by putting larger colors first, the jump probability of particle i is t i 1+t i .When t < 1, we see that if a particle falls behind the other it becomes discouraged and moves more slowly, while for t > 1 it becomes determined to catch up and moves more quickly.

Shuffling
In this section we show how the particle dynamics defined in Section 3 can be described as an algorithm which acts directly on the tilings by sliding, destroying, and creating dominos.In the one color case, this is known as the domino shuffling algorithm [11,22].

Interpretation in terms of Dominos: Schur Case
Domino shuffling is a sampling algorithm to generate a random rank-N domino tiling.Via the bijection between tilings and interlacing particle arrays, the particle dynamics defined in Section 3.1 coincides with the shuffling algorithm.Here we will review the shuffling algorithm, and refer the reader to Propp [22] for more details.
Recall that we assign weights c i to the horizontal dominos on the diagonal slice 2i−1, and b N −i+1 to the horizontal dominos on slice 2i, respectively, for i = 1, . . ., N .Take two additional numbers c N +1 , b N +1 .We now describe a way to randomly sample a tiling of rank N + 1 given one of rank N , such that if the original one is sampled with weights (c 1 , . . ., c N ), (b N , . . ., b 1 ), the one obtained from the algorithm will be sampled with weights (c 1 , . . ., c N +1 ), (b N +1 , . . ., b 1 ).
Given the checkerboard coloring of the squares of rank-N Aztec diamond, recall that there are four types of dominos which can appear in a tiling.They are shown in Figure 1.Label these four types of dominos as S (South), W (West), N (North), E (East), respectively.Given a tiling T of rank N , following three steps will lead to a tiling T of rank N + 1: 1. (Sliding) S dominos slide one unit South, W dominos slide one unit West, N dominos slide one unit North, and E dominos slide one unit East 2. (Destruction) If two dominos cross each other's path as they slide, both are destroyed and deleted from the tiling.
3. (Creation) What remains will be a partial tiling of rank N + 1.The untiled portion will be a disjoint union of 2 × 2 blocks (in a unique way).Fill in each 2×2 block independently with either a vertical pair or a horizontal pair of dominos, with probabilities vertical: where here 2i − 1 is the diagonal of the lower left square in the block, with respect to the indexing of diagonals on the rank-(N + 1) Aztec diamond.
The correspondence of the shuffling algorithm to the particle process is well known [3].The following facts, from which the equivalence of shuffling and the particle transition probabilities can be deduced, will be useful for us in the next section: • Slides west correspond to a particle being forced to stay.
• Slides north correspond to a particle being forced to jump.
• Creations correspond to a particle which can either stay or jump.
For an example, see the red tiling and red particles in Figure 6.

Interpretation in terms of Dominos: LLT Case
We now state the analogous shuffling algorithm corresponding to the LLT particle process.
Theorem 4.1.The following algorithm generates a random k-tiling of the rank-N Aztec diamond with probability proportional to its weight.
Algorithm: Start with a rank-0 Aztec diamond.To get from a rank-(T − 1) to rank-T k-tiling, 1. Slide and destroy as in the normal domino shuffle, independently for each color.
2. Fill in empty 2 × 2 squares according to the rule: (a) For the smallest color put a horizontal pair of dominoes with probability where where here 2i is the diagonal of the lower left square in the block, otherwise put a vertical pair.Do all of these first.
(b) Now do all the larger colors from smallest to largest.For color l > 1, put a horizontal pair of dominoes with probability where # 2 (l) = # colors m < l that locally have or and # 1 (l) is as in part (a), otherwise put a vertical pair.
This gives a k-tiling of the Aztec diamond whose rank has increased by one.
Repeat steps (2) and (3) until you get a rank-N Aztec diamond.
Proof.It suffices to show that the update step corresponds to the update of tuples of interlacing arrays of particles described in Proposition 3.9.For each color, the transition probabilities only differ from those of the usual domino shuffling algorithm through the creation probabilities.As in the one color domino shuffling, creations correspond to particles which can jump or stay, with creation of a horizontal pair corresponding to a particle jumping by 1 and creation of a vertical pair corresponding to staying.It is enough to show that the power of t in the creation step described above corresponds to the power of t in the jump probability of the particle process.Consider a color l, and a larger color m, which we represent by blue and red, respectively.Suppose that we are doing a creation for color l and in that 2 × 2 box the red configuration looks like one of the three configurations given in the definition of # 1 in the theorem statement.Given these domino configurations, we can give the possibilities for the corresponding red and blue particle configurations.These are shown in Figure 7.Note that in each case the particle configuration contributes a power of t to the blue particle's jump probability in (18).On the other hand, it is also clear from Figure 7 that these are all of the possible cases in which red particles contribute a power of t to the blue particle's jump probability in (18).

Case 1:
Case 2: Case 3: Figure 7: Three cases in which the dominos of a larger color m (in red) contribute a factor of t to the creation probabilities of a smaller color l (in blue).In each case we also show the relative positions of the blue particle, which is jumping, with the nearby red particles whose positions lead to the contribution of the power of t.We use the same notation as in Proposition 3.9; in particular, x (n,l) i and x (n,m) j denote the blue and red particles in the x (n) diagonal before the update.By the update rules for y particles (c.f.Proposition 3.9), these agree with the particle positions ỹ(n+1,l) i and ỹ(n+1,m) j in the y (n+1) diagonal after the update, which are shown in the figure.Note that the figures are still valid if j = n; in this case the bottom most red particle would be off of the Aztec diamond.The dashed blue circle indicates where the blue particle would be if it does not jump, the solid blue circle indicates where the particle would be if it jumped.Now consider a color l, and a smaller color m, which we represent by red and blue, respectively (so again blue is the smaller color, but now red is the color whose transition probability we are discussing).Suppose that we are doing a creation for the red dominos in a 2 × 2 box, and that the blue dominos there look like one of the two configurations in the definition of # 2 in the theorem statement.Then the possibilities for the corresponding red and blue particle configurations are shown in Figure 8. Again we see that one of these cases occurs if and only if the red particle's transition probability obtains a power of t from blue in (18).

Case 1:
Case 2: Case 3: Case 4: (n,m) j to denote particle positions of the x (n) diagonal at the previous time step, which agree with the particle positions ỹ(n+1) of the y (n+1) diagonal at the next time step.The dashed red circle indicates where the red particle would be if it does not jump, the solid red circle indicates where the particle would be if it jumped.
As a result, the powers of t that we pick up from # 1 (l) + # 2 (l) are exactly the powers of t described in Proposition 3.9.

Conclusion
In this article we generalize the domino shuffling algorithm for tilings of the Aztec diamond to shuffling algorithm for interacting k-tilings studied in [7].We present this algorithm in a variety of ways: 1.As a Markov chain on LLT processes, which generalize the well-studied Schur processes.
2. As a sequence of local moves on the dominos themselves.
3. And through a generalization of the 'spider move' on the underlying dimer models (see Appendix A).
It is intriguing that despite the increased complexity of the k-tilings compared to standard tilings of the Aztec diamond, the shuffling algorithm remains very familiar.In terms of moves on the dominos, the increase in complexity is confined only to how new pairs of dominos are created.This has interesting combinitorial consequences.In fact, it can be used to give a alternate proof that the partition function of the k-tilings of rank N is given by Z Using this shuffling algorithm allows fast sampling of the k-tilings.In simulations one finds that these k-tilings appear to show very interesting asymptotic features, including the presence of arctic curves and limit shapes.See Appendix B. One might hope that, just as for the standard tilings of the Aztec diamond, the shuffling algorithm may in the future be useful in the study of this asymptotic behavior.In what follows, when we show a small patch of a graph, dashed lines mean the patch is a portion of a larger graph.The labels indicate the weight assigned to a dimer that occupies the corresponding edge.When it is clear from the context, we sometimes also use the term 'cell' to refer to both the patch of graph and the dimers occupying edges of that face in a dimer configuration.As we are considering double dimer configurations, we will color the dimers red and blue to distinguish the two configurations.
Figure 10: Starting with a double dimer configuration on the Aztec diamond of rank 1, the shuffling algorithm produces one on an Aztec diamond of rank 2. The weight of the configuration before the shuffling can be computed as the product over cells of cell weights in the top right picture, using the local interactions listed in (20).This gives a weight of t.The weight of the configuration after the shuffling can be computed as the product of the local interactions listed in (21) over the cells in the bottom left picture; this also gives a total weight of t.One can check that for both the N = 1 and N = 2 domino tiling corresponding to the top left and bottom right configurations, the weight is indeed t, c.f. the domino interactions in Section 2.2.

Define interactions to be local configurations of the form (20)
A dimer occupying the dashed edge adjacent to a vertex v denotes that a dimer occupies one of the other edges adjacent v that are not part of the cell.
Note that if we consider the underlying dimer configuration of a 2-tiling of the Aztec diamond, the product over all cells of t # interactions agrees with the product of the domino interactions defined in Section 2.2.See Figure 10, top, for an example.Now consider applying the spider move to the cell above.This results in a new patch of graph with different edge weights a b c d where unlabelled edges have weight 1.After doing this local transformation to each cell in the Aztec diamond and then contracting all valence two vertices, one gets an Aztec diamond of one larger rank.We depict this process in Figure 9.We will also refer to these patches as cells.
After the spider moves, in the double dimer configuration we now count interactions of the form (21) since the product over cells of these interactions will agree with the domino interactions we get after the contraction of degree two vertices.For an example compare the interactions in the configurations in the bottom row of Figure 10.
Note that for the usual dimer model, if then for each choice of boundary condition for the patch, there is an equality of partition functions on the cell before and after the spider move, up to an overall factor of ∆ = ac + bd.To set some notation we list these equations below.Note that if x is a cell with a configuration of dimers, w(x) denotes the product of weights of occupied edges, and the weights are implicitly updated as described above after the spider move.
We label the boundary condition by the type of move that it corresponds in the shuffling when going from the domain on the LHS to the domain on the RHS.In the equations that follow, we will label a boundary condition for a cell by an arrow to indicate left/right/up/down, a 'c' for creation, or a 'd' for destruction.
For the double dimer model the situation is more complicated.A boundary condition (αβ) for a cell consists of a boundary condition α for the blue configuration and a boundary condition β for the red one.Define C as the set of boundary conditions (αβ) for a cell such that We have the following relation between partition functions before and after the spider move.which is consistent as c = a ac+bd .The following combinatorial lemma, whose proof we omit, will be useful.Suppose that in a double dimer configuration, we have a cell x with boundary condition (αβ).We define the transition probability P(x → x ) from x to a double dimer configuration x in the cell after performing the spider move as P(x → x ) = where Z αβ denotes the partition function in this cell after the spider move with the boundary conditions (αβ), as in Lemma A.1.Now we define the shuffling algorithm of Theorem 1.2 in terms of local moves; we will apply the spider move to each cell of the Aztec diamond and re-sample the double dimer configuration in each cell.In more detail, suppose that T (N +1) is a double dimer configuration of the rank-(N + 1) Aztec diamond.Note that T (N +1) corresponds uniquely to a double dimer configuration on the graph obtained from the rank-N Aztec diamond by decorating the boundary and performing spider moves at each cell (see Figure 9, bottom left).We denote by x (N +1) the local configuration in each cell of this graph.For a double dimer cover T (N ) of the rank-N Aztec diamond, we define P(T (N ) → T (N +1) ) = cells x P(x (N ) → x (N +1) ) where x (N ) is the configuration in cell x in T (N ) and x (N +1) is as described above.If all weights are constant, which is the uniform case, the transition probabilities defined above coincide with those of Theorem 1.2.
Alternate proof of Thm.1.2 for the case of two colors.Consider a random 2tiling T (N ) of the Aztec diamond of rank N .The tiling is made up of gluing together dimer configurations in each cell.The weight of a 2-tiling is given by product over cells of the weight in each cell.Let Z (N ) be the partition function for these 2-tilings.
By applying the spider move to each cell and then contracting all valence two vertices we get a 2-tiling of the rank-(N + 1) Aztec diamond (see Fig. 9).The probability to obtain a specific 2-tiling T (N +1) , whose configuration in cell x before the contraction step is denoted x (N +1) , is given by P(T (N +1) ) =

Figure 1 :
Figure 1: The four possible dominos and an example of a tiling of the Aztec diamond of rank 3.

Figure 3 :
Figure 3: Left and Center: Assigning partitions to slices of the Aztec diamond.The red line indicates the zero content line for the partitions.Right: The tiling corresponding to all partitions being empty for the Aztec diamond of rank 3.
numerator.Here the extra power of t in the numerator comes from the indicator 1( λ(m) j − j < λ(l) i − i + 1) in(11).Note that(12) does not contribute to the numerator as λ(m)j − j = λ(l) i − i + 1.The contribution from Eqn.(12) to the denominator can be split into two subcases:(a) If λ(m) j − j = λ (m) j − j + 1, it contributes a single power of t to the denominator.(b) Otherwise, it contributes nothing to the denominator.

Figure 6 :
Figure 6: Here we illustrate one possibility for the random update of the particle configuration of the k = 3-tiling of size N = 3 shown in Figure 4. First, shown is the initial configuration, augmented with an extra empty tuple of Maya diagrams.Second, the particles are shown after both step 1 and the forced jumps of step 2. Finally, one possible outcome of step 2 is shown.

Figure 8 :
Figure 8: Four cases in which the smaller color m (in blue) contributes a factor of t to the creation probability of a larger color l (in red).Similarly to Figure 7, we use x (n,l) i , x

Figure 9 :
Figure 9: Starting with an Aztec diamond of rank 2 the above local moves construct an Aztec diamond of rank 3. The first transformation is a boundary decoration, the second consists of a collection of spider moves, and the final transformation is the contraction of degree two vertices.The randomization of these local moves leads to the Markov chain on dimer configurations known as the shuffling algorithm.

Lemma A. 1 .=⇒ a 2
Let the weights after the spider move be For any pair of boundary conditions α, β ∈ {c, d, →, ←, ↑, ↓} for each color, denote by Z αβ the partition function of the domain with these boundary conditions before the spider move, and Z αβ that of the domain after the spider move.ThenZ αβ =∆ 2 ΓZ αβ , (αβ) ∈ C (22) Z αβ =∆ 2 Γ −1 Z αβ (αβ) ∈ D (23) Z αβ =∆ 2 Z αβ o.w.(24)where ∆ = ac + bd and Γ = ac+bd act+bd .Proof.As there are only 36 choices of boundary condition, one can check (by hand or via computer) that the required 36 equations are satisfied.As an example, Z d↓ means the partition function for the domain where the smaller color has the "Destruction" boundary condition while the larger color has the "Down" boundary condition.For boundary conditions of type (d, ↓) ct + abd = (ac + bd) 2 (act + bd) ac + bd c

Lemma A. 2 .
For any double dimer configuration on the Aztec diamond of rank N , along each SW-NE diagonal of cells the difference between the number of cells with local boundary condition of type (αβ) ∈ C and those of type (αβ) ∈ D is equal to 1.

T
(N )   P(T(N ) )P(T (N ) → T (N +1) ) =T (N ) w(T (N ) ) Z (N ) P(T (N ) → T (N +1) ) = 1 Z (N ) T (N ) cells x of T (N ) w(x (N ) )P(T (N ) → T (N +1) x consistent with T (N +1)w(x (N ) )P(x (N ) → x (N +1) ) line Γ(x) denotes the value of Γ for the weights a, b, c, d of the cell x, and we have implicitly updated the weights in the last line.In the last equality we use the relations from Lemma A.1.

Figure 13 :
Figure 13: A 2-tiling of rank N = 256 generated by the shuffling algorithm with t = 0.2 and c i = b i = 2 for all i = 1, . . ., N .

Figure 14 :
Figure 14: A 2-tiling of rank N = 256 generated by the shuffling algorithm with t = 0.2 and c i = b i = 0.5 for all i = 1, . . ., N .

Figure 15 :
Figure 15: A 2-tiling of rank N = 256 generated by the shuffling algorithm with t = 5 and c i = b i = 2 for all i = 1, . . ., N .

Figure 16 :
Figure 16: A 2-tiling of rank N = 256 generated by the shuffling algorithm with t = 0 and c i = b i = 5 for all i = 1, . . ., N .

Figure 17 :
Figure 17: A 2-tiling of rank N = 256 generated by the shuffling algorithm with t = 10000 and c i = b i = 1 for all i = 1, . . ., N .

Figure 18 :
Figure 18: A 3-tiling of rank N = 256 generated by the shuffling algorithm with t = 0.2 and c i = b i = 2 for all i = 1, . . ., N .