GUE corners limit of q-distributed lozenge tilings

We study asymptotics of $q$-distributed random lozenge tilings of sawtooth domains (equivalently, of random interlacing integer arrays with fixed top row). Under the distribution we consider each tiling is weighted proportionally to $q^{\mathsf{vol}}$, where $\mathsf{vol}$ is the volume under the corresponding 3D stepped surface. We prove the following Interlacing Central Limit Theorem: as $q\rightarrow1$, the domain gets large, and the fixed top row approximates a given nonrandom profile, the vertical lozenges are distributed as the eigenvalues of a GUE random matrix and of its successive principal corners. Our results extend the GUE corners asymptotics for tilings of bounded polygonal domains previously known in the uniform (i.e., $q=1$) case. Even though $q$ goes to $1$, the presence of the $q$-weighting affects non-universal constants in our Central Limit Theorem.


Introduction and main results
1.1. q-distributed lozenge tilings. We begin with defining our main object, a probability distribution P N q,ν on interlacing integer arrays of depth N with fixed top row ν = (ν 1 ≥ . . . ≥ ν N ), ν i ∈ Z. Here q > 0 is a parameter that will eventually tend to 1. By an interlacing array of depth N we mean a collection λ = λ k j ∈ Z : k = 1, . . . , N, j = 1, . . . , k satisfying the interlacing constraints λ k j ≤ λ k−1 j−1 ≤ λ k j−1 for all k, j. Interlacing integer arrays originated as index sets of basis vectors in irreducible representations of unitary groups, and consequently they are sometimes referred to as Gelfand-Tsetlin schemes or patterns. We will identify interlacing arrays with configurations of lozenges of three types as shown in Figure 1.
Define |λ k | := λ k 1 + . . . + λ k k , and note that this number is not necessarily nonnegative. We will stick to the convention that the volume of a configuration λ is  To interpret the numbers λ k j as distinct particles, we shift them to λ k j +k−j, and for each k place {λ k j +k−j} 1≤j≤k onto the kth level in the coordinate system as shown. Then we surround each particle by a vertical lozenge, and complete the tiling with lozenges of two other types to get a lozenge tiling of the corresponding sawtooth domain. For example, the top row here is a configuration (12, 10, 9, 6, 3, −1), so λ 6 = ν = (7, 6, 6, 4, 2, −1). This tiling has vol(λ) = 56.
One can see that for fixed N and λ N , up to an additive constant depending on λ N , vol(λ) can be interpreted as a volume behind (i.e., to the left of) the 3-dimensional surface as in Figure 1 under the agreement that the normal vector to the vertical lozenges points to the right.
Let P N q,ν be the probability distribution on the finite set of interlacing arrays λ of depth N and with fixed top row λ N = ν defined as where Z N q,ν is the normalization constant. We have Z N q,ν = s ν (1, q, . . . , q N −1 ) = where s ν is the Schur polynomial (see, e.g., [Pet14b] for details).
1.2. Asymptotic regime. Let us now describe the limit regime we consider. Define x ∈ [0, 1], see Figure 2. This is a weakly decreasing function.
Without loss of generality we can and will assume that f(1) = 0.
Remark 1.1. The simplest nontrivial case of the function f is where 0 < b < 1 and a > 0. Then the measures P N q,ν are supported by lozenge tilings of growing lattice hexagons, and after the rescaling the lattice hexagons approximate the hexagon with sides a, b, 1 − b, a, b, 1 − b, cf. In general, if f(s) is piecewise constant then the measures P N q,ν are supported on lozenge tilings of growing lattice polygons with fixed number of sides each of which grows linearly in N .
In Section 1.5 below we discuss other regimes of q 1 besides the one of Assumption 2. (1.2) We will denote the distribution of L K by GUE K (σ 2 ), and this will be referred to as the GUE eigenvalue distribution. Along with the eigenvalues L K of H, consider eigenvalues L r = (L r 1 ≥ . . . ≥ L r r ) ∈ R r of the principal r-corner [H ij ] r i,j=1 of H for each r = 1, . . . , K. These eigenvalues satisfy the interlacing constraints L r j ≤ L r−1 j−1 ≤ L r j−1 , r = 1, . . . , K, j = 2, . . . , r.
The marginal distribution of each L r is the same as (1.2) (with K replaced by r). Note that L 1 1 is simply a normal random variable with mean 0 and variance σ 2 . The joint distribution of the whole interlacing array will be referred to as the GUE corners process. We will denote it by GUE K×K (σ 2 ).
The GUE corners process (earlier also called the GUE minors process) is a classical object of random matrix theory. It appeared in the context of interacting particle systems in [Bar01] and also in connection with random domino and lozenge tilings in [Joh05], [OR06], and [JN06]. The latter three works contain formulas for the correlation kernel of the GUE corners process (which is a determinantal point process), but we do not need the kernel in the present paper. See also Section 1.6 below for an exposition of previous results on GUE corners limits in random lozenge tilings. as N → +∞: where λ r j = λ r j (N ) denote the locations of the vertical lozenges on the first K levels under the distribution P N q,ν . We call Theorem 1.3 on Gaussian asymptotics of vertical lozenges at finitely many bottommost rows of the tiling the Interlacing Central Limit Theorem (ICLT). The proof of Theorem 1.3 occupies Sections 2 to 5. Let us make two remarks in connection with this theorem.
1.4.1. Frozen boundary. Besides Theorem 1.3 describing the asymptotic behavior of finitely many rows of the interlacing array adjacent to the bottom boundary (cf. Figure 1) there are several other interesting limit regimes in various ensembles of random lozenge tilings. Let us only discuss the connection of ICLT and the law of large numbers established in [CLP98], [CKP01], and [KO07] both in the uniform and in the q vol -weighted case. This law of large numbers implies the existence of a frozen boundary curve separating the liquid phase (in which asymptotically all types of lozenges are present) from the frozen regions (consisting of lozenges of only one type), see Figure 4. The ICLT (as in Theorem 1.3) arises at a point where two frozen regions adjacent to the bottom boundary meet. Thus, the global location u (1.3) corresponds to the point where the frozen boundary curve is tangent to the bottom boundary.
To the best knowledge of the authors, the frozen boundary for q vol -weighted lozenge tilings is explicitly known only in the case of unbounded plane partitions [OR03], [OR07] and for the hexagon [BGR10]. 1 One can check that the tangent point of the frozen boundary from [BGR10] and the boundary of the hexagon coincides with u.
At generic points of the frozen boundary the locations of the vertical lozenges form a determinantal point process with the extended Airy kernel in the scaling limit [OR07], [DM17]. 1 In fact, this paper considers a more general weighting of lozenge tilings related to the q-Racah univariate orthogonal polynomials. In the uniform (q = 1) case the frozen boundary is known explicitly for a much wider family of boundary conditions, e.g., see [Pet14b], [BG15], [DM15].
We also note that [KO07] implies that for a bounded polygon (both in the uniform and in the q vol -weighted case) the frozen boundary solves an algebraic equation which can be found approximately via numeric homotopy. By explicitly knowing the frozen boundary we mean an explicit solution to such an equation in the form of, e.g., a parametrization like in [Pet14b] for sawtooth domains in the uniform case.

