Doubly periodic lozenge tilings of a hexagon and matrix valued orthogonal polynomials

Abstract We analyze a random lozenge tiling model of a large regular hexagon, whose underlying weight structure is periodic of period 2 in both the horizontal and vertical directions. This is a determinantal point process whose correlation kernel is expressed in terms of non‐Hermitian matrix valued orthogonal polynomials (OPs). This model belongs to a class of models for which the existing techniques for studying asymptotics cannot be applied. The novel part of our method consists of establishing a connection between matrix valued and scalar valued OPs. This allows to simplify the double contour formula for the kernel obtained by Duits and Kuijlaars by reducing the size of a Riemann–Hilbert problem. The proof relies on the fact that the matrix valued weight possesses eigenvalues that live on an underlying Riemann surface M of genus 0. We consider this connection of independent interest; it is natural to expect that similar ideas can be used for other matrix valued OPs, as long as the corresponding Riemann surface M is of genus 0. The rest of the method consists of two parts, and mainly follows the lines of a previous work of Charlier, Duits, Kuijlaars and Lenells. First, we perform a Deift–Zhou steepest descent analysis to obtain asymptotics for the scalar valued OPs. The main difficulty is the study of an equilibrium problem in the complex plane. Second, the asymptotics for the OPs are substituted in the double contour integral and the latter is analyzed using the saddle point method. Our main results are the limiting densities of the lozenges in the disordered flower‐shaped region. However, we stress that the method allows in principle to rigorously compute other meaningful probabilistic quantities in the model.


INTRODUCTION
A lozenge tiling of a hexagon is a collection of three different types of lozenges ( , , and ), which cover this hexagon without overlaps, see Figure 1 (left). There are finitely many such tilings; hence by assigning to each tiling  a nonnegative weight W( ), we define a probability measure on the tilings by where the sum is taken over all the tilings (and is assumed to be nonzero). Uniform random tilings of a hexagon (i.e., when W( ) = 1 for all  ) is a well-studied model. As the size of the hexagon tends to infinity (while the size of the lozenges is kept fixed), the local statistical properties of this model are described by universal processes. [1][2][3][4] We also refer to Refs. 5-7 for important early results and to Refs. 8, 9 for general references on tiling models. Uniform lozenge tilings of more complicated domains (nonnecessarily convex) have also been widely studied in recent years. [10][11][12][13][14] In this work, we consider the regular hexagon of (large) size F I G U R E 1 A tiling of a hexagon, and the associated non-intersecting paths F I G U R E 2 The graph  4 , and the periods of a 2 × 3 periodic weighting but we deviate from the uniform measure and study instead measures with periodic weightings. To explain what this means, we first briefly recall a well-known one-to-one correspondence between tilings of a hexagon and certain nonintersecting paths. This bijection can be written down explicitly, but is best understood informally. The paths are obtained by drawing lines on top of two types of lozenges and , as shown in Figure 1 (right). The paths associated to the tilings of  lie on a graph  , which depends only on the size of the hexagon, see Figure 2 (left). To each edge of  , we assign a nonnegative weight . The weight of a path is then defined as = ∏ ∈ , and the weight of a tiling  as W( ) = ∏ ∈ . Provided that at least one tiling has a positive weight, this defines a probability measure on the set of tilings by (1). If each edge is assigned the same weight, then we recover the uniform measure over tilings. We say that a lozenge tiling model has × periodic weightings if the weight structure on the edges is periodic of period in the vertical direction, and periodic of period in the horizontal direction, see Figure 2 (right) for an illustration with = 2 and = 3. Thus, a × periodic weighting is completely determined by 2 edge weights. Note that all paths share the same number of horizontal edges, and also the same number of oblique edges; hence lozenge tiling models with 1 × 1 periodic weightings are all equivalent to the uniform measure.
By putting points on the paths as shown in (3), each tiling of the hexagon gives rise to a point configuration, see also Figure 1 (right). Thus, the probability measure (1) on tilings can be viewed as a discrete point process. 15,16 For lozenge tiling models with × periodic weightings, it follows from the Lindström-Gessel-Viennot theorem 17,18 combined with the Eynard-Mehta theorem 19 that this point process is determinantal. Therefore, to understand the fine asymptotic structure (as → +∞) it suffices to analyze the asymptotic behavior of the correlation kernel. However, until recently, 20,21 the existing techniques were not appropriate for such analysis.
The main result of Ref. 20 is a double contour formula for the correlation kernels of various tiling models with periodic weightings (including lozenge tiling models of a hexagon as considered here). In this formula, the integrand is expressed in terms of the solution (denoted ) to a 2 × 2 Riemann-Hilbert (RH) problem. This RH problem is related to certain orthogonal polynomials (OPs), which are nonstandard in two aspects: • the OPs and the weight are × matrix valued, • the orthogonality conditions are non-Hermitian.
The size of the RH problem, the size of the weight, and the size of the OPs all depend on , but quite interestingly not on .
Lozenge tiling models of the hexagon with × periodic weightings are rather unexplored up to now. To the best of our knowledge, the model considered in Ref. 21 is the only one (other than the uniform measure) prior to the present work for which results on fine asymptotics exist. The model considered in Ref. 21 is 1 × 2 periodic and uses the formula of Ref. 20 as the starting point of the analysis. The techniques of Ref. 21 combine the Deift/Zhou steepest descent method 22 of (of size 2 × 2) with a nonstandard saddle point analysis of the double contour integral. However, because = 1, the associated OPs are scalar (this fact was extensively used in the proof) and it is not clear how to generalize these techniques to the case ≥ 2.
The aim of this paper is precisely to develop a method to handle a situation involving matrix valued OPs. We will implement this method on a particular lozenge tiling model with 2 × 2 periodic weightings, which presents one simply connected liquid region (which has the shape of a flower with six petals), six frozen regions, and six staircase regions (also called semifrozen regions); it will be presented in more detail in Section 2. The starting point of our analysis is the double contour formula from, 20 which expresses the kernel in terms of a 4 × 4 RH problem related to 2 × 2 matrix valued OPs. The method can be summarized in three steps as follows: 1. First, we establish a connection between matrix valued and scalar valued OPs. In particular, in Theorem 1 we obtain a new expression for the kernel in terms of the solution (denoted ) to 2 × 2 RH problem related to scalar OPs. This formula allows for a simpler analysis than the original formula from Ref. 20. 2. Second, we perform an asymptotic analysis of the RH problem for via the Deift-Zhou steepest descent method. The construction of the equilibrium measure and the associated -function is the main difficulty. The remaining part of the RH analysis is rather standard. 3. Third, the asymptotics for the OPs are substituted in the double contour integral and the latter is analyzed using the saddle point method.
The first step is the main novel part of the paper. The remaining two steps were first developed in Ref. 21 for a tiling model with 1 × 2 periodic weightings. Our main results, which are stated in Theorem 2, are the limiting densities of the different lozenges in the liquid region. However, we emphasize that the method also allows in principle to rigorously compute more sophisticated asymptotic behaviors in the model (such as the limiting process in the bulk).

