A Periodic Hexagon Tiling Model and Non-Hermitian Orthogonal Polynomials

We study a one-parameter family of probability measures on lozenge tilings of large regular hexagons that interpolates between the uniform measure on all possible tilings and a particular fully frozen tiling. The description of the asymptotic behavior can be separated into two regimes: the low and the high temperature regime. Our main results are the computations of the disordered regions in both regimes and the limiting densities of the different lozenges there. For low temperatures, the disordered region consists of two disjoint ellipses. In the high temperature regime the two ellipses merge into a single simply connected region. At the transition from the low to the high temperature a tacnode appears. The key to our asymptotic study is a recent approach introduced by Duits and Kuijlaars providing a double integral representation for the correlation kernel. One of the factors in the integrand is the Christoffel–Darboux kernel associated to polynomials that satisfy non-Hermitian orthogonality relations with respect to a complex-valued weight on a contour in the complex plane. We compute the asymptotic behavior of these orthogonal polynomials and the Christoffel–Darboux kernel by means of a Riemann–Hilbert analysis. After substituting the resulting asymptotic formulas into the double integral we prove our main results by classical steepest descent arguments.


Introduction
We study random lozenge tilings of large regular hexagons. We place the regular hexagon so that it has corners at (0, 0), (0, N ), (N , 2N ), (2N , 2N ), (2N , N ) and (N , 0) and consider tilings of the hexagon with the following three types of lozenges Type I Type I I and Type I I I , see also Fig. 1. The vertices of the lozenges are on the integer lattice and the vertical and horizontal edges have unit length. There are numerous ways of defining a probability measures on all possible tilings of the hexagon. In this paper, we will be interested in the case in which the probability of a tiling T is given by where W is a weight function on all possible tilings defined by for some fixed α ∈ (0, 1]. Note that if α = 1 all tilings occur with the same probability and the probability measure reduces to the uniform measure on all possible tilings. We exclude α = 0. In the limit α ↓ 0, there is only one possible tiling, see e.g. Fig. 3 below, and there is no randomness. The main results in this paper concern the asymptotic behavior of the random tilings as the size of the hexagon grows large, i.e., as N → ∞, and how this asymptotic behavior depends on the parameter α. Random tilings of planar domains have been extensively studied in the past decades and we refer to [6,[24][25][26]45,[50][51][52] for important early references, and to [13,46,48] for excellent introductions to the topic. When the domains are large, the statistical properties of the tilings are expected to be described by universal limiting processes. In various N N N Fig. 1. The hexagon (left) and an example of a tiling (right) of the hexagon by lozenges special classes, and especially in case the random measure is a determinantal point process, tools have been developed to compute the asymptotic behavior and verify the appearance of these universal processes. For instance, if the random measure is in the Schur class [62,64], then we have a double integral representation for the correlation kernel at our disposal to analyze the fine properties of the model. Random lozenge tilings of the hexagon are however typically not in the Schur class and asymptotic studies are often more complicated.
Although not being in the Schur class, the large N behavior of random lozenge tilings of the hexagon with the uniform measure (corresponding to α = 1 in our setup) has also been intensively studied by various authors. Based on a representation in terms of Hahn polynomials as found in [44] (see also [43]), the authors of [6] managed to perform a steepest descent analysis of the discrete Riemann-Hilbert (RH) problem for the Hahn polynomials and, consequently, describe the limiting disordered regions and the local universality laws. In [43] the local universality was obtained using methods developed in [16]. In a more general context, uniform lozenge tilings of more complicated domains were studied by means of double integral formulas [3,[35][36][37]66,67].
An important part of the recent literature on random tilings is concerned with proving the universality of the global fluctuations and the emergence of the Gaussian Free Field. For the uniform measure on all possible tilings of the hexagon there are now various techniques in the literature that prove this claim. In [67] the convergence of the global height fluctuations to the Gaussian Free Field was established using double integral formulas for the kernel. An alternative proof based on the recurrence coefficients of the Hahn polynomials was given in [33] extending the results on the fluctuations along vertical sections in [18]. Discrete loop equations can also be used [14] to compute the fluctuations along vertical sections. In [19,20], another approach is introduced using the notion of a Schur generating function. Each of these methods apply to their own general class of models and contain the uniform measure as a special case.
Measures on tilings of the (finite) hexagon that are not uniform are known to be difficult to analyze asymptotically and much less results are known. For instance, in [15] the authors introduced elliptic weights on the lozenge tilings, but a full asymptotic study of these models is still open. The situation 0 < α < 1, which is the topic of this paper, is a rather gentle way to break the uniform measure. Still, the above mentioned techniques do not apply. To study our model we will use a recently developed new approach [34] for studying determinantal point processes that are defined via products of minors of (scalar or block) Toeplitz minors. Although the original motivation of [34] was to analyze the socalled 2-periodic Aztec diamond (see also [7,22]), the methods apply to a much wider range of (tiling) models. The approach mainly consists of combining two important methods for asymptotic analysis: the classical steepest descent method for integrals and the Deift/Zhou steepest descent method for RH problems [27,29]. This opens up new possibilities for analyzing models that were thus far out of reach and the model studied in this paper is one such example.
It is possible to take the limit of our model in which the vertical sides of the hexagon tend to infinity (see, for example, [11] for an explanation that starts from the same setting as in the present paper). In that limit, our model is the same as a 2-periodic weighting of plane partitions against a linearly shaped back wall, as studied in [60] (see also [5] for a generalization to the setting of Macdonald processes). This model is then in the Schur class and thus double integral representations are available for asymptotic studies. It is important to note that the case of a finite hexagon does not only lead to technical challenges, but also more complicated phenomena occur. For instance, in our model a tacnode appears for α = 1/9.
In Fig. 2 we have plotted two sample tilings for large hexagons, one with 0 < α < 1 9 and the other with 1 9 < α < 1. We see that for 0 < α < 1 9 there appear two clouds in which the tiling shows randomness, while it is frozen outside. In the figure with 1 9 < α < 1, these two clouds seem to have merged. To understand why this phenomenon is happening, it is useful to view α as a temperature parameter. Indeed, after defining the energy of a tiling as we can write the weight of a tiling T as W (T ) = e (log α)E(T ) , and its probability as which is a Gibbs measure with inverse temperature β. Thus, T = − 1 log α may (and will) be viewed as the temperature parameter. The low temperature limit T ↓ 0 corresponds to α ↓ 0 and the high temperature limit T → ∞ to α ↑ 1.
For low temperatures, the number E(T ) is expected to be small. In fact, for T ↓ 0 the randomness disappears and the lozenge configurations freeze to the unique tiling with E(T ) = 0. This is the tiling that is shown in the left half of Fig. 3. It can be thought of as a staircase shaped wall where the floor and the ceiling only have tiles of type III. As the temperature increases, randomness starts appearing near the interfaces where the wall meets the ceiling and the floor. For T positive but small, we expect to observe two separate clouds that are far away from each other. When T increases further, the clouds meet and form one cloud. Eventually, as T → ∞, the model becomes the uniform measure on tilings and the cloud becomes the ellipse that is inscribed in the hexagon, as in the right part of Fig. 3. In other words, we expect that there is a critical point in the low to high temperature transition at which the topology of the disordered regime changes from being disconnected to being connected. As we will see, this transition indeed happens at α = 1 9 . We will therefore speak of 0 < α < 1 9 as the low temperature regime and of 1 9 < α ≤ 1 as the high temperature regime.
Our analysis follows a recent work [34]. The backbone of the approach in [34] is a connection to polynomials that satisfy an orthogonality relation (that could be matrix valued) on a contour in the complex plane. In the present paper we will be dealing with scalar orthogonality on a closed contour γ going once around the origin with counterclockwise orientation. Let p n be the monic polynomial of degree n such that 1 2πi γ p n (z)z j (z + 1) N (z + α) N z 2N dz = 0, j = 0, 1, . . . , n − 1. (1.2) It is important to note that (1.2) is an orthogonality condition with respect to a non-Hermitian bilinear form. It is therefore not evident that the polynomials p n are welldefined. We will prove that they are, provided that n ≤ 2N , see Proposition 5.1. The orthogonality (1.2) also changes with N , the size of the hexagon. It turns out that the random tilings naturally define a determinantal point process with a correlation kernel that can be expressed in terms of the polynomials p n . For the exact statement, we need to introduce a well-known correspondence between tilings of the hexagon and non-intersecting paths. For more background on determinantal point processes, random tilings and non-intersecting paths, we refer to [46].
We draw lines on two of the three types of lozenges as follows: and . The paths form a collection of non-intersecting paths π j : {0, . . . , 2N } → Z + 1 2 with initial points π j (0) = j + 1 2 and endpoints π j (2N ) = N + 1 2 + j for j = 0, . . . , N −1. It is well-known and easy to see that there is a one-to-one correspondence between tilings of the hexagon and non-intersecting up-right paths with these initial and end configurations.
The probability measure on the tilings defined in (1.1) induces a probability measure on such collections of non-intersecting paths. The Lindström-Gessel-Viennot lemma [41,55] tells us that the probability measure is proportional to where the T m are Z × Z matrices given by if m is even, and if m is odd. The probability (1.3) is a determinantal point process with a correlation kernel given by the Eynard-Metha formula [38].
In case the Z × Z matrices T m in (1.3) are (scalar or block) Toeplitz matrices, the paper [34] gives a double contour integral formula for the correlation kernel, which involves the (scalar or block) symbols of the Toeplitz matrices as well as a reproducing kernel for (scalar or matrix-valued) orthogonal polynomials, see also [8].
The matrices (1.4) and (1.5) are infinite Toeplitz matrices with only two non-zero diagonals. Their respective symbols are z + α and z + 1. Both Toeplitz matrices appear N times in the product (1.3) and this accounts for the orthogonality measure in (1.2). Then the general formula in [34] reduces to the following in the special situation of this paper. Proposition 1.1. Let α ∈ (0, 1] and let k ≥ 1 be an integer. Then for integers x 1 , . . . , x k , paths go through each of the points (x 1 , y 1 + 1 2 ), . . . , (x k , y k + 1 2 ) where the kernel K N is given by is the squared 'norm' of p n .
Proof. This is a special case of [34,Theorem 4.7], but for convenience of the reader we give more details on how to make the identification in the Appendix.
The above proposition is the starting point of our analysis. Clearly, to analyze the limiting behavior of the probabilities (2.27) it suffices to compute the asymptotic behavior of the kernel K N in (1.7) as N → ∞. To this end, we first compute the asymptotic behavior of the Christoffel-Daroux kernel R N corresponding to the orthogonal polynomials using Riemann-Hilbert techniques. After inserting the resulting asymptotics of R N into (1.7), we compute the asymptotic behavior of K N by a saddle point analysis. It should not come as a surprise to the experienced reader that there many possible fallpits and one may view the fact that this approach can indeed be carried out as the main result of our paper. With this approach one can, in principle, compute all fine asymptotic properties of the model. In an effort to limit the length of the paper, we restrict our main results to the description of the disordered region and the densities of the different types of lozenge there. We will though briefly comment on possible other limiting results that are within reach.

Statement of Results
In this section we state our main results. The proofs are postponed to later sections.

Preliminaries.
Our main result concerns the limiting densities of the lozenges as the size of the hexagon goes to infinity. We introduce the scaled variables (ξ, η) in the large N limit by where the point (ξ, η) belongs to the hexagon We will study the following probabilities , and P (x, y) .
Here (x, y) is the coordinate for the black dot. From simple geometric considerations, we note that these probabilities add up to 1. Our main result, Theorem 2.5 below, gives the limits of the probabilities (2.2) under the scaling (2.1) provided that (ξ, η) belongs to the liquid region. The result is stated in terms of a saddle point for the double contour integral in (1.7). The saddle points turn out to be solutions of an algebraic equation . See Definition 4.2 below for the precise definition of the branches of the logarithm and the square root in (2.10).
The following definition is central for the saddle point analysis of the double integral in (1.7).

Definition 2.2.
For each 0 < α ≤ 1 and (ξ, η) ∈ H, we define α (z) = α (z; ξ, η) as any solution of the equation In the low temperature regime 0 < α < 1 9 , we see from (2.7) that Q α is the square of a rational function. This means that (2.11) factorizes and α decouples into two rational functions with poles at −1, −α, 0 and a zero at ∞. This in turn implies that we obtain two well-defined rational functions α,± from (2.11): (2.12) α then is a meromorphic function defined on the Riemann surface R α associated with the equation w 2 = (z − z + (α))(z − z − (α)). It has two sheets R α,± , that are connected by a cut from z + (α) to z − (α) that we choose as where we recall from (2.8) that θ α = arg z + (α) = − arg z − (α). We take w = ((z − z + )(z − z − )) 1/2 with the branch of the square root that behaves like z as z → ∞ on the first sheet R α,+ and that behaves like −z as z → ∞ on the second sheet. 4. On the right, the two-sheeted Riemann surface for the high temperature case 1 9 < α ≤ 1 is displayed. The function α is meromorphic on the Riemann surface with simple poles at the indicated points − 1, − α, 0 on both sheets and a simple zero at both points at ∞. In the low temperature case 0 < α < 1 9 , the cuts from z + (α) to z − (α) disappear and the surface decouples, resulting in the picture that is displayed at the left Accordingly we have two branches of α , see also Fig. 4. The function α is meromorphic on the Riemann surface with simple poles at −1, −α, 0 on both sheets and a simple zero at both points at ∞. The four remaining zeros will be the saddle points for the double contour integral.

Saddle points and the liquid region.
We next describe the liquid region for general 0 < α ≤ 1. A reader acquainted with the asymptotic analysis of similar models for which the kernel can be represented in terms of double integral formulas, will recall that the liquid region in such cases is defined in terms of the saddle points of a phase function occurring in the integrand (see for example [12,32,63,66]). In the present situation, the function α from (2.12), (2.13) plays the role of the derivative of the phase function, which now turns out to be multivalued. The saddle points are the zeros of α . As was the case in previous works, we are interested in the particular saddle with strictly positive imaginary part (if it exists). For a given (ξ, η) ∈ L α with s = s(ξ, η; α), let T 1 and T α denote the triangles in C with vertex sets {−1, 0, s} and {−α, 0, s}, respectively. As indicated in Fig. 5, the angles of T 1 and T α are denoted by {φ 1 , φ 2 , φ 3 } and {ψ 1 , ψ 2 , ψ 3 }, respectively. Note that φ 3 = ψ 3 for any α, but φ j = ψ j for j = 1, 2 if and only if α = 1. The following is the main result of the paper.
π , x even. (2.14) π , x even, (2.15) and lim N →∞ (2.16) Theorem 2.5 follows from Proposition 7.7 below, and the proof of this proposition will be given in Sect. 7. Theorem 2.5 describes the situation in the liquid region L α , but it also explains the behavior at the boundary of L α . For each (ξ, η) ∈ L α , both s(ξ, η; α) and s(ξ, η; α) are simple zeros of α . When the point (ξ, η) approaches the boundary of L α , the saddle s(ξ, η; α) approaches the real line. Thus, at the boundary ∂L α , two zeros of α collide to form a double zero. Note also that when s(ξ, η; α) approaches the real line, the triangles T 1 and T α collapse with two of the angles approaching 0 and the third approaching π . In view of Theorem 2.5, this means that the tiling is frozen at the boundary of L α .
2.5. Structure in the low temperature regime. Let us now discuss the low temperature regime in more detail.
In the low temperature regime, each zero of α is a zero of one of the functions α,+ or α,− from (2.12). These zeros are easy to find since each of the functions α,± is as a rational function with a quadratic numerator. Setting the numerators equal to zero leads to the equations 2 ) . (2.17) with z ± = z ± (α). The equations (2.17) are quadratic in s with discriminants D ± = D ± (ξ, η) that depend on the coordinates ξ and η: (2.18) The equations D + (ξ, η) = 0, D − (ξ, η) = 0 represent two ellipses in the (ξ, η)-plane. The ellipses are inside the hexagon and each one of them is tangent to the boundary of the hexagon in four points. The two ellipses are disjoint for 0 < α < 1 9 , and they become tangent at the origin for α = 1 9 . Since a quadratic equation has two complex conjugate roots if and only if the discriminant is negative, we readily obtain the following proposition Proposition 2.6. For each 0 < α < 1 9 , the liquid region L α is the disjoint union of the two open ellipses L ± α defined by with D ± = D ± (ξ, η) given by (2.18). Moreover, the restrictions of (ξ, η) → s(ξ, η; α) to L ± α are diffeomorphisms onto C + .
See Sect. 3 for the proof, in particular of the statement about the diffeomorphisms. Let us now discuss the behavior of the ellipses near the boundary of the hexagon. The three poles z = 0, z = −α, z = −1 of α,± (z) together with the point at infinity correspond under the map s precisely to the points (ξ, η) where the ellipses touch the hexagon, see Fig. 6. A computation gives the following explicit expressions for the points of tangency: where the + and − signs correspond to the subscripts 1 and 2, respectively. Given two points P, Q on one of the ellipses ∂L ± α , we use the notation γ P Q ⊂ ∂L ± α to denote the counterclockwise subarc of the ellipse which starts at P and ends at Q. As (ξ, η) ∈ L α approaches a point in The liquid region (left) and the two disconnected sheets of R α (right) in the low temperature regime. The diffeomorphism (ξ, η) → s(ξ, η; α) maps the points A j , B j , C j , D j to − 1, − α, 0 and ∞, respectively The situation is more interesting on the arcs γ A 1 B 1 and γ A 2 B 2 . As (ξ, η) ∈ L α approaches one of these arcs, s(ξ, η; α) approaches the interval (−1, −α). In this limit we have φ 2 = π and ψ 1 = π , while all the other angles are zero. This means that at a point (x, y) near this part of the boundary of the liquid domain, we have i.e., there is an alternating pattern involving two different types of lozenges, as is clearly visible in Fig. 2. 2.6. Structure in the high temperature regime. In the high temperature regime 1 9 < α ≤ 1, the equation α (s; ξ, η) = 0 for the saddle points can be written after squaring as (2.23) The following proposition (which should be compared with Proposition 2.6) shows that s defines a diffeomorphism from the liquid region L α to the subset R + α of R α defined by The liquid region (left) and the two sheets of the Riemann surface R α (right) in the high temperature regime. The diffeomorphism (ξ, η) → s(ξ, η; α) maps the boundary points A j , B j , C j and D j to − 1, − α, 0, and ∞, respectively Proposition 2.7. For each 1 9 The boundary ∂L α of the liquid region is part of the zero set of the discriminant of the quadratic equation (2.23). Since the discriminant is invariant under the map (ξ, η) → (−ξ, −η), its zero set is symmetric with respect to the origin. Moreover, the zero set contains the line η = ξ/2, because (2.23) has a double zero at s = − √ α when η = ξ/2.
This line is however not part of the boundary of L α . The discriminant also vanishes at all points (ξ, η) which satisfy an algebraic equation of degree six. The real section of this algebraic curve is a curve inside the hexagon that touches the sides of the hexagon at the points (see Fig. 7) The liquid region is symmetric with respect to the line η = ξ/2. The cusp points are located at where η cusp = ξ cusp /2 and We also have η cusp = cos θ α 2 . Note that ξ cusp = 0 for α = 1/9 and ξ cusp = 1 for α = 1.
At points on the subarc of the boundary ∂L α between B j and C j we have (2.19), between C j and D j we have (2.20), and between D j and A j we have (2.21). This is a consequence of Theorem 2.5 and it is the same as in the low temperature regime. Finally, we have the alternating probabilities (2.22) between A 1 and B 2 , and between A 2 and B 1 .
A notable difference compared with the low temperature regime is that the liquid region in the high temperature regime is connected. As a result, the frozen region with the two types of tiles (sometimes called semi-frozen region) becomes disconnected into two disjoint components.
For α = 1, the Eq. (2.23) has a double root at s = −1 and two other roots that are the solutions of The latter two roots coincide if 4ξ 2 − 4ξη + 4η 2 = 3 and this is the equation for the ellipse that is tangent to all six sides of the hexagon. The semi-frozen region disappears for α = 1.