1.4.2.
Extension to all γ ∈ R. Sending γ → 0 in (1.3) and (1.4) recovers the parameters in the ICLT for the uniform case [GP15] (see also [Nov15]): Proposition 1.4 (Uniform case γ = 0). We have Proof. For shorter notation, letḡ denote the integral of g(s) over s ∈ [0, 1]. The first claim follows from a Taylor expansion: . This implies the first claim, and we will also need the next term in u to address the second claim. Namely, we have which implies the second claim.
The next statement allows us to extend Theorem 1.3 to the case q > 1, q 1, that is, to γ < 0. Observe that by reflecting the interlacing array with respect to zero and shifting it one can turn the measure P N q,ν , where q = e −γ/N < 1, γ > 0, into the measure P N 1/q,ν , whereν j = ν 1 − ν N +1−j for j = 1, . . . , N . If ν corresponds to a function f(s) as in Assumption 1, thenν corresponds to the functionf(s) = f(0) − f(1 − s). Clearly, if Theorem 1.3 holds for a sequence of measures of the form P N q,ν , then it also holds for the measures P N 1/q,ν , and the next proposition matches the limiting parameters: Remark 1.6. Propositions 1.4 and 1.5 imply that Theorem 1.3 holds in the regime q = e −γ/N for any fixed γ ∈ R. For simplicity in the rest of the paper we will continue to assume that γ is positive.
1.5. Slower q 1 regime. Theorem 1.3 establishes, in the regime q = e −γ/N 1 with γ > 0 fixed, a universal ICLT regardless of the (nonconstant) function f(s) governing the asymptotic behavior of the top row. Let us argue that when q = q(N ) converges to 1 slower, the presence of the GUE corners limit is not guaranteed and depends on the nature of the function f(s).
The law of large numbers [CLP98], [CKP01], [KO07] implies the existence of a nonempty liquid region at scale N for any fixed γ in the regime q = e −γ/N . The slower q 1 regime leads to a completely frozen situation at scale N , i.e., there is no liquid region, and the random tiling at this scale looks like the unique tiling of the minimal volume (see Figure 5 for examples of the latter). However, one can still ask if there is an ICLT at some fluctuation scale N . The slower convergence q 1 is equivalent to taking γ = γ(N ) depending on N such that γ(N ) → +∞ but γ(N ) N . Let us first formally send γ → +∞ in the limiting quantities (1.3), (1.4) in the following two cases of the top row: The function f 1 corresponds to the hexagon boundary conditions, and f 2 -to a top row having limiting density 1 2 (the latter can be modeled by taking ν j = N −j). Straightforward computations show that as γ → +∞: This suggests that if there could be an ICLT in the slower q 1 regime then its scale must depend on the form of the function f(s).
In the hexagon case there are, however, strong indications that an ICLT at a finite distance from the bottom row is not possible. Namely, in the slower q 1 regime the nontrivial asymptotic behavior (at a scale N ) occurs in a neighborhood of the unique three-lozenge combination in the minimal volume tiling (cf. Figure 5, left). Observe that the distance of from the bottom row is of order N . The configuration around should asymptotically coincide with the q vol -weighted random plane partition without boundary conditions (the latter model was studied in, e.g., [FS03], [OR03], [OR07]). The behavior of the bottommost lozenge of the type can then be guessed from the corresponding results on random plane partitions [Mut06], [VY06]. In particular, the fluctuations of the row number containing the bottommost lozenge should grow to infinity. Thus, even if we tune the speed at which q goes to 1 so that the plane partition behavior is visible (in any sense) at a finite distance from the bottom, still with positive probability the vertical lozenges at the bottommost rows are packed to the left (as in the minimal volume tiling). This suggests that for the hexagon (and similarly for any polygonal boundary conditions) the ICLT at a finite distance from the bottom does not occur. The situation in Figure 5, right, is different in an essential way. Namely, for this top row the interlacing constraints allow much more configurations with volume close to minimal, which suggests that the fluctuations of the vertical lozenges at finite distance from the bottom grow to infinity. This leaves open a possibility to have an ICLT at some fluctuation scale √ N . In terms of the weakly decreasing function f the difference between the situations in the left and the right pictures in Figure 5 is in the value of f (1) corresponding to the density of the vertical lozenges at the left end of the top row. This leads to the following conjecture: Conjecture 1.7. Let Assumption 1 hold, and replace Assumption 2 by N . If f (1) < 0, then an analogue of Theorem 1.3 holds under a suitable normalization determined by u(N ) and σ 2 (N ) with lim N →∞ u(N ) = lim N →∞ σ 2 (N ) = 0: 1.6. Background and methods. Let us discuss previous work on ICLT in random lozenge tilings and how our approach differs from the existing ones. The ICLT has been proven for various ensembles of random domino and lozenge tilings starting from [Joh05], [OR06], [JN06], [JN07]. In particular, [JN06] essentially establishes an ICLT for uniformly random lozenge tilings of the hexagon (see also [Nor09] for more details and for a discussion of ICLT for domino tilings), and [OR06] considers an ensemble of q vol -weighted random lozenge tilings of certain infinite regions. The latter paper also informally argues that the ICLT should hold for Gibbs (in particular, uniform) ensembles of random lozenge tilings at points where two frozen regions meet (such points are sometimes called turning points), see also [BD11] for an ICLT for a Gibbs measure on tilings of an infinite region. Here by a Gibbs property we mean that for each K ≥ 1 the distribution of levels 1, . . . , K − 1 of an interlacing array (discrete or continuous) conditioned on fixed configuration at level K is uniform over all configurations subject to interlacing. In non-Gibbs situation the limit at turning points might not be the universal GUE corners process but rather its certain splitting, cf. [Mkr14].
The ICLT results cited in the previous paragraph rely on the presence of explicit determinantal correlation kernels of the pre-limit models suitable for asymptotic analysis. Such kernels can be extracted from the formalism of Schur processes [OR03], [OR07], [BR05] for q vol -weighted or periodically weighted random tilings of infinite regions; via connections with orthogonal polynomials [JN06], [BGR10] for q vol -weighted tilings of the hexagon; or by an application of Eynard-Mehta type results [Pet14a], [DM15] for uniformly random tilings of general sawtooth domains. The latter method also produces a determinantal correlation kernel for the q vol -weighted tilings of general sawtooth domains [Pet14a], but it has a more complicated structure than in the uniform case, which so far has presented an impediment to its asymptotic analysis.
Recently other methods of proving ICLT (and other asymptotic results) for uniformly random lozenge tilings of general sawtooth domains as in Assumption 1 were developed in [GP15] and [Nov15]. The approach in the former paper is based on Schur generating functions and allows to also consider tilings with free boundary [Pan15] and to establish an ICLT for alternating sign matrices [Gor14]. The methods of [Nov15] are based on combinatorics of Hurwitz numbers. These more recent methods are related to free probability and go beyond the ICLT leading to limit shape results [BG15], [CNŚ16] (both approaches) and asymptotics of global fluctuations [BG16] (the Schur generating functions approach).
Our approach to the ICLT for q vol -weighted lozenge tilings of general sawtooth domains (Theorem 1.3) is closer to methods relying on determinantal kernels. However, while the determinantal structure of the measure q vol obtained in [Pet14a] is quite complicated, the distribution of the vertical lozenges at each K-th level of the interlacing array can be written as a much simpler K × K determinant [Pet14b] with matrix elements expressed as single contour integrals. This is the starting point of our work, and the asymptotic analysis of the single contour integrals yields our ICLT. 1.7. Outline. In Section 2 we recall the K × K determinantal formula from [Pet14b] for the distribution of the vertical lozenges at the Kth level of the interlacing array. In Section 3 we perform a preliminary asymptotic computation at the first (K = 1) level to determine the global location u (1.3) of the vertical lozenges. In Section 4 we analyze the critical points and the steepest descent contour of the integral entering the distribution of the vertical lozenges at the Kth level. In Section 5 we compute the asymptotics of the distribution of the vertical lozenges at any Kth level (which are dominated by the behavior in a neighborhood of a critical point), and complete the proof of Theorem 1.3.
1.8. Acknowledgments. We are grateful to conversations with Vadim Gorin, Jonathan Novak, and Greta Panova. SM was partially supported by the Simons Foundation Collaboration Grant No. 422190.
2. Projection of P N q,ν onto the Kth row The starting point of our analysis is a formula for the projection of our measure P N q,ν onto any fixed row K, i.e., a formula for the distribution of the particles λ K j , j = 1, . . . , K, under P N q,ν for any fixed 1 ≤ K < N . Throughout the paper we use the notation (with (a; q) 0 = 1), for the q-Pochhammer symbol.
Let us rewrite the functions A i (x) in a form more convenient for asymptotic analysis: . . , N } and the contour C is a positively oriented circle crossing the real line to the left of 0 and between q and 1.
Proof. Change the variables in (2.2) as z = q x w −1 . Then w −1 is integrated over a positively oriented circle encircling 1, q, q 2 , . . . and not q −1 , q −2 , . . .. That is, the integration contour for w −1 encircles the segment [0, 1]. Therefore, we can choose the contour C for w as in the claim, and note that there is an additional change of sign coming from the orientation of the contour. Thus, Rewriting the products under the integral, we get (2.3).