An expression for the kernel in terms of scalar OPs
The eigenvalues and eigenvectors of the 2 × 2 orthogonality weight play an important role in the first step of the analysis. They are naturally defined on a two-sheeted Riemann surface , which turns out to be of genus 0. This fact is crucial to obtain the new formula for the kernel in terms of scalar OPs. We expect that ideas similar to the ones presented here can be applied to other tiling models with periodic weightings, as long as the corresponding Riemann surface  is of genus 0. Lozenge tiling models of large hexagons with periodic weightings can feature all of the three possible types of phases known in random tiling models: the solid, liquid, and gas phases. A solid region (also called frozen region) is filled with one type of lozenges. In the liquid and gas phases, all three types of lozenges coexist. The difference between these two phases is reflected in the However, lozenge tiling models of the finite hexagon cannot be represented as models with infinitely many paths (as opposed to the Aztec diamond and the infinite hexagon). In particular, they do not belong to the class of models studied in Ref. 27 and thus the simplified formula from Ref. 27 cannot be used. This fact makes lozenge tiling models of the finite hexagon harder to analyze asymptotically (see also the comment in Ref. 27, beginning of section 6).

The figures
In addition to being in bijection with nonintersecting paths, lozenge tilings of the hexagon are also in bijection with dimer coverings, which are perfect matchings of a certain bipartite graph. We refer to Ref. 4 for more details on the correspondence with dimers (see also Ref. 10, fig. 1, for an illustration). The bijection with dimers is not used explicitly in this paper, but we do use it to generate the pictures via the shuffling algorithm. 32

MODEL AND BACKGROUND
In this section, we present a lozenge tiling model with 2 × 2 periodic weightings. We also introduce the necessary material to invoke the double contour formula from Ref. 20 for the kernel. In particular, we present the relevant 2 × 2 matrix valued OPs and the associated 4 × 4 RH problem.

Affine transformation for certain figures of lozenge tilings
For the presentation of the model and the results, it is convenient to define the hexagon and the lozenges as in (2) and (3). However, for the purpose of presenting certain figures of lozenge tilings, it is more pleasant to modify the hexagon and the lozenges by the following simple transformation: so that  is mapped by this transformation to a hexagon whose six sides are of equal length. Above the definition (2) of  , we used the standard terminology and called  "the regular hexagon"; note however that  becomes truly regular only after applying the transformation (4). In the figures, we will assign the colors red, green, and yellow for the three lozenges in (4), from left to right, respectively.