Local process in the bulk.
We chose to present Theorem 2.5 as our main result, but we stress that our method of proof allows us to compute much more complicated asymptotic behaviors (in this sense, our method of proof is the most important contribution of this paper). For instance, with a minor adaptation of the proof of Theorem 2.5 we compute the asymptotic behavior of local correlations in the bulk of the liquid region. (2.25) where ξ N and η N are such that and N ξ N and N η N are integers for every N ∈ N. We will additionally assume, without loss of generality, that N ξ N is even. The variables u 1 , u 2 , v 1 and v 2 are integer valued local variables independent of N . Then we have the limit (2.26) where s = s(ξ, η; α) and the integration path from s to s in (2.26 The proof of this theorem is given in Sect. 7. If u 1 = u 2 then the integral at the right-hand side of (2.26) can be computed explicitly to be the discrete sine kernel. For general u 1 and u 2 this is thus a kernel that is an extension of that discrete sine kernel. In fact, it falls into the class of extensions of the discrete sine kernel introduced in [10]. It is to note that the limiting kernel, and thus its associated point processes, depends on α. The periodicity in the horizontal direction is thus preserved in the limit. Theorem 2.8 gives the limiting correlation kernel for the point process of the paths. However, from the path picture one can compute the correlation functions for the different lozenges. For instance, the particle/hole duality tells us that the lozenges form a determinantal point processes with 1 − K N as correlation kernel. Under the same assumptions of Theorem 2.8 (but possibly with more than two points) we thus have , (2.27) whereK is the kernel at the right-hand side of (2.26).
2.8. Some comments on further asymptotic results. We end this section by commenting on further possible results on the asymptotic behavior of the random tilings.
Remark 2.9 (Frozen regions). The complement of the liquid region L α inside the hexagon, is called the frozen region. By definition, in the frozen region there are no solutions of α (z; ξ, η) = 0 in C + and all solutions are real. By using a saddle point analysis similar to the one we give in the proof of Theorem 2.5, one can show that this implies exponential decay of the fluctuations. Thus, in the frozen regions the randomness disappears rapidly and the tiling converges to deterministic patterns. In the corners of the hexagon the patterns are simple in the sense that we only have one type of lozenge in each corner. For α < 1 there are also other frozen regions near the centers of the vertical sides. Also here the randomness decays rapidly, but there are two types of lozenges forming a stair case pattern (as we also see in the degenerate siuation α = 0 as shown in the left picture in Fig. 3). Frozen regions that have different types of lozenges have appeared in other models. Some examples are [11,32] (after identifying Gelfand-Tsetlin patterns with lozenge tilings of the half plane). In fact, lozenge tilings of the infinite hexagon (or plane partitions) with an arbitrarily chosen back wall have been studied [17,60,61]. Part of this back wall can be a frozen region with more complicated patterns than the staircase pattern of the present paper.
Remark 2.10 (Edge universality). At the boundary of the liquid region (away from the points where the boundary touches the sides of the hexagon, and, in the high temperature regime, away from the cusp points) we expect Airy behavior. There is a vast amount of literature around this type of universality, and we only refer to [48] for an overview of results.
Remark 2.11 (Turning points). The turning points are the points where the boundary of the liquid region touches a side of the hexagon. Here we need to distinguish between the turning points that touch the hexagon at a vertical side from the other turning points. In both the low and high temperature regimes (assuming α < 1) there are four such points. They separate two frozen regions: one that contains two different types of lozenges, while the other has only one type of lozenges. We expect the local processes there to be the same as the processes that were found in (with a similar weight) in [60]. At the turning points that are not at the vertical sides of the hexagon we expect the GUE minor process [49] to appear.
Remark 2.12 (Cusp points). In the high temperature limit, the boundary of the liquid region has cusp points. Such cusp points have appeared before in the context of random tilings. It is known that the local limit process near such a cusp point is the Pearcey process [4,9,65,71].
We strongly believe that all the above universal behaviors can be verified using rather straightforward modifcations of the analysis that we present in this paper. More involved are the following remarks: Remark 2.13 (Tacnode). At the critical value α = 1 9 there is a transition from the low to high temperature regimes. The liquid region becomes a union of two ellipses that are tangent at the origin, and the origin is a tacnode. The tacnode process was first characterized in [1] and alternative characterizations were given shortly afterwards in [30,47]. See also [2,39]. Preliminary computations indicate that the same tacnode process appears, but we will return to this in a forthcoming paper.
Remark 2.14 (Height fluctuations). Another interesting feature of random tilings are the fluctuation of the height function. It was found in [51] that the limiting height function can be described by the complex Burgers equation. In [51] it is also conjectured that the fluctuations are described by the Gaussian Free Field. There is by now a long list of random tiling models where this conjecture has been verified, and we only mention [12,[19][20][21]32,33,67]. This turns out to be a very robust universality. Also in the model considered in this paper, we expect the Gaussian Free Field to appear, but with an interesting transition from the low to high temperature regimes. In the low temperature regime, the correlations between the different ellipses are expected to converge to zero exponentially and we expect to obtain two independent Gaussian Free Fields (in the appropriate coordinates), whilst we have only one Gaussian Free Field in the high temperature regime. It is natural to ask how these two fields merge to one in the transition from the low to high temperature regime. We plan to answer this question in a forthcoming paper.

2.9.
Overview of the rest of the paper. In the next section we first prove Propositions 2.3, 2.6 and 2.7.
The rest of the paper is devoted to the proof of Theorem 2.5. It is an asymptotic analysis of the double integral in (1.7) for K N (x, y, x, y) and for related double integrals that give the probabilities for each of the three lozenges. These double integrals are presented in Theorem 7.1 below.
The asymptotic analysis has two main parts. In the first part we analyze the orthogonal polynomials and their reproducing kernel R N (w, z) in the large N limit. The orthogonal polynomials are characterized by a RH problem that is essentially due to Fokas, Its and Kitaev [40]. This is recalled in Sect. 5.2. The reproducing kernel has a convenient formulation in terms of the solution of the RH problem, see Proposition 5.3. For the asymptotic analysis we use the Deift-Zhou steepest descent method for RH problems. A main ingredient for the analysis is the g-function, which in the present context is associated with an equilibrium measure on a contour in the complex plane.
This equilibrium measure is discussed in detail in Sect. 4. The transition at α = 1 9 is visible in the equilibrium measure since for 1 9 < α ≤ 1 the equilibrium measure is supported on a circular arc in the complex plane, while for 0 < α ≤ 1 9 it is supported on a full circle. We are able to give explicit formulas for the equilibrium measure, see Definition 4.2.
The steepest descent analysis of the RH problem is done in Sect. 5. We do not need strong asymptotics of the reproducing kernel R N , it suffices to have a uniform bound on R N (w, z)e N (g(w)−g(z)) (this is in Corollary 5.6) where R N (w, z) is a function related to the reproducing kernel, and which is given by (5.8).
The second part of the asymptotic analysis is a saddle point analysis of the double integrals like the one in (1.7). The saddle points depend on the asymptotic location (ξ, η) in the hexagon. We focus on the lower left part of the liquid region which corresponds to η ≤ ξ 2 ≤ 0. Then the saddle point s = s(ξ, η; α) is the zero of the derivative of a function α that is introduced in Sect. 6.1. We want to move the contours in the double integrals to contours γ z and γ w passing through the saddles s and s, and such that whenever w ∈ γ w \{s, s} and z ∈ γ z \{s, s}. To be able to do the deformation we need an analysis of the critical level set Re α (z) = Re α (s) of Re α passing through the saddle. This is done in Sect. 6.2.
The actual deformation and splitting of contours is done in Sect. 7. It turns out that the limiting probabilities in (2.14), (2.15), (2.16) come from residue contributions that arise from pole crossings during the deformations of contours. The remaining double contour integrals are then estimated and we only need they tend to zero as N → ∞. The details of the deformations are different for the low and high temperature regimes.

Proofs of Propositions 2.3, 2.6 and 2.7
In this section we prove Propositions 2.3, 2.6 and 2.7. We consider the low and high temperature regimes separately.
3.1. The low temperature regime. Since the saddle point equation α (s; ξ, η) = 0 reduces to the two quadratic equations (2.17) in the low temperature regime 0 < α < 1 9 , and also in the critical regime α = 1 9 , Proposition 2.3 is straightforward to prove in this regime. Proof of Proposition 2.3 for 0 < α ≤ 1 9 Any solution to α (s; ξ, η) = 0 is a solution to one of the quadratic equations in (2.17). The discriminants for these quadratic equations are given in (2.18). If, and only if, one of the discriminants is negative, then the corresponding quadratic equation has a zero in C + . Since the discriminants cannot be simultaneously negative, the statement follows.
Proof of Proposition 2.6. It is clear from the discussion preceding Proposition 2.6 that It is therefore enough to show that the restrictions of (ξ, η) → s(ξ, η; α) to L ± α are diffeomorphisms onto C + . We will show that for each s with Im s > 0, there are unique points Since ξ and η are real, we obtain the following two real equations by taking the real and imaginary parts of (3.1): .
We readily see that for s ∈ C + . Hence the 2 × 2 matrix on the left-hand side of (3.2) is invertible whenever Im s > 0. It follows that given s ∈ C + we can recover ξ ± and η ± uniquely by . (3.4) This proves that the restrictions of s to L ± α are bijections onto C + . The differentiability is also clear, and thus we have proved the statement.

The high temperature regime.
We now consider the high temperature regime and thus assume 1 9 < α ≤ 1. We start by defining the polynomial α by By (2.23), the zero set of α is the image of the zero set of α under the natural projection These are all simple calculations based on (3.5). The inequalities hold since Proof. If η = ξ/2 then, by parts (c), (d), and (e) of Lemma 3.1, α has a sign change, and therefore a zero, in each of the intervals (−1, − √ α) and (− √ α, −α). For η = ξ/2, α has a zero at − √ α by part (d), and in fact as can be checked from (3.5).
We now give the proof of Proposition 2.3 in the high temperature regime.
Proof of Proposition 2.3. for 1 9 < α ≤ 1. From Corollary 3.2 it follows in particular that there are at least two zeros of α in (−1, −α) in case α < 1. The remaining two zeros can also be real (frozen phase), or be a pair of complex conjugate non-real zeros (liquid phase). There is at most one complex conjugate pair of non-real zeros, and thus at most one zero with strictly positive imaginary part. By continuity this last fact also holds for α = 1. This proves Proposition 2.3 in the high temperature regime.
Proof of Proposition 2.7. The proof is similar to the proof of Proposition 2.6.
see (2.5) and (2.23). As in the proof of Proposition 2.6, we obtain two real equations by considering the real and imaginary parts. It follows that given s ∈ R + α , where R + α denotes the subset of R α defined in (2.24), we recover ξ and η from where the choice of square root in Q α (s) 1/2 is dictated by the location of s on the Riemann surface (different sign on different sheets). This shows that (ξ, η) → s(ξ, η; α) is a bijection from L α to R + α . It is clearly also differentiable (but not analytic!) and therefore it is a diffeomorphism. It also extends continuously to the boundary of L α mapping for example A 1,2 to −1, B 1,2 to −α, C 1,2 to 0, D 1,2 to ∞, and E 1,2 to − √ α, where the points with subscript 1 are mapped to the first sheet and points with subscript 2 to the second sheet, see also Fig. 7. We finally prove that the line segment For η = ξ/2, we see from (3.6) that α (s) has a double zero at − √ α while the two remaining zeros satisfy and with the aid of trigonometric identities one can show that s = √ αe iψ is a zero of α (s). If η increases from 0 to η cusp , then θ decreases from π to θ α , and ψ increases from θ α to π . It follows that s moves along the circle with radius √ α from z + to − √ α, that is, it moves along one side of the cut C on the Riemann surface. By symmetry, if η decreases from 0 to −η cusp then the saddle moves along the same circle but on the other side of C.

Preliminaries.
The orthogonality (1.2) does not depend on the specific choice of contour γ . By analyticity we can deform it to any other contour γ 0 that goes around 0 once in the positive direction. For the asymptotic analysis we need to select the 'correct' contour. The correct contour is typically (but not always...) the contour that attracts the zeros of the orthogonal polynomials as the degree tends to infinity. In (1.2) the orthogonality weight Such problems were studied in approximation theory where V is referred to as an external field [70]. Since the works of Stahl [69] and Gonchar-Rakhmanov [42] it is known that the zeros tend to a contour with a certain symmetry property for the logarithmic potential of its equilibrium measure. Such contours are now called S-contours. Later, Rakhmanov [68] made a systematic study of a max-min characterization of Scontours, and with Martínez-Finkelshtein [58] introduced the notion of a critical measure and identified the S-contours as trajectories of quadratic differentials. See [54,59] for further developments and historical remarks.
For α = 1 the external field (4.1) has only two logarithmic singularities and in such a case the orthogonal polynomials can be written in terms of classical Jacobi polynomials. Indeed, the nth degree polynomial p n is a multiple of the Jacobi polynomial in case α = 1. The Jacobi polynomial is non-standard, since one of the parameters is negative. The asymptotic zero distribution of Jacobi polynomials with varying nonstandard parameters was studied in [53,56,57]. The case (4.2) is contained in [57], see also [31], and it is known that the zeros of (4.2) tend to an arc on the unit circle as n, N → ∞ with n/N → 1.

Equilibrium measure.
In order to successfully apply the RH steepest descent analysis to the RH problem 5.2, we need a contour γ 0 going around 0 and a probability measure μ 0 on γ 0 with a corresponding g-function such that, for some constant ∈ C, Im with V as in (4.1). We call a probability measure μ 0 satisfying (4.3)-(4.5) an equilibrium measure in the external field V .
For a given γ we consider the probability measure μ on γ that minimizes the energy functional among all probability measures on γ . By classical results from logarithmic potential theory [70], there is a unique minimizer and it satisfies the conditions (4.4) on the real part of g + +g − −V . In order to be an equilibrium measure for V (as we defined it) we also need the condition (4.5) on the imaginary part. This condition characterizes S-contours. Indeed, by the Cauchy-Riemann equations the property (4.5) is equivalent to and ∂ ∂n ± denotes the normal derivatives on γ . This property is known as the S-property of 0 , and γ 0 is an S-contour.
We remark that the equilibrium measure is not necessarily unique. For example, if V (z) = log z then the normalized Lebesgue measure dμ = ds 2πis on any circle centered at the origin is an equilibrium measure for V . The radius is arbitrary and the equilibrium measure is not unique. This is a more general phenomenon in case the support is a full closed contour.

Construction of the equilibrium measure.
From conditions (4.4)-(4.5) it follows that we are looking for μ 0 such that g + + g − − V is piecewise constant on the support of μ 0 and therefore is analytic across the support of μ 0 . Thus Q is an analytic function in the complex plane with singularities determined by the singularities of V . We can furthermore recover μ 0 from Q. Indeed with an appropriate branch of the square root, and then by the Sokhotski Plemelj formula In our case of interest we have (4.1) and is rational with three simple poles. Therefore by (4.6) Q = Q α is a rational function with double poles at z = 0, z = −1, and z = −α. We can determine Q α explicitly, and it is given by the formulas in Definition 2.1, see also Sect. 4.6 below. We will prove that the associated measure (4.7) is indeed an equilibrium measure with external field V α .
while for 1 9 < α ≤ 1 the square root Q α (z) 1/2 was considered as a function on the first sheet of the Riemann surface R α shown in the right panel of Fig. 4. From now on it will be more convenient to change the branch cut of the Riemann surface from C to where θ α = arg z + = − arg z − . We also modify the definition of Q α (z) 1/2 so that now is defined and analytic for z ∈ C\ 0 with the square root such that Q α (z) 1/2 ∼ 1 z as z → ∞. The circular arc (4.10) will be the support of the equilibrium measure μ 0 .
We let γ 0 denote the circle of radius √ α centered at 0 oriented in the counterclockwise direction. With (4.9) and (4.11), we define the measure μ 0 , the associated g-function, and the variational constant as follows.
The associated g-function is defined by where z → log(z − s) is defined in the same way as in the high temperature regime.
(c) We define the variational constant ∈ C by (4.14) The definition (4.14) is such that equality holds in (4.4) at z = z + ∈ 0 for 1 9 < α ≤ 1 and at z = √ α ∈ 0 for 0 < α ≤ 1 9 . For the steepest descent analysis of the RH problem, it is convenient to introduce a function φ(z) which is a primitive function of Q α (z) 1/2 (with appropriate choices of the branch).
with Q 1/2 α given by (4.11), and the integration path from z + to z does not intersect with Q 1/2 α given by (4.9), and the integration path from The formulas (4.12) and (4.13) define μ 0 as a complex measure on 0 . The fact that it is a probability measure is part of the statement of the following proposition whose proof is given in Sect. 4.5.
for all z in the domain of φ. Moreover,  ∈ (a, b). A trajectory or an orthogonal trajectory is critical if it contains a zero or a pole of Q. Proof. Let z = z(t) = √ αe it , so that z = i z. For α ≥ 1 9 , we write z ± = √ αe ±iθ α with 0 < θ α ≤ π , and then by (2.5) This expression is indeed < 0 for −θ α < t < θ α and > 0 for θ α < t < π and −π < t < −θ α . For 0 < α < 1 9 , a similar computation using (2.7) and (2.6) gives (4.21) Since 0 < α < 1 9 we have 1+3α 2 > 2 √ α and therefore the numerator is always > 0. Thus Q α (z)(z ) 2 < 0 for every t ∈ [−π, π]. For α > 1 9 we recall that z ± are simple zeros of Q α . From the local structure of trajectories of a quadratic differential there are three critical trajectories emanating from each of the points z ± . One of these is an arc on the circle |z| = √ α, as we have seen. The other critical trajectories also connect z + with z − and a representative situation is shown in Fig. 8.
The trajectories of the quadratic differential Q α (z)dz 2 are level lines of Re φ, since φ is a primitive function of ±Q 1/2 α as follows from Definition 4.3. The orthogonal trajectories are level lines of Im φ.
Since √ α ∈ 0 we in fact have that Re φ = 0 on 0 as well as on the other critical trajectories (in the high temperature regime) that are shown in Fig. 8 for α = 0.3. The three critical trajectories are boundaries of three regions in the complex plane on which Re φ has a constant sign. Namely Re φ < 0 in the region containing −1, and Re φ > 0 in the region containing 0 and in the unbounded region.
To prove this we introduce Then 0 is contained in N φ , but N φ also contains other parts, see Figs. 9 and 10 for representative figures in the high and low temperature regimes. The first thing to observe is that Re φ extends to a continuous function on C away from −1, −α, and 0. Indeed, Q 1/2 α has simple poles at these three values, and therefore by integration as in definitions (4.15) and (4.16), we find that φ has logarithmic behavior. However, since the residues of Q 1/2 α are real, the real part of φ is single-valued. Thus Re φ is continuous on C\{−1, −α, 0} and harmonic on C\( 0 ∪ {−1, −α, 0}). We also note Fig. 9. The set N φ = {z ∈ C : Re φ(z) = 0} = −1 ∪ −α ∪ 0 is shown for α = 1 8 . This set divides C into three regions, and the sign of Re φ is shown in each of these regions In the high temperature regime the level set N φ consists of the critical trajectories of the quadratic differential Q α (z)dz 2 emanating from z + (α). Lemma 4.6. Let 1 9 < α ≤ 1. The set N φ consists of three analytic arcs connecting z + and z − which we denote by −1 , −α and 0 . The arc −1 intersects the real axis at x 1 ∈ (−∞, −1) and −α intersects the real axis at x 2 ∈ (−α, 0). The arc 0 is the support of the measure μ 0 and is part of the circle |z| = √ α.
Proof. Because of the local behavior of trajectories of a quadratic differential at a simple zero, there are three trajectories emanating from z + . One of these trajectories is 0 . The other two trajectories have to remain bounded and stay away from the poles −1, −α, 0 by (4.22). They have to come to the real axis. Indeed, if not, they would have to form a close loop in the upper haf plane and, since Re φ is harmonic inside this closed loop, we obtain a contradiction with the maximum/minimum principle for harmonic functions. Therefore, the trajectories come to the real axis and, by symmetry, they continue to the other simple zero z − = z + . The three trajectories enclose two bounded domains and Re φ = 0 on the boundary of these domains. Again, note that Re φ is harmonic in the interior, except at −1, −α, 0, where it tends to ±∞, see (4.22). By the maximum/minimum principle of harmonic functions each of the domains should contain at least one of the singularities. Again by (4.22) there are points x 1 ∈ (−∞, −1) and x 2 ∈ (−α, 0) with Re φ(x 1 ) = Re φ(x 2 ) = 0. Also Re φ( √ α) = 0 and we claim that x 1 , x 2 , √ α are the only points in To see this we recall that φ = Q 1/2 α , with a branch cut along 0 for the square root. From the formula (4.11) we then see that φ changes sign in the five values −1, − √ α, −α, 0, and √ α ∈ 0 . Thus φ > 0 (and Re φ is strictly increasing) on the intervals where the inequality holds because of the variational inequality (4.4) at − √ α ∈ γ 0 \ 0 , which in the high temperature regime is a strict inequality, see also (4.29). Therefore Re φ has no zeros in (−1, −α), and we proved the claim that We conclude that one critical trajectory comes to x 1 and another one to x 2 . This defines the contours −1 and −α .
It remains to prove there are no other parts in N φ . Any potential other part of N φ cannot intersect the real axis, as we already saw. Then such a part would be a closed contour in the upper or lower half plane and we arrive, again, at a contradiction because of the maximum/minimum principle for harmonic functions.
The structure of N φ is different in the low temperature regime, see Fig. 10. Lemma 4.7. Let 0 < α < 1 9 . The set N φ is the disjoint union of three analytic closed curves which we denote by −1 , −α and 0 . The closed curve 0 is the circle of radius √ α around 0, as before, and −1 , −α are two closed curves lying in the exterior/interior of 0 and going around −1 and −α, respectively.
Proof. Because of (4.22) the level set N φ is bounded and stays away from the poles −1, −α, and 0 of Q α . Since we already know from Lemma 4.5 that Re φ(− √ α) = 0, we infer from (4.16) that the zeros z ± of Q α are not on N φ . Therefore N φ does not contain any critical trajectories and hence consists of a finite union of disjoint closed curves. Because of the maximum/minimum principle for harmonic functions each component of C\N φ contains at least one of the singularities −1,−α, 0, or ∞.
A closer inspection of Re φ(z) for z ∈ R (also based on (4.9), (4.16) and (4.22) reveals that N φ has six intersection points with R. Two of them are the points ± √ α that belong to 0 . Then we have one point in each of the intervals (−∞, −1), (−1, − √ α), (− √ α, −α) and (−α, 0). This shows that there is a closed curve −α inside 0 and a closed curve −1 outside 0 as indicated in the statement.

Proof of Proposition 4.4.
We compute 0 dμ 0 by means of a residue calculation. Let us first consider the case 1 9 ≤ α ≤ 1. Then by (4.12) and the fact that Q α,+ (s) 1/2 = −Q α,− (s) 1/2 for s ∈ 0 , we have where C is a closed contour going around 0 once in the positive direction, and without enclosing any of the poles. Deforming the contour C to infinity, we pick up residue contributions at the poles. It is a straightforward calculation to show that The residues add up to zero, and since Q α (s) 1/2 = 1 s + O(s −2 ) as s → ∞, we thus find from (4.23) Let z(t) = √ αe it , −θ α < t < θ α , be a parametrization of 0 . Then the mapping has as its derivative which is real and non-zero for t ∈ (−θ α , θ α ) since Q α (z)(z ) 2 < 0 as 0 is a trajectory of the quadratic differential by Lemma 4.5 (a). Note also that the right-hand side of (4.26) vanishes for t = −θ α and equals 1 for t = θ α by (4.25). Therefore (4.26) is monotonically increasing from 0 to 1 as t goes from −θ α to θ α , and this is enough to conclude that μ 0 is a probability measure on 0 .
It now also follows (compare (4.15) and (4.26), and use Q 1/2 α,+ = −Q 1/2 α,− on 0 ) that φ − is purely imaginary along 0 and we have for z ∈ 0 . (4.27) Next we calculate g (z) = 0 dμ 0 (s) z−s . We write g as a contour integral with the same closed contour C as in (4.23), but we now also assume that z is in the exterior of C. We deform the contour to infinity where we now pick up a residue contribution from the pole at s = z as well, which is Q α (z) 1/2 . We use (4.24) to calculate the other residue contributions. There is no contribution from infinity and the result is that Integrating (4.28) from z + to z along a path that does not intersect which proves (4.17) for α ∈ [ 1 9 , 1] by the definition (4.14) of and the fact that φ(z + ) = 0.
We have also shown that φ − (z) ∈ iR for z ∈ 0 , and similarly φ(z) ∈ iR on the other critical trajectories that emanate from z + and z − , see Fig. 8. Moreover, Im φ is constant on orthogonal trajectories. We also saw that Im φ − (z) increases as z moves away from z − to z + along 0 . Then by the Cauchy-Riemann equations, we have Re φ > 0 in the domain on the minus side of 0 and by continuity it holds in the outer domain bounded by the critical trajectories. Then Re φ < 0 if we cross the critical trajectory going around −1, and in particular Re φ(z) < 0 for z on the critical orthogonal trajectory from z + to − √ α. In view of (4.17), this gives for z on this orthogonal trajectory, which is part of γ 0 \ 0 . This proves the inequality in (4.4). By symmetry the inequality also holds for z on the critical orthogonal trajectory from z − to − √ α. This completes the proof for the case α ≥ 1 9 . The proof for 0 < α < 1 9 is simpler. In this case (4.9) is a rational function with partial fraction decomposition The total integral of μ 0 defined by (4.13) is by a simple residue calculation with contributions only from the poles at s = 0 and s = −α. The total mass is 1 and as before it follows that μ 0 is a probability measure. We compute g (z) with another residue calculation Recalling the definition (4.16) of φ(z) and the expression (4.8) for V α (z), we conclude Integrating (4.30) from √ α to z along a path that does not intersect (−∞, 0] ∪ 0 , we Then (4.31) also holds for |z| < √ α, since g + ( √ α) = g − ( √ α) + πi, as can be verified from the definition of the branch of log(z − s) that was used in the definition of g. Thus (4.17) holds for 0 < α < 1 9 in the low temperature regime because of the definition of the constant . The identities (4.18) and (4.19) follow from (4.17) in the same way as in the case 1 9 < α ≤ 1.

4.6.
Calculations leading to Q α . The reader may wonder how to obtain the expressions (2.5) and (2.7). One clue is that we need the residues (4.24). This translates into the three conditions (which are consistent with (4.6)) lim z→0
To proceed, we make the one-cut assumption which says that Q α should have at least one multiple zero. It means that the discriminant of the numerator polynomial should be zero. The discriminant factors as (1 − α) 2 which leaves us with four possible choices for A, namely A 1 = 3 + α, A 2 = 3α + 1, For α = 1 we should recover (2.9) which means that we have to take A = A 4 for α = 1, and then by continuity also for α between 1 and a critical value of α. This leads to the formulas (2.5) and (2.4). The critical value is when z + (α) = z − (α), and this happens for α = 1/9.
For α = 1 9 , the two values A 2 and A 4 coincide, and for α < 1 9 we find that A 2 takes over. This leads to the formulas (2.7) and (2.6) with two double zeros of Q α .

Orthogonal Polynomials and Riemann-Hilbert Problem
We will now prove the existence of the orthogonal polynomials and pose a RH problem for the reproducing kernel R N (w, z) that appears in the double contour integral in the kernel (1.7).

Existence of the orthogonal polynomials.
Proposition 5.1. Let 0 < α ≤ 1 and N ∈ N. Then for every n = 0, 1, . . . , 2N there is a unique monic polynomial p n of degree n such that In addition, if n ≤ 2N − 1, then Proof. The orthogonality condition (5.1) is associated with the non-Hermitian bilinear form defined for polynomials f and g. The polynomial p n exists and is unique if and only if the n × n matrix of moments M n = z j , z k n−1 j,k=0 is invertible. We use the Lindström-Gessel-Viennot (LGV) lemma to prove that this is the case for n ≤ 2N . Consider the directed graph on Z 2 with an edge from (i, j) to (i , j ) if and only if i = i + 1 and j − j ∈ {0, 1}. The weights on the edges are w ((i, j), (i + 1, j) , j), (i + 1, j + 1)) = 1.
For two vertices A, B ∈ Z 2 we define where the sum is over all directed paths P on the graph from vertex A to vertex B. If there are no such paths then w(A, B) = 0. We assume 0 ≤ n ≤ 2N and we take vertices A j = (0, j) and B j = (2N , 2N −n+ j) for j = 0, 1, . . . , n − 1. The LGV lemma [41] states that det w(A j , B k ) n−1 j,k=0 is equal to the weighted sum of all non-intersecting path systems from A 0 , . . . A n−1 to B 0 , . . . , B n−1 . It is easy to verify that there exist such non-intersecting path systems (due to the fact that 0 ≤ n ≤ 2N ). Each non-intersecting path system has a positive weight since α > 0. Therefore det w(A j , B k ) n−1 j,k=0 > 0, which, in particular, implies that is an invertible matrix.
To calculate w(A j , B k ) we observe that any path from A j to B k is of length 2N with n − k + j horizontal edges. The weight of such a path is α l where l is the number of horizontal edges at an even level. We pick l out of the possible N even levels, and n − k + j − l out of the possible N odd levels, and we see that there are N l N n−k+ j−l paths from A j to B k with weight α l . Summing over l yields This sum over products of binomial coefficients is easily seen to be equal to the coefficient of z 2N −n+k− j in the product (z+1) N (z+α) N . Therefore, by Cauchy's theorem Comparing (5.3) and (5.4) we then see that M n is obtained from W n by reversing the order of the columns. Since W n is invertible, also M n is invertible, and it follows that p n uniquely exists.
To prove (5.2) let us assume that κ n = 0. Then by orthogonality we have p n , z j = 0 not only for j = 0, 1, . . . , n − 1 but also for j = n. It follows again by orthogonality of p n+1 in case n ≤ 2N − 1, that p n+1 + p n , z j = 0 for every j = 0, 1, . . . , n. However, we established that p n+1 is the only monic polynomial of degree n + 1 with these properties (if n ≤ 2N − 1). This contradiction shows that κ n = 0.

Riemann-Hilbert problem.
It is well-known that the orthogonal polynomials and the associated Christoffel-Darboux kernel can be characterized by a RH problem.

Riemann-Hilbert Problem 5.2. Let γ 0 be the circle of radius
√ α around 0 with positive direction. Find a function Y : C\γ 0 → C 2×2 with the following properties: The limits of Y (z) as z approaches γ 0 from inside and outside exist, are continuous on γ 0 and are denoted by Y + and Y − , respectively. Furthermore they are related by The RH problem 5.2 is due to Fokas, Its, and Kitaev [40]. Its solution contains the orthogonal polynomials of degrees N and N − 1 that uniquely exist by Proposition 5.1, for z ∈ C\γ 0 .

Proposition 5.3. (a) The kernel R N is given in terms of the solution Y of the RH problem 5.2 by
Proof. The formula (5.7) is a reformulation of the Christoffel-Darboux formula (1.8), as can be readily checked from (5.6) together with the fact that det Y ≡ 1. The formula (5.8) is obtained from (5.6) in a similar way.

First transformation of the RH problem.
The steepest descent analysis of the RH problem 5.2 for Y is fairly standard by now. It is modelled after the method developed by Deift et al. [28] for orthogonal polynomials on the real line. The extension to the complex plane is standard, once one has identified the correct contour γ 0 with the equilibrium measure μ 0 . In the high temperature regime we basically follow [28] including the construction of Airy parametrices for the local analysis at branch points z ± . The RH analysis in the low temperature regime is even simpler since we can separate contours and no local analysis is needed. The critical case α = 1/9 is more difficult, but can be handled with the construction of a local parametrix built out of Lax pair solutions associated with the Hastings-McLeod solution of Painlevé II. This is similar to the construction in [23] for orthogonal polynomials on the real line in cases where the equilibrium density vanishes quadratically at an interior point of its support. We will not give any details for this case.
In terms of the function V α defined in (4.1), the jump relation (5.5) for Y can be expressed as The first transformation Y → T uses the g-function from Definition 4.2 to normalize the RH problem at infinity. We define The jumps in the RH problem for T are conveniently expressed in terms of the function φ defined in (4.15) and (4.16). From the identities (4.17), (4.18), and (4.19) and the definition (5.9), we find the following RH problem.
(b) T has boundary values on γ 0 that satisfy Note that T depends on N , even though this is not indicated in the notation. What is important for us, is that T and T −1 remain bounded as N → ∞, provided we stay away from the branch points z ± (α) (only in the high temperature regime). We summarize what we need from the RH analysis in the following proposition.
The proposition is a result of the steepest descent analysis that we will perform next for the two regimes separately.
Because of (5.9) and the formula (5.8) for R N , we have (5.13) and before turning to the proof of Proposition 5.5 we note the following consequence. Proof. Parts (a) and (b) are immediate from (5.13) and Proposition 5.5. Because of (5.13) and the jump condition (5.11) for T along γ 0 \ 0 , the analytic continuation from part (c) is given by Since Re φ(w) < 0 for w in the region under consideration in part (c), see for example Fig. 9, part (c) follows from Proposition 5.5 as well.
We define , for z between γ 0 and γ + , , for z between γ 0 and γ − , Then S satisfies the following RH problem.

Riemann-Hilbert Problem
has boundary values on γ 0 , γ + and γ − that satisfy We remove the constant jump on γ 0 by defining Of course R should not be confused with the reproducing kernel R N , as these are totally different objects. The matrix valued function R satisfies the following RH problem.

Riemann-Hilbert
Since Re φ > > 0 for z ∈ γ + ∪γ − the jumps in the RH problem for R are exponentially close to the identity matrix as N → ∞. By standard estimates on small norm RH problems [27], we find R(z) = I + O(e − N ) as N → ∞, and in particular R and R −1 are uniformly bounded as N → ∞, uniformly on C. Tracing back the transformations (5.16) and (5.15) it then also follows that T and T −1 are uniformly bounded as N → ∞, uniformly on C, since Re φ ≥ 0 in the annular region bounded by γ + and γ − . This proves Proposition 5.5 for α < 1 9 . In case α = 1 9 , we are not able to choose γ + and γ − such that (5.14) holds on the full contours. Instead we let γ + and γ − go to γ 0 at the critical point − √ α = − 1 3 , and we can do it in such a way Re φ > 0 on (γ + ∪ γ − )\{− 1 3 }. Then we can proceed as in the case The jump contour γ 0 for the RH problem 5.2 for Y (black), the curves −1 and −α (red), and the points − 1, − α, 0 (black dots) in the high temperature regime α < 1 9 described above, except that we have to build a local parametrix at − 1 3 . This is done with the help of a special parametrix [23] that we will not describe here. We only need to know that it is uniformly bounded as N → ∞ and then Proposition 5.5 follows as before.

Proof of Proposition 5.5 (b) .
Proof. Suppose 1 9 < α ≤ 1 and let Y (z) denote the solution of the RH problem 5.2 with jump contour γ 0 . See Fig. 12 for γ 0 together with the contours −1 and −α that enclose the bounded domain where Re φ < 0 in the high temperature regime.
In small disks D z + and D z − around the endpoints of 0 we construct local parametrices P (z + ) and P (z − ) with the aid of Airy functions. This construction is standard by now and we do not give details. The only thing that concerns us is that the local parametrices depend on N and they slightly grow with N , namely uniformly for z ∈ D z ± . The third and final transformation S → R is defined by Then R is defined and analytic in with jump matrices that are I + O(N −1 ) as N → ∞. It follows that R(z) = I + O(N −1 ) uniformly for z ∈ C, and in particular R and R −1 remain bounded as N → ∞. Observe that in the construction of the local parametrics, the disks D z ± can be chosen arbitrarily small (but independent of N ), and we choose them with radii smaller than δ. Then following the transformations (5.17) and (5.20), and taking note of (5.19) we find that T and T −1 are uniformly of order N 1 6 as N → ∞. Outside the disks D z ± the global parametrix (5.18) applies, which does not change with N , and then T and T −1 remain uniformly bounded. Part (b) of Proposition 5.5 is now also proven.

Phase Functions α and α
6.1. Definitions. In the last two sections we analyzed the RH problem with the gfunction coming from the equilibrium measure as its main input. The outcome of this analysis is in Corollary 5.6 which states that R N (w, z)e N (g(w)−g(z)) remains uniformly bounded in certains regions, and actually (very roughly) We now turn to the asymptotic analysis of the double contour integrals coming from (1.7) and that give the probabilities for the three types of lozenges, see also Theorem 7.1 below.
After deforming contours and splitting up integrals, we are able to rewrite the integrals with an integrand containing as the main N -dependent entry, where F(z; x, y) = (z + 1) 3) see Propositions 7.8 and 7.9. Recall that x, y will be varying with N as in (2.1). Then in view of (6.1), (6.3) we see that (6.2) behaves roughly like e N ( α (z)− α (w)) with a certain function α that depends ons (ξ, η) ∈ H, and that is defined next, along with a companion function α . The equality leading to the third line in (6.4) follows from (4.17) and (4.1). Recall that φ = ±Q 1/2 α by Definition 4.3 and therefore by the definitions (6.4) and (6.5) we have that both α and α satisfy the algebraic equation (2.11) for α .
For each (ξ, η) ∈ L α , the saddle s(ξ, η; α) defined in Definition 2.4 is a zero of either α and α . Lemma 6.2. Let (ξ, η) ∈ L α and s = s(ξ, η; α). Then we have Let us consider the low temperature regime. From the formula (2.12) for α,± and (2.18) it follows that s is a zero of α,± if and only if D ± < 0, and we note that the regions D ± < 0 are contained in the regions η > ξ 2 and η < ξ 2 , respectively. Using (3.3) and (3.4) we see that, in the low temperature regime, ξ has the same sign as For the high temperature regime, we use Proposition 2.7 and the proof is analogous to the proof in the low temperature regime, but now (6.8) is replaced by ∓ Im s Q α (s) 1 2 , with the same choice of branch for the square root as in (3.7), i.e., Q α (s) 1 2 has a branch on C.