Global location u from first level asymptotics
According to the desired claim of Theorem 1.3, under the measures P N q,ν and for any fixed K the random locations of the vertical lozenges on the Kth row should behave as are the rescaled locations which we will later identify with the eigenvalues of the K × K GUE random matrix (see Definition 1.2). In this section we perform a preliminary computation for the simplest case K = 1 which will help determine the correct global asymptotic location u.
When K = 1, the projection formula of Section 2 becomes . (3.1) Since x should depend on N as x = uN + ξ √ N , we have from Assumption 2: Therefore, the exponential in N terms coming from the asymptotics of the contour integral should compensate e −N γ(u−1/2) . Let us now consider the exponential in N behavior of the terms under the integral in (3.1). The numerator behaves as 3) where in the second equality we used the approximation of the integral by the corresponding Riemann sum (which has error O(1/N )). Here and below we take log to be the standard branch having the cut along the negative real axis. Later in the course of our analysis it will become evident that the integral over the part of the contour C around the branch cut (i.e., e −γs −w < 0) is asymptotically negligible. This justifies our choice of the branch.
For the denominator in (3.1) we have Here and below for simplicity we omit the notation of the integer part, and write, e.g., The presence of the integer parts does not affect the asymptotics. The asymptotic behavior of such integrals can be analyzed using the steepest descent (also called Laplace) method. Namely, one finds a critical point w cr (in our case this indeed will be a single critical point) and a deformation of the integration contour C so that the deformed contour C passes through w cr and is a steepest descent contour for S(w) in the sense that Re (S(w) − S(w cr )) < 0, w ∈ C \ {w cr }. Then so the exponential in N terms in the integral have the form e N S(wcr) . Comparing this to (3.2) we see that one should find u and w cr such that The next lemma provides one such pair (u, w cr ) which will turn out to be the right one for asymptotics. . Throughout the rest of the paper we will fix the value of u as in (3.8) and will call it the global location.

