Elliptically Distributed Lozenge Tilings of a Hexagon

We present a detailed study of a four parameter family of elliptic weights on tilings of a hexagon introduced by Borodin, Gorin and Rains, generalizing some of their results. In the process, we connect the combinatorics of the model with the theory of elliptic special functions. Using canonical coordinates for the hexagon we show how the $n$-point distribution function and transitional probabilities connect to the theory of $BC_n$-symmetric multivariate elliptic special functions and of elliptic difference operators introduced by Rains. In particular, the difference operators intrinsically capture all of the combinatorics. Based on quasi-commutation relations between the elliptic difference operators, we construct certain natural measure-preserving Markov chains on such tilings which we immediately use to obtain an exact sampling algorithm for these elliptic distributions. We present some simulated random samples exhibiting interesting and probably new arctic boundary phenomena. Finally, we show that the particle process associated to such tilings is determinantal with correlation kernel given in terms of the univariate elliptic biorthogonal functions of Spiridonov and Zhedanov.


Introduction
This paper examines work began by Borodin, Gorin and Rains in [BGR10]. In op. cit., the authors examined q-distributed boxed plane partitions from several perspectives, but the q-distributions were obtained as limits of the elliptic distribution briefly appearing in the Appendix. The present paper takes the Appendix of [BGR10] and expands upon it, following the steps in [BGR10] and [BG09]. However, since we are working at the elliptic level (rather than a degeneration as in [BGR10]), new tools are needed to generalize the results of [BGR10]. These tools belong to the area of elliptic special functions, an active area of research in algebra and analysis generalizing, among other things, the Askey and q-Askey schemes of orthogonal polynomials (as described in [KS] for example). Thus, in some complementary sense, while being a generalization of [BGR10], the paper is an application of multivariate tools introduced by Rains in [Rai10] and [Rai06] (the first is more analytic, the second being more algebraic) which build upon univariate elliptic biorthogonal functions found by Spiridonov and Zhedanov a few years earlier in [SZ00]. Work in the area of elliptic special functions started with Frenkel and Turaev's discovery of elliptic (theta) hypergeometric series ( [FT97]) -the authors of op. cit. cite Baxter's work (see his book [Bax82]) as the genesis of the theory.
The history of the problem starts with random uniformly distributed boxed plane partitions. Much is known about these: asymptotics and frozen boundary behavior ( [CKP01], [CLP98], [KO07]); correlation kernel via Hahn orthogonal polynomials (see [Joh05], [BG09], [Gor08]); exact sampling algorithms ( [BG09]). Somewhat central to the subject is the topic of discrete Hahn orthogonal polynomials (which themselves are terminating generalized hypergeometric series). One level up and we arrive at qdistributions on boxed plane partitions ( [BGR10], and [KO07] for the variational problem used to derive the limit shape for the q ±V olume distributions). Almost as much as above is known about these, and central to the subject are certain discrete q-orthogonal polynomials (q-Racah, q-Hahn) from the q-Askey scheme, which themselves are terminating q-hypergeometric series (see [GR04] for a full description, [KS] for a distillation of the formulas).
The present work analyzes the elliptic level (the distribution was introduced in the Appendix of [BGR10], but also independently from a slightly different perspective in [Sch07]). We look at two aspects: exact sampling algorithm and correlation kernel. The third aspect in [BG09] and [BGR10] is obtaining asymptotics of the correlation kernel and through this obtaining the frozen boundary behavior in the large scale limit. While we indeed see a frozen boundary behavior in our case and can characterize it via variational techniques (and we present computer simulations of the results), we cannot yet analyze the asymptotics of elliptic biorthogonal functions (techniques used in previous works -e.g., in [BGR10] fail if we replace orthogonal polynomials by elliptic biorthogonal functions). More direct techniques like solving the variational problem described in [KO07] for the q-Hahn case and in Section 2.4 of [BGR10] seem computationally intractable so far. The reason is the associated complex Burgers equation one has to solve becomes considerably more complicated. Nevertheless, it is a (new) feature of the elliptic model that the frozen boundary can have 3 nodal points (as seen in the computer simulations).
From a different perspective, we try to create a bridge between elliptic special functions discussed in the references given above and combinatorics of tilings of hexagons (equivalently, dimer coverings of the appropriate graph etc. See Section 2 for interpretations). To wit, we give a combinatorial interpretation to several objects appearing in the theory of elliptic special functions: the (t = q case) multivariate elliptic difference operators discovered by Rains ([Rai10]), the ∆-symbols of [Rai06] and the (univariate) elliptic biorthogonal functions of Spiridonov and Zhedanov ([SZ00]).
This paper tries to emulate the organization of [BG09] and [BGR10], but with notation heavily influenced by [Rai06]. It is organized as follows: in the remainder of the introduction, we set up most the important notation and terminology.
We set up the combinatorial and probabilistic aspects in Section 2 (in the Appendix we consider a different way of assigning weights to rhombi that is manifestly more symmetric. We prefer to use the formulas from Section 2 though, at the cost of symmetry breaking since they lead to shorter computations and arguments). Also in this section we study positivity of our a priori complex measure and introduce various coordinate systems used throughout the paper, including the important canonical coordinates (which embed our model in a certain square of an elliptic curve).
In Section 3 we compute relevant distributions and transition probabilities. Sections 2 and 3 are an expansion and in depth analysis of the Appendix in [BGR10].
Section 4 recalls some definitions and properties of elliptic tools introduced by Rains ( [Rai10], [Rai06] -we refer the reader to these works for the proofs we omit) and then connects these with the probability and combinatorics being studied. We show that the constraints of the model are intrinsically captured by the elliptic difference operators under discussion.
Section 5 describes a perfect sampling algorithm for such elliptic distributed boxed plane partitions. It is based on the idea of forming a new measure-preserving Markov chain out of two old quasi-commuting ones (as in [BF10]; see also [DF90]). The algorithm starts from a deterministic parallelogram shape and samples relatively easy distributions to successively transform the parallelogram into a hexagon accordingly distributed (by increasing one side by 1, and decreasing another side by 1; a parallelogram can be seen as a hexagon with two sides of length 0). We use the quasi-commutation relations for the elliptic difference operators of Section 4 to construct this algorithm.
Section 6 deals with the correlation kernel. We start by recalling facts about univariate elliptic biorthogonal functions and show that the time increasing (decreasing) Markov process is indeed determinantal, with correlation kernel given as a determinant of (a matrix composed of) elliptic biorthogonal functions. These replace the orthogonal polynomials discussed above.
In Section 7 we present some computer simulations obtained from the algorithm described in Section 5. Choice of parameters for obtaining the trinodal cases (surfaces where the arctic circle has 3 nodes at 3 vertices of the hexagon) are also explained.
We end with the Appendix, which provides a highly symmetric view of the entire picture For the remainder of the section, we will set the notation that will appear in the rest of the paper.
We define the theta and elliptic Gamma functions ( [Rui97]) as follows: If f (x 1 , ..., x n ) is a function of n variables defined on (C * ) n , we call it BC n -symmetric if it is symmetric (does not change under permutation of the variables) and invariant under x k → 1/x k for all k. We will call it a BC n -symmetric theta function of degree m if in addition, it satisfies the following: f (px 1 , ..., x n ) = ( 1 px 2 1 ) m f (x 1 , ..., x n ).
The prototypical example of a BC n -symmetric theta function of degree 1 is: We now define the following function (which will play an important subsequent role): Note ϕ is BC 2 -antisymmetric (ϕ(z, w) = −ϕ(w, x)) of degree 1. It is a consequence of the addition formula for Riemann theta functions that ϕ(x, y) = ϕ(z, x) ϕ(w, x) − ϕ(z, y) ϕ(w, y) ϕ(w, x)ϕ(w, y) ϕ(z, w) for arbitrary z, w. We observe that the expression in parentheses appearing above is a Vandermonde-like factor in transcendental coordinates X = ϕ(z,x) ϕ(w,x) , Y = ϕ(z,y) ϕ(w,y) , so ϕ(z i , z j ) is an "elliptic analogue" of the (Vandermonde) difference z i − z j . This is indeed the case if one takes the right limit and then a product over i < j. To wit: Notationally, for a function f of n variables, we will use the abbreviation f (...x k ...) to stand for f (x 1 , ..., x n ).
We will make reference to the delta symbols defined in [Rai10] (see also [Rai06] -we are in the case t = q in the notation from both references), which we define here. We fix λ ∈ m n a partition (that is, a partition with at most n parts all bounded by m). Define the partition 2λ 2 by (2λ 2 ) i = 2(λ i/2 ).
Of interest will be the ∆-symbol with six parameters t 0 , t 1 , t 2 , t 3 , u 0 , u 1 satisfying the balancing condition q 2n−2 t 0 t 1 t 2 t 3 u 0 u 1 = q. Because the usual balancing condition has pq on the right hand side (the reader should consult the Appendix of [Rai10] for more on why this is necessary), we multiply u 1 by p (this choice is arbitrary, so a priori some symmetry is broken, but this will not affect our results). We define the discrete elliptic Selberg density as: where l i = n − i + λ i , z i = q li t 0 and the constant is independent of λ and present in the formula to make the ∆-symbol elliptic in all of its arguments (the value of the constant can be computed, but such constants will be immaterial for the rest of the paper). This discrete elliptic Selberg density is the weight function for the discrete elliptic multivariate biorthogonal functions defined in [Rai06]. Notice it can be written more symmetrically in terms of the z i 's and the elliptic Gamma functions as We will denote by E the elliptic (Tate) curve C * / p for some |p| < 1. An elliptic function f (of 1 variable) will just be a function defined on E (that is, f (px) = f (x)).
Throughout the remainder, constants (by which we mean factors independent of the variables usually denoted by x k , y k , z k ) will largely be ignored (and we will write const wherever this appears), but they are there to make measures into probability measures (i.e., normalizing factors) or to make certain functions elliptic (i.e., invariant under p-shifts). Their values can often be recovered, and we comment on how to recover them whenever possible.
Finally, throughout this paper we will freely use two different systems of coordinates for our model (related by a simple affine transformation -see the next section). While this may seem redundant, coordinatizing in two different ways will more aptly reveal different features of the elliptic special functions and difference operators under study. details are set out in Section 2.2. We will find it more convenient to encode the hexagon via the following three numbers: Equivalently, these tilings can be thought of as dimer matchings on the dual honeycomb lattice (every rhombus in a tiling is a line matching two vertices in the dual lattice), stepped surfaces, boxed plane partitions (b × c rectangles with positive integers ≤ a filled in that decrease weakly along rows and columns starting from the top left corner box) or 3D Young diagrams (any section parallel to one of the three bounding walls is a Young diagram).
A yet different way of viewing such tilings, important hereinafter, is as collections of non-intersecting paths in the square lattice. The paths start at N consecutive points on the vertical axis (counting from the origin upwards) and end at N consecutive points on the vertical line with coordinate T . Each path is composed of horizontal segments or diagonal (Southwest to Northeast, slope 1) segments, and the paths are required not to intersect. Figure 2 explains this, and also introduces the coordinate frame (t, x) that will be used for computational convenience in various sections to follow: (i, j) = (t, x − t/2).
Following the notation in [BG09], let Ω(N, S, T ) denote the set of N non-intersecting paths in the lattice N 2 starting from positions (0, 0), ..., (0, N − 1) and ending at positions (T, S), ..., (T, S + N − 1). Each path has segments of slope 0 or 1 (paths go either horizontally or diagonally upwards from left to right). Set T is the set of all possible particle positions in a section vertical section of our hexagon with horizontal coordinate t (in (t, x) coordinates). X S,t N,T is the set of all possible N -tuples of particles in the same vertical section.
For X ∈ Ω(N, S, T ), we have X = (X(t)) 0≤t≤T and each X(t) ∈ X S,t N,T . X is a discrete time Markov chain as it will be shown.