Definition of the model
The regular hexagon  has corners located at (0,0), (0, ), ( , 2 ), (2 , 2 ), (2 , ), and ( , 0). We normalize the lozenges such that they cover each a surface of area 1, and the vertices of the lozenges have integer coordinates. We recall that each lozenge tiling of  gives rise, through (3), to a system of nonintersecting paths. These paths live on the graph  , which is illustrated in ). We denote the paths by and they satisfy the initial positions (0) = + . The particular 2 × 2 periodic lozenge tiling model that we consider depends on a parameter ∈ (0, 1]. The weightings are defined on the 2 × 2 bottom left block of the lattice as shown in Figure 3 (left), and is then extended by periodicity as shown in Figure 3 (right). More formally, if = (( 1 , 1 + otherwise. For any values of ∈ (0, 1], the weightings (6) are such that W( ) > 0 for all  , and thus we have a well-defined probability measure via (1). On the other hand, if = 0, then several edges have weights 0, and it is easy to see (e.g., from Figure 3 (right)) that W( ) = 0 for all  . So in this case, (6) does not induce a probability measure, and this explains why we excluded = 0 in the definition of the model. If = 1, all tilings have the same weight, and we recover the uniform distribution. Proposition 1 states that for < 1, there is a particular tiling  max that is more likely to appear than any other tiling.  max is illustrated in Figure 4 (left) for = 60. F I G U R E 5 Three tilings taken at random accordingly to the measure induced by (6) with = 100 and = 0.01 (left), = 0.05 (middle), = 0.2 (right) Proposition 1. Let ∈ (0, 1) and let ≥ 1 be an integer. There exists a unique tiling if is even, Proof. The proof of Proposition 1 is based from a careful inspection of  , and we omit the details. ■ It follows from Proposition 1 that, as → 0, the randomness disappears because the tiling  max becomes significantly more likely than any other tiling. Therefore, our model interpolates between the uniform measure over the tilings (for = 1) and a particular totally frozen tiling  max (as → 0), see Figures 4 and 5. Intriguingly, these figures show similarities with the rectangle-triangle tiling of the hexagon obtained by Keating and Sridhar in Ref. 33, fig. 18.
Several tiling models in the literature (e.g., those considered in Refs. 34 and 21) are defined by weightings on the lozenges, instead of weightings on the edges. To ease possible comparisons with these models, we give an alternative definition of our model. The weight W( ) of a tiling  can alternatively be defined as where is the weight function over the lozenges given by where ∈ (0, 1]. The above lozenge weightings depend only on the parity of and , and thus are periodic of period 2 in both directions. By using the correspondence (3) between lozenge tilings and nonintersecting paths, it is straightforward to verify that the weightings (7) and (8) define the same measure as the weightings (6).
More explicitly, this gives , if is even, , if is odd, and we can retrieve the entries of , +1 from its symbol by ( , +1 (2 1 , 2 2 ) , +1 (2 1 , 2 2 + 1) , +1 (2 1 + 1, 2 2 ) , +1 (2 1 + 1, 2 2 + 1) where is any closed contour going around 0 once in the positive direction. The symbol associated to  is then obtained by taking the following product (see Ref. 20, eq. (4.9)): To limit the length and technicalities of the paper, from now we take the size of the hexagon even, i.e., = 2 , where is a positive integer. This is made for convenience; the case of odd integer could also be analyzed in a similar way, but then a discussion on the partity of is needed. Because = 2 , following Ref. 20, eq. (4.12), the relevant orthogonality weight to consider is We consider two families { } ≥0 and { } ≥0 of 2 × 2 matrix valued OPs defined by where 0 2 denotes the 2 × 2 zero matrix, 2 is the identity matrix, and is, as before, a closed contour surrounding 0 once in the positive direction. Because the weight (13) The matrix is characterized as the unique solution to the following RH problem.

RH problem for
The limits of ( ) as approaches from inside and outside exist, are continuous on and are denoted by + and − , respectively. Furthermore, they are related by (c) As → ∞, we have ( ) = ( 4 + ( −1 ))

Double contour formula from Ref. 20 for the kernel
As mentioned in the Introduction, the point process obtained by putting points on the paths, as shown in (3), is determinantal. We let denote the associated kernel. By definition of determinantal point processes, for integers ≥ 1, and 1 , … , , 1 where , is defined by and  is given by As particular cases of the above, we obtain the following formulas.

STATEMENT OF RESULTS
The new double contour formula for the kernel in terms of scalar OPs is stated in Theorem 1. In this formula, the large behavior of the integrand is roughly Ξ , for a certain phase function Ξ, which in our case is defined on a two-sheeted Riemann surface  . The restriction of Ξ on the first and second sheet are denoted by Φ and Ψ, respectively. The saddle points are the solutions ∈ ℂ for which either Φ ′ ( ) = 0 or Ψ ′ ( ) = 0. In the liquid region, Proposition 3 states that there is a unique saddle, denoted , lying in the upper half plane. This saddle plays an important role in our analysis, and some of its properties are stated in Propositions 4 and 5. The limiting densities for the lozenges in the liquid region are given explicitly in terms of in Theorem 2.
Remark 1. If = 1, our model reduces to the uniform measure and the kernel can be expressed in terms of scalar-valued OPs. However, our approach is based on the formulas (21) and (22), and even though these formulas are still valid for = 1, this case requires a special attention (because of a different branch cut structure in the analysis). Because the limiting densities for the lozenges in this case are already well-known, 39 from now we will assume that ∈ (0, 1) to avoid unnecessary discussions.

New formula for the kernel in terms of scalar OPs
We define the scalar weight by and consider the following 2 × 2 RH problem.

RH problem for
(a) ∶ ℂ ⧵ ℂ → ℂ 2×2 is analytic, where ℂ is a closed curve surrounding and −1 once in the positive direction, but not surrounding 0. (b) The limits of ( ) as approaches ℂ from inside and outside exist, are continuous on ℂ and are denoted by + and − , respectively. Furthermore, they are related by (c) As → ∞, we have ( ) = ( 2 + ( −1 )) Remark 2. Theorem 1 is proved in Section 6. It is based on an unpublished idea of A. Kuijlaars that matrix valued OPs in a genus zero situation can be reduced to scalar orthogonality. In our case, the scalar orthogonality appears in (26) and (27) and a main part of the proof of Theorem 1 consists of relating the matrix valued reproducing kernel  from (20) to the scalar reproducing kernel  from (28).

The rational function 
The function  is a meromorphic function that appears in the equilibrium problem associated to the varying weight . Its explicit expression is obtained after solving a nonlinear system of five equations with five unknowns. Here, we just state the formula for , and refer to Section 8 for a more constructive approach. We define  as follows: where is given by (24), 1 , 2 , and 3 are given by and + and − are given by The zero + of  lies in the upper half plane, − = + , and the other zeros and poles of  are real. Furthermore, for all ∈ (0, 1), they are ordered as follows:

Saddle points and the liquid region
For each ( , ) ∈ , there are in total eight saddles for the double contour integral (29), which are the solutions to the algebraic equation where ( ) is given by (32). Following the previous works, 10,21,[40][41][42] we define the liquid region as the subset of  for which there exists a saddle lying in the upper half-plane ℂ + = { ∈ ℂ ∶ Im > 0}. Proposition 3 states that there is a unique such saddle (whenever it exists), which is denoted by ( , ; ). This saddle plays a particular role in the analysis of Section 11 and appears in the final formulas for the limiting densities (40).
Remark 3. For a given set , the notation stands for the closure of .
By definition, the saddles lie in the complex plane. We show here that they can be naturally projected on a Riemann surface. Define ( ) 1∕2 with a branch cut joining − to + along Σ 1 , such that ( ) 1∕2 ∼ 1 2 as → ∞, and denote the associated Riemann surface by  : This is a two-sheeted covering of the -plane glued along Σ 1 , and the sheets are ordered such that = ( ) 1∕2 on the first sheet and = −( ) 1∕2 on the second sheet. For each solution to (41), there exists a satisfying 2 = ( ), and such that

Limiting densities in the liquid region
Theorem 2 states that the limits (40) are expressed in terms of the angles shown in Figure 8. (38) with ( , ) ∈  . We obtain the following limits: These limits can equivalently be stated as follows: ) , By combining (51) with Theorem 2, we obtain the following immediate corollary: (38) with ( , ) ∈  . We have ) , where the three matrices inside each brackets correspond, from left to right, to = 1, 2, 3.

Outline of the rest of the paper
The proofs of Propositions 3, 4, and 5 are rather direct and are presented in Section 4. In Section 5, we follow an idea of Ref. 20 and perform an eigendecomposition of the matrix valued weight. The eigenvalues and eigenvectors are naturally related to a two-sheeted Riemann surface . The proof of Theorem 1 is given in Section 6, and relies on the fact that  is of genus 0. The proof of Theorem 2 is done via a saddle point analysis in Section 11, after considerable preparations have been carried out in Sections 7-10: • In Section 7, we use Theorem 1 to find double contour formulas for the lozenges in terms of scalar OPs. We also use the symmetry in our model to conclude that it is sufficient to prove Theorem 2 for the lower left quadrant of the liquid region. • In Section 9, we will perform a Deift/Zhou 22 steepest descent analysis on the RH problem for . This analysis goes via a series of transformations ↦ ↦ ↦ . The first transformation ↦ uses a so-called -function, which is obtained in Section 8. • The functions Φ and Ψ denote the restrictions of the phase function Ξ to the first and second sheets of  and play a central role in the large analysis of the kernel. In Section 10, we study the level set  Φ = { ∈ ℂ ∶ Re Φ( ) = Re Φ( )}, which is of crucial importance to find the contour deformations that we need to consider for the saddle point analysis.
As mentioned in Remark 1, we will always assume that ∈ (0, 1), even though it will not be written explicitly.

PROOFS OF PROPOSITIONS 3, 4(a), and 5 4.1 Proof of Proposition 3
By (50), the saddles are the zeros of the polynomial given by Because the coefficients of are real, Proposition 3 follows if has at least six zeros on the real line. This can be proved by a direct inspection of the values of ( ) at = −∞, 1 , 0, , 2 , , , 3 , −1 , +∞: the leading coefficients of are 1 − ( − ) 2 > 0. We conclude the following: 1. if ≠ − , has at least one simple root on (−∞, 1 ) and at least one simple root on , has at least one simple root on ( , 2 ) and at least one simple root on ( 2 , ), 3. if ≠ 2 , has at least one simple root on ( , 3 ) and at least one simple root on ( 3 , 1 ).

Proof of Propositions 4 and 5
We start with the proof of Proposition 5. By rearranging the terms in (47), we see that the saddles are the solutions to where satisfies 2 = ( ). This can be rewritten as Taking the real and imaginary parts of (55), and recalling that , ∈ ℝ, we get Thus, the 2 × 2 matrix at the left-hand side of (56) is invertible, and we get (49). This shows that ( , ) ↦ ( ( , ; ), ( , ; )) is a bijection from  to  + . This mapping is clearly differentiable, and therefore it is a diffeomorphism. Replacing ( , ) ↦ ( , − ) in the right-hand side of (49), we see that the left-hand side becomes ( , ) ↦ (− , − ). This implies the symmetry ( , ; ) = (− , − ; ). It remains to prove that ( , ) ∈  is mapped to a point ( ( , ; ), ( , ; )) lying in the upper half plane of the first sheet. The proof of this claim is splitted in the next two lemmas.
Proof. Consider the function defined by By the fundamental theorem, for each ∈ [0, +∞), there are eight solutions ∈ ℂ to ( ) = .
By Lemma 1, (58) is 0 if and only if ∈ Σ 1 , which implies that the sign of (58) is constant for where > 0, we conclude that (58) is negative for all sufficiently large and lying in ℂ + , and the claim follows. ■ This finishes the proof of Proposition 5 and Proposition 4(a). In principle, it is also possible to use (49) to prove parts (b)-(f) of Proposition 4, but it leads to more involved analysis. However, by rearranging the terms in (55), we can find other expressions than (49) for the mapping ( , ) ↦ ( , ) that lead to simpler proofs of (b)-(f). We only sketch the proof of (e). First, we rewrite (55) as Next, we verify that ) > 0, for Im > 0, which implies that = ( , ; ) has the same sign as ) .
Finally, in a similar way as in Lemma 1, we show that this quantity is 0 if and only if ∈ ℝ ∪ Σ , which proves part (e). We omit the proofs of parts (b), (c), (d), and (f).