Critical points and contour deformation
In this section we study in detail the critical points of the function S(w) as well as the contour plot of Re S(w). One of the key ingredients is an approximation of the weakly decreasing function f(x), x ∈ [0, 1], of Assumption 1 by piecewise constant weakly decreasing functions f m (x) (we also require that each f m has a finite range). The approximating functions f m (x) do not depend on N , and the whole approximation is done independently of the main limit N → +∞. Note that the case of f(x) piecewise constant corresponds to considering random tilings of polygons with fixed number of proportionally growing sides. Recall that we always assume that f(1) = 0. Note that {s j } and {α j } will depend on m, but we suppress this dependence in the notation.

Proof. A straightforward computation.
The expression for S m (w) given in (4.4) is valid with standard branches of the logarithms for w sufficiently close to 0. Instead of looking at (4.4) directly we will focus on the algebraic equation (4.5) implied by S m (w) = 0 in part because (4.5) does not involve choosing the branches.
Clearly, w = 0 satisfies (4.5) regardless of the value of u (but depending on u the solution w = 0 might not lead to a critical point of S m (w)). Let us examine other possible critical points: Proposition 4.2. All solutions to (4.5) are real. Moreover, besides w = 0, all solutions to (4.5) belong to the segment e −γ(u+1) , e −γ(u−α 1 ) .
Proof of Proposition 4.2. Let N m (w) be the numerator and D m (w) the denominator of (4.5). Both are polynomials in w of degree m + 1. We will count the number of real solutions to the equation N m (w) = D m (w).
First, note that the ratio of the leading coefficients in N m (w) and D m (w) is e γ m j=1 e γ(u−α j +s j−1 ) m j=1 e γ(u−α j +s j ) = e γ(1+s 0 −sm) = 1 by (4.1). Thus, N m (w) − D m (w) is a polynomial of degree at most m, so the equation N m (w) = D m (w) has at most m solutions.
Observe that in our piecewise constant case the expression (3.8) for u simplifies to .
(4.6) This readily implies that N m (w) and D m (w) have the same linear in w terms. Since, in addition, the constant terms are both 1, we see that w = 0 is a double solution to N m (w) = D m (w). Consider separately the roots of N m (w) and D m (w), and denote n j := e −γ(u−α j +s j−1 ) and d j := e −γ(u−α j +s j ) , j = 1, . . . , m. Note that N m (w) in addition has the root w = e −γ , and D m (w) in addition has the root w = 1. From (4.1) we see that the other roots interlace as d m < n m < d m−1 < n m−1 < . . . < n 2 < d 1 < n 1 .