Probabilistic model
We will now define the probability measure on Ω(N, S, T ) that will be the object of study. For a tiling T corresponding to an X ∈ Ω(N, S, T ) we define its weight to be: where by a horizontal lozenge we mean a lozenge whose diagonals are parallel to the i and j axes. The probability of such a tiling would then simply be: The weight function w on horizontal lozenges l is defined by w(l) = (u 1 u 2 ) 1/2 q j−1/2 θ p (q 2j−1 u 1 u 2 ) θ p (q j−3i/2−1 u 1 , q j−3i/2 u 1 , q j+3i/2−1 u 2 , q j+3i/2 u 2 ) where (i, j) is the coordinate of the top vertex of the horizontal lozenge l, u 1 , u 2 , q, p are complex parameters (|p| < 1) and u 1 = q −S v 1 , u 2 = v 2 (the reason for this break in symmetry is that it will make other formulas throughout the paper more symmetric).
Remark 2.1. Only considering weights of horizontal lozenges for a tiling of a hexagon is equivalent to considering all types of lozenges but assigning the other two types weight 1 (i.e., each lozenge that is not horizontal has weight 1). This is a break in symmetry that can easily be fixed. However, for the remainder of the paper we prefer this non-symmetric weight assignment system as it makes computations easier. Nevertheless, we show in Appendix 8 that we can assign weights to the 3 types of lozenges in a S 3 -invariant way (i.e., invariant under permuting the 3 types of lozenges or equivalently the 3 spatial directions). This weight on dimer coverings of a hexagon was derived in [BGR10] (see also [Sch07] for elliptic enumeration of lattice paths).
The connection with elliptic functions will now be explained. Fix a horizontal coordinate i, denote by w(i, j) the weight of the horizontal lozenge with top vertex coordinates (i, j), and observe that for two consecutive vertical positions we have (u 1 u 2 u 3 = 1): In 3-dimensional coordinates (x, y, z) pictured in Figure 3 (note we only consider surfaces in 3 dimensions that are stepped, meaning there is a 1-1 correspondence between the 2D tiling picture and the 3D surface picture) with i = x − y, j = z − (x + y)/2, the weight ratio looks like whereũ 1 = q y+z−2x u 1 ,ũ 2 = q x+z−2y u 2 ,ũ 3 = q x+y−2z u 3 , u 1 u 2 u 3 = 1 and (x, y, z) is the 3-dimensional centroid of the 1 × 1 × 1 full cube (on the left in Figure 4) with top lid the horizontal lozenge with top vertex coordinate (i, j).
The word elliptic now becomes clear as r in (8) is an elliptic function of q (that is, defined on Esee the Introduction for details). Moreover, r is the unique elliptic function of q with zeros atũ 1 ,ũ 2 ,ũ 3 and poles at 1/ũ 1 , 1/ũ 2 , 1/ũ 3 normalized such that r(1) = 1. Of interest is also that r is elliptic inũ k for k = 1, 2, 3 subject to the condition that 3 k=1ũ k = 1. Remark 2.2. r is invariant under the natural action of S 3 permuting theũ k 's (and of course the 3 axes: x, y, z).
We can view our tilings as stepped surfaces composed of 1 × 1 × 1 cubes bounded by the 6 planes x = 0, y = 0, z = 0, x = b, y = c, z = a. Then the 2 dimensional picture in Figure 1 can be viewed as a projection of the 3 dimensional stepped surface onto the plane x + y + z = 0.
For T a tiling, we have where (i, j) are the coordinates of the top vertex of a horizontal lozenge. Grouping all 1 × 1 × 1 cubes into columns in the z direction with fixed (x, y) coordinates (see Figure 3), we obtain: where the product is taken over all cubes (visible and hidden) of the boxed plane partition and (i, j) is the top coordinate of the bounding hexagon of a 1 × 1 × 1 cube. Note to get to this equality we have merely observed that wt(empty box) is a constant independent of i and j. We can further refine this (rearranging the terms in the product and gauging away more constants -see Section 2.3 of [BGR10] for more details) as: where v = (x 0 , y 0 , z 0 ) ranges over all vertices on the border (but not on the bounding hexagon) of the stepped surface with x 0 , y 0 , z 0 integers (equivalently, v ranges over all vertices of the triangular lattice inside the hexagon, but we view v in 3 dimensions). h(v) is the distance from v to the plane x + y + z = 0 divided by