ANALYSIS OF THE RH PROBLEM FOR
To describe the behavior of as → +∞, one needs to control the 2 × 2 upper right block of the jumps, which is ( ) 2 −2 . To do this, we follow an idea of Duits and Kuijlaars 20 and proceed with the eigendecomposition of . Then, we use this factorization to perform a first transformation ↦ on the RH problem.

Eigendecomposition of
The matrix ( ) defined in (12) has the following eigenvalues: where the + and − signs read for 1 and 2 , respectively, and + and − are given by and satisfy − < −1 < + < 0 and + − = 1. We define the square root √ ( − + )( − − ) such that it is analytic in ℂ ⧵ [ − , + ], with an asymptotic behavior at ∞ given by The eigenvectors of are in the columns of the following matrix: , and we have the factorization where Λ( ) = diag( 1 ( ), 2 ( )) is the matrix of eigenvalues. The matrix ( ) is analytic for ∈ ℂ ⧵ [ − , + ], and satisfies

First transformation ↦
The first transformation of the RH problem diagonalizes the 2 × 2 upper right block of the jumps, and is defined by Remark 4. By standard arguments, 43 we have det ≡ 1. Note however that the ↦ transformation does not preserve the unit determinant. Indeed, because det ( Using the jumps for given by (62), we verify that satisfies the following RH problem.

RH problem for
is analytic, where we recall that is a closed contour surrounding 0 once in the positive direction. (b) The jumps for are given by , where  ∶= ∩ [ − , + ]. Depending on ,  can be the empty set, a finite set, or an infinite set. If  contains one or several intervals, then on these intervals the jumps are .

PROOF OF THEOREM 1
First, we use the factorization of obtained in (61) together with the transformation ↦ given by (64), to rewrite the formulas (21) and (22) as follows: where  is given by The property (23) of  implies the following reproducing property for  : for every 2 × 2 matrix valued polynomial of degree ≤ − 1. Now, we introduce the Riemann surface  associated to the eigenvalues and eigenvectors of . This Riemann surface is of genus 0 and therefore there is a one-to-one map between it and the Riemann sphere (called the -plane).

The Riemann surface  and the -plane
The Riemann surface  is defined by and has genus zero. We represent it as a two-sheeted covering of the -plane glued along [ − , + ].
On the first sheet we require = + (1) as → ∞, and on the second sheet we require = − + (1) as → ∞. To shorten the notations, a point ( , ) lying on the Riemann surface will simply be denoted by when there is no confusion, that is, we will omit the -coordinate. If we want to emphasize that the point ( , ) is on the -th sheet, ∈ {1, 2}, then we will use the notation ( ) . With this convention, the two points at infinity are denoted by ∞ (1) and ∞ (2) . The function (0 (2) ) = −1, (0 (1) ) = 1.
The following identities will be useful later, and can be verified by direct computations: We define by From straightforward calculations using (73), we have

The reproducing kernel  
For ( ) on the -th sheet of  and ( ) on the -th sheet, we define   ( ( ) , ( ) ) by It is convenient for us to consider formal sums of points on , which are called divisors in the literature. More precisely, a divisor can be written in the form and we say that ≥ 0 if 1 , … , ≥ 0. The divisor of a nonzero meromorphic function on  is defined by where 1 , … , 1 are the zeros of of multiplicities 1 , … , 1 , respectively, and 1 +1 , … , 2 are the poles of of order 1 +1 , … , 2 , respectively. Given a divisor , we define (− ) as the vector space of meromorphic functions on  given by The following divisors will play an important role: Thus, ∶= (− ) is the vector space of meromorphic functions on , with poles at ∞ (1) and ∞ (2) only, such that the pole at ∞ (1) is of order at most − 1, and the pole at ∞ (2) is of order at most . Similarly we define * = (− * ). From the Riemann-Roch theorem, we have because there is no holomorphic differential (other than the zero differential) on a Riemann surface of genus 0.

Lemma 3.
We have ↦   ( , ) ∈ * for every ∈  * , (c)   is a reproducing kernel for in the sense that for every ∈ , where  is a closed contour surrounding once 0 (1) and 0 (2) on the Riemann surface  in the positive direction (in particular  visits both sheets).
Proof. Using the definitions of   and  given by (81) and (68), we can rewrite   as where From properties (a) and (b) of the RH problem for , the functions are analytic in  * . By combining the large asymptotics of ( ) (given by (63)) with property (c) of the RH problem for , we obtain ( ) from which we conclude that the functions 's have poles of order at most at ∞ (1) and at most + 1 at ∞ (2) . Therefore, we have shown that The numerator in (83) is, for each fixed ∈  * a linear combination of the functions , so belong to (−( + ∞ (1) + ∞ (2) )) as a function of . By definitions of   and  , the numerator vanishes for = (1) and for = (2) . Thus, the division by − in (83) does not introduce any poles, but it reduces the order of the poles at ∞ (1) and ∞ (2) by one, and therefore ↦   ( , ) ∈ as claimed in part (a). Now, we turn to the proof of part (b). First, we note that Therefore, because det ≡ 1, by using condition (c) of the RH problem for , we have ) as → * ∈ { + , − }, and we conclude from (84) that the functions are also analytic in  * . On the other hand, by using the asymptotics ( ) = 4 + ( −1 ) as → ∞ together with the fact that det ≡ 1, we can obtain asymptotics for −1 ( ) as → ∞ using (64). After some simple computations, we get from which it follows that ∈ (−( * + ∞ (1) + ∞ (2) )), = 1, 2, 3, 4.
We conclude the proof of part (b) as in part (a), by noting that   ( , ) in (83) has no pole at = (on any sheet). Finally, let us take ( ) = ( ) 1 = ( ) By multiplying the above from the right by , we obtain We denote (1) and (2) for the projections of on the first and second sheets of , respectively. Using (81), the above can be rewritten as for every ∈ ℂ and for any ∈ {1, 2}. The two integrals combine to one integral over a contour  surrounding both 0 (1) and 0 (2) on  in the positive direction, and thus for every ∈  * . Let us now take ( ) = ( ) 2 = ( ) Let ∶= { ∶ ( ) = 1 ( ) + 2 ( )( ( ) − 2 − ) with 1 , 2 two polynomials of degree ≤ − 1}. Because ↦ − 2 − has a simple pole at ∞ (2) (and no other poles), we conclude that ⊆ . Note also that dim = dim = 2 , and thus we have = . This finishes the proof. ■

The reproducing kernel 
To ease the notations, we define = ( ) and = ( ) by with the same convention as in (74), that is, (resp. ) is on the first sheet if | | > 1 (resp. | | > 1), and on the second sheet if | | < 1 (resp. | | < 1). We define  in terms of   as follows: Proposition 6. Let and be defined as in (24).  is a bivariate polynomial of degree ≤ 2 − 1 in both and . It satisfies for every scalar polynomial of degree ≤ 2 − 1, where ℂ is a closed curve in the complex plane going around and −1 once in the positive direction, but not going around 0.

Proposition 7.
The reproducing kernel  defined by (88) can be rewritten in terms of as follows: Proof. By Ref. 20, lemma 4.6(c), there is a unique bivariate polynomial  of degree ≤ 2 − 1 in both and which satisfies (89). Therefore, it suffices to first replace  in the left-hand side of (89) by ) , and then to verify that (89) still holds with this replacement. The rest of the proof goes exactly as in the proof of Ref. 20, proposition 4.9, so we omit it. ■