Critical level set of Re α .
In what follows we focus on the case (a) of Lemma 6.2, namely (ξ, η) ∈ L α with η < ξ 2 < 0, and its extension η = ξ 2 < 0, which is the lower left part of the liquid region. The corresponding saddle s = s(ξ, η; α) satisfies α (s) = 0 and |s| < √ α if η < ξ 2 . For η = ξ 2 < 0 (which is only relevant in the high temperature regime) we have |s| = √ α with θ α < arg s < π, and we still have α (s) = 0. We are interested in the level set of Re α that contains s, We emphasize that α has a branch cut along 0 . However Re α is well-defined and continuous, also on 0 . Typical behaviors of N are shown in Figs. 14, 15 and 16. The level set N makes a cross locally at s since it is a simple saddle. Four curves emanate from s that are denoted by 1 , …, 4 in the figures.
It is important for us that three of these curves stay inside 0 (in low temperature regime) or inside 0 ∪ −1 and connect s with s. Only one of them (denoted by 4 in the figures) meets with either 0 or 0 ∪ −1 .
To be able to prove this we need information on the behavior of the two functions z → log |z| and z → log (z+1)(z+α) z on −1 ∪ 0 . We start with a lemma. Lemma 6.3. We have the following for 0 < α ≤ 1, if and only if z = √ α and z ∈ (γ 0 \ 0 ) ∪ R.
Proof. (a) We consider the case 0 < α < 1. Observe that z 2 Q α (z) tends to 1 as z → ∞, and there are no sign changes on the real line. Thus z 2 Q α (z) ≥ 0 for real values of z, with double poles at z = −1 and z = −α, and a local minimum at z = √ α. There is a minimum at z = − √ α in case α ≥ 1 9 , and a local maximum at z = − √ α in case α < 1 9 . In the latter case there are local minima at z = z ± . It can be verified that From an inspection of the graph, it follows that for any x > αQ α ( √ α), x = 1, there are four real solutions to the equation For x = 1 there are three real solutions and a solution at infinity, while for α Q α (− √ α) < x < αQ α ( √ α) there are two real solutions. If α ≤ 1 9 , there are again four real solutions (counting multiplicities) for each 0 ≤ x ≤ α Q α (− √ α).
(6.11) and in the cases (6.11) there are only two real solutions.
Then, with an argument similar to the one we used to prove part (a), we check that these are the only z for which (6.13) is in (0, ∞). This proves part (c).  Proof. Indeed, from the definition (6.4) and the fact that Re φ = 0 on −1 and 0 , we obtain for z ∈ −1 ∪ 0 , and by Lemma 6.4 the sum at the right-hand-side of (6.17) is strictly decreasing since ξ < 0 and ξ 2 − η ≥ 0. Due to Corollary 6.5, we see that the level set (6.9) has at most one point of intersection with ( −1 ∪ 0 ) ∩ C + , because Re α is strictly decreasing there. Therefore at least three of the j 's, say 1 , 2 , 3 , do not intersect ( −1 ∪ 0 ) ∩ C + , which means that they have to go to the real line inside the domain enclosed by −1 ∪ 0 (or inside the disk bounded by 0 in the low temperature regime), and then by symmetry end at s inside that domain. Taking p j ∈ j ∩ R for j = 1, 2, 3, we choose the ordering of the j 's such that p 1 < p 2 < p 3 . The contours 1 , 2 , 3 enclose two bounded domains for which Re α is constant on the boundaries and harmonic inside, except at the singularities −1, −α, 0, where Re α is unbounded by (6.7). By the maximum principle for harmonic functions, each of the two domains has to contain at least one of the singularities. Also Re( α − α (s)) has opposite signs on the two bounded domains. Then again by (6.7) one domain contains 0 and the other domain contains −α, and possibly also −1, since at both these points Re α tends to −∞. Thus Fig. 15. The level set N (blue) and the contours −1 ∪ 0 in the high temperature regime (here α = 1 8 ) in case −1 < p 1 < −α. The set N divides the plane into five regions and the sign of Re( α − α (s)) is indicated in each of these five regions by + or − If 4 would remain inside −1 ∪ 0 as well, then it would also go to the real line, say at a point p 4 , and continue to s inside this domain. If p 3 < p 4 < √ α then 4 and 3 would enclose a domain with Re α is constant on the boundary, and harmonic inside, and we would have a contradiction with the maximum principle. If p 4 < p 1 then 4 and 1 enclose a bounded domain within and we find a contradiction in the same way.
Thus 4 comes to ( −1 ∪ 0 ) ∩ C + and continues into the outer domain of C\N φ . It cannot go to infinity because of (6.7) and so it has to go to the real line at a point p 4 and by symmetry it continues in the lower half plane where it crosses −1 ∪ 0 again and ends at s.
As Re decreases along ( −1 ∪ 0 ) ∩ C + from left to right, we find Re α ( √ α) < Re α (s). Since Re α (z) → +∞ as z → ∞, the level set N intersects the real line at a point > √ α. This can only be at p 4 . Thus 4 and 3 enclose a domain where Re α < Re α (s) and that contains (part of) 0 where α has its branch cut, and where Re α is not harmonic. Hence there is no contradiction with the maximum principle.
To summarize, we have a situation as in Fig. 14 in case p 1 < −1, or as in Fig. 15 if −1 < p 1 < −α. In the latter case, there is also a separate part 5 of N that goes around −1. Figures 14 and 15 are for the high temperature regime. In the low temperature regime we have that 0 is the full circle of radius √ α. Then in the above discussion we can replace −1 ∪ 0 by 0 . It follows that 1 , 2 , 3 stay inside the disk of radius √ α, and so 1 does not go around −1. There is always a part 5 going around −1 in the low temperature regime, as shown in Fig. 16.
It is now clear that we can find contours as described next. See also Figs. 17 and 18 below. Corollary 6.6. Let (ξ, η) ∈ L α with η ≤ ξ 2 < 0. (a) In the low temperature regime there are closed contours γ z and γ w,in , γ w,out such that • γ w,out lies outside the circle γ 0 , does not go around −1, and is such that The level set N (blue) and the contours −1 and 0 in the low temperature regime (here α = 1 10 ). The set N α divides the plane into five regions and the sign of Re( α − α (s)) is indicated in each of these five regions by + or − • γ w,in lies inside the circle γ 0 , goes around −α, and it passes through s and s in such a way that Re α (w) > Re α (s), w ∈ γ w,in \{s, s}.
• γ z lies inside the circle γ 0 , goes around 0, and it passes through s and s in such a way that Re α (z) < Re α (s), z ∈ γ z \{s, s}.
(b) In the high temperature regime there exist contours γ z and γ w such that • γ w lies in the domain bounded by 0 ∪ −1 , it goes around −1, and it passes through s and s in such a way that Re α (w) > Re α (s), w ∈ γ w \{s, s}, • γ z lies inside the circle γ 0 , goes around 0, and it passes through s and s in such a way that In the low temperature regime we will also use γ w = γ w,in ∪ γ w,out .