Positivity of the weight
The content of the previous subsection shows that in order to make the whole model well defined as a probabilistic model, it suffices to establish positivity of the elliptic weight ratio r(i, j) = w(i, j)/w(i, j −1) defined in (7) (where (i, j) is the location of a given horizontal tiling and ranges over all possible horizontal tilings inside the hexagon). Recall that r(i, j) = q 3 θ p (ũ 1 /q,ũ 2 /q,ũ 3 /q) θ p (qũ 1 , qũ 2 , qũ 3 ) whereũ 1 = q j−3i/2 u 1 ,ũ 2 = q j+3i/2 u 2 ,ũ 3 = q −2j u 3 and u 1 u 2 u 3 = 1. We recall that r is elliptic inũ k for k = 1, 2, 3 as well as in q. In order to make r positive, we will first restrict ourselves to the case where r is real valued. This means r is defined over a real elliptic curve, and we have −1 < p = 0 < 1 (a priori, p is complex of modulus less than 1; p ∈ (−1, 1) − {0} is equivalent to E being defined over R -for more on real elliptic curves, we refer the reader to Chapter 5 of [Sil94]). We then ensure positivity of r by an explicit computation. We will of course have two cases: p < 0 and p > 0. We deal with the case p > 0 throughout (and make remarks when necessary for p < 0). Now that we have restricted ourselves to real elliptic curves E, we first note that q ∈ E (i.e., r is elliptic as a function of q). For a chosen 0 < p < 1 there are two non-isomorphic elliptic curves defined over R (since Gal(C/R) = Z/2Z), both homeomorphic to a disjoint union of two circles (every real elliptic curve is topologically homeomorphic to a circle if p < 0 or with a disjoint union of two circles if p > 0 -one can just see this by plotting the Weierstrass equation in R 2 and compactifying): We will call the first case real and the second trigonometric (abusing terminology, since both are real elliptic curves). We will analyze the trigonometric case, but the real case is similar. In the trigonometric case, the curve has two connected components (circles): the identity component (it contains the points 1 and -1) and another component that contains the other 2-torsion points: ± √ p. There will be 3 cases to be analyzed which we list now and motivate after (if p < 0 there is only one component so the 3 cases coalesce to only one -Case 2.): • Case 1. q lies on the non identity component (|q| = √ p).
• Case 2. q and all the u k 's (and so theũ k 's) lie on the identity component (|q| = |u 1 | = |u 2 | = |u 3 | = 1) • Case 3. q and one of the u k 's lies on the identity component, the other two u k 's lie on the non identity component To analyze positivity at a fixed site (i, j) inside the hexagon, we note that r(q) has zeros atũ k and poles at 1/ũ k (k = 1, 2, 3). We note r = ±1 at q = ±1 so at least one u k (along with its reciprocal/complex conjugate 1/u k ) needs to be on the identity component (so that r can change signs on the identity component). Since r = −1 at q = ± √ p and u 1 u 2 u 3 = 1, either exactly one or all three of the u's need to be on the identity component. This motivates the three choices above. Case 1. will never lead to positivity for all four admissible sites (i, j) inside a 1 × 2 × 2 hexagon (see Figure 5), so we can eliminate it (if a 1 × 2 × 2 hexagon is never positive, much larger ones which are of interest to us will also never be as they contain the 1 × 2 × 2 case). For a proof, we suppose that u 1 is on the identity component, and u 2 , u 3 are (along with q) on the non identity component (the case where all three u's are on the identity component is handled similarly). Theũ's differ from the u's by integer powers of q given in the last three columns of the following table (for the four admissible (i, j) pairs in the 1 × 2 × 2 hexagon): -3 3 0 1/2 3 -4 5 -1 Figure 5: The admissible sites (i, j) inside a 1 × 2 × 2 hexagon Notice mod 2 (and we only care about mod 2 as q 2 is on the identity component), the four vectors (from the last 3 columns of the table) above are (1, 0, 1), (0, 0, 0), (1, 1, 0), (0, 1, 1). The corresponding u k 's we get by multiplying each u k by q to the power coming from the vector (0, 1, 1) -(ũ 1 ,ũ 2 ,ũ 3 ) = (q −4 v 1 , q 5 v 2 , q −1 v 3 ) -will all be on the identity component, which means the elliptic weight ratio will be negative at the site (i, j) = (1/2, 3) as q is on the non identity component. This is a contradiction. The other cases are handled similarly, leading to contradictions. This proves q must be on the identity component, so only cases 2. and 3. above can lead to positive hexagons. Figure 6: The identity component of E ∼ =R {u ∈ C * /p Z : |u| 2 ∈ {1, p}}. For positivity of r throughout the hexagon (i.e., for all admissibleũ k 's), q must always be closer to 1 than anyũ ±1 k as depicted.
I will next discuss the case where q and all u k are on the identity component (case 2. above; for case 3. the reasoning is similar). For a fixed site (i, j) inside the hexagon, the 3ũ k 's and their reciprocals (complex conjugates) break down the unit circle into 6 arcs (see Figure 6) and q must be on one of the three arcs where r is positive (as depicted in the figure). If we want to ensure positivity of the ratio for all 4 admissible sites (i, j) within a given 1 × 2 × 2 hexagon ( Figure 5), we first observe that for |x| = 1: So we reduce to positivity of the corresponding four functions i=1,2,3 1−ũi/q 1−ũiq . Through standard trigono-metric manipulations we thus want positivity of each of the following functions: where 2πα i = arg u i , α 1 +α 2 +α 3 ∈ {0, 1, 2}, 2πα = arg q and (α, α 1 , α 2 ) are defined on the 3-dimensional unit torus R 3 /Z 3 . If we restrict to the fundamental domain [0, 1] 3 and look at all the regions (polytopes) cut out by the planes (linear functions) in the arguments of the sines above (divided by π), we find (by solving the appropriate linear programs via Mathematica) that there exists only one region of positivity for all 4 functions. We can characterize the region best in terms of Figure 6. That is, as (i, j) range over all 4 sites inside a 1 × 2 × 2 hexagon, there should not be anyũ k (k = 1, 2, 3) or anyũ −1 k on the arc subtended by 1 and q (and that does not contain -1).
Remark 2.3. In view of the above, for any choice of a reasonably large hexagon (say one that contains a 1 × 2 × 2 hexagon) and parameters u 1 , u 2 , u 3 (satisfying the balancing condition), the set of q giving rise to nonnegative weights is a symmetric closed arc containing 1.

Degenerations of the weight
Certain degenerations of the weight have been studied before (among the relevant sources for our purposes are [BG09], [BGR10], [Joh05], [KO07], [Gor08]) from many angles. For example, when q = 1 the weight in (6) becomes a constant independent of the position of the horizontal lozenges, and so we are looking at uniformly distributed tilings of the appropriate hexagon. An exact sampling algorithm to sample such a tiling was constructed in [BG09] and the theory behind this (as well as behind other results for such tilings) is closely connected to the theory of discrete Hahn orthogonal polynomials (see [Joh05], [BG09], [Gor08]). The frozen boundary phenomenon (the shape of a "typical boxed plane partition") was first proven in [CLP98] and then via alternate techniques (and generalized) in [CKP01] and [KO07].
A more general limit than the above is the following: in (6) we let v 1 = v 2 = κ √ p and then let p → 0. This is the q-Racah limit (named after the discrete orthogonal polynomials that appear in the analysis). This limit is the most general limit that can be analyzed by orthogonal polynomials (as q-Racah polynomials sit at the top of the q-Askey scheme -see [KS]). Up to gauge equivalence, we obtain the weight of a horizontal lozenge with top corner (i, j) as: This weight was studied in [BGR10]. If we take κ to 0 or ∞, we see the q-Racah weight is an interpolation between two types of weights: A direct alternative limit from the elliptic level is given by v 1 = v 2 = p 1/3 , p → 0 (and then replace q 2 by q or 1/q). These two weights give rise to tilings weighted proportional to q Volume or q -Volume (where Volume = number of 1 × 1 × 1 cubes in the stepped surface representing a tiling). This is the q-Hahn weight (as the analysis leads to q-Hahn orthogonal polynomials). The frozen boundary phenomenon for this type of weight was first studied in [KO07], and then via alternative methods in [BGR10].
Finally, the Racah weight is the limit q → 1 in (9) (we denote k = log q (κ) and need κ → 1 as q → 1). The weight function becomes Notice in all these limits the weight of a horizontal lozenge is independent of the horizontal coordinate of its top vertex.
Taking these limits corresponds to the hypergeometric hierarchy of special functions involved in the analysis via the orthogonal polynomial (OP) or biorthogonal elliptic functions (down arrows denote limits): Elliptic hypergeometric (elliptic weights; elliptic biorthogonal ensembles) ⇓ q-hypergeometric (q-weights; q-OP ensembles) ⇓ Hypergeometric (uniform/Racah weight; Hahn/Racah OP ensemble) As a side final note, the most general degeneration of the weight is the top level trigonometric limit p → 0, which gives rise to a 3 parameter family of weights (the use of the word trigonometric here should not be confused with its usage in Section 2.3). Being more general (more parameters) than the q-Racah limit, its analysis requires q rational biorthogonal functions rather than orthogonal polynomials. We will not use this limit hereinafter, as we can approximate the trigonometric level by choosing p really small at the elliptic level.