Proof of formula (29)
Now, using the results of Sections 6.1, 6.2, and 6.3, we give a proof for formula (29). From (21) and (22), for ∈ {1, … , 2 − 1}, ∈ ℤ, and ∈ {0, 1}, we have where is a closed contour surrounding 0 once in the positive direction. The proof consists of using the successive transformations  ↦  ↦   ↦  . We first use the eigendecomposition 61 of and the  ↦  transformation given in (68) to rewrite (91) as Using (60), we can write ( ) and ( ) −1 as where we have also used the relation to obtain (94). The identities (93) and (94) allow to rewrite the integrand of (92) by noting that ( 1 Therefore, using also the  ↦   transformation given by (81), we obtain [ (2 + , 2 + , 2 + , 2 + )] where  is a closed contour surrounding once 0 (1) and 0 (2) on  in the positive direction. By performing the change of variables = ( ) and = ( ) as in (86) and (87), using the factorization (75) and (77), the identity (76), and also the   ↦  transformation given by (88), we get where ℂ is a closed curve surrounding and −1 once in the positive direction, such that it does not surround 0. On the other hand, using again (75) and (77), we have By using the definition (24) Using the identities it is a simple computation to verify that (97) can be rewritten as (30) and (31). This finishes the proof.

LOZENGE PROBABILITIES
This section is about the lozenge probabilities ( , ), = 1, 2, 3, defined in (37). In Section 7.1, we use Theorem 1 to find double contour formulas for ( , ), = 1, 2, 3, in terms of  . In the rest of this section, we follow Ref. 21, section 7, and use the symmetries in our model to restrict our attention to the lower left part ≤ 2 ≤ 0 of the liquid region for the proof of Theorem 2.

Double contour formulas
Formula (29) for the kernel can be rewritten as where and̃are given by Proof. This is an immediate consequence of (3) and (100). ■ The next lemma establishes a double integral formula for the expectation value of the height function.

Symmetries
Let ( , ) be a 2 × 2 meromorphic function in both and , whose only possible poles in each variable are at 0, , −1 , , and −1 . Furthermore, we assume that all the poles of are of order 1 and that ( , ) is bounded as and/or tend to ∞. For ∈ {1, 2, … , 2 − 1} and ∈ ℤ, we define Because the poles of are of order at most 1, recalling (24), the only poles of the integrand are at 0, , and −1 , in both the and variables. The following star-operation will play an important role for a symmetry property of : Let 1 be the circle centered at −1 of radius 1 . The star-operation maps 1 into itself, but reverses the orientation. Furthermore, it satisfies ( * ) * = for all ∈ ℂ ∪ {∞}. We start by proving some symmetries for  .
from which we deduce (118). Now we prove (b). Note that the first column of only contains polynomials, which are independent of the choice of the contour ℂ that appears in the formulation of the RH problem for . Therefore,  is independent of the choice of ℂ as well by (120). Because 1 encloses both and −1 , and does not enclose 0, 1 is a valid choice of contour. We use the freedom we have in the choice of ℂ by letting be the solution to the RH problem for associated to the contour 1 . We can verify by direct computations that Recalling also (119) and (121), (124) follows by deforming back 1 to the original contour ℂ (in each variable). ■ We recall that ( , ; ) is defined for ( , ) ∈  as the unique solution of (41) lying in the upper half-plane, and that  is defined by (32). These quantities will appear naturally in the analysis of the next sections. For now, we simply note the following symmetries for ( , ; ).
where * denotes the star-operation defined in (117).
Proof. The symmetry (128) is part of Proposition 4 and has already been proved in Section 4. It remains to prove (129). We define the function as follows: so that (41) can be rewritten as Note that both and  depend on , even though this is not indicated in the notation. It is a long but direct computation to verify that Because the star operation maps the upper half-plane to the lower half-plane, ( , ; ) * lies in the lower half-plane. Therefore, applying the conjugate operation in (132), and noting that ( ) = ( ) and ( ) = ( ), we infer that ( , − ) ∈  if and only if ( , ) ∈  , and that (129) holds. ■

-FUNCTION
In Section 9, we will perform a Deift/Zhou 22 steepest descent analysis on the RH problem for . The first transformation ↦ consists of normalizing the RH problem and requires considerable preparation. This transformation uses a so-called -function, 43 which is of the form where is a probability measure, is its density, and supp is its (bounded and oriented) support. For any choice of , the -function satisfies −1)), in the sense that ( ) −2 ( ) 3 = 2 + ( −1 ) as → ∞. Also, we note that in the definition of , the contour ℂ can be chosen arbitrarily, as long as it is a closed curve surrounding −1 and once in the positive direction, which does not surround 0. However, to successfully perform an asymptotic analysis on the RH problem for , we need to choose and ℂ appropriately so that the jumps for have "good properties." In this section, we find the key ingredients for the ↦ transformation of Section 9, that is, we find a -function (built in terms of ) and a relevant contour ℂ . Let us rewrite as follows: where the potential is given by and we take the principal branch for the logarithms. We require and ℂ to satisfy the following criteria (we define supp( ) as an oriented open set for convenience): (a) ℂ is a closed curve surrounding −1 and once in the positive direction, but not surrounding 0.
In approximation theory, the equality (143) together with the inequality (144) are usually referred to as the Euler-Lagrange variational conditions, 45 and is the Euler-Lagrange constant. A measure satisfying (143) and (144) is called the equilibrium measure 45 in the external field , because it is the unique minimizer of among all probability measures̃supported on supp( ). Here, we require in addition that (145) is satisfied. This extra condition characterizes supp( ) as a so-called -curve. 46-51