Moreover, by (4.1) and Lemma 3.2 we have
Thus, all roots of N m (w) and D m (w) (and not just the n j 's and the d j 's) "almost interlace". Namely, there is exactly one consecutive pair of roots of N m (w), one of them is e −γ , surrounded by a pair of roots of D m (w), and similarly there is only one consecutive pair of roots of D m (w), one of them is 1, surrounded by a pair of roots of N m (w). See Figure 6. Figure 6. Positioning of the roots of N m (w) (hollow) and D m (w) (solid) for m = 10, as well as schematic plots of both polynomials. One readily sees that these plots must intersect at least m − 2 = 8 times.
This almost interlacing implies that there are at least m − 2 real solutions to N m (w) = D m (w) lying from e −γ(u+1) to e −γ(u−α 1 ) . Adding to this the double solution w = 0 we see that the claim holds. Proof. If f(x) = f m (x) is piecewise constant, this follows from Proposition 4.2 and the fact that the equation S m (w) = 0 implies (4.5).
If f(x) is arbitrary, let us approximate it by a sequence of piecewise constant decreasing functions f m (x), m → +∞, such that the corresponding derivatives S m (w) converge to S (w) uniformly in w belonging to compact subsets of C \ R. If there is a nonreal critical point of S(w), then by Rouche's theorem one of the functions S m (w) also would have a nonreal critical point, which is impossible by Proposition 4.2.
To control locations of real critical points note that S (w) is well defined (by differentiating under the integrals in (3.5)) for real w outside [e −γ(u+1) , e −γ(u−f(0)) ], and so can be uniformly approximated there by the S m (w)'s. This completes the proof. 4.2. Steepest descent contour. After describing the critical points of the function S(w) let us focus on the steepest descent contour C := {w ∈ C : Im S(w) = 0} passing through the critical point w cr = 0 and along which Re S(w) attains its only maximum at w cr (recall that S(0) = γ(u − 1/2) is real). The goal of this subsection is to justify that the contour C and the region where Re S(w) < S(0) look as in Figure 7. Our first observation is that the direction at which the steepest descent contour C passes through zero is vertical: We will prove that S m (0) > 0 by induction on m. When m = 1, using (4.1) we have S 1 (0) = 0. For general m, S m (0) depends on α 1 in the following way: 1 − e γ = 1 + e γ − (e γ − 1) where C 1,2 and D 1,2 do not depend on α 1 . Differentiating this with respect to A, we obtain The denominator is equal to the cube of which is positive by (4.1). Therefore, the sign of the left-hand side of (4.7) is the same as the sign of e −γα j (e γs j − e γs j−1 ) e −γα j (e γs j + e γs j−1 ) − e −γα 1 (e γs 1 + 1) .
From (4.1) we see that each summand above is positive. Therefore, S m (0) decreases in A = e −γα 1 , hence it increases in α 1 . Since our goal is to show that S m (0) > 0 and we have α 1 > α 2 , it's enough to prove S m (0) ≥ 0 when α 1 = α 2 . However, the latter corresponds to reducing m by one. Hence, by induction, we are done with the case of a piecewise constant f. It follows that if m is large, all nonzero critical points of S m are larger than 1 2 e −γ(f(0)+1) , so S m cannot have nonzero critical points in the 1 2 e −γ(f(0)+1) neighborhood of 0 as m → ∞. This is a contradiction which completes the proof of the lemma.
The steepest descent contour C passing through the origin in the vertical direction could in principle escape to infinity. Let us show that this is not the case: Therefore, the steepest descent contour C emanating from 0 must intersect the real line again. The next lemma specifies where this point of intersection is: Lemma 4.6. The contour C = {w : Im S(w) = 0} passing through 0 intersects the real line again at a point between e −γ and 1.
Proof. We will prove the lemma by looking at the behavior of Im S m (x + i ) for x ∈ R and small > 0. First, observe that Im log(x + i ) = arg(x + i ) ∼ π1 x<0 .
We see that there must be at least one point in (e −γ , 1) for which Im S m (x + i ) = 0.
There are two cases. First, assume that C intersects R at a critical point of S m . Then the limit of Im S m (x + i ) as 0 is zero in a neighborhood (within R) of this point of intersection. From Proposition 4.2 we know that there cannot be such an intersection of C with (e −γ(u−fm(0)) , +∞) or with (−∞, e −γ(u+1) ) (besides the original critical point w cr = 0). From the above description of monotonicity of Im S m (x+i ) it readily follows that a neighborhood within (e −γ(u+1) , e −γ(u−fm(0)) ) where Im S m (x) = 0 should be inside (e −γ , 1), and therefore the point of intersection of C with R is also inside (e −γ , 1). In the other case C intersects R at a point which is not a critical point of S m , and this corresponds to a transversal intersection of the graph of Im S m (x + i ) with the x-axis (this case is shown in Figure 8). The monotonicity of Im S m (x + i ) implies that such a transversal intersection can occur only inside (e −γ , 1), which completes the proof in the case of a piecewise constant function f = f m .
In the case of an arbitrary f we argue similarly to the proof of Proposition 4.3 and approximate f by piecewise constant functions f m as m → +∞. The segments on which Im S(x + i ) is weakly monotone for small > 0 are the same as in the piecewise constant case, and so there is x ∈ (e −γ , 1) where Im S(x + i ) = 0. This x corresponds to the intersection of the contour C with R.
Lemmas 4.4 to 4.6 imply that the steepest descent contour C looks as in Figure 7. Note that the numerator of the integrand in A i of Proposition 2.2 has zeros at w = q r , r = 1, . . . , N −K−1, while the denominator in this integrand does not have double zeros. This implies that the integrand in A i has no poles between q N −K−1 = e −γ+γ K+1 N ∼ e −γ and q = e − γ N ∼ 1. Moreover, the integrand in A i also has no poles on the negative real line. Therefore, the deformation of the integration contour C in formula (2.3) of Proposition 2.2 to the steepest descent contour C does not give rise to any residues. In the sequel we will thus assume that the integration contour in (2.3) is already C.