Canonical coordinates
It will be convenient for various computations to express the geometry of an elliptic lozenge tiling in terms of coordinates on a certain product of elliptic curves. First we will introduce 6 parameters A, B, C, D, E, F depending on q, t, S, T, N, v 1 , v 2 (note we have listed, other than q, 6 parameters, of which 4 are discrete and dictate the geometry: t, S, T, N ). t here is a (discrete) time parameter and ranges from 0 to T . It will be explained better in Section 3. It corresponds to the fact that we will be interested in distributions of particles on a certain vertical line: that is, tilings of hexagons that have prescribed positions of particles (or holes) on the vertical line with horizontal coordinate t. The set of parameters is: Observe that Recall that the weight function (to be more precise, the ratio of weights of a full to an empty 1 × 1 × 1 box -see (8)) depends on the geometry of the hexagon via the three parametersũ 1 ,ũ 2 ,ũ 3 ( ũ k = 1) which in the (i, j) coordinates are:ũ Remark 2.4. We can see from above that there exists a bijection between the six bounding edges of our hexagon and the 6 parameters A, B, C, D, E, F . That is, to an edge we assign the parameter that appears to the power ±3 above. The 6 parameters are not independent (they satisfy one balancing condition ABCDEF = q 3−2N ), but then neither are the 6 edges (they must satisfy the condition that the hexagon they form is tillable by the three types of rhombi, which in this case tautologically means the edges form the 6 visible frame-edges of a rectangular parallelepiped). See Figure 7. With (12) in mind we have a (local) map R 2 → E 2 (where E 2 is isomorphic to the subvariety of E 3 with coordinates (ũ 1 ,ũ 2 ,ũ 3 ) and relation ũ i = 1) which embeds our hexagon in E 2 : Note that E 2 is the square of a real elliptic curve if parameters are chosen so that the weight ratio (of full to empty 1 × 1 × 1 box) is real positive. Hence as E is homeomorphic to a circle or a disjoint union of two circles, the above embeds our hexagon in a 2-dimensional real torus (base field = R).

Distributions and transition probabilities
In this section we compute the N -point correlation function and transitional probabilities for the model under study. We refer the reader to the Appendix of [BGR10] for the relevant application of Kasteleyn's theorem which makes these computations easy and to Kasteleyn's original paper for the theory itself Take a collection of N non-intersecting lattice paths in Ω(N, S, T ). Fix a vertical line inside the hexagon with horizontal integer coordinate t (0 ≤ t ≤ T ). This vertical line will contain N particles Depending on the geometry of our hexagon, there are four ways in which we can fix a vertical line with horizontal coordinate t inside a collection of N non-intersecting paths in Ω(N, S, T ). They are described below (see also Figure 8 in which the four cases are depicted -we only depict the outside bounding hexagon and the middle vertical line which is the desired particle line): We make the following notations: L t (X) = sum of products of weights corresponding to holes (horizontal lozenges) to the left of the vertical line with coordinate t. The sum is taken over all possible ways of tiling the region to the left of this line. Equivalently, it is taken over all families of paths starting at ((0, 0), ..., (0, N − 1)) and ending at ((t, x 1 ), ..., (t, x N )).
R t (X) = sum of products of weights corresponding to holes to the right of the vertical line with coordinate t. The sum is taken over all possible ways of tiling the region to the right of this line. Equivalently, it is taken over all families of paths starting at ((t, x 1 ), ..., (t, x N )) and ending at ((T, S), ..., (T, S + N − 1)).
C t (X) = product of weights corresponding to the holes on this vertical line. Let Remark 3.1. As observed in the introduction, ϕ t,S (x, y) = −ϕ t,S (y, x) so the product k<l ϕ t,S (x k , x l ) is the "elliptic" analogue of the Vandermonde product k<l (x k − x l ) (to which it tends in the limit p → 0, q → 1 as explained in the Introduction).
Proposition 3.2. We have Proof. This follows from an elaborate calculation and Lemma 10.2 in Appendix A of [BGR10] (which is in essence an application of Kasteleyn's theorem and the computation of the appropriate inverse Kasteleyn matrix). First, as is in the case of the aforementioned lemma, we restrict ourselves to the case S < t < T − S (Case 2. in (13); Computations are similar for the other 3 cases). Note in such a case we have N particles and S holes on the line with abscissa t. We then apply the particle-hole involution (as the weight in Lemma 10.2 in Appendix A of [BGR10] is given in terms of the positions of the holes = horizontal lozenges on the t-line). There are two types of products appearing in the total weight in question: a univariate one over the holes and a bivariate Vandermonde-like (again over the holes). For the first product, we just reciprocate to turn it into a product over particles (as the total product over holes and particles of the functions involved is a constant dependent only on t, S, T, N, q, p, v 1 , v 2 ). For the Vandermonde-like product, we note for a function f satisfying f (y i , y j ) = −f (y j , y i ) we have (up to a possible sign not depending on holes or particles): where y's represent locations of holes (top vertices of horizontal lozenges) and x's locations of particles. We take f = ϕ t,S as defined in (14). Finally, in Appendix A of [BGR10], the convention is that particles and holes are counted from the top going down. This is opposite to the convention in this paper, so we substitute x l → S + N − 1 − x l . After standard manipulations with theta-Pochhammer symbols we arrive at the desired result.
Proof. Similar to the previous proof except we use Lemma 10.3 in Appendix A of [BGR10].
Proposition 3.4. We have Proof. This weight is (up to a constant not depending on holes or particles) the reciprocal of the total weight of the S holes (horizontal lozenges) on the t line and the latter is readily computed from the definition (6).
Remark 3.6. The above distribution is what was called in the Introduction the discrete elliptic Selberg density. That is to say, where λ ∈ m n (m = S + N − 1, n = N ) and λ i + N − i = x N +1−i (to account for the fact that x 1 < x 2 < ... < x N whereas partitions are always listed in non-increasing order). The particle-hole involution invoked in Proposition 3.2 then takes the following form. If λ p is the partition associated to the particle positions (at time t) via the above equation and λ h is the partition associated to the whole positions at the same time (in the case above, there are S holes), then where λ c denotes the complemented partition corresponding to λ ∈ m n (λ c i = m − λ n+1−i ) and λ denotes the dual (transposed) partition (λ i = number of parts of λ that are ≥ i). The fact that both probabilities (in terms of holes and in terms of particles) are ∆-symbols can be observed directly as shown in Proposition 3.2 or using the following relations appearing in [Rai06]: We will for brevity denote the measure described in Theorem (3.5) by ρ S,t (note it also depends on N, T, v 1 , v 2 , p, q, but it is the dependence on S and t that will be of most interest to us). Observe we can transform the factor appearing in the univariate product of the above probability into something proportional to and absorbing into the initial constant anything independent of x (of the particle positions x k ). After using (10) our probability distribution becomes where w is the weight function for the discrete elliptic univariate biorthogonal functions discovered by Spiridonov and Zhedanov (see [SZ00], [SZ01]). It is of course also the discrete elliptic Selberg density for N = 1 (hence a ∆-symbol in n = 1 variable as seen in (4)). Notice in (18) above B and E play a special role, as does F . This will become more transparent in Section 6. w is elliptic in q, v 1 , v 2 and q {t,S,T,N } (or, analogously, in A, B, C, D, E, F, q).
Remark 3.7. Note that in the definition of w above, the first line is given in terms of the geometry of the hexagon and the choice of the particular particle line (Case 2. in (13) as previously discussed), while the second line is intrinsic and the geometry of the hexagon only comes in after using (10). We can also define the equivalent of (10) in the other 3 cases described in (13) (and the three other choices of 6 parameters differ from (10) by (a): interchanging S ant t, (b): shifting the 6 parameters in (10) by q ±(t+S−T ) ), or (c): a combination of both (a) and (b)). We will not use this any further, as all calculations will be done in Case 2. from (13).
Remark 3.8. The limit v 1 = v 2 = κ √ p, p → 0 gives the distributions present in [BGR10] at the q-Racah level. Also, as will be seen in Section 6, such probabilities are structurally a product of a "Vandermondelike" determinant squared (the first two products in (18)) and a product over the particles of univariate weights of elliptic biorthogonal functions. Indeed, under the appropriate limits, one can arrive from (18) to a much simpler (prototypical) such N -point function: the joint density of the N eigenvalues of a GUE N × N random matrix.
The transition and co-transition probabilities for the Markov chain X(t) are given by the next two statements.
We are now in a position to define six stochastic matrices (Markov chains) needed in what will follow. Their stochasticity along with other properties will be proven in Section 4, although we know the first two are stochastic as they represent the transition probabilities obtained in this section. To condense notation, we denote z k = F q x k . Let: be defined by: The normalizing constants are independent of the x k 's and the y k 's. They will become explicit in Section 4.
Note that t− P S,t S− , under interchanging t and S, becomes P S,t t− . Under the same procedure t+ P S,t S+ becomes P S,t t+ . We can think of P S,t t+ (P S,t t− ) as a Markov chain that increases (decreases) t, while t± P S,t S+ ( t± P S,t S− ) increases (decreases) S.