Definition of  and related computations
By taking the derivative in (143), we have and by condition (b), ′ is analytic in ℂ ⧵ supp( ). Therefore, the function is meromorphic on ℂ. By (142), we get from which we conclude that  has a double zero at ∞, and double poles at 0, , −1 , , and −1 . Because a meromorphic function on the Riemann sphere (genus 0) has as many poles as zeros,  has eight other zeros. As → ∞, we have ′ ( ) = −1 + ( −2 ), from which we get ( ) = 2 −2 −2 + ( −3 ). Therefore,  can be written in the form where Π is a monic polynomial of degree 8 which remains to be determined. If we assume that ′ ( ) remains bounded for ∈ ℂ, then we can deduce from (147) and (148) the leading order term for ( ) as → * ∈ {0, , −1 , , −1 }: By combining these asymptotics with (149), we get Π(0) = 4 , Π ( ) = ( 1 − ) 8  This gives five linear equations for the eight unknown coefficients of Π, which is not enough to determine Π (and hence, ). Therefore, one needs to make a further assumption: we assume that we can find Π in the form As we will see, Assumption (156) implies that supp( ) consists of a single curve ("one-cut regime"). This assumption is justified if we can: (1) find 1 , 2 , 3 , + , − so that (155) holds and (2) construct a -function via (147), which satisfies the properties (a), (b), and (c). Substituting (156) in (155), we obtain five nonlinear equations for the five unknowns 1 , 2 , 3 , + , − . This system turns out to have quite a few solutions-we need to select "the correct one." Let us define 1 , 2 , 3 , + , − by (33) and (34). It is a simple computation to verify that indeed (155) holds in this case. We will show in Section 8.4 that this definition of 1 , 2 , 3 , + , − is "the correct solution" to (155), in the sense that it allows to construct a -function satisfying the properties (a), (b), and (c).
Remark 5. Let us briefly comment on how to find (33) and (34). Unfortunately, we were not able to solve analytically the nonlinear system obtained after substituting (156) into (155). Instead, we have solved numerically (using the Newton-Raphson method) this system for a large number of values of ∈ (0, 1). As already mentioned, the system (155) possesses several solutions. To ensure numerical convergence to "the correct solution," we choose starting values of 1 , 2 , and 3 so that (35) holds. The expressions (33) and (34)