4.3.
First level CLT. Before proving Theorem 1.3 in full generality in Section 5 below, let us use the results of Sections 4.1 and 4.2 to establish the simpler particular case K = 1 of this theorem. That is, we are interested in establishing a Central Limit Theorem for the probability measure on Z given in (3.1), where the integration contour is C, the steepest descent one. In the limit at N → +∞ the contribution to the integral outside a neighborhood of w cr = 0 of size, say, N − 1 6 , is negligible. Indeed, outside this neighborhood N Re (S(w) − S(0)) is bounded from above by a negative constant times N 2 3 , which dominates the O( √ N ) terms in the exponent under the integral (in, e.g., (3.6)). Therefore, we can focus on the behavior of the integral over the part of the contour C in a N − 1 6 neighborhood of w cr = 0. In this neighborhood of zero we change the variables as w = −iN − 1 2 t, t ∈ R. The minus sign comes from the orientation of the contour C. With this change of variables and with scaling x = uN + ξN 1 2 , ξ ∈ R, (3.1) becomes . (4.9) Here we took the factors out of the products in order to Taylor expand in the small variable N − 1 2 t for fixed t ∈ R (note that r/N and ν r /N stay bounded as N → +∞): In the denominator we also need to keep track of the term e γξN −1/2 = 1 + N − 1 2 γξ + O(N −1 ), and so we have Every sum in the exponent in the previous two formulas is a Riemann sum for the corresponding integral, and combining all exponents in (4.9) together we get an exponent of (which immediately follows from (3.5)). Therefore, we can continue (4.9) and evaluate the Gaussian integral which leads to a Gaussian density in the variable ξ ∈ R (the factor N − 1 2 corresponds to the rescaling of the space from discrete to continuous). Thus, we have established the particular case K = 1 of Theorem 1.3: Proposition 4.7 (First level CLT). Under Assumptions 1 and 2, as N → +∞, we have the following convergence in distribution of the location λ 1 1 = λ 1 1 (N ) of the vertical lozenge on the first level: Here the global location u and the limiting variance σ 2 are given by (1.3) and (1.4), respectively.