Elliptic difference operators
In this section we explain how recent results on elliptic special functions and elliptic difference operators intrinsically capture the model we described thus far. The main two references are [Rai10] and [Rai06] and we will state results from these without going into the proofs (with a few exceptions where the proofs are short and revealing of common techniques employed in the area). The focus will be on certain elliptic difference operators satisfying normalization, quasi-commutation and quasi-adjointness relations. We define them abstractly in the first subsection. We then turn to motivating the definitions and interpreting the operators probabilistically.

Definitions and some properties
In [Rai10] (see also [Rai06] for an algebraic description) Rains has introduced a family of difference operators acting nicely on various classes of BC n -symmetric functions. To define it, we let r 0 , r 1 , r 2 , r 3 ∈ C * satisfy r 0 r 1 r 2 r 3 = pq 1−n . Then define D(r 0 , r 1 , r 2 , r 3 ) (also depending on q, p, n) by: Remark 4.1. The difference operator above described is the special case t = q of the more general elliptic (q, t) difference operator mentioned in the references.
In view of r 0 r 1 r 2 r 3 = pq 1−n we will break symmetry and denote the difference operator by D(r 0 , r 1 , r 2 ), the fourth parameter being implicit from the balancing condition.
Remark 4.2. D takes BC n -symmetric functions to BC n -symmetric functions.
By letting D act on the function f ≡ 1, we obtain the following important lemma, whose proof we sketch following [Rai10]: Proof. By direct computation the LHS above is invariant under z k → pz k for all k (this is insured by the fact r 0 r 1 r 2 r 3 = pq 1−n ). It is also BC n -symmetric (invariant under permutations of z k 's and under z k → 1/z k ). Finally, by multiplying LHS by R = k z −1 k θ p (z 2 k ) k<l ϕ(z k , z l ) we will have cleared potential poles of the LHS. Because R is BC n -antisymmetric the result will end up being a multiple of R: R · LHS = const · R showing LHS has no singularities in the variables and is thus independent of the z i 's. Evaluating then at z i = r 0 q n−i yields the result. Observe the main point here was to prove the LHS is elliptic and has no poles in the variables, and indeed any analysis that shows this will prove the result.
The difference operators described above satisfy a number of identities, including a series of commutation relations. For an elegant proof which relies on the action of these operators on a suitably large space of functions (more precisely, on the action of the difference operators on BC n -symmetric interpolation abelian functions), see [Rai10] or [Rai06].
Remark 4.5. d λ is a special version of the interpolation theta functions P * (m,n) λ (...x k ...; a, b; q; p) defined in [Rai06] (matching the notation in the reference with ours, a = u, b = q −m−n+1 /a)). They are defined (up to normalization) by two properties: being BC n -symmetric of degree m (which happens for d λ 's) and vanishing at µ = λ (which trivially happens in our case).
so in a precise way, d λ is an interpolation Kronecker-delta theta-function. We then immediately have the following proposition: Proposition 4.6. Fix τ ∈ {±1} n . Let z k = uq λ k +n−k . Then Proof. Immediate by substituting into the definition of the difference operator (25). For any σ = τ , q σ k /2−τ k /2 z k will be of the form uq µ k +n−k with µ = λ and the corresponding summand will be 0.
A useful final property of the difference operators is their quasi-adjointness. It was shown in [Rai10] that the D's satisfy a certain "adjointness" relation that we will need in the next section. We start with 6 parameters t 0 , t 1 , t 2 , t 3 , u 0 , u 1 satisfying the balancing condition q 2n−2 t 0 t 1 t 2 t 3 u 0 u 1 = pq.
We fix the number of variables at n and λ will be a partition in m n . As in the introduction, we denote l i = λ i + n − i. We define the discrete Selberg inner product , (depending on p, q and the 6 parameters) by f, g = 1 Z λ⊆m n f (...t 0 q li ...)g(...t 0 q li ...)∆ λ (q 2n−2 t 2 0 |q n , q n−1 t 0 t 1 , q n−1 t 0 t 2 , q n−1 t 0 t 3 , q n−1 t 0 u 0 , q n−1 t 0 u 1 ; q) where f, g belong to some sufficiently nice set of functions (we will assume they are BC n -symmetric) and Z is an explicit constant that makes 1, 1 = 1. This is a discrete analogue of the continuous inner product introduced in [Rai10] and can be obtained from that by residue calculus. If the above conditions are satisfied, then ([Rai10]): where (t 0 , t 1 , t 2 , t 3 , u 0 , u 1 ) = (q 1/2 t 0 , q 1/2 t 1 , q −1/2 t 2 , q −1/2 t 3 , q 1/2 u 0 , q −1/2 u 1 ) and , is the inner product defined in (27) with primed parameters inserted throughout.

Interpretation of difference operators and their properties
We now show how the difference operators and their properties discussed in the previous section can be given probabilistic interpretations.
In what follows h k (h k ) is the location of the k-th particle on the vertical line i = t (i = t + 1) in the (i, j) frame (note according to the t → t + 1 dynamics the particles move either up or down by 1/2). We can prove the following proposition: If σ k = 1, the corresponding k-th particle at vertical position h k moves up to h k = h k + 1/2 (and if σ k = −1, the k-th particle moves down). Next observe that in the univariate product appearing in any term of (D(A, B, C)1)(...z k ...), we can change θ p (uz −b i ) (b = 1, 2) to θ p (z b i /u) by the reflection formula for theta functions and it will now match with the univariate product appearing in P S,t t+ . The product k:y k =x k +1 (...) y k =x k (...) now indeed is identical (modulo constants independent of the particle positions) to k:h k =h k +1/2 (...) k:h k =h k −1/2 (...) which is nothing more than k:σ k =1 (...) k:σ k =−1 (...) in (25). The elliptic Vandermonde product k<l appearing in (19) is the same product (modulo constants independent of the particles) as the Vandermonde-like product in any term of (D(A, B, C)1)(...z k ...) once we've transformed (in the latter product) θ p (z l /z k ) into θ p (z k /z l ) and θ p (1/z k z l ) into θ p (z k z l ) (picking up appropriate multipliers in front that will be powers of q appearing the Vandermonde-like product in (19)). The extra powers of q appearing in (19) will also surface in the difference operator once we've performed the aforementioned transformations. Finally observe that the ratio reduces (modulo the power of q up front already accounted for) to a ratio of only 2 theta functions (of the 4 initially present) because either h k − h l = h k − h l or h k + h l = h k + h l (depending whether particles k and l moved both in the same or in different directions).
Remark 4.9. We describe how the difference operators capture the particle interpretation of the model intrinsically. In their definition specialized appropriately as in the statement of the above proposition, if two consecutive particles k, k + 1 are 1 unit apart (h k+1 − h k = 1), the bottom one cannot move up and the top one down to collide because the summand in the difference operator is 0 (indeed θ p (qz k z −1 k+1 ) = θ p (1) = 0 in the cross terms). Thus, the non-intersecting condition on the paths is intrinsically built into the difference operator. A similar reasoning shows that top-most and bottom-most particles are not allowed to leave the bounding hexagon either. To exemplify, for the difference operator D(A, B, C) corresponding to the t → t + 1 transition (particles moving from left most vertical line to the right), we observe that the restriction on top (bottom) particle is not to cross the NE (SE) edge labeled C (A) in Figure 7 (or indeed not to "walk too far" to the right by crossing the B edge). However A and C are two of the parameters of the difference operator, and the corresponding terms in the univariate product in the appropriate summand in (25) become 0 once the top (bottom) particle tries to leave the hexagon. Same reasoning applies to the particles not being able to "walk too far right". Hence the difference operators intrinsically capture the boundary constraints of our model.
Remark 4.10. Proposition 4.8 is even more general, as we obtain 6 3 = 20 different stochastic matrices (Markov chains) from the 20 different difference operators (6 of them already described).
We are now in a position to prove that the 6 matrices defined in section 3 are indeed stochastic and measure preserving.
Proof. There is one way to prove these statements which works for 4 of the 6 matrices. Observe that the results for t± follow from Theorems 3.9 and 3.10, and then to observe that under t ↔ S, we have X S,t = X t,S , and ρ S,t = ρ t,S and then under interchanging S and t, P S,t t+ becomes t+ P S,t S+ (and P S,t t− becomes t− P S,t S− , respectively). This idea worked both at the q-Racah level and Hahn level (see [BGR10] and [BG09]).
Alternatively we can observe that the first two equalities are, by using (10) and Proposition 4.8, restatements of Lemma 4.3 for difference operators corresponding to parameters (A, B, C) (for P S,t t+ ), (D, E, F ) (for P S,t t− ), (A, B, D) (for t+ P S,t S+ ), (B, C, F ) (for t+ P S,t S− ), (D, E, A) (for t− P S,t S+ ), (E, F, C) (for t− P S,t S− ). Moreover, the normalizing constants that we omitted in defining the transition matrices can be recovered easily from Proposition 4.8.
The last two statements are a special case of the adjointness relation. We will prove the third statement for the t+ operator. Similar results exist for the other 5 operators. We recall that ρ S,t (X) is nothing more than the discrete elliptic Selberg density ∆ λ X (q 2N −2 F 2 |q N , q N −1 AF, q N −1 (pB)F, q N −1 CF, q N −1 DF, q N −1 EF ) defined in the introduction, with λ X,k + n − k = x n+1−k . We also define the partition λ Y to be the one corresponding to vertical line t + 1 and particle positions given by Y : λ Y,k + n − k = y n+1−k . Then one sees ρ S,t+1 (Y ) = X P S,t t+ (X, Y ) · ρ S,t (X) is equivalent to: where the prime parameters and , are defined in the previous section. The above equality (29) is only "morally correct" as we encounter the following issue: the (summands in the) difference operators D correspond to transitional probabilities in the (i, j) coordinates where particles move up or down by 1/2 from the t vertical line to the t + 1 vertical line (from Proposition 4.8). d λ , ∆ λ as well as the definition of the inner product (27) correspond to coordinates (x, t) where particles either move horizontally 1 step to the right or diagonally up by 1 from vertical line t to vertical line t + 1 (see the previous subsection and recall λ k + n − k = x n+1−k ). But this can be easily fixed since (i, j) = (t, x − t/2). However, writing an "approximate" version conveys the meaning of the quasi-adjointness of the difference operators in a notationally uncluttered way.
We finish this section with a graphical description of the 6 Markov processes described thus far. The key is to look at the domain and codomain of the difference operators in canonical coordinates. We will exemplify with the difference operator D (A, B, D), corresponding to Markov chain t+ P S,t S+ . Recall this Markov chain quasi-commutes with the t → t + 1 chain. The key is the following relation (a restatement of Theorem 4.11):  In this section, which follows closely the notation and proofs of [BG09] (see also [BGR10]), we will define a stochastic matrix that is measure preserving: it preserves the elliptic measure µ(N, S, T ) -the total mass of a hexagon tiling (collection of N non-intersecting lattice paths) in Ω(N, S, T ). Recall µ is defined as a product of the weights of the individual horizontal lozenges inside the hexagon. Viewed as a Markov chain, the input for P S S →S+1 is a hexagon of size a × b × c and the output a hexagon of size a × (b − 1) × (c + 1). Both the input and the output will turn out to be distributed according to µ(N, S, T ) and µ(N, S + 1, T ) respectively.
We define the matrix P S →S+1 : Ω(N, S, T ) × Ω(N, S + 1, T ) → [0, 1] by Theorem 5.1. The matrix P S →S+1 is stochastic and µ-measure preserving, in the sense that (32) Proof. (following [BG09]) We want to show that where the sum is taken over all Y = (Y (0), ..., Y (T )) ∈ Ω(N, S + 1, T ) such that We first sum over Y (T ) and because Y (T ) is distributed according to a singleton measure, the respective sum is 1. Next we deal with the sum , X(T )) > 0 (because of (33)). Because of the quasi-commutation relations from Theorem 4.4, we have We are summing over Y (T − 1) such that the LHS above is non-vanishing, but if it vanishes, then by the above inequality so does t+ P S+1,T −1 S− (Y (T − 1), X(T )) (one of the numerator terms in the sum over Y (T − 1) considered). This means we can drop the condition that (P S+1,T −1 t+ · t+ P S+1,T S− )(Y (T − 1), X(T )) > 0 and sum over all Y (T − 1). We obtain 1 for this sum (the denominator is independent of the summation variable, and summing the numerator over Y (T − 1) we obtain the denominator). We next sum inductively over Y (T − 2) and so on until we are left over with a sum over Y (0). This sum only has 1 term, so we obtain the desired result.