Critical trajectories of 
In this subsection, we study the critical trajectories of , which are relevant to define thefunction and study its properties. Let ↦ ( ), ∈ [ , ] be a smooth parameterization of a curve , satisfying ′ ( ) ≠ 0 for all ∈ ( , ). is a trajectory of the quadratic differential ( ) 2 if ( ( )) ′ ( ) 2 < 0 for every ∈ ( , ), and an orthogonal trajectory if ( ( )) ′ ( ) 2 > 0 for every ∈ ( , ). is critical if it contains a zero or a pole of . Note that these definitions are independent of the choice of the parameterization.
Because + and − are simple zeros of , there are three critical trajectories (and also three orthogonal critical trajectories) emanating from each of the points ± . Recall the definitions of 0 , , 1 , Σ 0 , Σ , and Σ 1 given in Section 3.4.
We verify by direct computations that and thus the right-hand side of (158) is negative for ∈ (− 1 , 1 ), positive for ∈ (− , − 1 ) ∪ ( 1 , ), and zero for = − , − 1 , 1 , . We conclude that Σ 1 is a critical trajectory and that 1 ⧵ Σ 1 is the union of two orthogonal critical trajectories. The statement about Σ , ⧵ Σ , Σ 0 , 0 ⧵ Σ 0 can be proved a similar way, and we provide less details. For = ( ) = 0 , ∈ [− , ], after long but straightforward computations, we obtain Because Therefore, we deduce from an inspection of (160) that Σ is a critical trajectory and that ⧵ Σ is the union of two orthogonal critical trajectories. This finishes the proof. ■

Branch cut structure and the zero set of Re
As can be seen in (147), ′ can be expressed as for a certain branch of ( ) 1∕2 . To obtain a -function with the desired properties (a), (b), and (c), it turns out that the branch cut needs to be taken along the critical trajectory Σ 1 (as in Section 3.4).

