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 coo...
Saved in:
| Published in: | Symmetry, Integrability and Geometry: Methods and Applications |
|---|---|
| Date: | 2018 |
| Main Author: | |
| Format: | Article |
| Language: | English |
| Published: |
Інститут математики НАН України
2018
|
| Online Access: | https://nasplib.isofts.kiev.ua/handle/123456789/209540 |
| Tags: |
Add Tag
No Tags, Be the first to tag this record!
|
| Journal Title: | Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| Cite this: | Elliptically Distributed Lozenge Tilings of a Hexagon / D. Betea // Symmetry, Integrability and Geometry: Methods and Applications. — 2018. — Т. 14. — Бібліогр.: 28 назв. — англ. |
Institution
Digital Library of Periodicals of National Academy of Sciences of Ukraine| _version_ | 1860270649892667392 |
|---|---|
| author | Betea, D. |
| author_facet | Betea, D. |
| citation_txt | Elliptically Distributed Lozenge Tilings of a Hexagon / D. Betea // Symmetry, Integrability and Geometry: Methods and Applications. — 2018. — Т. 14. — Бібліогр.: 28 назв. — англ. |
| collection | DSpace DC |
| container_title | Symmetry, Integrability and Geometry: Methods and Applications |
| description | 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 BCn-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 with such tilings is determinantal with a correlation kernel given in terms of the univariate elliptic biorthogonal functions of Spiridonov and Zhedanov.
|
| first_indexed | 2025-12-07T19:06:03Z |
| format | Article |
| fulltext |
Symmetry, Integrability and Geometry: Methods and Applications SIGMA 14 (2018), 032, 39 pages
Elliptically Distributed Lozenge Tilings of a Hexagon
Dan BETEA
Paris, France
E-mail: dan.betea@gmail.com
Received October 27, 2017, in final form April 06, 2018; Published online April 12, 2018
https://doi.org/10.3842/SIGMA.2018.032
Abstract. 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 BCn-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 cer-
tain 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 simu-
lated 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 Spiri-
donov and Zhedanov.
Key words: boxed plane partitions; elliptic biorthogonal functions; particle systems; exact
sampling
2010 Mathematics Subject Classification: 33E05; 60C05; 05E05
1 Introduction
This paper examines work begun by Borodin, Gorin and Rains in [6]. 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 their Appendix. The
present paper takes the Appendix of [6] and expands upon it, following the steps in [5, 6]. How-
ever, since we are working at the elliptic (hypergeometric) level (rather than a degeneration as
in [6]), new tools are needed to generalize the results of [6]. 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 [18]
for example). In some complementary sense, while being a generalization of [6], the paper is an
application of multivariate tools introduced by Rains in [20, 21] (the first is more analytic, the
second being more algebraic). They build upon the univariate elliptic biorthogonal functions of
Spiridonov and Zhedanov from a few years earlier [26]. Work in the field of elliptic special func-
tions started with Frenkel and Turaev’s discovery of elliptic (theta) hypergeometric series [12] –
the authors of op. cit. cite Baxter’s work (see for example [1]) 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 [8, 9, 17]; correlation
kernels via orthogonal polynomials (see [5, 14, 15]); exact sampling algorithms [5]. Somewhat
central to the subject is the topic of discrete Hahn orthogonal polynomials (which themselves are
This paper is a contribution to the Special Issue on Elliptic Hypergeometric Functions and Their Applications.
The full collection is available at https://www.emis.de/journals/SIGMA/EHF2017.html
mailto:dan.betea@gmail.com
https://doi.org/10.3842/SIGMA.2018.032
https://www.emis.de/journals/SIGMA/EHF2017.html
2 D. Betea
terminating generalized hypergeometric series). One level up and we arrive at the q-distributions
on boxed plane partitions in [6] (see also [17] for the variational problem used to derive the limit
shape for the q±Volume distributions). Central to this subject are certain discrete q-orthogonal
polynomials (q-Racah, q-Hahn) from the q-Askey scheme, which themselves are terminating
q-hypergeometric series (see [13] for a full description or [18] for a distillation of the results).
The present work analyzes the elliptic level. The elliptic distribution was introduced in the
Appendix of [6], but also independently from a slightly different perspective in [24]. We look
at two aspects: exact sampling algorithms and correlation kernels. The third aspect in [5, 6] 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 [6] – fail if we replace orthogonal polynomials by elliptic
biorthogonal functions. More direct techniques like solving the variational problem described
in [17] for the q-Hahn case and in [6, Section 2.4] 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 apparent frozen
boundary can have three 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 above and combinatorics of tilings of hexagons (equivalently, dimer
coverings of the appropriate graph). We give a combinatorial interpretation to several objects
appearing in the theory of elliptic special functions: the (t = q case) multivariate elliptic dif-
ference operators discovered by Rains [21], the ∆-symbols of [20] and the (univariate) elliptic
biorthogonal functions of Spiridonov and Zhedanov [26].
This paper tries to emulate the organization of [5] and [6], but with notation heavily influenced
by [20]. It is organized as follows: in the remainder of the Introduction, we set up most of the
important notation and terminology. We set up the combinatorial and probabilistic aspects in
Section 2. 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 in depth expansion of the Appendix in [6].
Section 4 recalls some definitions and properties of elliptic tools introduced by Rains [20, 21]
(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 elliptically 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 [4]; see also [10]). The algorithm starts from a deterministic
parallelogram shape and samples relatively easy distributions to successively transform the par-
allelogram into a hexagon accordingly distributed by increasing one side by one, and decreasing
another by one; a parallelogram can be seen as a hexagon with two sides of length zero. We use
the quasi-commutation relations for the elliptic difference operators of Section 4 to construct
this algorithm.
Section 6 deals with correlations in the model. We start by recalling facts about univariate
elliptic biorthogonal functions and show that the time increasing (decreasing) Markov process is
determinantal, with correlation kernel given as a determinant of elliptic biorthogonal functions.
These replace the orthogonal polynomials discussed above.
We end with two appendices. Appendix A provides a highly symmetric view of the entire
picture. In Appendix B we present some computer simulations obtained from the algorithm
described in Section 5.
Elliptically Distributed Lozenge Tilings of a Hexagon 3
For the remainder of the section, we will set the notation that will appear in the rest of the
paper. We define the theta function and elliptic gamma function [23] as follows
θp(x) =
∏
k≥0
(
1− pkx
)(
1− pk+1/x
)
, Γp,q(x) =
∏
k,l≥0
1− pk+1ql+1/x
1− pkqlx
.
Note the elliptic gamma function is symmetric in p and q. The theta-Pochhammer symbol
(a generalization of the q-Pochhammer symbol) is defined, for m ≥ 0, as
θp(x; q)m =
∏
0≤i<m
θp
(
qix
)
.
As is customary in this area, presence of multiple arguments before the semicolon (inside theta
or elliptic gamma functions) will mean multiplication. To wit
θp
(
uz±1; q
)
m
= θp(uz; q)mθp(u/z; q)m, Γp,q(a, b) = Γp,q(a)Γp,q(b).
We have the following important, if simple, identities (for n ≥ 0 an integer)
θp(x) = θp(p/x), θp(px) = θp(1/x) = −(1/x)θp(x),
Γp,q
(
qnx
)
= θp(x; q)nΓp,q(x). (1.1)
The last identity in (1.1) can be extended for n < 0 or even for non integer n to provide
a generalization of the theta-Pochhammer symbol for negative or even non-integer lengths.
Theta-Pochhammer symbols satisfy various simple identities, a few of which we list and later
use without explicit mention
θp(a; q)n+k = θp(a; q)nθp
(
aqn; q
)
k
,
θp(a; q)n = θp
(
q1−n/a; q
)
n
(−a)nq(
n
2),
θp(a; q)n−k =
θp(a; q)n
θp
(
q1−n/a; q
)
k
(
−q
a
)k
q(
k
2)−nk,
θp
(
aq−n; q
)
k
=
θp(a; q)kθp(q/a; q)n
θp
(
q1−k/a; q
)
n
q−nk,
θp(a; q)−n =
1
θp
(
aq−n; q
)
n
=
1
θp(q/a; q)n
(
−q
a
)n
q(
n
2),
θp
(
aqn; q
)
k
=
θp(a; q)kθp
(
aqk; q
)
n
θp(a; q)n
=
θp(a; q)n+k
θp(a; q)n
.
If f(x1, . . . , xn) is a function of n variables defined on (C∗)n, we call it BCn-symmetric if it is
symmetric (does not change under permutation of the variables) and invariant under xk → 1/xk
for all k. It is called a BCn-symmetric theta function of degree m if in addition, it satisfies the
following
f(px1, . . . , xn) =
(
1
px2
1
)m
f(x1, . . . , xn).
The prototypical example of a BCn-symmetric theta function of degree one is
∏
1≤k≤n
θp
(
ux±1
k
)
.
The function ϕ(z, w) = z−1θp(zw, z/w) plays an important role. It is BC2-skewsymmetric
(symmetric under reciprocation, skewsymmetric under permutation: ϕ(z, w) = −ϕ(w, x)) of
degree one. The Weierstrass addition formula for theta functions
θp
(
xw±1
)
θp
(
yz±1
)
− θp
(
xz±1
)
θp
(
yw±1
)
= yw−1θp
(
xy±1
)
θp
(
wz±1
)
4 D. Betea
has as consequence that
ϕ(x, y) =
(
ϕ(z, x)
ϕ(w, x)
− ϕ(z, y)
ϕ(w, y)
)
ϕ(w, x)ϕ(w, y)
ϕ(z, w)
for arbitrary z, w. We note 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 ϕ(zk, zl) is an “elliptic analogue”
of the (Vandermonde) difference zk − zl. This is indeed the case if one takes the right limit
lim
q→1
lim
p→0
ϕ
(
iqxk , iqxl
)
i
(
q − q−1
) = xk − xl.
Notationally, for a function f of n variables, we will use the abbreviation f(. . . xk . . . ) to
stand for f(x1, . . . , xn).
We will make reference to the delta symbols defined in [20, 21] (we are in the case t = q in
the notation of both references). We fix λ ∈ mn 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 (λdi/2e). Then
C0
λ(x; q) =
∏
1≤i
θp
(
q1−ix; q
)
λi
,
C+
λ (x; q) =
∏
1≤i≤j
θp
(
q2−i−jx; q
)
λi+λj
θp
(
q2−i−jx; q
)
λi+λj+1
=
∏
i<j
θp
(
q2−i−jx
)
θp
(
q2−i−j+λi+λjx
) ∏
1≤i
θp
(
q2−2ix; q
)
2λi
θp
(
q2−i−nx; q
)
λi
,
C−λ (x; q) =
∏
1≤i≤j
θp
(
qj−ix; q
)
λi−λj+1
θp
(
qj−ix; q
)
λi−λj
=
∏
i<j
θp
(
qj−i−1x
)
θp
(
qj−i+λi−λj−1x
) ∏
1≤i
θp
(
qn−ix; q
)
λi
,
∆λ(a | . . . bi . . . ; q) =
C0(. . . bi . . . ; q)
C0
(
. . . pqabi . . . ; q
) · C0
2λ2(pqa; q)
C−λ (pq, q; q)C+
λ (pa, a; q)
.
Of interest will be the ∆-symbol with six parameters t0, t1, t2, t3, u0, u1 satisfying the
balancing condition q2n−2t0t1t2t3u0u1 = q. Because the usual balancing condition has pq on
the right-hand side (the reader should consult the Appendix of [21] for more on why this is
necessary), we multiply u1 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
∆λ
(
q2n−2t20 | qn, qn−1t0t1, q
n−1t0t2, q
n−1t0t3, q
n−1t0u0, q
n−1t0(pu1); q
)
= const ·
∏
i<j
(ϕ(zi, zj))
2
∏
1≤i
qli(2n−1)θp
(
z2
i
)θp(t20, t0t1, t0t2, t0t3, t0u0, t0u1; q
)
li
θp
(
q, q t0t1 , q
t0
t2
, q t0t3 , q
t0
u0
, q t0u1 ; q
)
li
= const′ ·
∏
i<j
(ϕ(zi, zj))
2 ·
∏
i
z2n−1
i θp
(
z2
i
) Γp,q(t0zi, t1zi, t2zi, t3zi, u0zi, u1zi)
Γp,q
( q
t0
zi,
q
t1
zi,
q
t2
zi,
q
t3
zi,
q
u0
zi,
q
u1
zi
) , (1.2)
where li = n − i + λi, zi = qlit0 and the constants are independent of λ and present to make
the ∆-symbol elliptic in all of its arguments. Their values are explicit [20]. This discrete elliptic
Selberg density is the weight function for the discrete elliptic multivariate biorthogonal functions
defined in [20].
We will denote by E the elliptic curve C∗/〈p〉 for some complex |p| < 1. An elliptic function f
(of one 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 xk, yk, zk) will largely be ignored and we will write const wherever this
appears; they are there to make measures into probability measures (i.e., normalizing factors)
Elliptically Distributed Lozenge Tilings of a Hexagon 5
Figure 1. A tiling of a 3× 2× 3 hexagon and the associated stepped surface.
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 as can be seen in 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.
2 The model
2.1 Interpretations
We consider random tilings of an a × b × c regular hexagon embedded in the triangular lattice
(with Cartesian coordinates (i, j)) by tiles of three types, as can be seen in the Fig. 1. The
probabilistic details are set out in Section 2.2. We will find it more convenient to encode the
hexagon via the following three numbers
N = a, T = b+ c, S = c.
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.
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 one) segments, and the paths are required not to intersect. Fig. 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 [5], let Ω(N,S, T ) denote the set of N non-intersecting paths in the
lattice N2 starting from positions (0, 0), . . . , (0, N−1) and ending at positions (T, S), . . . , (T, S+
N − 1). Each path has segments of slope zero or one (paths go either horizontally or diagonally
upwards from left to right). Set
XS,tN,T = {x ∈ Z : max(0, t+ S − T ) ≤ x ≤ min(t+N − 1, S +N − 1)},
XS,t
N,T =
{
X = (x1, . . . , xN ) ∈
(
XS,tN,T
)N
: x1 < x2 < · · · < xN
}
.
6 D. Betea
Figure 2. Duality between tilings and non-intersecting paths.
XS,tN,T is the set of all possible particle positions in a vertical section of our hexagon with horizontal
coordinate t (in (t, x) coordinates). XS,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) ∈ XS,t
N,T . X is a discrete time
Markov chain as it will be shown.
2.2 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
w(T) =
∏
l∈{horizontal lozenges}
w(l),
where by a horizontal lozenge we mean a lozenge whose diagonals are parallel to the i and j
axes respectively. The probability of such a tiling is
P(T) =
w(T)∑
S∈Ω(N,S,T )
w(S)
.
The weight function w on horizontal lozenges l is defined by
w(l) =
(u1u2)1/2qj−1/2θp
(
q2j−1u1u2
)
θp
(
qj−3i/2−1u1, qj−3i/2u1, qj+3i/2−1u2, qj+3i/2u2
)
=
(v1v2)1/2qj−S/2−1/2θp
(
q2j−S−1v1v2
)
θp
(
qj−3i/2−S−1v1, qj−3i/2−Sv1, qj+3i/2−1v2, qj+3i/2v2
) , (2.1)
where (i, j) is the coordinate of the top vertex of the horizontal lozenge l, u1, u2, q, p are complex
parameters, |p| < 1 and u1 = q−Sv1, u2 = v2 – 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 equiv-
alent to considering all types of lozenges but assigning the other two types weight one. This
is a break in symmetry that can easily be fixed – see Appendix A. However, for the remainder
of the paper we prefer this non-symmetric weight assignment system as it makes computations
easier.
This weight on lozenge tilings of a hexagon was introduced in [6] (see also [24] for an equivalent
weight on lattice paths).
Elliptically Distributed Lozenge Tilings of a Hexagon 7
Figure 3. From three dimensions to two dimensions.
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, for u1u2u3 = 1, the following weight
ration
r(i, j) =
w(i, j)
w(i, j − 1)
=
q3θp
(
qj−3i/2−1u1, q
j+3i/2−1u2, q
−2j−1u3
)
θp
(
qj−3i/2+1u1, qj+3i/2+1u2, q−2j+1u3
)
=
q3θp
(
qj−3i/2−S−1v1, q
j+3i/2−1v2, q
−2j+S−1/v1v2
)
θp
(
qj−3i/2−S+1v1, qj+3i/2+1v2, q−2j+S+1/v1v2
) . (2.2)
In three-dimensional coordinates (x, y, z) pictured in Fig. 3 with i = x− y, j = z− (x+ y)/2,
the weight ratio looks like
r(x, y, z) =
w(full box)
w(empty box)
=
q3θp(ũ1/q, ũ2/q, ũ3/q)
θp(ũ1q, ũ2q, ũ3q)
, (2.3)
where
ũ1 = qy+z−2xu1, ũ2 = qx+z−2yu2, ũ3 = qx+y−2zu3, u1u2u3 = 1,
and (x, y, z) is the three-dimensional centroid of the 1 × 1 × 1 full cube with top lid the
horizontal lozenge with top vertex coordinate (i, j).
The word elliptic now becomes clear as r in (2.3) is an elliptic function of q. 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 S3 permuting the ũk’s (and of course
the three axes: x, y, z).
We can view our tilings as stepped surfaces composed of 1× 1× 1 cubes bounded by the six
planes x = 0, y = 0, z = 0, x = b, y = c, z = a. Then the two-dimensional picture in Fig. 1 can
be viewed as a projection of the three-dimensional stepped surface onto the plane x+ y+ z = 0.
For T a tiling, we have
wt(T) =
∏
∈ T
w(i, j),
8 D. Betea
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 Fig. 3), we obtain
wt(T) = const ·
∏ w(i, j)
w(i, j − 1)
,
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 as
wt(T) = const ·
∏
v
(
w(i, j)
w(i, j − 1)
)h(v)
= const ·
∏
v
r(i, j)h(v),
where v = (x0, y0, z0) ranges over all vertices on the border (but not on the bounding hexagon)
of the stepped surface with x0, y0, z0 integers (equivalently, v ranges over all vertices of the
triangular lattice inside the hexagon, but we view v in three dimensions). h(v) is the distance
from v to the plane x+ y + z = 0 divided by
√
3.
2.3 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 (2.2), 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) =
q3θp(ũ1/q, ũ2/q, ũ3/q)
θp(qũ1, qũ2, qũ3)
,
where ũ1 = qj−3i/2u1, ũ2 = qj+3i/2u2, ũ3 = q−2ju3 and u1u2u3 = 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 6= 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, see of [25, Chapter 5]). We can 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 R2 and compactifying)
E ∼=R R∗/pZ and E ∼=R
{
u ∈ C∗/pZ : |u|2 ∈ {1, p}
}
.
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 three cases to be analyzed which we list now and motivate after (if
p < 0 there is only one component so the three cases coalesce to only one – Case 2):
• Case 1: q lies on the non-identity component: |q| = √p;
Elliptically Distributed Lozenge Tilings of a Hexagon 9
Figure 4. The admissible sites (i, j) inside a 1× 2× 2 hexagon.
• Case 2: q and all the uk’s (and so all the ũk’s) lie on the identity component (|q| = |u1| =
|u2| = |u3| = 1);
• Case 3: q and one of the uk’s lies on the identity component, the other two uk’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 uk (along with
its reciprocal/complex conjugate 1/uk) needs to be on the identity component (so that r can
change signs on the identity component). Since r = −1 at q = ±√p and u1u2u3 = 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
depicted in Fig. 4). It can thus be eliminated (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 u1 is on the identity component, and u2, u3 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 (listed are the four admissible (i, j) pairs in the 1× 2× 2 hexagon):
j i j − 3i
2 j + 3i
2 −2j
1/2 1 −1 2 −1
1 2 −2 4 −2
0 2 −3 3 0
1/2 3 −4 5 −1
Notice mod 2 (and we only care about mod 2 as q2 is on the identity component), the four
vectors (from the last three columns of the table) above are (1, 0, 1), (0, 0, 0), (1, 1, 0), (0, 1, 1).
The corresponding ũk’s we get by multiplying each uk by q to the power coming from the vector
(0, 1, 1), that is (ũ1, ũ2, ũ3) =
(
q−4v1, q
5v2, q
−1v3
)
, 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.
We will next discuss the case where q and all uk 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 three ũk’s and
their reciprocals (complex conjugates) break down the unit circle into six arcs (see Fig. 5) 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 four admissible sites (i, j) within a given 1× 2× 2 hexagon
(Fig. 4), we first observe that for |x| = 1 we have
θp(x) = (1− x)
∏
i≥1
∣∣1− pix∣∣2,
10 D. Betea
Figure 5. The identity component of E ∼=R {u ∈ C∗/pZ : |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.
so we reduce to positivity of the corresponding four functions
∏ 1−ũi/q
1−ũiq . Through standard
trigonometric manipulations we thus want positivity of each of the following functions
sinπ(α1 − α)
sinπ(α1 + α)
· sinπ(α2 − α)
sinπ(α2 + α)
· sinπ(α3 − α)
sinπ(α3 + α)
,
sinπ(α1)
sinπ(α1 + 2α)
· sinπ(α2 − 3α)
sinπ(α2 − α)
· sinπ(α3)
sinπ(α3 + 2α)
,
sinπ(α1 − 3α)
sinπ(α1 − α)
· sinπ(α2)
sinπ(α2 + 2α)
· sinπ(α3 + α)
sinπ(α3 + 3α)
,
sinπ(α1 + α)
sinπ(α1 + 3α)
· sinπ(α2 − 2α)
sinπ(α2)
· sinπ(α3 − 2α)
sinπ(α3)
,
where 2παi = arg ui, α1 +α2 +α3 ∈ {0, 1, 2}, 2πα = arg q and (α, α1, α2) ∈ R3/Z3. One way to
make all of these positive, checked by direct calculation, is depicted in Fig. 5. That is, as (i, j)
range over all four 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 not containing −1. Furthermore, numerical simulations
in Mathematica suggest this is the only way.
Remark 2.3. In view of the above, for any reasonably large hexagon (i.e., containing a 1×2×2
hexagon) and parameters u1, u2, u3 satisfying the balancing condition
∏
ui = 1, the set of q’s
giving rise to nonnegative weights is conjecturally a symmetric closed arc containing 1. This is
the only case we shall consider in what follows.
2.4 Degenerations of the weight
Certain degenerations of the weight have been studied before (among the relevant sources for
our purposes are [5, 6, 14, 15, 17]) from many angles. For example, when q = 1 the weight
in (2.1) 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 [5] and the theory behind this is closely
connected to the theory of discrete Hahn orthogonal polynomials (see [5, 14, 15]). The frozen
boundary phenomenon (the shape of a “typical boxed plane partition”) was first proven in [9]
and then via alternate techniques in [8, 17].
A more general limit than the above is the following: in (2.1) we let v1 = v2 = κ
√
p and
then let p → 0. This is the q-Racah limit (named so the discrete orthogonal polynomials that
appear in the analysis). This limit is the most general limit that can be analyzed by orthogonal
Elliptically Distributed Lozenge Tilings of a Hexagon 11
Elliptic hypergeometric (elliptic weights; elliptic biorthogonal ensembles)
↓
q-hypergeometric (q-weights; q-orthogonal polynomial ensembles)
↓
Hypergeometric (uniform/Racah weight; Hahn/Racah orthogonal polynomial ensembles)
Figure 6. A barebones schematic view of the Askey hierarchy of hypergeometric functions.
polynomials (as q-Racah polynomials sit atop the q-Askey scheme – see [18]). Up to gauge
equivalence, we obtain the weight of a horizontal lozenge with top corner (i, j) as
w(i, j) = κqj − 1
κqj
. (2.4)
This weight was studied in [6]. Upon renormalizing, if we take κ to 0 or ∞, we see the
q-Racah weight is an interpolation between two types of weights
w(i, j) = qj and w(i, j) = q−j .
A direct alternative limit from the elliptic level is given by v1 = v2 = p1/3, p → 0 (and then
replace q2 by q or 1/q). These two weights give rise to tilings weighted proportional to qVolume
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 q-Hahn orthogonal polynomials appear in its analysis.
The frozen boundary phenomenon for this type of weight was first studied in [17], and then via
alternative methods in [6].
Finally, the Racah weight is the limit q → 1 in (2.4) (we denote k = logq(κ) and need κ→ 1
as q → 1). The weight function becomes
w(i, j) = k + j.
Notice in all these limits the weight of a horizontal lozenge is independent of the horizon-
tal coordinate of its top vertex. They correspond to the hypergeometric hierarchy of special
functions involved in the algebra and analysis, depicted in Fig. 6 (down arrows are limits).
As a final sidenote, the most general degeneration of the weight is the top level trigonometric
limit p→ 0, which gives rise to a three parameter family of weights (the use of the word trigono-
metric 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.
2.5 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 six
parameters A, B, C, D, E, F depending on q, t, S, T , N , v1, v2. Note we have listed, other
than q, six parameters, of which four 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 (absence of
rhombi) on a certain vertical line: that is, tilings of hexagons that have prescribed positions of
12 D. Betea
particles (or holes) on the vertical line with horizontal coordinate t. The set of parameters is
A = qt/2+S/2−T+1/2√v1v2, B = qt/2+S/2+T+1/2
√
v2
v1
,
C = qt/2−S/2−N+1/2 1
√
v1v2
, D = q−t/2+S/2−N+1/2 1
√
v1v2
,
E = q−t/2−S/2+1/2
√
v1
v2
, F = q−t/2−S/2+1/2√v1v2. (2.5)
Observe that q2N−2ABCDEF = q.
Recall that the weight function (to be more precise, the ratio of weights of a full unit box to
an empty one in (2.3)) depends on the geometry of the hexagon via the three parameters ũ1,
ũ2, ũ3 (
∏
ũk = 1) which in the (i, j) coordinates are:
ũ1 = qj−3i/2−Sv1, ũ2 = qj+3i/2v2, ũ3 = q−2j+S/v1v2.
We want to change coordinates from (i, j) (two-dimensional) or (x, y, z) (three-dimensional)
to (ũ1, ũ2, ũ3) via the above formulas. We call these new coordinates canonical. In practice each
line of interest in the geometry has an equation in the (i, j) plane which can then be translated
in terms of the ũk’s by solving in (2.5) for t, S, T , N , v1, v2 in terms of A, B, C, D, E, F . We
thus find the following equations for the relevant edges of our hexagon:
left vertical edge (corresp. eq.: i = 0):
ũ1
ũ2
= q−Sv1/v2 =
(
ABC
DEF
)1/2
E3q−3/2,
right vertical edge (corresp. eq.: i = T ) :
ũ1
ũ2
= q−3T−Sv1/v2 =
(
ABC
DEF
)1/2
B−3q3/2,
NW edge (corresp. eq.: j = i/2 +N) :
ũ3
ũ1
= q2S−3N1/v2
1v2 =
(
ABC
DEF
)1/2
D3q−3/2,
SE edge (corresp. eq.: j = i/2− (T − S)) :
ũ3
ũ1
= q3T−S1/v2
1v2 =
(
ABC
DEF
)1/2
A−3q3/2,
NE edge (corresp. eq.: j = −i/2 + S +N) :
ũ2
ũ3
= q2S+3Nv1v
2
2 =
(
ABC
DEF
)1/2
C−3q3/2,
SW edge (corresp. eq.: j = −i/2) :
ũ2
ũ3
= q−Sv1v
2
2 =
(
ABC
DEF
)1/2
F 3q−3/2,
vertical particle line (corresp. eq.: i = t) :
ũ1
ũ2
= q−3t−Sv1/v2 =
DEF
ABC
. (2.6)
Remark 2.4. We can see from above that there exists a bijection, depicted in Fig. 7, between
the six bounding edges of our hexagon and the six parameters A, B, C, D, E, F : to an edge
we assign the parameter that appears to the power ±3 above. The six parameters are not
independent: they satisfy one balancing condition ABCDEF = q3−2N , but then neither are the
six edges: they must satisfy the condition that the hexagon they form is tilable by the three
types of rhombi.
With (2.6) in mind we have a (local) map R2 → E2, where E2 is isomorphic to the subvariety
of E3 with coordinates (ũ1, ũ2, ũ3) and relation
∏
ũi = 1, which embeds our hexagon in E2
(i, j) 7→ (ũ1, ũ2, ũ3).
Note that E2 is the square of a real elliptic curve if parameters are chosen so that the weight
ratio (of full to empty unit boxes) 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 two-dimensional real torus.
Elliptically Distributed Lozenge Tilings of a Hexagon 13
Figure 7. Correspondence between edges and the six parameters.
Figure 8. The four ways of choosing a vertical particle line (dashed) inside a hexagon. In all cases
N = 5 particles, T = 8, S ∈ {3, 5}.
3 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 [6] for the relevant application
of Kasteleyn’s theorem which makes these computations easy and to Kasteleyn’s original paper
for the theory itself [16].
Take a collection of N non-intersecting lattice paths in Ω(N,S, T ). Fix a vertical line inside
the hexagon with integer abscissa t (0 ≤ t ≤ T ). This vertical line will contain N particles
X = (x1 < · · · < xN ) ∈ XS,t
N,T . Depending on the geometry of our hexagon, there are four ways
in which we can fix this vertical line. They are described below and in Fig. 8:
Case 1: t < S, t < T − S, 0 ≤ xk ≤ t+N − 1;
Case 2: S ≤ t ≤ T − S, 0 ≤ xk ≤ S +N − 1;
Case 3: T − S ≤ t < S, t+ S − T ≤ xk ≤ t+N − 1;
Case 4: t ≥ T − S, t ≥ S, t+ S − T ≤ xk ≤ S +N − 1. (3.1)
We make use of the following notations:
Lt(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, x1), . . . , (t, xN )).
Rt(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
14 D. Betea
this line. Equivalently, it is taken over all families of paths starting at ((t, x1), . . . , (t, xN ))
and ending at ((T, S), . . . , (T, S +N − 1)).
Ct(X) = product of weights corresponding to the holes on this vertical line.
Furthermore let
ϕt,S(xk, xl) = q−xkθp
(
qxk−xl , qxk+xl+1−t−Sv1v2
)
. (3.2)
Remark 3.1. The product
∏
k<l
ϕt,S(xk, xl) is an elliptic analogue of the Vandermonde product∏
k<l
(xk − xl).
Proposition 3.2. We have
Lt(X = (x1, . . . , xN )) = const ·
∏
k<l
ϕt,S(xk, xl)
×
∏
1≤k≤N
qNxlθp
(
q2xl+1−t−Sv1v2
)θp(q1−N−t, q1−t−Sv1, q
tv2, q
1−t−Sv1v2; q
)
xl
θp
(
q, q2−2t−Sv1, qv2, q1+N−Sv1v2; q
)
xl
.
Proof. This follows from an elaborate calculation and Lemma 10.2 in Appendix A of [6] which
itself follows from Kasteleyn’s theorem.
First, we restrict ourselves to the case S < t < T − S (Case 2 in (3.1); computations are
similar for the other three cases). Note in such a case we have N particles and S holes on the line
with abscissa t. We then need to apply a particle–hole involution, as the weight in Lemma 10.2
in Appendix A of [6] 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, v1, v2). For
the Vandermonde-like product, we note for a function f satisfying f(yi, yj) = −f(yj , yi) we have∏
1≤i<j≤S
f(yi, yj) =
∏
1≤i<j≤N
f(xi, xj)
∏
0≤u<v≤S+N−1
f(u, v)
×
∏
1≤i≤N
1∏
0≤u<xi
f(xi, u)
∏
xi<u≤S+N−1
f(u, xi)
,
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 (3.2). Finally, in Appendix A of [6], 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 xl 7→ S +N − 1− xl. After standard manipulations with theta-
Pochhammer symbols we arrive at the desired result. �
Proposition 3.3. We have
Rt(X = (x1, . . . , xN )) = const ·
∏
k<l
ϕt,S(xk, xl)
×
∏
1≤k≤N
qNxlθp
(
q2xl+1−t−Sv1v2
) θp
(
q1−N−S , q−2t−Sv1, q
1+T v2, q
1−T v1v2; q
)
xl
θp
(
q1−S−t+T , q1−t−S−T v1, q2+tv2, q1+N−tv1v2; q
)
xl
.
Proof. Similar to the previous proof except we use Lemma 10.3 in Appendix A of [6]. �
Elliptically Distributed Lozenge Tilings of a Hexagon 15
Proposition 3.4. We have
Ct(X = (x1, . . . , xN )) = const ·
∏
1≤k≤N
θp
(
qxl−2t−Sv1, q
xl−2t−S+1v1, q
xl+tv2, q
xl+t+1v2
)
qxlθp
(
q2xl+1−t−Sv1v2
) .
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 (2.1). �
Theorem 3.5. We have
P(X(t) = (x1, . . . , xN )) = const ·
∏
k<l
(ϕt,S(xk, xl))
2
∏
1≤k≤N
q(2N−1)xkθp
(
q2xk+1−t−Sv1v2
)
×
∏
1≤k≤N
θp
(
q1−N−t, q1−N−S , q1−t−Sv1, q
1+T v2, q
1−T v1v2, q
1−t−Sv1v2; q
)
xk
θp
(
q, q1−S−t+T , q1−t−T−Sv1, qv2, q1+N−Sv1v2, q1+N−tv1v2; q
)
xk
= const ·
∏
k<l
(ϕt,S(xk, xl))
2
∏
1≤k≤N
q(2N−1)xkθp
(
q2xkF 2
)θp(AF,BF,CF,DF,EF, F 2; q
)
xk
θp
(
q, qAF , q
B
F , q
C
F , q
D
F , q
E
F ; q
)
xk
.
Proof. It follows from P(X(t) = (x1, . . . , xN )) ∝ Lt(X)Ct(X)Rt(X). �
Remark 3.6. The above distribution is what was called in the Introduction the discrete elliptic
Selberg density. That is to say,
P(X(t) = (x1, . . . , xN ))
= ∆λ
(
q2N−2F 2 | qN , qN−1AF, qN−1(pB)F, qN−1CF, qN−1DF, qN−1EF
)
,
where λ ∈ mn (m = S + N − 1, n = N) and λi + N − i = xN+1−i (to account for the fact
that x1 < x2 < · · · < xN whereas parts of partitions are always 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
λh =
(
λcp
)′
,
where λc denotes the complemented partition corresponding to λ ∈ mn (λci = 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 [20]
∆λ′(a | . . . bi . . . ; 1/q) = ∆λ
(
a/q2 | . . . bi . . . ; q
)
,
∆λc(a | . . . bi . . . ; q)
∆mn(a | . . . bi . . . ; q)
= ∆λ
(
q2m−2
q2na
∣∣∣ . . . qn−1bi
qma
. . . , qn, pqn, q−m, pq−m; q
)
.
We will for brevity denote the measure described in Theorem (3.5) by ρS,t (note it also
depends on N , T , v1, v2, 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
qxq(2N−2)x θp
(
q1−t−Sv1, q
1+T v2
)
θp
(
q1−t−S−T v1, qv2
)
16 D. Betea
appearing in the univariate product of the above probability into something proportional to
qx
θp
(
qN−t−Sv1, q
N+T v2
)
θp
(
q2−N−t−S−T v1, q2−Nv2
) 1
θp
(
qx+1−t−Sv1, q−x+t+S+T /v1, qx+1+T v2, q−x/v2
)
N−1
,
by using
θp
(
AqN−1; q
)
x
=
θp(A; q)xθp
(
Aqx; q
)
N−1
θp(A; q)N−1
,
θp
(
Aq1−N ; q
)
x
=
q(1−N)xθp(A; q)xθp(q/A; q)N−1
θp
(
q1−x/A; q
)
N−1
and absorbing into the initial constant anything independent of x (of the particle positions xk).
After using (2.5) our probability distribution becomes
P(X(t) = (x1, . . . , xN )) = const ·
∏
k<l
(ϕt,S(xk, xl))
2
×
∏
1≤k≤N
1
θp
(
B(Fqxk)±1, E(Fqxk)±1; q
)
N−1
∏
1≤k≤N
w(xk), (3.3)
where
w(x) =
qxθp
(
q2x+1−t−Sv1v2
)
θp
(
q1−N−t, q1−N−S, qN−t−Sv1, q
N+Tv2, q
1−Tv1v2, q
1−t−Sv1v2; q
)
x
θp
(
q1−t−Sv1v2
)
θp
(
q, q1−S−t+T, q2−N−t−T−Sv1, q2−Nv2, q1+N−Sv1v2, q1+N−tv1v2; q
)
x
=
qxθp
(
F 2q2x
)
θp
(
AF,BF
( q
ABCDEF
) 1
2 , CF,DF,EF
( q
ABCDEF
) 1
2 , F 2; q
)
x
θp
(
F 2
)
θp
(
F
Aq,
F
B q
(
ABCDEF
q
) 1
2 , FC q,
F
Dq,
F
E q
(
ABCDEF
q
) 1
2 , q; q
)
x
.
Here we note w is the weight function for the discrete elliptic univariate biorthogonal functions
of Spiridonov and Zhedanov, see [26, 27]. It is of course also the discrete elliptic Selberg density
for N = 1, hence a ∆-symbol in one variable as seen in (1.2). Notice in (3.3) above B and E
play a special role, as does F . This will become more transparent in Section 6. The weight w is
elliptic in q, v1, v2 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 geom-
etry of the hexagon and the choice of the particular particle line (Case 2 in (3.1) as previously
discussed), while the second line is intrinsic and the geometry of the hexagon only comes in
after using (2.5). We can also define the equivalent of (2.5) in the other three cases described
in (3.1) (and the three other choices of six parameters differ from (2.5) by (a): interchanging S
ant t, (b): shifting the six parameters in (2.5) 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 (3.1).
Remark 3.8. The limit v1 = v2 = κ
√
p, p → 0 gives the distributions present in [6] at the
q-Racah level. Such probabilities are also structurally a product of a Vandermonde-like deter-
minant squared (the first two products in (3.3)) and a product over the particles of univariate
weights. Indeed, under the appropriate limits, one can arrive from (3.3) 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.
Elliptically Distributed Lozenge Tilings of a Hexagon 17
Theorem 3.9. If Y = (y1, . . . , yN ) and X = (x1, . . . , xN ) such that yk − xk ∈ {0, 1} ∀ k, then
P(X(t+ 1) = Y |X(t) = X) = const ·
∏
k<l
ϕt+1,S(yk, yl)
ϕt,S(xk, xl)
∏
k : yk=xk+1
w1(xk)
∏
k : yk=xk
w0(xk),
where
w0(x) =
q−x−N+1θp
(
qx+T−t−S , qx−T−t−Sv1, q
x+t+1v2, q
x+N−tv1v2
)
θp
(
q2x+1−t−Sv1v2
) ,
w1(x) = −
q−xθp
(
qx+1−N−S , qx−2t−Sv1, q
x+T+1v2, q
x−T+1v1v2
)
θp
(
q2x+1−t−Sv1v2
) .
Proof. The formula
P(X(t+ 1) = Y |X(t) = X) =
Lt(X)Ct(X)Ct+1(Y )Rt+1(Y )
Lt(X)Ct(X)Rt(X)
=
Ct+1(Y )Rt+1(Y )
Rt(X)
along with the formulas for L, R and C yield the result. �
Theorem 3.10. If Y = (y1, . . . , yN ) and X = (x1, . . . , xN ) such that yk − xk ∈ {0,−1} ∀ k,
then
P(X(t− 1) = Y |X(t) = X) = const ·
∏
k<l
ϕt−1,S(yk, yl)
ϕt,S(xk, xl)
∏
k : yk=xk−1
w′1(xk)
∏
k : yk=xk
w′0(xk),
where
w′0(x) = −
q−xθp
(
qx−N−t+1, qx−t−S+1v1, q
x+tv2, q
x−t−S+1v1v2
)
θp
(
q2x+1−t−Sv1v2
) ,
w′1(x) =
q−x−N+1θp
(
qx, qx−2t−S+1v1, q
xv2, q
x+N−Sv1v2
)
θp
(
q2x+1−t−Sv1v2
) .
Proof. As before
P(X(t− 1) = Y |X(t) = X) =
Lt−1(X)Ct−1(X)Ct(Y )Rt(Y )
Lt(X)Ct(X)Rt(X)
=
Lt−1(Y )Ct−1(Y )
Lt(X)
. �
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 zk = Fqxk . Let
PS,tt+ : XS,t × XS,t+1 → [0, 1], PS,tt− : XS,t × XS,t−1 → [0, 1],
t+P
S,t
S+ : XS,t × XS+1,t → [0, 1], t+P
S,t
S− : XS,t × XS−1,t → [0, 1],
t−P
S,t
S+ : XS,t × XS+1,t → [0, 1], t−P
S,t
S− : XS,t × XS−1,t → [0, 1]
be defined by
PS,tt+ (X,Y ) =
const ·
∏
k<l
ϕt+1,S(yk, yl)
ϕt,S(xk, xl)
×
∏
k : yk=xk+1
−
q−xkθp
(
Azk, Bzk, Czk, q
1−Nzk/ABC
)
θp
(
z2
k
)
×
∏
k : yk=xk
q−xk−N+1θp
(
zk/A, zk/B, zk/C, q
N−1zkABC
)
θp
(
z2
k
) ,
if yk − xk ∈ {0, 1} ∀ k,
0, otherwise;
(3.4)
18 D. Betea
PS,tt− (X,Y ) =
const ·
∏
k<l
ϕt−1,S(yk, yl)
ϕt,S(xk, xl)
×
∏
k : yk=xk−1
q−xk−N+1θp
(
zk/D, zk/E, zk/F, q
N−1zkDEF
)
θp
(
z2
k
)
×
∏
k : yk=xk
−
q−xkθp
(
Dzk, Ezk, Fzk, q
1−Nzk/DEF
)
θp
(
z2
k
) ,
if yk − xk ∈ {0,−1} ∀ k,
0, otherwise;
(3.5)
t+P
S,t
S+(X,Y ) =
const ·
∏
k<l
ϕt,S+1(yk, yl)
ϕt,S(xk, xl)
×
∏
k : yk=xk+1
−
q−xkθp
(
Azk, Bzk, Dzk, q
1−Nzk/ABD
)
θp
(
z2
k
)
×
∏
k : yk=xk
q−xk−N+1θp
(
zk/A, zk/B, zk/D, q
N−1zkABD
)
θp
(
z2
k
) ,
if yk − xk ∈ {0, 1} ∀ k,
0, otherwise;
(3.6)
t+P
S,t
S−(X,Y ) =
const ·
∏
k<l
ϕt,S−1(yk, yl)
ϕt,S(xk, xl)
×
∏
k : yk=xk+1
−
q−xkθp
(
Bzk, Czk, F zk, q
1−Nzk/BCF
)
θp
(
z2
k
)
×
∏
k : yk=xk
q−xk−N+1θp
(
zk/B, zk/C, zk/F, q
N−1zkBCF
)
θp
(
z2
k
) ,
if yk − xk ∈ {0,−1} ∀ k,
0, otherwise;
(3.7)
t−P
S,t
S+(X,Y ) =
const ·
∏
k<l
ϕt,S+1(yk, yl)
ϕt,S(xk, xl)
×
∏
k : yk=xk−1
q−xk−N+1θp
(
zk/D, zk/E, zk/A, q
N−1zkDEA
)
θp
(
z2
k
)
×
∏
k : yk=xk
−
q−xkθp
(
Dzk, Ezk, Azk, q
1−Nzk/DEA
)
θp
(
z2
k
) ,
if yk − xk ∈ {0, 1} ∀ k,
0, otherwise;
(3.8)
t−P
S,t
S−(X,Y ) =
const ·
∏
k<l
ϕt,S−1(yk, yl)
ϕt,S(xk, xl)
×
∏
k : yk=xk−1
q−xk−N+1θp
(
zk/E, zk/F, zk/C, q
N−1zkEFC
)
θp
(
z2
k
)
×
∏
k : yk=xk
−
q−xkθp
(
Ezk, Fzk, Czk, q
1−Nzk/EFC
)
θp
(
z2
k
) ,
if yk − xk ∈ {0,−1} ∀ k,
0, otherwise.
(3.9)
Elliptically Distributed Lozenge Tilings of a Hexagon 19
The normalizing constants are independent of the xk’s and the yk’s. They will become explicit
in Section 4. Note that t−P
S,t
S−, under interchanging t and S, becomes PS,tt− . Under the same
procedure t+P
S,t
S+ becomes PS,tt+ . We can think of PS,tt+ (PS,tt− ) as a Markov chain that increases
(decreases) t, while t±P
S,t
S+ (t±P
S,t
S−) increases (decreases) S.
Remark 3.11. In the q-Racah limit v1 = v2 = κ
√
p, p→ 0, the chains t±P
S,t
S+ coalesce into one
(PS,tS+ in [6]). Likewise for t±P
S,t
S−.
4 Elliptic difference operators
In this section we explain how recent constructions in the field of elliptic special functions
intrinsically capture the model we described thus far. The main two references are [20, 21]
and we will state results from these without proofs (with a few exceptions where the proofs are
short). 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.
4.1 Definitions and some properties
In [20, 21] Rains has introduced a family of difference operators acting on various classes of
BCn-symmetric functions. To define them, we let r0, r1, r2, r3 ∈ C∗ satisfy r0r1r2r3 = pq1−n.
Then define D(r0, r1, r2, r3) (also depending on q, p, n) by
(D(r0, r1, r2, r3)f)(. . . zk . . . )
=
∑
σ∈{±1}n
∏
1≤k≤n
∏
0≤s≤3
θp
(
rsz
σk
k
)
θp
(
z2σk
k
) ∏
1≤k<l≤n
θp
(
qzσkk zσll
)
θp
(
zσkk zσll
) f( . . . qσk/2zk . . . ). (4.1)
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 – note in this remark alone t
has the meaning of Macdonald’s [19, Chapter VI] parameter t and not of time.
In view of r0r1r2r3 = pq1−n we will break symmetry and denote the difference operator by
D(r0, r1, r2), the fourth parameter being implied by the balancing condition.
Remark 4.2. D takes BCn-symmetric functions to BCn-symmetric functions.
By letting D act on the function f ≡ 1, we obtain the following important lemma, whose
proof we sketch following [21].
Lemma 4.3. For r0r1r2r3 = pq1−n we have
∑
σ∈{±1}n
∏
1≤k≤n
∏
0≤s≤3
θp
(
rsz
σk
k
)
θp
(
z2σk
k
) ∏
1≤k<l≤n
θp
(
qzσkk zσll
)
θp
(
zσkk zσll
) =
∏
0≤k<n
θp
(
qkr0r1, q
kr0r2, q
kr1r2
)
.
Proof. By direct computation the left-hand side above is invariant under zk → pzk for all k
(this is insured by the fact r0r1r2r3 = pq1−n). It is also BCn-symmetric. Finally, by multiplying
the left-hand side by R =
∏
k z
−1
k θp
(
z2
k
) ∏
k<l
ϕ(zk, zl) we will have cleared its potential poles.
Because R is BCn-skewsymmetric the result will end up being a multiple of R: R·LHS = const·R
showing the left-hand side has no singularities in the variables and is thus independent of the zi’s.
Evaluating at zi = r0q
n−i yields the result. �
20 D. Betea
Hereinafter we will use D for the normalized difference operator, so that D(r0, r1, r2)1 = 1,
following Lemma 4.3.
The difference operators described above satisfy a number of identities, including a series of
quasi-commutation relations. For an elegant proof which relies on the action of these operators
on a suitably large space of functions see [21] or [20].
Lemma 4.4. If U , V , W , Z are four parameters, then
D(U, V,W )D
(
q1/2U, q1/2V, q−1/2Z
)
= D(U, V, Z)D
(
q1/2U, q1/2V, q−1/2W
)
.
Next we look at the action of the difference operators on special classes of functions. For
λ ∈ mn a partition, let
dλ(. . . xk . . . ) =
∏
1≤k≤n
∏
1≤l≤m+n
θp
(
uql−1x±1
k
)
∏
1≤l≤n
θp
(
uqλl+n−lx±1
k
) .
By direct computation, we see that dλ
(
. . . uqµk+n−k . . .
)
= δλ,µcλ.
Remark 4.5. dλ is a special version of the interpolation theta functions
P
∗(m,n)
λ (. . . xk . . . ; a, b; q; p)
defined in [20] (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 BCn-symmetric of degree m
and vanishing at µ 6= λ.
If we now define dλ = dλ
cλ
, we see that
dλ
(
. . . uqµk+n−k . . .
)
= δλ,µ,
so that in a precise way, dλ is an interpolation Kronecker delta theta function. We then imme-
diately have the following proposition.
Proposition 4.6. Fix τ ∈ {±1}n. Let zk = uqλk+n−k. Then
(D(r0, r1, r2)dλ)
(
. . . q−τk/2zk . . .
)
=
∏
k
θp
(
r0z
τk
k , r1z
τk
k , r2z
τk
k ,
(
pq1−n/r0r1r2
)
zτkk
)
θp
(
z2τk
k
) ∏
k<l
θp
(
qzτkk z
τl
l
)
θp
(
zτkk z
τl
l
) .
Proof. Immediate by substituting into the definition of the difference operator (4.1). For any
σ 6= τ , qσk/2−τk/2zk will be of the form uqµk+n−k with µ 6= λ and the corresponding summand
will be zero. �
A useful final property of the difference operators is their quasi-adjointness. It was shown
in [21] that the D’s satisfy a certain adjointness relation that we will need in the next section.
We start with six parameters t0, t1, t2, t3, u0, u1 satisfying the balancing condition
q2n−2t0t1t2t3u0u1 = pq.
We fix the number of variables at n and λ will be a partition in mn. As in the introduction,
we denote li = λi + n− i. We define the discrete Selberg inner product 〈·, ·〉 (depending on p, q
and the six parameters) by
〈f, g〉 =
1
Z
∑
λ⊆mn
f
(
. . . t0q
li . . .
)
g
(
. . . t0q
li . . .
)
×∆λ
(
q2n−2t20 | qn, qn−1t0t1, q
n−1t0t2, q
n−1t0t3, q
n−1t0u0, q
n−1t0u1; q
)
, (4.2)
Elliptically Distributed Lozenge Tilings of a Hexagon 21
where f , g belong to some sufficiently nice set of functions (we will assume they are BCn-
symmetric) and Z is an explicit constant that makes 〈1, 1〉 = 1. This is a discrete analogue
of the continuous inner product introduced in [21] and can be obtained from that by residue
calculus.
If the above conditions are satisfied, then [21]
〈D(u0, t0, t1)f, g〉 =
〈
f,D(u′1, t
′
2, t
′
3)g
〉′
,
where
(t′0, t
′
1, t
′
2, t
′
3, u
′
0, u
′
1) =
(
q1/2t0, q
1/2t1, q
−1/2t2, q
−1/2t3, q
1/2u0, q
−1/2u1
)
and 〈 , 〉′ is the inner product defined in (4.2) with primed parameters inserted throughout.
4.2 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. First, observe from (2.5) that q2n−3ABCDEF = 1.
In what follows hk (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). The following proposition links difference operators with combinatorics.
Proposition 4.7. For A, B, C, D, E, F and zk = Fqhk given by (2.5), the summands in
(D(A,B,C)1)(. . . zk . . . )
(see equation (4.1)), appropriately normalized using (4.3), are equal to the transition probabilities
(entries in the stochastic matrix) PS,tt+ (H,H ′) defined in (3.4), after switching coordinates from
(x, y) back to (i, j). This statement also holds for
D(D,E, F ) and PS,tt− ,
D(A,B,D) and t+P
S,t
S+,
D(B,C, F ) and t+P
S,t
S−,
D(D,E,A) and t−P
S,t
S+,
D(E,F,C) and t−P
S,t
S−.
Proof. We will only prove the statement for D(A,B,C) and t+ (the equivalent statement for
D(D,E, F ) and t− is proved much the same way). The proof is immediate in view of (2.5),
the change of variables (X,Y ) 7→ (H,H ′) in (3.4) (to the (i, j) coordinates) and the following
observations.
First, a choice of σk ∈ {±1} for all k in the definition of D(A,B,C) is equivalent to
a choice of which particles move up/down from the position vector H (at vertical line t)
to the position vector H ′ (at vertical line t + 1). If σk = 1, the corresponding k-th par-
ticle at vertical position hk moves up to h′k = hk + 1/2 (and if σk = −1, the k-th par-
ticle moves down). Next observe that in the univariate product appearing in any term of
(D(A,B,C)1)(. . . zk . . . ), we can change θp
(
uz−bi
)
(b = 1, 2) to θp
(
zbi /u
)
by the reflection formula
for theta functions and it will now match with the univariate product appearing in PS,tt+ . The
product
∏
k : yk=xk+1
(. . . )
∏
yk=xk
(. . . ) now indeed is identical (modulo constants independent of the
particle positions) to
∏
k : h′k=hk+1/2
(. . . )
∏
k : h′k=hk−1/2
(. . . ) which is nothing more than∏
k : σk=1
(. . . )
∏
k : σk=−1
(. . . ) in (4.1).
22 D. Betea
The elliptic Vandermonde product
∏
k<l
appearing in (3.4) is the same product (modulo con-
stants) as the Vandermonde-like product in any term of (D(A,B,C)1)(. . . zk . . . ) once we have
transformed (in the latter product) θp(zl/zk) into θp(zk/zl) and θp(1/zkzl) into θp(zkzl), picking
up appropriate multipliers in front that will be powers of q appearing the Vandermonde-like
product in (3.4). The extra powers of q appearing in (3.4) will also surface in the difference
operator once we have performed the aforementioned transformations. Finally observe that the
ratio
ϕt+1,S(h′k,h
′
l)
ϕt,S(hk,hl)
reduces (modulo the power of q up front already accounted for) to a ratio
of only two theta functions (of the four initially present) because either h′k − h′l = hk − hl or
h′k + h′l = hk + hl (depending whether particles k and l moved both in the same or in different
directions). �
Remark 4.8. 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 one unit apart (hk+1 − hk = 1),
the bottom one cannot move up and the top one down to collide because the summand in
the difference operator is zero (indeed θp
(
qzkz
−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
Fig. 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 (4.1) become zero 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.9. Proposition 4.7 is even more general, as we obtain
(
6
3
)
= 20 different stochastic
matrices (Markov chains) from the twenty different difference operators, six of which we have
already described.
We are now in a position to prove that the six matrices defined in Section 3 are indeed
stochastic and measure preserving.
Theorem 4.10. We have∑
Y
PS,tt± (X,Y ) = 1,∑
Y
t±P
S,t
S±(X,Y ) = 1,
ρS,t±1(Y ) =
∑
X
PS,tt± (X,Y ) · ρS,t(X),
ρS±1,t(Y ) =
∑
X
t±P
S,t
S±(X,Y ) · ρS,t(X).
Proof. There is one way to prove these statements which works for four of the six 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
XS,t = X t,S and ρS,t = ρt,S ,
Elliptically Distributed Lozenge Tilings of a Hexagon 23
and then under interchanging S and t, PS,tt+ becomes t+P
S,t
S+ (and PS,tt− becomes t−P
S,t
S−, respec-
tively). This idea also worked in the q-Racah and Hahn limits (see [5, 6]).
Alternatively we can observe that the first two equalities are, by using (2.5) and Propo-
sition 4.7, restatements of Lemma 4.3 for difference operators corresponding to parameters
(A,B,C) (for PS,tt+ ), (D,E, F ) (for PS,tt− ), (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.7.
The last two statements are special cases of the adjointness relation. We will prove the third
statement for the t+ operator. Similar results exist for the other five operators. We recall that
ρS,t(X) is nothing more than the discrete elliptic Selberg density
∆λX
(
q2N−2F 2 | qN , qN−1AF, qN−1(pB)F, qN−1CF, qN−1DF, qN−1EF
)
defined in the introduction, with λX,k + n − k = xn+1−k. We also define the partition λY to
be the one corresponding to vertical line with abscissa t+ 1 and particle positions given by Y :
λY,k + n− k = yn+1−k. Then one sees ρS,t+1(Y ) =
∑
X P
S,t
t+ (X,Y ) · ρS,t(X) is equivalent to〈
D(A,B,C)dλY , 1
〉
=
〈
dλY ,D(D′, E′, F ′)1
〉′
, (4.3)
where the prime parameters and 〈·, ·〉′ are defined in the previous section.
The right-hand side in (4.3) equals
∑
µ dλY (. . . F qµk+n−k . . . )∆′µ = ∆′λY = ρS,t+1(Y ) (observe
∆′ = ∆ with prime parameters corresponds to the distribution of particles at the line t + 1)
while the left-hand side equals
∑
λX
P(λY |λX) · ∆λX =
∑
X P
S,t
t+ (Y |X) · ρS,t(X). The result
follows. �
We now give a graphical description of the six 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.10:∑
X
P(Y |X;A,B,D)P(X;A,B,C,D,E, F ) = P(Y ;A′, B′, C ′, D′, E′, F ′),
where
(A′, B′, C ′, D′, E′, F ′) =
(
q
1
2A, q
1
2B, q−
1
2C, q
1
2D, q−
1
2E, q−
1
2F
)
.
We note t+P
S,t
S+ corresponding to the difference operator D(A,B,D) maps marked random
tilings of hexagons determined by parameters (A,B,C,D,E, F ) to random tilings of hexagons
determined by parameters (A′, B′, C ′, D′, E′, F ′) (marked here refers to the particle line cor-
responding to parameter t). We figure what happens to the edges of such hexagons when
parameters get shifted by q±1/2 by using equations (2.6). Fig. 9 is a graphical description. In
particular, we observe t+P
S,t
S+ increases S by one. Similarly for the other difference operators:
they increase (decrease) S or t by one while leaving the other constant.
Remark 4.11. We finish by returning to the original (q, t) difference operators of [21], where
again in this remark alone the parameter t is Macdonald’s t. Using these operators one can
construct a non-determinantal elliptic process generalizing the celebrated Macdonald processes
of [3]. A dual approach to the same construction would be to use Rains’ Pieri, branching and
Cauchy identities of [20, 22]. These deserve further study. Indeed the situation is not clear
even in the Hall–Littlewood limit (essentially p, q → 0), as even this limit involves bounded
Cauchy/Littlewood identities, in contrast to the situation in [3]. Moreover, in [3], it was the
q-Whittaker limit that gave the authors a lot of traction. We do not know how to access
a q-Whittaker-like limit from the elliptic level as degenerations lead to principally specialized
Macdonald polynomials, which make setting t = 0 trivial.
24 D. Betea
Figure 9. Action of the difference operator D(A,B,D) on a tiling of a N = 2, S = 4, T = 7 hexagon
drawn in canonical coordinates. The source is marked 1 and the destination 2. Only edges relevant to
the model are considered: the six bordering edges and the particle line at horizontal displacement t from
the leftmost vertical edge. Note the slight shifting, the increase in S by one, and the fact that the particle
line’s displacement from the left vertical edge (= t) is kept constant (though particle positions are shifted
by a third step).
5 Perfect Markov chain sampling algorithm
5.1 The S 7→ S + 1 step
In this section, which follows closely the notation and proofs of [5, 6], we define a stochastic
matrix
PSS 7→S+1 : Ω(N,S, T )× Ω(N,S + 1, T )→ [0, 1],
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 ). Viewed as
a Markov chain, the input for PSS 7→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.
Given a collection of non-intersecting paths X = (X(0), . . . , X(T )) ∈ Ω(N,S, T ), we will con-
struct a (random) new collection Y = (Y (0), . . . , Y (T )) ∈ Ω(N,S+1, T ) by defining a stochastic
transition matrix PSS 7→S+1(X,Y ). Observe that Y (0) ∈ XS+1,0 = (0, . . . , N−1) is unambiguously
defined. Next we perform a sequential (inductive) update. That is, we describe how to obtain
Y (t + 1) given knowledge of Y (0), . . . , Y (t) and X. Y (t + 1) will be defined according to the
distribution
P(Y (t+ 1) = Z) =
PS+1,t
t+ (Y (t), Z) · t+PS+1,t+1
S− (Z,X(t+ 1))(
PS+1,t
t+ · t+PS+1,t+1
S−
)
(Y (t), X(t+ 1))
=
t−P
S,t+1
S+ (X(t+ 1), Z) · PS+1,t+1
t− (Z, Y (t))(
t−P
S,t+1
S+ · PS+1,t+1
t−
)
(X(t+ 1), Y (t))
,
where the last equality follows from the fact that
ρS+1,t+1(A)PS+1,t+1
t− (A,B) = ρS+1,t(B)PS+1,t
t+ (B,A)
(this is nothing more than the equality P(A ∩B) = P(A)P(B |A) = P(B)P(A |B)).
Elliptically Distributed Lozenge Tilings of a Hexagon 25
We define the matrix PS 7→S+1 : Ω(N,S, T )× Ω(N,S + 1, T )→ [0, 1] by
PS 7→S+1(X,Y ) =
T−1∏
t=0
PS+1,t
t+ (Y (t), Y (t+ 1)) · t+PS+1,t+1
S− (Y (t+ 1), X(t+ 1))(
PS+1,t
t+ · t+PS+1,t+1
S−
)
(Y (t), X(t+ 1))
,
if
T−1∏
t=0
(
PS+1,t
t+ · t+PS+1,t+1
S−
)
(Y (t), X(t+ 1)) > 0,
0, otherwise.
Theorem 5.1. The matrix PS 7→S+1 is stochastic and measure preserving, in the sense that
µ(N,S + 1, T )(Y ) =
∑
X∈Ω(N,S,T )
PS 7→S+1(X,Y )µ(N,S, T )(X). (5.1)
Proof. (following [5]) We want to show that
∑
Y
PS 7→S+1(X,Y ) =
∑
Y
T−1∏
t=0
PS+1,t
t+ (Y (t), Y (t+ 1)) · t+PS+1,t+1
S− (Y (t+ 1), X(t+ 1))(
PS+1,t
t+ · t+PS+1,t+1
S−
)
(Y (t), X(t+ 1))
= 1,
where the sum is taken over all Y = (Y (0), . . . , Y (T )) ∈ Ω(N,S + 1, T ) such that
T−1∏
t=0
(
PS+1,t
t+ · t+PS+1,t+1
S−
)
(Y (t), X(t+ 1)) > 0. (5.2)
We first sum over Y (T ) and because Y (T ) is distributed according to a singleton measure,
the respective sum is one. Next we deal with the sum
∑
Y (T−1)
PS+1,T−2
t+ (Y (T − 2), Y (T − 1)) · t+PS+1,T−1
S− (Y (T − 1), X(T − 1))(
PS+1,T−2
t+ · t+PS+1,T−1
S−
)
(Y (T − 2), X(T − 1))
over Y (T −1) satisfying
(
PS+1,T−1
t+ · t+PS+1,T
S−
)
(Y (T −1), X(T )) > 0 (because of (5.2)). Because
of the quasi-commutation relations from Theorem 4.4, we have(
PS+1,T−1
t+ · t+PS+1,T
S−
)
(Y (T − 1), X(T ))
=
(
t+P
S+1,T−1
S− · t+PS,T−1
S−
)
(Y (T − 1), X(T ))
≥ t+P
S+1,T−1
S− (Y (T − 1), X(T − 1)) · PS,T−1
t+ (X(T − 1), X(T )).
We are summing over Y (T − 1) such that the left-hand side 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 )). This means we can
drop the condition that
(
PS+1,T−1
t+ · t+PS+1,T
S−
)
(Y (T − 1), X(T )) > 0 and sum over all Y (T − 1).
We obtain one 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 one
term, so we obtain the desired result.
To show PS 7→S+1 preserves the measure µ, observe first that
µ(N,S, T )(X) = m0(X(0)) · PS,0t+ (X(0), X(1)) · · ·PS,T−1
t+ (X(T − 1), X(T )),
26 D. Betea
where m0 is the unique probability measure on any singleton set (in this case XS,0). Then the
right-hand side of (5.1) becomes
∑
X
m0(X(0))
T−1∏
t=0
PS,tt+ (X(t), X(t+ 1))
×
T−1∏
t=0
PS+1,t
t+ (Y (t), Y (t+ 1)) · t+PS+1,t+1
S− (Y (t+ 1), X(t+ 1))(
PS+1,t
t+ · t+PS+1,t+1
S−
)
(Y (t), X(t+ 1))
. (5.3)
Pulling out factors independent of the summation variables, replacing 1 = m0(X(0)) with 1 =
m0(Y (0)), using t+P
S+1,T
S− (Y (T ), X(T )) = t+P
S+1,0
S− (Y (0), X(0)) = 1 and PS,tt+ · t+P
S,t+1
S− =
t+P
S,t
S− · P
S−1,t
t+ , we transform (5.3) into
m0(Y (0))
T−1∏
t=0
PS+1,t
t+ (Y (t), Y (t+ 1))
∑
X
T−1∏
t=0
t+P
S+1,t
S− (Y (t), X(t)) · PS,tt+ (X(t), X(t+ 1))(
t+P
S,t
S− · P
S+1,t
t+
)
(Y (t), X(t+ 1))
.
Now we sum first over X(T ), then over X(T − 1) and so on like in the previous argument to
finally obtain on the left-hand side the desired result
m0(Y (0))
T−1∏
t=0
PS+1,t
t+ (Y (t), Y (t+ 1)) = µ(N,S + 1, T )(Y ). �
5.2 Algorithmic description of the S 7→ S + 1 step
As before, whenever possible, we try to keep the notation similar to [5]. For x ∈ N we define
p(x) =
qθp
(
qx−t−S+T−1, qx−t−T−1v1, q
x+t+1v2, q
x−t−S−1v1v2
)
θp
(
qx+1, qx−2t−S−1v1, qx−S+T+1v2, qx−T+1v1v2
) θp
(
q2x−t−S+1v1v2
)
θp
(
q2x−t−S−1v1v2
) .
Note p also depends on S, T , v1, v2, q, p, but we will omit these for simplicity of notation.
Also note p is an elliptic function of q, qS , qT , qt, v1, v2, qx. Consider (again omitting most
parameter dependence)
P (x; s) =
s∏
i=1
p(x+ i− 1).
P is just a ratio of five length-s theta-Pochhammer symbols over five others (multiplied by qs−1 to
make everything elliptic). We define the following probability distribution on the set {0, 1, . . ., n}:
P(s) = D(x;n)(s) =
P (x; s)
n∑
j=0
P (x; j)
. (5.4)
For the exact sampling algorithm, given X = (X(0), . . . , X(T )) ∈ Ω(N,S, T ), we will con-
struct Y = (Y (0), . . . , Y (T )) ∈ Ω(N,S + 1, T ) by first observing that Y (0) = (0, . . . , N − 1) is
uniquely defined. We then perform T sequential updates. At step t+1 we obtain Y (t+1) based
on Y (t) and X(t + 1). Suppose X(t + 1) = (x1, . . . , xN ) ∈ XS,t+1 and Y (t) = (y1, . . . , yN ) ∈
XS+1,t. We want to define/sample Y (t+ 1) = (z1, . . . , zN ) ∈ XS+1,t+1. Y (t) and X(t+ 1) satisfy
xi− yi ∈ {0,−1, 1} (follows by construction from (PS+1,t
t+ · t+PS+1,t+1
S− )(Y (t), X(t+ 1)) > 0). We
thus have three cases, and in each case we describe how to choose zi.
• Case 1: Consider all i such that xi − yi = 1. Then zi = xi is forced;
Elliptically Distributed Lozenge Tilings of a Hexagon 27
•
◦
◦
•
◦
•
•
•
◦
•
◦
•
◦
◦
◦
•
•
•
•
◦
◦
•
◦
◦
•
◦
◦
•
•
•
•
◦
•
◦
◦
•
◦
◦
•
•
◦
•
•
◦
•
◦
S S + 1
X
(t
)
X
(t
+
1)
Y
(t
)
Y
(t
+
1)
could not jump (first case)
split point determined by third case
was forced to jump (second case)
2
(k
,l
)-
b
lo
ck
Figure 10. Sample block split.
• Case 2: Consider all i such that xi − yi = −1. Then zi = yi 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 yi−1 < k − 1, yi+l > k + l and the block consists
of
xi = yi = k, xi+1 = yi+1 = k + 1, . . . , xi+l−1 = yi+l−1 = k + l − 1.
For each such block independently, we sample a random variable ξ according to the dis-
tribution D(k; l). We set zi = xi for the first ξ consecutive positions in the block, and we
set zi = xi + 1 for the remainder of the l− ξ positions. We provide an example in Fig. 10
below.
Theorem 5.2. By constructing Y this way, we have simulated a S 7→ S + 1 step of the Markov
chain PS 7→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 xi = yi in the case of interest):
P(Y (t+ 1)=Z) =
PS+1,t
t+ (Y (t), Z) · t+PS+1,t+1
S− (Z,X(t+ 1))(
PS+1,t
t+ · t+PS+1,t+1
S−
)
(Y (t), X(t+ 1))
=(factors independent of Z)
×
∏
i : zi=yi
q−yi−N+1 θp
(
qyi+T−S−t−1, qyi−T−S−t−1v1, q
yi+t+1v2, q
yi+N−tv1v2
)
θp
(
q2yi−t−Sv1v2
)
×
∏
i : zi=yi+1
q−yi
θp
(
qyi−S−N , qyi−2t−S−1v1, q
yi+T+1v2, q
yi−T+1v1v2
)
θp
(
q2yi−t−Sv1v2
)
×
∏
i : zi=xi
q−xi
θp
(
qxi−S−N , qxi−t−T−1v1, q
xi+T+1v2, q
xi−t−S−1v1v2
)
θp
(
q2xi−t−S−1v1v2
)
×
∏
i : zi=xi+1
q−xi−N
θp
(
qxi+1, qxi−T−S−t−1v1, q
xi−S+T+1v2, q
xi+N−tv1v2
)
θp
(
q2xi−t−S+1v1v2
) . (5.5)
28 D. Betea
We thus see the blocks split independently due to the evident product structure. The prob-
ability 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 one) is, by using the above formula
j−1∏
i=0
qθp
(
qk+i−t−S+T−1, qk+i−t−T−1v1, q
k+i+t+1v2, q
k+i−t−S−1v1v2
)
θp
(
q2k+2i−t−S−1v1v2
)
×
l−1∏
i=j
θp
(
qk+i+1, qk+i−2t−S−1v1, q
k+i−S+T+1v2, q
k+i−T+1v1v2
)
θp
(
q2k+2i−t−S+1v1v2
) ×(factors independent of j),
where in (5.5) we have gauged away everything independent of the split position j. This prob-
ability is nothing more than the distribution D we defined in (5.4). This finishes the proof. �
5.3 Algorithmic description of the S 7→ S − 1 step
Similar to the PS 7→S+1 matrix described in the previous two sections, we can construct a PS 7→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. Given X ∈ Ω(N,S, T ) and
Y (0), Y (1), . . . , Y (t) already defined inductively, we choose Y (t+ 1) from the distribution:
P(Y (t+ 1) = Z) =
PS−1,t
t+ (Y (t), Z) · t+PS−1,t+1
S+ (Z,X(t+ 1))(
PS−1,t
t+ · t+PS−1,t+1
S+
)
(Y (t), X(t+ 1))
.
We define
PS 7→S−1(X,Y ) =
T−1∏
t=0
PS−1,t
t+ (Y (t), Y (t+ 1)) · t+PS−1,t+1
S+ (Y (t+ 1), X(t+ 1))(
PS−1,t
t+ · t+PS−1,t+1
S+
)
(Y (t), X(t+ 1))
,
if
T−1∏
t=0
(
PS−1,t
t+ · t+PS−1,t+1
S+
)
(Y (t), X(t+ 1)) > 0,
0, otherwise.
We will also sketch the algorithm for sampling using PS 7→S−1. We need to define the equivalent
for p from the previous section. For x ∈ N we define
p′(x) =
qθp
(
qx−t−N−1, qx−t−2Sv1, q
x+tv2, q
x−t+N−1v1v2
)
θp
(
qx−S−N+1, qx−2t−Sv1, qx+Sv2, qx−S+N+1v1v2
) θp(q2x−t−S+1v1v2
)
θp
(
q2x−t−S−1v1v2
) .
As before, p′ is an elliptic in q, qS , qN , qt, v1, v2, qx. We also define P ′(x; s) =
s∏
i=1
p′(x+ i−1)
and the following distribution on {0, 1, . . . , n}:
P(s) = D′(x;n)(s) =
P ′(x; s)
n∑
j=0
P ′(x; j)
.
Assuming we have X ∈ Ω(N,S, T ) with X(t + 1) = (x1 < · · · < xN ) and inductively
Y (0), . . . , Y (t) = (y1 < · · · < yN ), we sample Y (t + 1) = (z1 < · · · < zN ) by first observing
that xi − yi ∈ {0, 1, 2} (because (PS−1,t
t+ · t+PS−1,t+1
S+ )(Y (t), X(t+ 1)) > 0) and then performing
appropriate updates for the following three simple cases:
Elliptically Distributed Lozenge Tilings of a Hexagon 29
• Case 1: For all i with xi − yi = 0 we set zi = xi;
• Case 2: For all i with xi − yi = 2 we set zi = yi + 1;
• Case 3: For the remaining indices (for which xi − yi = 1), group them in blocks and
consider one such called a (k, l)-block. That is, we have yi−1 < k− 1, yi+l > k+ l and the
block consists of
xi = yi + 1 = k, xi+1 = yi+1 + 1 = k + 1, . . . , xi+l−1 = yi+l−1 + 1 = k + l − 1.
For each such block independently, we sample a random variable ξ according to the dis-
tribution D′(k; l). We set zi = yi for the first ξ consecutive positions in the block, and we
set zi = yi + 1 for the remainder of the l − ξ positions. See Fig. 10.
An analogous of Theorem 5.2 exists and is proved in a similar way to show the above three
steps are all that is necessary to simulate the Markov chain PS 7→S−1.
6 Correlation kernel and determinantal representations
In this section we will show the process X(t) corresponding to a tiling of the hexagon is determi-
nantal with correlation kernel given in terms of the elliptic biorthogonal functions of Spiridonov
and Zhedanov [20, 26]. We start by a brief overview of the necessary facts about biorthogo-
nal functions, and continue with the heart of the proof: an application of the Eynard–Mehta
theorem.
6.1 A brief overview of elliptic biorthogonal functions
We will first gather together a few results about univariate discrete elliptic biorthogonal func-
tions. The notation and exposition will mostly be following [20]. We will need to make brief use
of univariate interpolation abelian functions. They were introduced in [20, 21] and are, for a fixed
integer l, BC1-symmetric ratios of BC1-symmetric theta functions of degree l with prescribed
poles and zeros. To wit
R∗l (x; a, b) =
θp
(
ax±1; q
)
l
θp
(
bq−lx±1; q
)
l
.
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 univariate biorthogonal functions Rl(x; t0 : t1, t2, t3;
u0, u1) of [26] can be defined in terms of the interpolation functions following [20]. Fix |p| < 1, q
as well as six parameters t0, t1, t2, t3, u0, u1 such that t0t1t2t3u0u1 = pq. Then (dependence
on p, q implied but not written)
Rl(x; t0 : t1, t2, t3;u0, u1) =
∑
0≤k≤l
dkR
∗
k(x; t0, u0) = dlR
∗
l + lower order terms,
where the formula for the dk’s is explicitly given in [20] and is independent of x (but of course
depends on t0, t1, t2, t3, u0, u1, q, p and k). These functions have poles at shifts of u±1
0 (we will
say u0 controls the poles of Rl(x; t0 : t1, t2, t3;u0, u1)). They are elliptic in the six 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
t0t1 = q−m
30 D. Betea
for some m > 0 an integer, the functions with poles controlled by u0 and those with poles
controlled by u1 satisfy the following discrete biorthogonality relation on {0, . . . ,m}∑
0≤s≤m
Rl
(
t0q
s; t0 : t1, t2, t3;u0, u1
)
Rk
(
t0q
s; t0 : t1, t2, t3;u1, u0
)
×∆s
(
t20 | q, t0t1, t0t2, t0t3, t0u0, t0u1
)
= δl,kcl,
where ∆s is the univariate delta symbol defined in the Introduction and
cl = const ·∆l
(
1/u0u1 | q, t0t1, t0t2, t0t3, 1/t0u0, 1/t0u1
)−1
= const ·∆
(
t̂20 | q, t̂0t̂1, t̂0t̂2, t̂0t̂3, t̂0û0, t̂0û1
)−1
.
The “hat” parameters are defined by the relations
t̂0 =
√
t0t1t2t3
pq
, t̂0t̂i = t0ti,
ûj
t̂0
=
uj
t0
for i = 1, 2, 3 and j = 0, 1. The “hat” is an involution and the hat parameters satisfy the same
balancing conditions as the original parameters. They are important because by hatting we can
exchange the variable and the index of the biorthogonal functions as follows:
Rl
(
t0q
s; t0 : t1, t2, t3;u0, u1
)
= Rs
(
t̂0q
l; t̂0 : t̂1, t̂2, t̂3; û0, û1
)
.
The biorthogonal functions described above have t0 as a special normalization parameter,
distinguished among the ti’s. That is, Rl(t0; t0 : t1, t2, t3;u0, u1) = 1. The normalized difference
operators of Section 4 act on the biorthogonal functions as follows
D(u0, t0, t1)Rl
((
q1/2t0
)
qs; q1/2t0 : q1/2t1, q
−1/2t2, q
−1/2t3; q1/2u0, q
−1/2u1
)
= Rl
(
t0q
s; t0 : t1, t2, t3;u0, u1
)
. (6.1)
Finally, we can exchange t0 with another tj at the choice of picking up a factor (this is in
essence a renormalization so that R takes value one at tj rather than at t0):
Rl(x; t1 : t0, t2, t3;u0, u1) =
Rl(x; t0 : t1, t2, t3;u0, u1)
Rl(t1; t0 : t1, t2, t3;u0, u1)
. (6.2)
6.2 Determinantal representations
We now show the processes t 7→ t ± 1 are determinantal point processes. For a review of such
processes we direct the reader to [2]. We will do the calculation for the t 7→ t−1 Markov process
as it leads to less complicated formulas, but analogous results hold for t 7→ t+ 1.
For the remainder, it is convenient to change the set of parameters A, B, C, D, E, F to the
set t0, t1, t2, t3, u0, u1 in order for certain symmetries to become more prominent using the
following:
A = t2, qN−1B = u1, C = t3, D = t1, qN−1E = u0, F = t0. (6.3)
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 u0u1t0t1t2t3 =q.
Since the balancing condition for the biorthogonal functions requires a pq on the right-hand side,
we will again multiply u1 by p.
We state the Eynard–Mehta theorem, in a “decreasing-time” form convenient for our com-
putations (see [2, 11] for a review and [7] for an elementary proof):
Elliptically Distributed Lozenge Tilings of a Hexagon 31
Theorem 6.1. Assume we are given the following:
• a discrete biorthonormal system
(
f tl , g
t
l
)
l≥0
on l2({0, 1, . . . , L}) for each time t = 0, . . . , T ;
• a matrix
vt→t−1(x, y) =
∑
l≥0
f t−1
l
(
tt−1
0 qx
)
gtl
(
tt0q
y
)
,
for n ≥ 0, t = 1, . . . , T and a parameter t0 changing with time;
• a discrete time Markov chain X(t) (with time decreasing from T to 0) taking values in state
spaces X t (set of possible particle positions at time t) with one-dimensional distributions
proportional to
det
1≤k,l≤N
(
f tk−1
(
tt0q
xl
))
det
1≤k,l≤N
(
gtk−1
(
tt0q
xl
))
and transition probabilities proportional to
det1≤k,l≤N (vt→t−1(xk, yl)) det1≤k,l≤N
(
f t−1
k−1
(
tt−1
0 qyl
))
det1≤k,l≤N
(
f tk−1
(
tt0q
xl
)) .
Then
P(x1 ∈ X(τ1), . . . , xs ∈ X(τs)) = det
1≤k,l≤s
(K(τk, xk; τl, xl))
where
K(τ1, x1; τ2, x2) =
∑
s≥0
f τ1s
(
tτ10 q
x1
)
gτ2s
(
tτ20 q
x2
)
, if τ1 ≥ τ2,
−
∑
s≥N
f τ1s
(
tτ10 q
x1
)
gτ2s
(
tτ20 q
x2
)
, if τ1 < τ2.
The first step in showing the required determinantal formulas needed to apply the Eynard–
Mehta theorem is the following determinantal formula due to Warnaar [28]:
Lemma 6.2. We have
det
1≤k,l≤n
Rl−1(zk; t0 : t1, t2, t3;u0, pu1) = const ·
∏
k<l
ϕ(zk, zl)
∏
k
1
θp
(
q1−nu0z
±1
k ; q
)
n−1
,
where zk = t0q
xk , the constant is independent of the zk’s and nonzero.
Proof. This proof is essentially the same as that of Lemma 5.3 in [28], 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 zero). 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 zk’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 zk, 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 u0 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
zk = zl, zk = 1/zl for l 6= k (and p shifts thereof). Plugging in zk = zl into L makes two
32 D. Betea
columns the same, so L vanishes. Since univariate biorthogonal functions are BC1-symmetric
in the variable, L also vanishes if zkzl = 1 for some l 6= 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. We arrived at the above formula noticing the right-hand side appears in Corol-
lary 5.4 of [28]. What appear in the determinant on the left are the abelian interpolation
functions R∗l discussed in the previous section
det
1≤k,l≤n
θp
(
az±1
k ; q
)
n−l
θp
(
bz±1
k ; q
)
n−l
= a(n2)q(
n
3)
∏
k<l
ϕ(zk, zl)
∏
k
θp
(
b/a, abq2n−2k; q
)
k−1
θp
(
bz±1
k ; q
)
n−1
.
The above formula in fact allows us to compute the constant explicitly by expanding the
biorthogonal functions in terms of abelian interpolation functions. Only the leading coefficient
is of interest for the determinant, and it is explicitly given in [20].
To simplify notation hereinafter we let
Φt
l
(
t0q
s
)
:= Rl
(
t0q
s; t0 : t1, t2, t3;u0, pu1
)
, Ψt
l
(
t0q
s
)
:= Rl
(
t0q
s; t0 : t1, t2, t3; pu1, u0
)
.
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 (ti and uj are implicit for tti, u
t
j respectively; see (6.3) and (2.5)).
We will also denote
Ψ̃l
(
t0q
s
)
= Ψl
(
t0q
s
)
∆s
(
t20 | q, t0t1, t0t2, t0t3, t0u0, pt0u1
)
/cl, (6.4)
so that∑
s≥0
Φk
(
t0q
s
)
Ψ̃l
(
t0q
s
)
= δk,l.
Thus Lemma 6.2 along with (3.3) and (6.3) yields the following.
Proposition 6.4. We have
P(X(t) = (x1, . . . , xN )) = const · det
1≤k,l≤n
Φt
l−1
(
t0q
xk
)
· det
1≤k,l≤n
Ψt
l−1
(
t0q
xk
)
·
∏
k
∆xk
= const · det
1≤k,l≤n
Φt
l−1
(
t0q
xk
)
· det
1≤k,l≤n
Ψ̃t
l−1
(
t0q
xk
)
.
Proposition 6.5. We have
vt→t−1(k, l) :=
∑
s≥0
Φt−1
s
(
tt−1
0 qk
)
Ψ̃t
s
(
tt0q
l
)
=
1
Z
(
w′0δk,l + w′1δk+1,l
)
with w′0 and w′1 as in Theorem 3.10 and Z = 1
θp
(
ut−1
0 tt−1
0 ,ut−1
0 tt−1
1 ,tt−1
0 tt−1
1
) .
Proof. We observe that∑
s≥0
Φt
s
(
tt0q
k
)
Ψ̃t
s
(
tt0q
l
)
= δk,l,
Elliptically Distributed Lozenge Tilings of a Hexagon 33
which expresses the relation BA = 1 where A(k, l) = Φt
k
(
t0q
l
)
, B(k, l) = Ψ̃t
l
(
t0q
k
)
and we
know AB = 1 by definition (see (6.4)). We now apply the difference operator D
(
ut−1
0 , tt−1
0 , tt−1
1
)
corresponding to the Markov transition t 7→ 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 (6.1)). 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 [5, 6] formulas as in Proposition 6.5 were proven via the three term recur-
rence relation satisfied by the orthogonal polynomial ensembles considered (q-Racah and Hahn
respectively). Such a relation exists for biorthogonal functions as well [26] and in conjunction
with arguments from say [5] provides an alternative proof of Proposition 6.5.
Remark 6.7. A similar result holds if we apply the transition t 7→ t+ 1 which corresponds to
the operator D(u1, t2, t3). For that though, we have to renormalize the biorthogonal functions
at either t2 or t3 (see (6.1) and (6.2)), so the bidiagonal matrix that will appear on the right-
hand side 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). Using the notation from Theorem 3.10 for w′0, w′1, X, Y )
det
1≤k,l≤N
(vt→t−1(xk, yl)) = const ·
∏
k : yk=xk−1
w′1(xk)
∏
k : yk=xk
w′0(xk).
We thus obtain the following.
Proposition 6.8. We have
P(X(t− 1) = Y |X(t) = X) = const · det
1≤k,l≤N
(vt→t−1(xk, yl))
det1≤k,l≤N
(
Φt−1
k−1
(
tt−1
0 qyl
))
det1≤k,l≤N
(
Φt
k−1
(
tt0q
xl
)) .
It finally leads to
Theorem 6.9. The Markov processes t 7→ 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. �
Remark 6.10. For obtaining quantitative arctic boundary-type results about our measures, we
can try to look at the asymptotics of the diagonal of the correlation kernel of the process (the
probability that a particle is present at that site)
K(x, x) =
S+N−1∑
i=0
Rti
(
t0q
x | t0 : t1, t2, t3;u0, pu1
)
Rti
(
t0q
x | t0 : t1, t2, t3; pu1, u0
)
×∆x
(
t20 | q, t0t1, t0t2, t0t3, t0u0, pt0u1
)
×∆i
(
1/(pu0u1) | q, t0t1, t0t2, t0t3, 1/(t0u0), 1/(pt0u1)
)
,
but said asymptotics appear complicated and we do not pursue them here.
34 D. Betea
Figure 11. Shifting lozenges in the triangular lattice, we multiply the labels by q or q−1 as depicted.
Figure 12. An α × β × γ hexagon with canonical coordinates of the edges on the outside and edge
lengths on the inside.
A Symmetric lozenge weights
In this appendix we show how to assign S3-invariant weights to the three types of rhombi
(lozenges) that make up a tiling of a hexagon. We start with the 2 × 2 × 2 triangle that
contains an overlap of the three types of rhombi considered. To the three different types of
rhombi in this triangle we assign labels ũ1, ũ2, ũ3 that multiply to one ũ1ũ2ũ3 = 1 using the
convention depictued here: .
Each ũi will eventually be a power of q times ui, see Section 2.2. First, we can obviously
shift any such rhombus along the directions given by its edges, either upwards or downwards.
If we shift the horizontal lozenge labeled ũ3 upwards-right or upwards-left, the label of the new
lozenge 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 (a composite
of an upwards-right and upwards-left shift). A similar rule is used for lozenges with labels ũ2
and ũ3. The process is depicted in Fig. 11. Translating any lozenge along its long diagonal does
not change its label.
To a lozenge with label ũi (i = 1, 2, 3) we assign the following weight
wt(lozenge with label ũi) = ũ
−1/2
i θp(ũi), i = 1, 2, 3,
where
ũ1 = qy+z−2xu1, ũ2 = qx+z−2yu2, ũ3 = qx+y−2zu3, u1u2u3 = 1,
u1, u2, u3 are three complex numbers that multiply to one and (x, y, z) is the three-dimensional
coordinate of the center of a lozenge. At this point we need to fix a choice of square roots:
√
q,√
u1,
√
u2,
√
u3 such that
√
u1
√
u2
√
u3 = 1. Further note the three-dimensional coordinates are
only defined up to the diagonal action of Z. The three lozenges with labels ui (x = y = z = 0)
have their centers at the hidden corner of the hexagon (the origin in Fig. 12).
This way of assigning weights is manifestly S3-invariant. The weight of a tiling of the hexagon
is the product of weights of all lozenges comprising the tiling. Furthermore, as a probability
measure, we recover the same probability distribution as in Section 2.2. To see this, one can
Elliptically Distributed Lozenge Tilings of a Hexagon 35
simply check the weight ratio of a full 1× 1× 1 box to an empty 1× 1× 1 box and observe the
result is the same as in (2.2).
The S3-invariance can be viewed at the level of the partition function (the sum of weights of
all tilings in a hexagon written in this gauge) as follows. We start with an α × β × γ hexagon.
The origin is at the hidden corner of the 3D box. In canonical coordinates
(ũ1, ũ2, ũ3) =
(
qy+z−2xu1, q
x+z−2yu2, q
x+y−2zu3
)
the six bounding edges have the following equations (see Fig. 12 for correspondence between
edges and Li’s)
L0 = ũ1/ũ2 = q3βu1/u2, L1 = ũ3/ũ1 = q−3γu3/u1,
L2 = ũ2/ũ3 = q3γu2/u3, L3 = ũ1/ũ2 = q−3αu1/u2,
L4 = ũ3/ũ1 = q3αu3/u1, L5 = ũ2/ũ3 = q−3βu2/u3.
With these weights we have the following.
Proposition A.1. The partition function for an α× β × γ hexagon is equal to
P × lim
ρ→1
Γp,q,q
(
q1+α+β+γρ, q1+αρ, q1+βρ, q1+γρ
)
Γp,q,q
(
qρ, q1+α+βρ, q1+α+γρ, q1+β+γρ
)
×
Γp,q,q
(
q1−α+β+γu1, q
1−αu1, q
1−β+α+γu2, q
1−βu2, q
1−γ+α+βu3, q
1−γu3
)
Γp,q,q
(
q1−α+βu1, q1−α+γu1, q1−β+αu2, q1−β+γu2, q1−γ+αu3, q1−γ+βu3
)
= P × lim
ρ→1
Γp,q,q
(
q(L0L2L4)
1
3 ρ, q(L0L4L5)
1
3 ρ, q(L0L1L2)
1
3 ρ, q(L2L3L4)
1
3 ρ
)
Γp,q,q
(
qρ, q(L0/L3)
1
3 ρ, q(L4/L1)
1
3 ρ, q(L2/L5)
1
3 ρ
)
×
Γp,q,q
(
q(L0L2L3)
1
3, q(L0L3L5)
1
3, q(L2L4L5)
1
3, q(L1L2L5)
1
3, q(L0L1L4)
1
3, q(L1L3L4)
1
3
)
Γp,q,q
(
q(L0/L4)
1
3, q(L3/L1)
1
3, q(L5/L3)
1
3, q(L2/L0)
1
3, q(L4/L2)
1
3 , q(L1/L5)
1
3
) ,
where
P = qαβγ−
αβ2+βα2+αγ2+γα2+βγ2+γβ2
4 u
−βγ
2
1 u
−αγ
2
2 u
−αβ
2
3 ,
Γp,q,t(x) =
∏
i,j,k≥0
(
1− pi+1qj+1tk+1/x
)(
1− piqjtkx
)
.
It is left invariant by S3 permuting the coordinates ũi. Furthermore, this invariance can be
expanded to the group W (G2) = S3 oZ2 = Dih6 (the symmetry group of a regular hexagon) with
the missing involution being the transformation
(u1, u2, u3) 7→
(
1
qAu1
,
1
qBu2
,
1
qCu3
)
,
where A = −2α+ β + γ, B = α− 2β + γ, C = α+ β − 2γ.
Proof. We start with the elliptic MacMahon identity derived in the Appendix of [6]∑
tilings T wt(T,G)
wt(0, G)
= qαβγ
∏
1≤x≤α
1≤y≤β
1≤z≤γ
θp
(
qx+y+z−1, qy+z−x−1u1, q
x+z−y−1u2, q
x+y−z−1u3
)
θp
(
qx+y+z−2, qy+z−xu1, qx+z−yu2, qx+y−zu3
) ,
where 0 denotes the empty tiling (box) and G is any gauge equivalent to the ones used in this
paper (that is to say, both sides are gauge-independent). For G the S3-invariant gauge herein
36 D. Betea
discussed, the formula for the empty tiling multiplied by the right-hand side above simplifies
the partition function via straightforward computations. We arrive at the desired result using
the following transformations for gamma functions
Γp,q(qx) = θp(x)Γp,q(x), Γp,q,t(tx) = Γp,q(x)Γp,q,t(x).
The limit ρ→ 1 is needed for technical reasons to avoid zeros of triple gamma functions.
For S3-invariance, it suffices to show how edges transform under the 3-cycle (ũ1, ũ2, ũ3) →
(ũ2, ũ3, ũ1) (a 120◦ clockwise rotation) and the transposition ũ1 ↔ ũ2 (a reflection in the z axis).
For the 3-cycle, the new edges (denoted with primes) have equations
L′i = Li+2,
where +2 is taken mod 6, while for the transposition we have
L′0 = 1/L3, L′1 = 1/L2, L′2 = 1/L1, L′3 = 1/L0, L′4 = 1/L5, L′5 = 1/L4.
Both these transformations leave the partition function invariant. The extra involution giving
the group W (G2) is a reflection through the centroid of the hexagon having coordinates(
qA/2u1, q
B/2u2, q
C/2u3
)
,
so that the edges transform as
L′i = 1/Li+3,
where +3 is taken 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
(
qlqmx, x
)
Γq,q
(
qlx, qmx
) = (−x)mlq−l(
m
2 )−m(l2). �
B 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 one, until it
reaches the desired size – after 200 iterations in our case). Under each figure we list the values of
the four parameters p, q, v1, v2. Computations and simulations are done using double precision,
the S 7→ S + 1 algorithm polynomial algorithm described above, and a custom program written
in Java that can handle large hexagons (in excess of N = 1000 particles) fast enough on modern
CPUs.
In Fig. 13 we observe that the sample looks like one from the uniform measure with the arctic
ellipse theoretically predicted in [9] clearly visible. Figs. 14 and 15 exhibit a new behavior for
the arctic circle: the curve seems to acquire three conjectural nodes at the three 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 (2.2) is “barely positive” as described in Section 2.3. To be more precise, we have q = e
2πi
T−1 ,
v1 = q2T−1, v2 = 1/q. This fixes three of the four parameters of the measure and we have the
extra degree of freedom p and so we obtain a one parameter family of conjecturally trinodal
Elliptically Distributed Lozenge Tilings of a Hexagon 37
Figure 13. p = 10−7, q = 0.999999995, v1 = 0.0000214, v2 = 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).
Figure 14. An instance of a trinodal arctic boundary. p = 0.00186743, arg q = 0.000835422, arg v1 =
0.667502, arg v2 = −0.000835422.
arctic boundaries. All simulations are taken from the trigonometric positivity case (q, v1, v2
are of unit modulus – see Section 2.3). While the first arctic boundary looks like an equilateral
“flat” triangle, the second looks “thinner”. The change from Fig. 14 to 15 is an increase in p.
Indeed if we increase p further the triangle will get thinner and thinner, until it will degenerate
into a union of the three coordinate axes as p → 1. The limit p → 0 yields the same “thinning
behavior” in the real positivity case. Finally in Fig. 16 we exhibit a trinodal case in the top
level trigonometric case p = 0 when q, v1, v2 are of unit modulus.
Acknowledgements
The author would like to thank Alexei Borodin, Fokko van de Bult, Vadim Gorin, and Eric
Rains for their help through numerous conversations. He is also indebted to Igor Pak and Greta
Panova for putting the tiling picture herein described into perspective, and to three anonymous
38 D. Betea
Figure 15. Another instance of a trinodal arctic boundary. p = 0.2, arg q = 0.000835422, arg v1 =
0.667502, arg v2 = −0.000835422. Note p is larger in this case than in the previous.
Figure 16. Top level trigonometric p = 0 case. As above, arg q = 0.000835422, arg v1 = 0.667502,
arg v2 = −0.000835422.
referees for improving the clarity of the manuscript. This article was written while the author was
a graduate student in the Department of Mathematics at the California Institute of Technology
to which many remerciements are due for all its support during the five years the author spent
there.
References
[1] Baxter R.J., Exactly solved models in statistical mechanics, Academic Press, Inc., London, 1982.
[2] Borodin A., Determinantal point processes, in The Oxford Handbook of Random Matrix Theory, Oxford
University Press, Oxford, 2011, 231–249, arXiv:0911.1153.
[3] Borodin A., Corwin I., Macdonald processes, Probab. Theory Related Fields 158 (2014), 225–400,
arXiv:1111.4408.
https://arxiv.org/abs/0911.1153
https://doi.org/10.1007/s00440-013-0482-3
https://arxiv.org/abs/1111.4408
Elliptically Distributed Lozenge Tilings of a Hexagon 39
[4] Borodin A., Ferrari P.L., Anisotropic growth of random surfaces in 2 + 1 dimensions, Comm. Math. Phys.
325 (2014), 603–684, arXiv:0804.3035.
[5] Borodin A., Gorin V., Shuffling algorithm for boxed plane partitions, Adv. Math. 220 (2009), 1739–1770,
arXiv:0804.3071.
[6] Borodin A., Gorin V., Rains E.M., q-distributions on boxed plane partitions, Selecta Math. (N.S.) 16 (2010),
731–789, arXiv:0905.0679.
[7] Borodin A., Rains E.M., Eynard–Mehta theorem, Schur process, and their Pfaffian analogs, J. Stat. Phys.
121 (2005), 291–317, math-ph/0409059.
[8] Cohn H., Kenyon R., Propp J., A variational principle for domino tilings, J. Amer. Math. Soc. 14 (2001),
297–346, math.CO/0008220.
[9] Cohn H., Larsen M., Propp J., The shape of a typical boxed plane partition, New York J. Math. 4 (1998),
137–165, math.CO/9801059.
[10] Diaconis P., Fill J.A., Strong stationary times via a new form of duality, Ann. Probab. 18 (1990), 1483–1522.
[11] Eynard B., Mehta M.L., Matrices coupled in a chain. I. Eigenvalue correlations, J. Phys. A: Math. Gen. 31
(1998), 4449–4456, cond-mat/9710230.
[12] Frenkel I.B., Turaev V.G., Elliptic solutions of the Yang–Baxter equation and modular hypergeometric
functions, in The Arnold–Gelfand Mathematical Seminars, Birkhäuser Boston, Boston, MA, 1997, 171–204.
[13] Gasper G., Rahman M., Basic hypergeometric series, Encyclopedia of Mathematics and its Applications,
Vol. 96, 2nd ed., Cambridge University Press, Cambridge, 2004.
[14] Gorin V.E., Nonintersecting paths and the Hahn orthogonal polynomial ensemble, Funct. Anal. Appl. 42
(2008), 180–197, arXiv:0708.2349.
[15] Johansson K., Non-intersecting, simple, symmetric random walks and the extended Hahn kernel, Ann. Inst.
Fourier (Grenoble) 55 (2005), 2129–2145, math.PR/0409013.
[16] Kasteleyn P.W., Graph theory and crystal physics, in Graph Theory and Theoretical Physics, Academic
Press, London, 1967, 43–110.
[17] Kenyon R., Okounkov A., Limit shapes and the complex Burgers equation, Acta Math. 199 (2007), 263–302,
math-ph/0507007.
[18] Koekoek R., Lesky P.A., Swarttouw R.F., Hypergeometric orthogonal polynomials and their q-analogues,
Springer Monographs in Mathematics, Springer-Verlag, Berlin, 2010.
[19] Macdonald I.G., Symmetric functions and Hall polynomials, 2nd ed., Oxford Mathematical Monographs,
Oxford Science Publications, The Clarendon Press, Oxford University Press, New York, 1995.
[20] Rains E.M., BCn-symmetric Abelian functions, Duke Math. J. 135 (2006), 99–180, math.CO/0402113.
[21] Rains E.M., Transformations of elliptic hypergeometric integrals, Ann. of Math. 171 (2010), 169–243,
math.QA/0309252.
[22] Rains E.M., Elliptic Littlewood identities, J. Combin. Theory Ser. A 119 (2012), 1558–1609,
arXiv:0806.0871.
[23] Ruijsenaars S.N.M., First order analytic difference equations and integrable quantum systems, J. Math.
Phys. 38 (1997), 1069–1146.
[24] Schlosser M., Elliptic enumeration of nonintersecting lattice paths, J. Combin. Theory Ser. A 114 (2007),
505–521, math.CO/0602260.
[25] Silverman J.H., Advanced topics in the arithmetic of elliptic curves, Graduate Texts in Mathematics,
Vol. 151, Springer-Verlag, New York, 1994.
[26] Spiridonov V., Zhedanov A., Spectral transformation chains and some new biorthogonal rational functions,
Comm. Math. Phys. 210 (2000), 49–83.
[27] Spiridonov V., Zhedanov A., Generalized eigenvalue problem and a new family of rational functions biorthog-
onal on elliptic grids, in Special Functions 2000: Current Perspective and Future Directions (Tempe, AZ),
NATO Sci. Ser. II Math. Phys. Chem., Vol. 30, Kluwer Acad. Publ., Dordrecht, 2001, 365–388.
[28] Warnaar S.O., Summation and transformation formulas for elliptic hypergeometric series, Constr. Approx.
18 (2002), 479–502, math.QA/0001006.
https://doi.org/10.1007/s00220-013-1823-x
https://arxiv.org/abs/0804.3035
https://doi.org/10.1016/j.aim.2008.11.008
https://arxiv.org/abs/0804.3071
https://doi.org/10.1007/s00029-010-0034-y
https://arxiv.org/abs/0905.0679
https://doi.org/10.1007/s10955-005-7583-z
https://arxiv.org/abs/math-ph/0409059
https://doi.org/10.1090/S0894-0347-00-00355-6
https://arxiv.org/abs/math.CO/0008220
https://arxiv.org/abs/math.CO/9801059
https://doi.org/10.1214/aop/1176990628
https://doi.org/10.1088/0305-4470/31/19/010
https://arxiv.org/abs/cond-mat/9710230
https://doi.org/10.1017/CBO9780511526251
https://doi.org/10.1007/s10688-008-0027-1
https://arxiv.org/abs/0708.2349
https://doi.org/10.5802/aif.2155
https://doi.org/10.5802/aif.2155
https://arxiv.org/abs/math.PR/0409013
https://doi.org/10.1007/s11511-007-0021-0
https://arxiv.org/abs/math-ph/0507007
https://doi.org/10.1007/978-3-642-05014-5
https://doi.org/10.1215/S0012-7094-06-13513-5
https://arxiv.org/abs/math.CO/0402113
https://doi.org/10.4007/annals.2010.171.169
https://arxiv.org/abs/math.QA/0309252
https://doi.org/10.1016/j.jcta.2012.03.001
https://arxiv.org/abs/0806.0871
https://doi.org/10.1063/1.531809
https://doi.org/10.1063/1.531809
https://doi.org/10.1016/j.jcta.2006.07.002
https://arxiv.org/abs/math.CO/0602260
https://doi.org/10.1007/978-1-4612-0851-8
https://doi.org/10.1007/s002200050772
https://doi.org/10.1007/978-94-010-0818-1_14
https://doi.org/10.1007/s00365-002-0501-6
https://arxiv.org/abs/math.QA/0001006
1 Introduction
2 The model
2.1 Interpretations
2.2 Probabilistic model
2.3 Positivity of the weight
2.4 Degenerations of the weight
2.5 Canonical coordinates
3 Distributions and transition probabilities
4 Elliptic difference operators
4.1 Definitions and some properties
4.2 Interpretation of difference operators and their properties
5 Perfect Markov chain sampling algorithm
5.1 The S S+1 step
5.2 Algorithmic description of the S S+1 step
5.3 Algorithmic description of the S S-1 step
6 Correlation kernel and determinantal representations
6.1 A brief overview of elliptic biorthogonal functions
6.2 Determinantal representations
A Symmetric lozenge weights
B Computer simulations
References
|
| id | nasplib_isofts_kiev_ua-123456789-209540 |
| institution | Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| issn | 1815-0659 |
| language | English |
| last_indexed | 2025-12-07T19:06:03Z |
| publishDate | 2018 |
| publisher | Інститут математики НАН України |
| record_format | dspace |
| spelling | Betea, D. 2025-11-24T10:49:46Z 2018 Elliptically Distributed Lozenge Tilings of a Hexagon / D. Betea // Symmetry, Integrability and Geometry: Methods and Applications. — 2018. — Т. 14. — Бібліогр.: 28 назв. — англ. 1815-0659 2010 Mathematics Subject Classification: 33E05; 60C05; 05E05 arXiv: 1110.4176 https://nasplib.isofts.kiev.ua/handle/123456789/209540 https://doi.org/10.3842/SIGMA.2018.032 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 BCn-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 with such tilings is determinantal with a correlation kernel given in terms of the univariate elliptic biorthogonal functions of Spiridonov and Zhedanov. The author would like to thank Alexei Borodin, Fokko van de Bult, Vadim Gorin, and Eric Rains for their help through numerous conversations. He is also indebted to Igor Pak and Greta Panova for putting the tiling picture herein described into perspective, and to three anonymous referees for improving the clarity of the manuscript. This article was written while the author was a graduate student in the Department of Mathematics at the California Institute of Technology, to which many remerciements are due for all its support during the five years the author spent there. en Інститут математики НАН України Symmetry, Integrability and Geometry: Methods and Applications Elliptically Distributed Lozenge Tilings of a Hexagon Article published earlier |
| spellingShingle | Elliptically Distributed Lozenge Tilings of a Hexagon Betea, D. |
| title | Elliptically Distributed Lozenge Tilings of a Hexagon |
| title_full | Elliptically Distributed Lozenge Tilings of a Hexagon |
| title_fullStr | Elliptically Distributed Lozenge Tilings of a Hexagon |
| title_full_unstemmed | Elliptically Distributed Lozenge Tilings of a Hexagon |
| title_short | Elliptically Distributed Lozenge Tilings of a Hexagon |
| title_sort | elliptically distributed lozenge tilings of a hexagon |
| url | https://nasplib.isofts.kiev.ua/handle/123456789/209540 |
| work_keys_str_mv | AT betead ellipticallydistributedlozengetilingsofahexagon |