Algorithmic description of the S → S + 1 step
As before, whenever possible, we try to keep the notation similar to [BG09]. For x ∈ N we define Note p also depends on S, T, v 1 , v 2 , q, p, but we will omit these for simplicity of notation. Also note p is an elliptic function of q, q S , q T , q t , v 1 , v 2 , q x . Consider (again omitting most parameter dependence) P is just a ratio of 5 length s theta-Pochhammer symbols over 5 others (multiplied by q s−1 to make everything elliptic). We define the following probability distribution on the set {0, 1, ..., n}. P rob(s) = D(x; n)(s) = P (x; s) n j=0 P (x; j) .
• Case 2. Consider all i such that x i − y i = −1. Then z i = y i is forced.
• Case 3. For the remaining indices, group them in blocks and consider one such called a (k, l) block (where k is the smallest particle location in the block, and l is the number of particles in the block).
That is, we have y i−1 < k − 1, y i+l > k + l and the block consists of For each such block independently, we sample a random variable ξ according to the distribution D(k; l). We set z i = x i for the first ξ consecutive positions in the block, and we set z i = x i + 1 for the remainder of the l − ξ positions. We provide an example in Figure 10 below: Figure 10: Sample block split. Same picture appears in [BG09] for uniformly distributed tilings.
Theorem 5.2. By constructing Y this way, we have simulated a S → S + 1 step of the Markov chain P S →S+1 .
Proof. We perform the following computation (and are interested in Case 3. described above, that is on how to split a (k, l) block; note x i = y i in the case of interest): We thus see the blocks split independently. The probability that the first j particles in a (k, l) block stay put from Y (t) to Y (t + 1) (and the rest of l − j jump by 1) is, by using the above formula: where in (36) we have gauged away everything independent of the split position j. This probability is nothing more than the distribution D we defined in (35). This finishes the proof.

Algorithmic description of the S → S − 1 step
Similar to the P S →S+1 matrix described in the previous two sections, we can construct a P S →S−1 measure preserving Markov chain that takes random tilings in Ω(N, S, T ) and maps them to random tilings in Ω(N, S − 1, T ). We proceed exactly as in Section 5.1 and will omit most details and theorems as they transfer verbatim from Section 5.1 and the previous section (we refer the reader to [BG09] for more details on this algorithm). Given X ∈ Ω(N, S, T ) and Y (0), Y (1), ..., Y (t) already defined inductively, we choose Y (t + 1) from the distribution: We define We will also sketch the algorithm for sampling using P S →S−1 . We need to define the equivalent for p from the previous section. For x ∈ N we define As before, p is an elliptic in q, q S , q N , q t , v 1 , v 2 , q x . We also have P (x; s) = s i=1 p (x + i − 1) and the following distribution on {0, 1, ..., n}: Assuming we have X ∈ Ω(N, S, T ) with X(t + 1) = (x 1 < ... < x N ) and inductively Y (0), ..., Y (t) = (y 1 < ... < y N ), we sample Y (t + 1) = (z 1 < ... < z N ) by first observing that x i − y i ∈ {0, 1, 2} (because (P S−1,t t+ · t+ P S−1,t+1 S+ )(Y (t), X(t + 1)) > 0) and then performing appropriate updates for the following three simple cases: • Case 2. For all i with x i − y i = 2 we set z i = y i + 1.
• Case 3. For the remaining indices (for which x i − y i = 1), group them in blocks and consider one such called a (k, l) block (where k is the smallest particle location in the block, and l is the number of particles in the block). That is, we have y i−1 < k − 1, y i+l > k + l and the block consists of x i = y i + 1 = k, x i+1 = y i+1 + 1 = k + 1, ..., x i+l−1 = y i+l−1 + 1 = k + l − 1.
For each such block independently, we sample a random variable ξ according to the distribution D (k; l). We set z i = y i for the first ξ consecutive positions in the block, and we set z i = y i + 1 for the remainder of the l − ξ positions. See Figure 10.
An analogous of Theorem 5.2 exists and is proved in a similar way to show the above 3 steps are all that is necessary to simulate the Markov chain P S →S−1 .