Definition 3.
We define  1∕2 as where the branch cut for √ ( − + )( − − ) is taken on Σ 1 such that It will also be convenient to define a primitive of  1∕2 .
where the path of integration does not intersect (−∞, −1 We first state some basic properties of . By (150)-(154),  1∕2 has simple poles at 0, , −1 , , and −1 , and the residues are real. Also, because Σ 1 is a critical trajectory of , we have ± ( ) ∈ ℝ for ∈ Σ 1 . Therefore, Re is single-valued and continuous in ℂ ⧵ In the rest of this subsection, we determine the zero set  of Re . This will be useful in Section 8.4 to establish the (a), (b), and (c) properties of the -function. Let us define (165) Lemma 9. We have In particular,  divides the complex plane in three regions. The sign of Re in these regions is as shown in Figure 10.
Thus,  divides the complex plane in three regions in which Re does not change sign. The signs in each of these regions is then determined immediately by (164) (or equivalently, by (169)). ■

Definition and properties of
Definition 5. We define the measure by where Σ 1 = supp( ) is given by (44), and is oriented from − to + ; so  − ( ) 1∕2 denotes the limit of ( ) 1∕2 as → ∈ Σ 1 with in the exterior of the circle 1 .
The next proposition shows, among other things, that Definition 6 for is consistent with (161), and that satisfies (143).

Proposition 14. The functions and are related by
and we have Proof. We first prove (179). For a fixed ∈ ℂ ⧵ Σ 1 , we have where  is a closed curve surrounding Σ 1 once in the positive direction, but not surrounding any of the poles of , and not surrounding . By deforming  into another contour surrounding 0, , −1 , , −1 , and , we obtain where  = {0, , −1 , , −1 }. By deforming to ∞, noting that ( ) 1∕2 = ( −1 ) as → ∞, the integral on the right-hand side of (182) is 0. The sum can be evaluated using the residues (173), and we get Using (148) and ′ =  1∕2 , the above can be rewritten as Integrating this identity from + to along a path that does not intersect (−∞, − which is strictly decreasing as goes from − to + . So (145) holds as well, and hence (c), which finishes the proof. ■

STEEPEST DESCENT FOR
In this section, we will perform an asymptotic analysis of the RH problem for as → +∞, by means of the Deift/Zhou steepest descent method. 22 As mentioned in Section 8, the relevant contour to consider for the RH problem for is ℂ = 1 . The analysis is split in a series of transformations ↦ ↦ ↦ . The first transformation ↦ of Section 9.1 uses the -function obtained in Section 8 to normalize the RH problem at ∞. The opening of the lenses ↦ is realized in Section 9.2. The last step ↦ requires some preparations that are done in Section 9.3: it consists of constructing approximations (called "parametrices") for in different regions of the complex plane. Finally, the ↦ transformation is carried out in Section 9.4.

First transformation: ↦
We normalize the RH problem with the following transformation: where and are defined in Definition 6. Using (179), we can write the jumps for in terms of the function of Definition 4. From (180) and (181), we find that satisfies the following RH problem.
The following estimates for will be important for the saddle point analysis of Section 11.

Parametrices
In this subsection, we find good approximations to in different regions of the complex plane. By Lemma 9, Re ( ) > 0 for ∈ + ∪ − , Re ( ) < 0 for ∈ 1 ⧵ Σ 1 , and Re ( ) = 0 for ∈ Σ 1 . So the jumps for on + ∪ − ∪ ( 1 ⧵ Σ 1 ) are exponentially close to the identity matrix matrix as → ∞, uniformly outside fixed neighborhoods of − and + . By ignoring these jumps, we are left with the following RH problem, whose solution is denoted as (∞) . We will show in Section 9.4 that (∞) is a good approximation to away from + and − .
The condition on the behavior of (∞) ( ) as → * ∈ { + , − } has been added to ensure existence of a solution. This RH problem is independent of , and its unique solution is given by where ( ) ∶= ( − + − − ) 1∕4 is analytic in ℂ ⧵ Σ 1 and such that ( ) ∼ 1 as → ∞.
Note that (∞) is not a good approximation to in small neighborhoods of + , − ; this can be seen from the behaviors Let > 0 in Proposition 15 be fixed, and let  + and  − be small open disks of radius ∕2 centered at + and − , respectively. We now construct local approximations ( + ) and ( − ) (called "local parametrices") to in  + and  − , respectively. We require ( ± ) to satisfy the same jumps as inside  ± , to remain bounded as → ± , and to satisfy the matching condition uniformly for ∈  ± . The density of vanishes like a square root at the endpoints + and − , and therefore ( ± ) can be built in terms of Airy functions. 52 These constructions are well known and standard, so we do not give the details. What is important for us is that uniformly for ∈  ± .

PHASE FUNCTIONS AND
In Section 11, we will prove Proposition 11 via a saddle point analysis of the double contour integral (116). As it will turn out, the dominant part of the integrand as → +∞ will be in the form 2 (Φ( ; , )−Φ( ; , )) , for a certain function Φ which is described below. The analytic continuation of Φ to the second sheet of  is denoted Ψ-it will also play a role in the saddle point analysis and is presented below. The content of this section is a preparation for the saddle point analysis of Section 11. We will study the level set and also find the relevant contour deformations to consider.

Preliminaries
We start with a definition.
In the formulas that will be used in Section 11, Φ and Ψ will always appear in the form and similarly   Proof. This is an immediate consequence of Propositions 4 and 5. ■

The level set 
We study the set in case ≤ 2 < 0. We have represented  Φ for different values of ( , , ) in Figures 12, 13, and 14.
We show with the next two lemmas that the set  Φ ∩ (Σ ∪ Σ 1 ) ∩ ℂ + is either the empty set or a singleton.
For ∈ ℂ ⧵ {0, , −1 , , −1 }, we define the following functions: where 1 is a closed curve surrounding either or −1 , such that int( 1 ) ∩  Φ = ∅, and 2 is a closed curve surrounding either or −1 , such that int( 2 ) ∩  Φ = ∅. Because  Φ intersects both (0, ) and ( −1 , ) exactly once, 1 surrounds and 2 surrounds −1 . Then, the max/min principle implies that Re Φ is constant on int(Γ 3 ∪ Γ 4 ) ⧵ int( 2 ), which is a contradiction. The four other cases than (206) can be treated similarly, so we omit the proofs. ■ The set  Φ is represented in blue, and Σ ∪ Σ 1 in red. The parameters are ( , , ) = (0.4, −0.12, −0.86) as in Figure 12 (left). The contour * is represented in green, and * in black. The black dots represent 0, , −1 , , and −1 and the blue dots are and Lemma 12 states that Γ 4 crosses Σ ∪ Σ 1 exactly once. We know from (203) that Re Φ( ) → +∞ as → ∞, so Γ 4 intersects the real line, and then by symmetry ends at . So each of the Γ 's intersects ℝ. We denote for the intersection point of Γ with ℝ, and choose the ordering such that , we know from Proposition 4 (b) that lies on 1 ⧵ Σ 1 . For the saddle point analysis, we will need * lying inside 1 (not necessarily strictly inside). To prove existence of such a contour * , we need to know that Re Φ( ) − Re Φ( ) is strictly negative for ∈ 1 ⧵ Σ 1 (at least in small neighborhoods of and ).

SADDLE POINT ANALYSIS
In this section, we prove Proposition 11 by means of a saddle point analysis that mainly follows the lines of Ref. 21. This analysis relies mostly on Sections 9 and 10 and is only valid for ( , ) in the lower left part of the liquid region, that is for ( , ) ∈  ∩ { ≤ 2 ≤ 0}. We divide the proof in three subcases: ≤ 2 < 0, < 2 = 0, and = = 0.
Remark 6. By adapting the analysis of this section and of Section 10, it is possible to carry out similar saddle point analysis when ( , ) lies in the other quadrants of the liquid region. Note however that this is not needed, thanks to the symmetries of Section 7.2 (see also Proposition 11).

The case ≤ <
The double integral  is defined in (116). The associated two contours of integration can be chosen freely, as long as they are closed curves surrounding and −1 once in the positive direction, and not surrounding 0. From now, it will be convenient to take different contours in the and variables, so we indicate this in the notation by rewriting (116) Only the first column of appears in (210), which is independent of the choice of the contour ℂ associated to the RH problem for . However, by using the jumps for , we will find (just below) another formula for  in terms of the second column of . Therefore, the choice of ℂ will matter. To be able to use the steepest descent of Section 9, we assume from now that ℂ = 1 . Recall that is expressed in terms of via (183), and definẽ By Proposition 15, ( , ) is uniformly bounded as and stay bounded away from + and − . We will need the analytic continuation in of ( , ) from the interior of 1 to the bounded region delimited by Σ 1 ∪ Σ (see Figure 10). We denote it , ( , ), and by (185)  ) , ∈ int(( 1 ⧵ Σ 1 ) ∪ Σ ), ∈ ℂ ⧵ 1 . (212) By Lemma 9, Re ( ) < 0 for ∈ int(( 1 ⧵ Σ 1 ) ∪ Σ ), so , ( , ) remains bounded as → +∞, uniformly for and bounded away from + and − , as long as ∈ int(Σ 1 ∪ Σ ). Our next goal is to prove the following. . We do not indicate whether we take the + or − boundary values in the integrand of (213). This is without ambiguity, because has no jumps on 1 (this can be verified using (184) and (185)).
For ∈ Σ 1 , we use (221) to split the integrand into two parts, and again we deform the integral associated to the first term slightly inward, and the other one slightly outward. As a result, both deformed integrals have exponentially decaying integrands.
For ∈ Σ , , ( , ) is given by the second line of (212), and thus the dominant -part in the integrand is For the first term, we deform Σ outward so that Re ( ) > 0, and for the first term, we deform Σ inward so that Re ( ) < 0.
On the deformed contours, the integrand is uniformly exponentially small, as long as and are bounded away from + and − . For and close to ± , by Proposition 15 we have ( ) = ( 1∕6 ) and −1 ( ) = ( 1∕6 ). The contribution to (223) when and are close to + and − is thus bounded by for certain 1 , 2 , 3 > 0 and for all large enough . This finishes the proof of Proposition 11.

A C K N O W L E D G M E N T S
The first step of the method is based on an unpublished idea of Arno Kuijlaars. I am very grateful to him for allowing me and even encouraging me to use his idea. I also thank Jonatan Lenells, Arno