Lozenge probabilities.
In the final part of the analysis we are going to deform contours in the double contour integral to the ones from Corollary 6.6, which leads to the proof of Theorem 2.5. We start by expressing the probabilities for the three types of lozenges as double contour integrals.
We use F(z; x, y) as in (6.3) and for a function (w, z) → H (w, z), We will use (7.1) only for functions (w, z) → H (w, z) that are products of a rational function in w and a rational function in z, both with poles at −1, −α, and 0 only. In addition, the integrand in (7.1) will have singularities for w = 0 and z = 0 only, and the contour γ 0 can be deformed to an arbitrary closed contour around 0, and we can take different contours for the two integrals.
Proof. The proof is straightforward.
The next step is a double integral formula for the expectation value of the height function.
Proof. By the determinantal structure of the correlations (see Proposition 1.
We insert the double contour integral formula of Lemma 7.3 and combine the two integrals by subtracting the two integrands. Since , if x is odd, which we can check from (6.3) separately for x even or odd, the formula (7.2) follows. Note also that the pole at z = w disappeared when we took the difference, and thereforẽ γ can be moved back to γ in (7.2).
The proof of (7.3) is similar, and (7.4) is immediate from the structure of the determinantal point process, as already noted after the statement of Theorem 7.1.

Contour deformation in the low temperature regime
We start the analysis of the double integral (7.1) with a contour deformation. There are several ways to deform the contours, and the ones we are going to present will be useful for the lower left part of the liquid region, that is for (ξ, η) ∈ L α with η ≤ ξ/2 < 0 as in Corollary 6.6. The deformations will be different for the low and high temperature regimes. Proposition 7.8. Let 0 < α ≤ 1 9 and (ξ, η) ∈ L α with η < ξ 2 < 0. Let γ z , γ w,in and γ w,out be closed contours as in Corollary 6.6 (a), (see also Fig. 17). Then (7.1) is equal to where R N is given by (5.8) and F is given by (6.3).
Proof. In (7.1) we use γ z for the integral with respect to the z variable, and γ 0 (initially) for the w variable. By the conditions in Corollary 6.6 (a), the contour γ z lies inside γ 0 . By Sokhotskii-Plemelj formula and (5.8) we have for w ∈ γ 0 , where the ± boundary values are with respect to the w variable. This we substitute into the double integral (7.1) to obtain the difference of two double integrals, We deform γ 0 inwards to γ w,in in the first double integral and outwards to γ w,out in the second double integral. (Recall that +-side refers to the interior of γ 0 and −-side to its exterior.) We do not encounter any singularites of the integrand if we do the deformation into the exterior domain, since by assumption γ w,out does not go around −1. Thus by Cauchy's theorem we obtain the last term in (7.15).
In the deformation of the first integral we pick up residue contributions for those z ∈ γ z that are in the exterior of γ w,in . This is due to the pole at w = z that we encounter when deforming γ 0 into γ w,in . Since R N (z, z) = 1, the contribution of the poles leads to the first term in (7.15). The remaining double integral is the second term in (7.15).

Contour deformation in the high temperature regime
In the second proposition (relevant for the high temperature case) we modify the definition (5.8). We use a large circle γ ρ centered at the origin of radius ρ > 10 and define Note that (7.16) coincides with (5.8) for w inside γ 0 , and it is the analytic continuation (in the w variable) of (5.8) with |w| < α to the disk |w| < ρ. Because of (5.13) and the jump (5.11) of T , we have Proposition 7.9. Let 1 9 < α < 1 and (ξ, η) ∈ L α with η ≤ ξ 2 < 0. Suppose γ z and γ w are closed contours as in Corollary 6.6 (b), (see also Fig. 18). Let (x, y) be coordinates inside the hexagon. Then the double contour integral (7.1) is equal to (7.18) where R N is given by (7.16) and F is given by (6.3).
Proof. As in the proof of Proposition 7.8 we have (but now we use (7.16)) Fig. 18. The contours γ z (green) and γ w (black) in the high temperature regime. The contours satisfy the conditions of Corollary 6.6 (b) and Proposition 7.9 with w ∈ γ ρ , and the ± boundary values are for w ∈ γ ρ .
We choose γ ρ for the contour in the w integral in (7.1) and γ z for the z integral. Then the double contour integral is a difference of two double integrals with γ z inside γ ρ . The integrand in the second double integral has no singularities for |w| > ρ, since the poles are at w = z, w = −1, w = −α, and they are all inside. For |w| > ρ we have R N (w, z) = R(w, z). From the asymptotic behavior in the RH problem 5.2 for Y we get as w → ∞, and thus by (5.8) Also by the definition of F, see (6.3), we have (F(w; x 2 , y 2 )) −1 = O(w y 2 −x 2 ) as w → ∞. By combining with (7.5), we see that the full integrand in (7.19) is therefore O w −N +y 2 −x 2 −1 as w → ∞. Since (x, y) is a point inside the hexagon, we have inequalities −N < y 2 −x 2 < N . Thus, since we are dealing with integers, the integrand is O(w −2 ) as w → ∞. Therefore the second double integral in (7.19) vanishes identically.
In the first double integral we deform γ ρ to γ w as in the statement of the proposition. We pick up a residue contribution at the pole w = z for those z ∈ γ z that lie in the exterior of γ w . This gives the first term in (7.18). The remaining double integral is the second term in (7.18). 7.5. Proof of Proposition 7.7. We are now ready for the proof of Proposition 7.7 which, as already noted leads to the proof of Theorem 2.5. We also noted that it suffices to prove the proposition for (ξ, η) ∈ L α with η ≤ ξ 2 ≤ 0. We first assume ξ < 0 and later deal with the modifications that are necessary for ξ = 0.
We write x = x N = (1 + ξ N )N , y = y N = (1 + η N )N , and we are in the situation where For N large enough, we then also have (ξ N , η N ) ∈ L α with ξ N 2 < 0. We may also assume that η N ≤ ξ N 2 < 0, because of symmetries as in Proposition 7.4 (b) and Proposition 7.6.
7.5.1. Low temperature regime with η < ξ 2 < 0 Let γ (N ) z and γ (N ) w,in , γ (N ) w,out be contours as in Corollary 6.6 (a) and Proposition 7.8 but corresponding to the parameters (ξ N , η N ) and s = s N . Then by (7.15) and in view of (7.20) it is enough to show that the two double integrals in (7.21) tend to 0 as N → ∞. By Corollary 5.6 (a) there exists a constant C 1 > 0 such that , if x N is odd.
The contours stay away from −α and −1, therefore the extra factor in case x N is odd remains bounded and bounded away from 0. Combining this with (7.22) we obtain for some constant C 2 > 0, Fig. 19. The sets N (left) and N (right) in the high temperature regime for ξ = 0 and η < 0. The signs of Re( α − α (s)) (left) and Re( α − α (s)) (right) are indicated with ± tends to 0 as N → ∞. We recall that w → R N (w, z) is the analytic continuation of w → R N (w, z) from the disk |w| < √ α into the large disk |w| < ρ. It then follows from Corollary 5.6 (b) and (c) that R N (w, z) ≤ C 1 e N (g(z)−g(w)) (7.26) whenever w is in the domain bounded by 0 ∪ −1 and z ∈ C with w, z bounded away from the branch points z ± . This is the estimate that is analogous to (7.22) in the low temperature regime. By Corollary 6.6 (b) the contour γ (N ) w is inside 0 ∪ −1 , and we can apply (7.26) in the estimation of (7.25). The rest of the proof is the same as in the low temperature regime with ξ < 0.
For ξ = 0, the saddle is on the branch cut 0 for the functions φ and α . We need additional deformation of contours to handle this case. For definiteness we focus on the high temperature regime, but the low temperature regime can be done similarly.
We treat the case (0, η) ∈ L α with η < 0 as a limit of the case (ξ, η) with η < ξ 2 < 0 that we considered before. In this limit the contours from Corollary 6.6 (b) can be chosen in such a way that they tend to contours γ z and γ w that partly overlap with 0 , such that the following hold (see Fig. 20 together with Fig. 19, left) • γ w contains the subarcs Re α (z) < Re α (s), z ∈ γ z \ 0 , Re α,+ (z) = Re α (s), z ∈ γ z ∩ 0 . (7.28) We want to estimate the double integral in (7.18) with x = x N = (1 + o(1))N and y = y N = (1 + η + o(1))N as N → ∞. To avoid the use of N dependent contours as in the proofs above (which can be handled but would obscure the exposition) we assume x N = N + O(1) and y N = (1 + η)N + O(1) as N → ∞. Then by combining (6.3), (6.4) with (7.17) we find that R N (w, z) F(z;x N ,y N ) F(w;x N ,y N ) (which is the main part of the integrand in (7.18)) is equal to 1 −e 2N φ(w) T −1 (w)T (z) 1 0 , w ∈ γ w , |w| > √ α (7.29) times a factor that remains bounded as N → ∞. In (7.29) we take + boundary values for α and T whenever w and/or z are on 0 .
Because of (7.27) and (7.28) we see that (7.29) becomes exponentially small as N → ∞ unless w ∈ γ w ∩ 0 and z ∈ γ z ∩ 0 . Here we also use that Re φ(w) < 0 for w ∈ γ w , |w| > √ α, and that T and T −1 remain bounded as N → ∞ if we stay away from the branch points, see Proposition 5.5 (b).
On γ z ∩ 0 we use the identity (7.30) which follows from the jump (5.10) of T across 0 . Using (7.30) in (7.29) we split the integral over γ z ∩ 0 into a sum of two integrals, and deform both integrals away from 0 . The integral with the first term of the right-hand side of (7.30) is deformed to the interior, that is to a contour from s to s lying inside the disk |z| = √ α. The dominant part of the integrand is e N ( α (z)−2φ(z)) and Re α (z) > Re α (s) and Re φ(z) > 0 for z on the deformed contour. Fortunately, Re( α (z) − 2φ(z)) < Re α (s), and this can be seen as follows. By (6.4) and (6.5) we have α − 2φ = α . Since ξ = 0 we also find from (6.4) and (6.5) that α + α = −2η log z. Thus indeed Re α (z) =− Re α (z)−2η log |z| < − Re α (s)−2η log |z| < Re α (s) =−η log √ α for z on the deformed contour, since Re α (z) > Re α (s)|z| < √ α < 1 there. We also use η < 0. Thus the deformed integral coming from the first term of (7.30) becomes small as N → ∞.
The integral with the second term is moved outwards, again to a contour from s to s but now lying in |z| > √ α. Since α,+ = α,− the deformed integral has the exponentially varying factor e N α . The contour can be taken such that Re α (z) < 0 on the contour (see Fig. 19, right), and again the contribution becomes small as N → ∞.
The integral (in the w-variable) over γ w ∩ 0 can be dealt with analytic continuation only. We note that by (5.10) which remains bounded if we analytically continue it to the exterior of 0 . We deform γ w ∩ 0 to a contour from s to z + (α) lying in the exterior of γ 0 together with its mirror image in the real, which is a contour from z − (α) to s. Since α,+ (w) = α,− (w) on 0 , the main term in the analytic continuation of (7.29) across γ w ∩ 0 becomes e −N α (w) . We are able to deform contours such that Re α (w) > 0 on the deformed contour (from Fig. 19, right), where we also take note of the local behavior near the saddle points s and s. The result is that the integral over the deformed contour becomes small as N → ∞. What remains are local contributions near the saddles s and s and also near the branch points z ± (α), since we cannot move γ w away from the branch points. The contributions from the saddles can be estimated as was done in detail for the low temperature regime with η < ξ 2 < 0. The contributions from the branch points are estimated similarly, but we have to note that T −1 (w) = O(N 1/6 ) for w close to the branch points, see Proposition 5.5 (b). This slight increase however still leads to a decay in the estimate and the conclusion is that all contributions vanish as N → ∞. for ξ = η = 0, and Re α (s) = 0 where s = s(0, 0; α) = z + (α). We approach this case as a limit of (ξ, η) ∈ L α with η ≤ ξ 2 < 0. In this limit the contours from Corollary 6.6 (b) tend to contours γ w and γ z that we may take as follows • γ w contains −1 and its analytic continuation (which is a critical orthogonal trajectory, see Fig. 8) such that Re α (w) > 0, w ∈ γ w \ −1 .
The integrand of the double integral in (7.18) behaves like (7.29) as N → ∞. With the above choice of contours the integrand is exponentially small unless w ∈ −1 and z ∈ 0 . The case z ∈ 0 is handled using the identity (7.30) that we also used in the case ξ = 0, η < 0. It allows us to split the integral into two integrals, deform one of them outwards and the other one inwards, and both deformed integrals have exponentially decaying integrands.
For w ∈ −1 we use the second line of (7.29) which tells us that the main wdependent part is which naturally splits into a sum (recall also α = φ) e −N φ(w) 1 0 T −1 (w) − e N φ(w) 0 1 T −1 (w) and a corresponding splitting and deformation of the w-integral. Namely the integral with the first term is deformed from −1 to a contour from z + (α) to z − (α) lying outside −1 (where Re φ > 0) and the integral with the second term is deformed inwards (where Re φ < 0).
Then there is exponentially decay on the deformed contours as N → ∞, except for w and z near the branch points z ± (α). T and T −1 have moderate growth there, both of O(N 1/6 ). They combine to give an increase in T −1 (w)T (z) of O(N 1/3 ). Local estimates still lead to a decay in the integrals, as required.
This completes the proof of Proposition 7.7 in all cases.

Proof of Theorem 2.8.
Proof. With the coordinates in (2.25) (and the fact that N ξ N is assumed to be even) we can rewrite the kernel K N in (1.7) as K N (x 1 , y 1 , x 2 , y 2 ) = − χ u 1 >v 2 2πi γ H K (z, z; u 1 , v 1 , u 2 , v 2 )dz + I N (N ξ N , N η N ; H K ) (7.31) where I N is as in (7.1) with H K (w, z; u 1 , v 1 , u 2 , v 2 ) = (z + 1) The first integral in (7.31) is independent of N . The asymptotic behavior of I N (N ξ N , N η N ; H K ) as N → ∞ is already computed in Proposition 7.7. The first integral and the limit from Proposition 7.7 can be combined naturally into one single integral, which is the right-hand side of (2.26). This finishes the proof.
Acknowledgement. Open access funding provided by Royal Institute of Technology.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

A. Proof of Proposition 1.1
Proof of Proposition 1.1. This is a special case of Theorem 4.7 in [34]. To identify the formula in [34] with (1.7), we first of all note that p = 1 and K N is a scalar kernel. We have to identify (m, x, m , y) and (N , M, L) in [34] with (x 1 , y 1 , x 2 , y 2 ) and (N , N , 2N ) in the setting of our paper.
Furthermore, for 0 ≤ i < j ≤ 2N , the notation A i, j (z) in [34] stands for A i, j (z) =  N (w, z) is a bivariate polynomial of degree ≤ N − 1 in both variables that is uniquely characterized by the property that for every polynomial q of degree ≤ N − 1, see Lemma 4.6 in [34]. Since all orthogonal polynomials p n of degrees n ≤ 2N exist (we prove this in Proposition 5.1), the sum in (1.8) is well-defined, and by orthogonality using (1.9) it defines a kernel with the required reproducing property (A.1). The expression in the second line of (1.8) is known as the Christoffel-Darboux formula, and it continues to hold for non-Hermitian orthogonality on a contour, with the same proof as for usual orthogonal polynomials on the real line.