Correlation kernel and determinantal representations
In this section we will show the process X(t) corresponding to a tiling of the hexagon is determinantal with correlation kernel given in terms of elliptic biorthogonal functions due to Spiridonov and Zhedanov ([SZ00], [Rai06]). We start by a brief overview of the necessary facts about biorthogonal functions, and continue with the heart of the proof: an application of the Eynard-Mehta theorem.

A brief overview of elliptic biorthogonal functions
We will first gather together a few results about univariate discrete elliptic biorthogonal functions. The notation and exposition will mostly be following [Rai06]. We will need to make brief use of univariate interpolation abelian functions. They were introduced in [Rai10] (see also [Rai06] for a description closer to our purposes). They are, for a fixed integer l, BC 1 -symmetric (i.e., symmetric under x → 1/x) ratios of BC 1 -symmetric theta functions of degree l with prescribed poles and zeros. To wit: Observe R * l has zeros at finitely many q-shifts of a and poles at finitely many q-shifts of b (up to taking reciprocals and shifting by p). The biorthogonal functions R l (x; t 0 : t 1 , t 2 , t 3 ; u 0 , u 1 ) discovered by Spiridonov and Zhedanov ([SZ00] in the univariate case; see [Rai10] for continuous and [Rai06] for discrete multivariate analogs) can be defined in terms of the interpolation functions as follows ( [Rai06]). Fix |p| < 1, q as well as six parameters t 0 , t 1 , t 2 , t 3 , u 0 , u 1 such that t 0 t 1 t 2 t 3 u 0 u 1 = pq. Then (dependence on p, q implied but not written) where the formula for the d k 's is explicitly given in [Rai06] and is independent of x (but of course depends on t 0 , t 1 , t 2 , t 3 ; u 0 , u 1 , q, p and k). These functions have poles at shifts of u ±1 0 (we will say u 0 controls the poles of R l (x; t 0 : t 1 , t 2 , t 3 ; u 0 , u 1 )). They are elliptic in the 6 parameters (provided the balancing condition is satisfied) as well as in the variable x. Furthermore, if in addition to the balancing condition, one also has (for some m > 0 an integer), the functions with poles controlled by u 0 and those with poles controlled by u 1 satisfy the following discrete biorthogonality relation on {0, ..., m}: 0≤s≤m R l (t 0 q s ; t 0 : t 1 , t 2 , t 3 ; u 0 , u 1 )R k (t 0 q s ; t 0 : t 1 , t 2 , t 3 ; u 1 , u 0 )∆ s (t 2 0 |q, t 0 t 1 , t 0 t 2 , t 0 t 3 , t 0 u 0 , t 0 u 1 ) = δ l,k c l where ∆ s is the univariate weight discussed in the Introduction (also appearing in section 3) and c l = const · ∆ l (1/u 0 u 1 |q, t 0 t 1 , t 0 t 2 , t 0 t 3 , 1/t 0 u 0 , 1/t 0 u 1 ) −1 = const · ∆(t 2 0 |q,t 0t1 ,t 0t2 ,t 0t3 ,t 0û0 ,t 0û1 ) −1 .
(45) The "hat" parameters are defined by the relationŝ for i = 1, 2, 3 and j = 0, 1. The "hat" is an involution. Also observe the hat parameters satisfy the same balancing conditions the original parameters satisfy. They are important because by hatting we can exchange the variable and the index of the biorthogonal functions as follows: R l (t 0 q s ; t 0 : t 1 , t 2 , t 3 ; u 0 , u 1 ) = R s (t 0 q l ;t 0 :t 1 ,t 2 ,t 3 ;û 0 ,û 1 ).
Finally, we can exchange t 0 with another t j at the choice of picking up a factor (this is in essence a renormalization so that R takes value 1 at t j rather than t 0 ):

Determinantal representations
In this section we will show the processes t → t ± 1 are determinantal point processes. For a review of such processes we direct the reader to [Bor11]. We will do the calculation for the t → t − 1 Markov process as it leads to less complicated formulas, but analogous results hold for t → t + 1. For the remainder, it is now convenient to relabel and rescale the parameter set {A, B, C, D, E, F } as {t 0 , t 1 , t 2 , t 3 , u 0 , u 1 } in order for certain symmetries to become more prominent (and in doing so, we will use the notation set forth in the previous section). To wit: Note these parameters depend on t (the time parameter), and such dependence will be made more explicit when it becomes important. Notation is as in the previous section. Note u 0 u 1 t 0 t 1 t 2 t 3 = q. Since the balancing condition for the biorthogonal functions requires a pq on the right hand side, we will again multiply u 1 by p. These are the parameters of the univariate biorthogonal functions discussed in the previous section. u 0 and u 1 control the poles of the pair of biorthogonal functions.
At the core of the computations will be the Eynard-Mehta theorem, which we now state in a "decreasing-time" form convenient for our computations (see [EM98], [Bor11] for a review and [BR05] for an elementary proof): Theorem 6.1. Assume we are given the following: • a discrete biorthonormal system (f t l , g t l ) l≥0 on l 2 ({0, 1, ..., L}) for each time t = 0, ..., T and transition probabilities proportional to Then P rob(x 1 ∈ X(τ 1 ), ..., x s ∈ X(τ s )) = det 1≤k,l≤s The first step in showing the required determinantal formulas needed to apply the Eynard-Mehta theorem is the following determinantal formula, a version of which was discovered by Warnaar (see [War02] Lemma 5.3 and Corollary 5.4 for comparison): where z k = t 0 q x k , the constant is independent of the z k 's and nonzero.
Proof. This proof is essentially the same as that of Lemma 5.3 in [War02], but is reproduced here for clarity. A first observation is that the constant in front of the right hand side will not matter much, and because it is ignored, the proof is somewhat simpler (of course, something has to be said about it not being 0). If we denote the left hand side by L and the right hand side by R, we notice both L and R are elliptic in the z k 's (for R this is a direct calculation, and for L the biorthogonal functions inside the determinant are elliptic as mentioned in the previous section though one can just see this from the definition in terms of abelian interpolation functions). Fixing a variable z k , we see poles for L/R come from the zeros of R or the poles of L. For the latter, the poles are controlled by u 0 but are exactly canceled by the zeros of 1/R appearing in the univariate product (one can see this from the definition of biorthogonal functions in terms of abelian interpolation functions). For the former the zeros of R possibly leading to poles are z k = z l , z k = 1/z l for l = k (and p shifts thereof). Plugging in z k = z l into L makes two columns the same, so L vanishes. Since univariate biorthogonal functions are BC 1 -symmetric in the variable (a fact made explicit in the previous section in the definition in terms of abelian interpolation functions; see also [Rai06]), L also vanishes if z k z l = 1 for some l = k. Hence all the poles of L/R are removable, and since L/R is elliptic, it must be constant. To show the constant is nonzero, we notice that the functions inside the determinant are linearly independent, so the columns of the determinant are linearly independent. This concludes the proof.
Remark 6.3. A more convoluted way to arrive at such determinantal representations (but the way that nevertheless suggested the formula above) would be to take the right hand side of the above formula and observe it appears in Corollary 5.4 of [War02]. What appear in the determinant on the left are the abelian interpolation functions R * l discussed in the previous section: The above formula in fact allows us to compute the constant explicitly by expanding the biorthogonal functions in terms of abelian interpolation functions (only leading coefficient is of interest for the determinant, and it is explicitly given in [Rai06]).
To simplify notation hereinafter we let Φ t l (t 0 q s ) := R l (t 0 q s ; t 0 : t 1 , t 2 , t 3 ; u 0 , pu 1 ) Ψ t l (t 0 q s ) := R l (t 0 q s ; t 0 : t 1 , t 2 , t 3 ; pu 1 , u 0 ). The t superscript for these functions stands for the fact their arguments, as it will become apparent in the next proposition, are essentially locations of the particles at time t. Likewise the parameters depend on t (t i and u j are implicit for t t i , u t j respectively; see (50) and (10)). We'll also denotẽ Thus Lemma 6.2 along with (18) and (50) yields: Proposition 6.5. We have with w 0 and w 1 as in Theorem 3.10 and Z = 1 θp(u t−1 which expresses the relation BA = 1 where A(k, l) = Φ t k (t 0 q l ), B(k, l) =Ψ t l (t 0 q k ) and we know AB = 1 by definition (see (51)). We now apply the difference operator D(u t−1 0 , t t−1 0 , t t−1 1 ) (corresponding to the Markov transition t → t − 1) to both sides and observe the parameters at time t are the required q shifts of the parameters at time t − 1 (see (48)). Finally on the right hand side we have a delta function which is acted upon by the difference operator to produce the desired result (see Proposition 4.6).
Remark 6.6. In [BGR10] and [BG09] formulas as in the above proposition involved discrete orthogonal polynomials (q-Racah and Hahn respectively) and were proven via the three term recurrence relation satisfied by these polynomials (which is an identity between hypergeometric or q-hypergeometric series). Such a relation exists for biorthogonal functions as well (we refer the reader to [SZ00] for an explicit form, though with different notation) and can be used to prove the above proposition, but the computations are more involved.
Remark 6.7. A similar result holds if we apply the transition t → t + 1 which corresponds to the operator D(u 1 , t 2 , t 3 ). For that though, we have to renormalize the biorthogonal functions at either t 2 or t 3 (see (48) and (49)), so the bidiagonal matrix that will appear on the RHS will be of the above form conjugated by two diagonal matrices (coming from the renormalization coefficients). This is an artifact of our choice of coordinates (we are counting particles going up from the bottom left edge of the hexagon).
Finally, in applying Theorem 6.1 to the t → t − 1 Markov chain X(t) we need to check that the transition probabilities have the required determinantal form. This is a consequence of Theorem 3.10, Lemma 6.2 and the following computation (the proof of which is immediate from Theorem 3.10 and Proposition 6.5; we use the notation from 3.10 for w 0 , w 1 , X, Y ): We thus obtain: Proposition 6.8.
Theorem 6.9. The Markov processes t → t ± 1 discussed in Section 3 meet the assumptions of Theorem 6.1 and are therefore determinantal.
Proof. This follows from all the results gathered in this Section for the t− Markov chain with f = Φ and g =Ψ in the notation of Theorem 6.1. For t+ see Remark 6.7.