Completing the proof of Theorem 1.3
In this section we finish the proof of our main result, Theorem 1.3. We first consider the asymptotics of the distribution at each level K = 1, 2, . . . (using the formulas of Theorem 2.1 and Proposition 2.2 together with our results from Sections 3 and 4), and then explain how one obtains a joint ICLT for all K(K + 1)/2 lozenges λ = λ k j : k = 1, . . . , K, j = 1, . . . , k .
5.1. CLT at level K. Let us begin with the behavior of the prefactor in (2.1): Lemma 5.1. Fix K ≥ 1. Let κ = (κ 1 ≥ . . . ≥ κ K ) depend on N as κ j = uN + ξ j N 1 2 , where u is given by (3.8), and ξ 1 ≥ . . . ≥ ξ K are fixed real numbers. Then as N → +∞ the prefactor in front of the determinant in (2.1) behaves as Proof. We have by (1.1): and, moreover, |κ| = uKN + N 1 2 K i=1 ξ i . This implies the claim. Let us now address the asymptotic behavior of the individual entries A i (κ j − j) of the K × K determinant in (2.1). The analysis is similar to Section 4.3, but some care is required to pass from the asymptotics of the individual entries to the asymptotics of the determinant. This dictates our formulation of Lemma 5.2 below.
Lemma 5.2. Fix K ≥ 1 and let x = uN + ξN 1 2 . Then A i (x) given by (2.3) has the following asymptotics as N → +∞: where the remainder o(1) in front does not depend on i while each i,l,ξ (N ) = o(1) may depend on i, l, and ξ.
Proof. Arguing as in Section 4.3, we can restrict the integration to a small neighborhood of w cr = 0, and change the variables as w = −itN − 1 2 , t ∈ R. Then we have The prefactor scales as 2) The part of the integrand not depending on i behaves as where we employed Taylor expansions and approximation of sums by the corresponding Riemann integrals similarly to Section 4.3 (the difference by finitely many factors in the numerator does not affect the approximation). The product over I i,N,K contains only K − 1 factors, which is finite and thus cannot be approximated by the exponent of an integral. We have for r ≤ i − 1: Similarly, for r ≥ N − K + i + 1: Therefore, we have Lemma 5.2 implies that the determinant det[A i (κ j − j)] of Theorem 2.1 is, up to a prefactor, asymptotically equal to a product of two simpler determinants, one of the c i,l 's, and the other one involving the G l (ξ j )'s (where κ j = uN + ξ j N 1 2 ), as long as both of these simpler determinants are nonzero. Indeed, under the latter assumption in the expansion of det K−1 l=0 c i,l G l (ξ j )(1 + i,l,ξ j (N )) as a sum over permutations the asymptotically dominating terms do not contain any of the i,l,ξ j 's, and lead to a product of two nonzero determinants. Let us now compute these two simpler determinants separately: Proof. Introduce variables y 1 , . . . , y K . Then det 1≤i≤K, 0≤l≤K−1 (1 + y j ) i−1 (1 + e γ y j ) K−i by the very definition of c i,l . The latter determinant can be rewritten as (1 + y j ) i−1 (1 + e γ y j ) K−i = K j=1 (1 + e γ y j ) K−1 K det i,j=1 1 + y j 1 + e γ y j i−1 = K j=1 (1 + e γ y j ) K−1 i<j 1 + y i 1 + e γ y i − 1 + y j 1 + e γ y j = (1 − e γ ) K(K−1) 2 i>j (y i − y j ), as desired.
Lemmas 5.3 and 5.4 show that both determinants whose product enters the determinant of A i (κ j − j) in (2.1) are indeed nonzero as long as the ξ j 's are distinct. When some of the ξ j 's are equal to each other, the GUE eigenvalue density also vanishes, and one can readily check that the scaling limit of (2.1) is zero, as it should be. Therefore, we arrive at the following result: Proposition 5.5. (Level K CLT) Under Assumptions 1 and 2, as N → +∞, for the locations λ K j = λ K j (N ), j = 1, . . . , K, of the vertical lozenges at level K we have the following convergence in distribution on R K : Here u and σ 2 are given by (1.3) and (1.4), respectively, and GUE K is the eigenvalue distribution of the K × K Gaussian Unitary Ensemble (cf. Definition 1.2).
Because the prefactor N − K 2 corresponds to the rescaling of the space from discrete to continuous, we see that the right-hand side of (5.7) leads to the eigenvalue density GUE K (σ 2 ). This completes the proof. 5.2. Gibbs property. The last step required to complete the proof of Theorem 1.3 is to show that the joint distribution of all K(K + 1)/2 locations λ = {λ k j : k = 1, . . . , K, j = 1, . . . , k} of the vertical lozenges at levels 1, . . . , K is approximated by the GUE corners process L = {L k j } ∼ GUE K×K (σ 2 ) (see Definition 1.2 for the latter). We employ the Gibbs property of GUE K×K (σ 2 ): conditioned on the fixed Kth row L K , the distribution of the rest of the array {L k j : 1 ≤ k ≤ K − 1, 1 ≤ j ≤ k} is uniform over the polytope 4 described by the interlacing conditions L k j ≤ L k−1 j−1 ≤ L k j−1 for all k = 1, . . . , K and j = 2, . . . , k. At the same time, the pre-limit joint distribution of the locations of the vertical lozenges λ = {λ k j } satisfies a q-deformation of the Gibbs property. Under the latter, the uniform conditional probabilities of the rows λ 1 , . . . , λ K−1 conditioned on the fixed row λ K are replaced by the probabilities proportional to q vol . Let us explain why in our limit regime the q-Gibbs property leads to the Gibbs one.
In Proposition 5.5 we showed that the Kth row λ K of the discrete interlacing array corresponding to the vertical lozenges in our q vol -distributed lozenge tiling converges after rescaling to the Kth row L K of the GUE corners process. This rescaling in particular implies that the distances between the λ K j 's are of order √ N . Therefore, conditioned on the fixed configuration λ K on the Kth row with λ K 1 −λ K K of order √ N , the difference between the maximal and the minimal volume of the tiling is also of order √ N . Thus, since q = e −γ/N , we have q max(vol) = (1 + o(1))q min(vol) . This means that the conditional distribution of the lower rows λ 1 , . . . , λ K−1 becomes uniform (up to the interlacing constraints) in the N → +∞ limit. This immediately leads to the ICLT and thus completes the proof of Theorem 1.3.