Computer simulations
In this section we present computer simulations of the exact sampling algorithm from Section 5. We are (with one exception) looking at 200 × 200 × 200 hexagons, and parameters are chosen so the elliptic measure sampled is positive throughout the range of the algorithm (recall that the algorithm starts with a 200 × 400 × 0 box and increases c while decreasing b by 1, until it reaches the desired sizeafter 200 iterations in our case). Under each figure we list the values of the four parameters p, q, v 1 , v 2 . Computations and simulations are done using double precision, the S → S + 1 algorithm polynomial algorithm described above, and a custom program written in Java and that can handle large hexagons (in excess of N = 1000 particles) fast enough on modern CPUs. Figure 11: p = 10 −7 , q = 0.999999995, v 1 = 0.0000214, v 2 = 1.00675. 400 × 400 × 400. Because q is very close to 1, the limit shape looks uniform (recall that q = 1 gives rise to the uniform measure).
In Figure 11 we observe that the sample looks like one from the uniform measure with the arctic ellipse theoretically predicted in [CLP98] clearly visible. Figures 12 and 13 exhibit a new behavior for the arctic circle: the curve seems to acquire 3 nodes at the 3 vertices of the hexagon seen in the pictures. To obtain these shapes, the parameters have been tweaked so that the elliptic weight ratio vanishes (or = ∞) at the respective corners (in other words, the weight ratio (7) is "barely positive" as described in Section 2.3). To be more precise, we have: This fixes 3 of the 4 parameters of the measure and we have the extra degree of freedom p and so we obtain a 1-parameter family of trinodal arctic boundaries. All simulations are taken from the trigonometric positivity case (q, v 1 , v 2 are of unit modulus -see Section 2.3). While the first arctic boundary looks like an equilateral "flat" triangle, the second looks like an equilateral "thin/hyperbolic" triangle. The change from Figure 12 to 13 is an increase in p (and indeed if we increase p further the triangle will get thinner and thinner, until it will degenerate into a union of the 3 coordinate axes as p → 1). The limit p → 0 yields the same "thinning behavior" in the real positivity case.
Finally in Figure 14 we exhibit a trinodal case in the top level trigonometric case p = 0 when q, v 1 , v 2 are of unit modulus (in the case q and v i are real, arctic boundary is the union of the coordinate axes as stated above).

Appendix
In this Appendix we show how one can assign S 3 -invariant weights to the three types of rhombi (lozenges) that make up a tiling of a hexagon in the triangular lattice. We start with the 2×2×2 triangle (inside the Figure 12: An instance of a trinodal arctic boundary. p = 0.00186743, arg q = 0.000835422, arg v 1 = 0.667502, arg v 2 = −0.000835422. Figure 13: Another instance of a trinodal arctic boundary. p = 0.2, arg q = 0.000835422, arg v 1 = 0.667502, arg v 2 = −0.000835422. Note p is larger in this case than in the previous. triangular lattice) depicted in Figure 15 that contains an overlapping of the 3 types of rhombi considered for our tilings.
To each such type of rhombus we assign a label from the set {ũ 1 ,ũ 2 ,ũ 3 } (see Figure 16) such that if the rhombi are as described overlapping inside a 2 × 2 × 2 triangle we havẽ u 1ũ2ũ3 = 1.
Eachũ i will eventually be a power of q times u i (see Section 2.2). First, we can obviously shift any of such rhombi along the directions given by their edges, either upwards or downwards. If we shift the horizontal lozenge labeledũ 3 upwards-right or upwards-left, the label of the new lozenges will be   multiplied by q −1 . If we shift it downwards-right/left, the label will get multiplied by q. Naturally, if we shift directly upwards, the label will be multiplied by q −2 as a composite of an upwards-right and and upwards-left shift. A similar rule is used for lozenges with labelsũ 2 andũ 3 . The process is depicted in Figure 17, with the caveat that for labelsũ 1 andũ 2 we only show the directions in which the label gets multiplied by q (it gets multiplied by q −1 in the opposite two directions than the ones depicted). Clearly translating any lozenge along its long diagonal does not change its label.
whereũ 1 = q y+z−2x u 1 ,ũ 2 = q x+z−2y u 2 ,ũ 3 = q x+y−2z u 3 , u 1 u 2 u 3 = 1, u 1 , u 2 , u 3 are three complex numbers that multiply to 1 and (x, y, z) is the 3-dimensional coordinate of the center (intersection of the diagonals) of a lozenge. At this point we need to fix a choice of square roots: √ q, √ u 1 , √ u 2 , √ u 3 such that √ u 1 √ u 2 √ u 3 = 1. Note the 3-dimensional coordinates are only defined up to the diagonal action of Z. Figure 18 depicts the 3 lozenges with labels u i (x = y = z = 0) in the chosen coordinate system. This way of assigning weights is manifestly S 3 -invariant. To recover the same probability distribution as in Section 2.2 (i.e., a gauge-equivalent weight for tilings) we again require that the weight of a tiling of a hexagon is the product of weights of lozenges inside it. To check this, one can simply check the weight ratio of a full 1 × 1 × 1 box to an empty 1 × 1 × 1 box (this is a gauge-invariant quantity) under the present assumptions and observe the result is the same as in (7). Figure 18: Choice of a coordinate system and the 3 special parameters (lozenges centered at the origin) needed.
We then have the following proposition. Throughout, the S 3 -invariant weight is assumed.
Both these transformations leave the partition function invariant. The extra involution giving the group W (G 2 ) is a reflection through the centroid of the hexagon having coordinates: (q A/2 u 1 , q B/2 u 2 , q C/2 u 3 ) The edges transform as: where addition is mod 6. We look at the first form of the partition function written in the statement. We use the following two difference equations to simplify the calculations and arrive at the original form: Γ p,q,q (q/x) = Γ p,q,q (pqx) = Γ q,q (qx)Γ p,q,q (qx) Γ q,q (q l q m x, x) Γ q,q (q l x, q m x) = (−x) ml q −(l( m 2 )+m( l 2 ))