The two periodic Aztec diamond and matrix valued orthogonal polynomials

We analyze domino tilings of the two-periodic Aztec diamond by means of matrix valued orthogonal polynomials that we obtain from a reformulation of the Aztec diamond as a non-intersecting path model with periodic transition matrices. In a more general framework we express the correlation kernel for the underlying determinantal point process as a double contour integral that contains the reproducing kernel of matrix valued orthogonal polynomials. We use the Riemann-Hilbert problem to simplify this formula for the case of the two-periodic Aztec diamond. In the large size limit we recover the three phases of the model known as solid, liquid and gas. We describe fine asymptotics for the gas phase and at the cusp points of the liquid-gas boundary, thereby complementing and extending results of Chhita and Johansson.


Introduction
We study domino tilings of the Aztec diamond with a two periodic weighting. This model falls into a class of models for which existing techniques for studying fine asymptotics are not adequate and only recently first important progress has been made [9,20]. We introduce a new approach based on matrix valued orthogonal polynomials that allows us to compute the determinantal correlations at finite size and their asymptotics as the size of the diamond gets large in a rather orderly way. We strongly believe that this approach will also prove to be a good starting point for other tiling models with a periodic weighting.
Random tilings of planar domains have been studied intensively in the past decade. Such models exhibit a rich structure including a limit shape and fluctuations that are expected to fall in various important universality classes (see [8,22,23,43,45,46,47,48] for general references on the topic and [18,58,59] for recent contributions). When the correlation structure is determinantal, there is hope to understand the fine asymptotic structure by studying the asymptotic behavior of the correlation kernel. An important source of examples of such models is the Schur process [56]. For these models the correlation kernel can be explicitly computed in terms of a double integral representation, opening up to the possibility of performing an asymptotic analysis by means of classical steepest decent (or stationary phase) techniques.
Of course, the Schur process is rather special and many models of interest fall outside this class. In particular this is true for random tilings or dimer models with doubly periodic weightings. Yet, these models have exciting new features and have therefore been discussed in the physics and mathematics literature [9,20,21,36,48,55]. An important feature is the appearance of a so-called gas phase. For instance, in the two periodic weighting for domino tiling of the Aztec diamond, the diamond can be partitioned into three regions: the solid, liquid and gas region [56] (as we will see in Figure 7 below). The gas region has not been observed in models that are in the Schur class. The 2-point correlations (for an associated particle process) in the gas region behave differently when compared to the liquid regions. Indeed, the correlation kernel decays exponentially with the distance d between the points, instead of ∼ 1/d in the liquid region. At the liquid-solid boundary one expects the Airy process to appear, but the situation at the gas-liquid boundary is far more complicated [9,20].
To the best of our knowledge, the two periodic Aztec diamond is the only model with periodic weightings for which rigorous results on fine asymptotics exist [9,20]. Inspired by a formula for the Kasteleyn matrix found by Chhita and Young [21], Chhita and Johansson [20] found a way to compute the asymptotic behavior of the Kasteleyn matrix as the size of the Aztec diamond goes to infinity. We will follow a different approach to studying such models with periodic weightings.
As we will recall in Section 3, the Aztec diamond can be described by non-intersecting paths (we refer to [45] and the references therein for more background on the relation between dimers, tilings, non-intersecting paths and all that). For a general class of discrete non-intersecting paths with p-periodic transition matrices (which includes p-periodic weightings for domino tilings of the Aztec diamond and p-periodic weightings for lozenge tilings of the hexagon), we show in Section 4 how the correlation kernel can be written as a double integral formula involving matrix valued polynomials that satisfy a non-hermitian orthogonality.
We believe that this general setup has a high potential for a rigorous asymptotic analysis. The key fact is that these matrix valued orthogonal polynomials can be characterized in terms of the solution of a 2p × 2p matrix valued Riemann-Hilbert problem. With the highly developed Riemann-Hilbert toolkit at hand, we may thus hope to compute the asymptotic behavior of the polynomials, and more importantly the correlation kernel. The formalism will be worked out in Section 4. It provides a new perspective even on the classical examples of uniform domino tilings of the Aztec diamond and lozenge tilings of a hexagon, as we will discuss briefly in Sections 4.7 and 4.8.
The main focus of the paper is to show how the Riemann-Hilbert approach can be exploited to find an asymptotic analysis for the two periodic Aztec diamond. Remarkably, in this case the result of the Riemann-Hilbert analysis is a surprisingly simple double integral formula for the correlation kernel. It is not an asymptotic result, but an exact formula valid for fixed finite N . This representation also appears to be more elementary than the one given in Chhita and Johansson [20]. The Riemann-Hilbert analysis is given in Section 5. We analyze the double integral formula for the kernel asymptotically using classical steepest descent techniques in Section 6.
The model of the two periodic Aztec diamond is explained and the main results are summarized in the next section.

Statement of results
In this section we will introduce the two periodic Aztec diamond and state our main results for this model.

Definition of the model
The Aztec diamond is a region on the square lattice with a sawtooth boundary that can be covered by 2 × 1 and 1 × 2 rectangles, called dominos. The squares have a black/white checkerboard coloring and a possible tiling of the Aztec diamond of size 4 is shown in Figure 1 dominos, namely North, West, East, and South, that are also shown in the figure. The Aztec diamond model was first introduced in [33].
In the two periodic Aztec diamond we assign a weight to each domino in a tiling, depending on its shape (horizontal or vertical) and its location in the Aztec diamond. We assume the Aztec diamond is of even size.
To describe the two periodic weighting we introduce a coordinate system where (0, 0) is at the center of the Aztec diamond. The center of a horizontal domino has coordinates (x, y + 1/2) with x, y ∈ Z. We then say that the horizontal domino is in column x. The center of a vertical domino has coordinates (x + 1/2, y) with x, y ∈ Z, and we say that the vertical domino is in row y. The row and column numbers run from −N + 1 to N − 1, where 2N is the size of the Aztec diamond.
We fix two positive numbers a and b and define the weights as follows.
Definition 2.1. The weight of a domino D in a tiling T of the Aztec diamond is and the probability for T is where Z N = T w(T ) (sum over all possible tilings T of an Aztec diamond of size 2N ) is the partition function.
In the example from Figure 1 the weights are shown in Figure 2. The weight of the tiling is w(T ) = a 10 b 10 .
The model is homogeneous in the sense that the probabilities (2.3) do not change if we multiply a and b by a common factor. We may and do assume ab = 1. In what follows it will be more convenient to work with α = a 2 , β = b 2 (2.4) instead of a and b. We have αβ = 1, and without loss of generality we assume α ≥ 1. If α = β = 1 then the model reduces to the uniform weighting on domino tilings, and so the true interest is in the case α > 1, and this is what we assume from now on.

Particle system and determinantal point process
By putting a particle in the black square of the West and South dominos, we obtain a random particle system. In our running example the particle systems is shown in Figure 3. Any possible tiling of the Aztec diamond gives rise to a subset X = X (T ) of B N containing the squares that are occupied by a particle. We use (m, n) ∈ B N to denote an element B N and we will refer to m as the level in B N . Any X that comes from a tiling will have 2N − m particles at level m for each m = 0, 1, . . . , 2N . Therefore the cardinality is |X | = N (2N + 1).
There are also interlacing conditions that are satisfied when comparing the particles at level m with those at level m + 1.
The probability measure (2.3) on tilings gives rise to a probability measure on subsets X , that turns out to be determinantal. This means that there exists a kernel K N : with the property that for any subset S ⊂ B N This is a discrete determinantal point process [15]. We found an explicit double contour integral formula for the kernel K N . We take (m, n), (m , n ) ∈ B N , and instead of K N ((m, n), (m , n )) we write K N (m, n; m , n ). We collect K N (m, n; m , n ) with some of its neighbors in a 2 × 2 matrix K N (m, n; m , n ) = K N (m, n; m , n ) K N (m, n + 1; m , n ) K N (m, n; m , n + 1) K N (m, n + 1; m , n + 1) (2.7) and this matrix appears in our formula (2.8).
Remark 2.3. The square root factor in (2.10) is defined and analytic for z ∈ C \ ((−∞, −α 2 ] ∪ [−β 2 , 0]) with the branch that is positive for real z > 0. This is one sheet of the Riemann surface R associated with the cubic equation y 2 = z(z + α 2 )(z + β 2 ), (2.11) that will play an important role in what follows. It is a two sheeted surface consisting of two sheets C \ ((−∞, −α 2 ] ∪ [−β 2 , 0]) glued together along the two cuts (−∞, −α 2 ] and [−β 2 , 0] in the usual crosswise manner. The surface has genus 1 unless α = β = 1 in which case the genus drops to 0. The matrix valued function F (z) from (2.10) is considered on the first sheet. Its analytic continuation to the second sheet is given by I 2 − F (z). Remark 2.4. The eigenvalues of A(z), see (2.9), are equal to are the eigenvalues of 2αz α(z + 1) βz(z + 1) 2βz . Thus (2.12) are the two branches of the meromorphic function ρ = (α + β)z + y on the Riemann surface R associated with (2.11), with ρ 1 on the first sheet and ρ 2 on the second sheet. The eigenvectors can also be considered on the Riemann surface. There is a matrix E(z) whose columns are the eigenvectors such that See (5.6) below for the precise formula for E(z). It turns out that (see (2.10) for the definition of F (z)) (2.14) Thus F (z) has eigenvalues 0 and 1, and F (z) commutes with A(z). We will also work with which in view of (2.13) has eigenvalue decomposition (2.17) These eigenvalues are the two branches of a meromorphic function λ on the Riemann surface R, with λ 1 defined on the first sheet and λ 2 on the second sheet.
Remark 2.5. The only singularity for the w integral in (2.8) that is inside the contour γ 1 is the pole at w = 1. Because of (2.9) we see that A N −m (w) has a pole of order N −m (it is a zero if m > N ) and therefore the integrand in the double integral of (2.8) has a pole of order N − m + N = 2N − m at w = 1. There is no pole if m = 2N , and thus it follows that (2.8) vanishes identically for m = 2N . This is in agreement with the fact that there are no particles at level 2N . Level 0 is full of particles, which means that K N (0, n; 0, n) = 1 for all n = 0, . . . , 2N − 1. This can be seen from formula (2.8) as follows. For n even, the formula gives for K N (0, n; 0, n), In view of (2.14), (2.15), and (2.16) we have We will see in Lemma 5.2 (b) that λ 1 (w) has a double pole and λ 2 (w) has a double zero at w = 1. There are no other zeros and poles. Thus (2.19) has a pole of order N at w = 1, and (2.20) has a zero of order N at w = 1. It follows that the w-integrand in (2.18) has a pole of order 2N at w = 1. However, if we replace F (w) by I 2 − F (w), then the integrand does not have a pole at w = 1 anymore, and thus by Cauchy's theorem the double integral is zero. It follows that the value of the double integral (2.18) remains the same if we remove the factor F (w). Then the integrand that remains is rational in w, with a simple pole at w = z and a decay O(w −1−N +n/2 ) as w → ∞. We apply the residue theorem to the exterior of γ 1 and the only contribution comes from the pole at w = z. The z-integral that remains in (2.18) reduces to 1 2πi γ 0,1 dz z I 2 = I 2 and so by looking at the diagonal entries, see (2.7), we get K N (0, n; 0, n) = K N (0, n + 1; 0, n + 1) = 1, as claimed.
Remark 2.6. Observe that the particle system on the black squares is not equivalent to the tiling of the Aztec diamond. However, it can be extended to a larger system where we also also put a particle in the white square of Figure 5: Extended particle system that is equivalent to a domino tiling of the Aztec diamond.
the West and South dominos, as shown in Figure 5. This extended process is equivalent to the tiling. It is also determinantal and we will now give a double contour integral formula for the kernel. We chose to work with the particle system determined by the black squares only since the formulas and their asymptotics stated in the next paragraphs take an easier form. The analogous asymptotics for the extended processes are straightforward from these results and the formula (2.25) below. So in addition to the set of black squares (2.5) we also consider the white squares Then the extended particle system has a correlation kernel (that we continue to denote by K N ) To write down the extended kernel we use the notation that will also be used later in (3.3) and (3.5), namely, for an integer m ∈ Z, and for m, m ∈ 1 2 Z, Then the extended kernel (2.22) is as follows. We assume N is even (m, n) ∈ B N ∪ W N and (m , n ) ∈ B N ∪ W N with both m + n and m + n even. Then K N (m, n; m , n ) K N (m, n + 1; m , n ) K N (m, n; m , n + 1) K N (m, n + 1; m , n + 1) Observe that by (2.9) and (2.22) and so by the definition (2.24), we have

Matrix valued orthogonal polynomials
The starting point of our approach to Theorem 2.2 is the non-intersecting path reformulation of the Aztec diamond and the Lindström-Gessel-Viennot lemma. This will be developed in Section 3. The novel ingredient in the further analysis is the use of matrix valued orthogonal polynomials (MVOP). A matrix valued polynomial of degree k and size d is a function matrix on a set Γ in the complex plane.
Definition 2.7. Suppose P N is a matrix valued polynomial of degree N with an invertible leading coefficient. Then P N is a matrix valued orthogonal polynomial with weight matrix w on Γ if for all matrix polynomials Q of degree ≤ N −1, where Q t denotes the matrix transpose.
The integral in (2.26) is to be taken entrywise, and 0 d denotes the d × d zero matrix. We note that the order of the factors in the integrand in (2.26) is important since we are dealing with matrices.
For us the weight matrix will be 1 2πi W N (z) on the closed contour γ 1 around 1. Thus d = 2, and W N denotes the N th power of W , so that the weight matrix is varying with N . Recall that W is defined in (2.15), and explicitly we have The existing literature on MVOP mostly deals with the case of orthogonality on an interval of the real line, with a positive definite weight matrix w with all existing moments. In such a case the MVOP exists for every degree n, and they can be normalized in such a way that However, it is interesting to note that MVOP first appeared in connection with prediction theory where the orthogonality is on the unit circle, see [11] for a recent survey. The interest in MVOP on the real line has been steadily growing since the early 1990s. The analytic theory of MVOP on the real line is surveyed in [25] with [6] as one of the pioneering works. MVOP satisfy recurrence relations [32] and special cases satisfy differential equations [31].
Interesting examples of MVOP come from matrix valued spherical functions, see [39,49] as well as many other papers. We deviate from the usual set-up of MVOP in several ways • Γ is a closed contour in the complex plane, • the weight matrix w(z) = 1 2πi W N (z) is complex-valued on Γ and varies with N , • the weight matrix is not symmetric or Hermitian (let alone positive definite), or have any other property that would imply existence and uniqueness of the MVOP.
Since there is no complex conjugation in (2.26), we are thus dealing with non-Hermitian matrix valued orthogonality with varying weights on a closed contour in the plane. As already noted, existence and uniqueness of the MVOP are not guaranteed in this general setting. However, for the weight 1 2πi W N we can show that the monic MVOP up to degrees N all exist and are unique. However, since the weight matrix is not symmetric, we cannot normalize to obtain orthonormal MVOP as in (2.28). In our case the MVOP of degrees > N do not exist.
In Section 4 we consider a situation that is more general than the two periodic Aztec diamond. It deals with a multi-level particle system that is determinantal, and transitions between the levels are periodic. See Assumptions 4.1 and 4.2 for the precise assumptions. In this general setting we make a connection with matrix valued (bi)orthogonal polynomials and our main result in Section 4 is Theorem 4.7 that expresses the correlation kernel as a double contour integral containing a reproducing kernel for the matrix polynomials.
In the special situation of the two periodic Aztec diamond it gives that the matrix K N with the correlation kernels as in (2.8) is given by where R N (w, z) is the reproducing kernel associated with the matrix polynomials of degrees ≤ N − 1. That is, R N (w, z) is a bivariate matrix valued polynomial of degree N − 1 in both w and z, such that holds for every matrix valued polynomial Q of degree ≤ N − 1. The MVOP of degree N is characterized by a Riemann-Hilbert problem and the reproducing kernel R N (w, z) can be expressed in terms of the solution of the Riemann-Hilbert problem. This is known from work of Delvaux [29] and we recall it in Section 4.6. Then we perform an analysis of the Riemann-Hilbert problem, and quite remarkably this produces the exact formula (2.8).

Classification of phases
The explicit formula (2.8) in Theorem 2.2 is suitable for asymptotic analysis as N → ∞. See Figure 6 for a sampling of a large 2-periodic Aztec diamond. In this figure three regions emerge where the tiling appears to have different statistical behavior. We first describe how we can distinguish these three phases (solid, liquid and gas) in the model. This classification will depend on the location of saddle points for the double integral in (2.8).
We fix coordinates −1 < ξ 1 < 1 and −1 < ξ 2 < 1 and choose m, m ≈ (1 + ξ 1 )N and n, n ≈ (1 + ξ 2 )N . Then from the formula (2.8), we see that the z integral of the double contour integral is dominated as N → ∞ by the expression Figure 6: A sample of a two-periodic Aztec diamond of size 500. We colored the West and South dominos blue. The East and North dominons are colored yellow. The gas phase is visible in the middle. This figure is generated by a code that was kindly provided to us by Sunil Chhita.
where W is given by (2.15). In view of the eigenvalue decomposition (2.16) this is with Φ j (z) = 2 log(z − 1) − (1 + ξ 2 ) log z + ξ 1 log λ j (z), j = 1, 2. Hence we are led to consider as a function on the Riemann surface R, depending on parameters ξ 1 and ξ 2 . It is multi-valued, but its differential is a single-valued meromorphic differential with simple poles at z = 1 (1) , z = 1 (2) , z = 0 and z = ∞, see also Section 6.2. (For j = 1, 2, we use 1 (j) to denote the value 1 on the jth sheet of the Riemann surface.) There are also four zeros, counting multiplicities, since the genus is 1.
Definition 2.8. The saddle points are the zeros of Φ (z)dz.
The real part R r of the Riemann surface consists of all real tuples (z, y) satisfying the algebraic equation (2.11) together with the point at infinity. The real part is the union of two cycles, where C 1 is the union of the intervals [−α 2 , −β 2 ] on the two sheets, and C 2 is the union of the two intervals [0, ∞] on both sheets. It turns out that there are always at least two distinct saddle points on the cycle C 1 , see Proposition 6.4 below. The location of the other two saddle points determines the phase. Definition 2.9. Let −1 < ξ 1 , ξ 2 < 1.
(b) If two saddles are outside of the real part of the Riemann surface, then (ξ 1 , ξ 2 ) is in the liquid phase, and we write (ξ 1 , ξ 2 ) ∈ L .
(c) If all four saddles are simple and belong to C 1 , then (ξ 1 , ξ 2 ) is in the gas phase, and we write (ξ 1 , ξ 2 ) ∈ G.
Transitions between phases take place when two or more saddle points coalesce.
(d) If there is a double saddle point on C 2 , then (ξ 1 , ξ 2 ) is on the solidliquid transition.
(e) If there is a double or triple saddle point on C 1 then (ξ 1 , ξ 2 ) is on the liquid-gas transition.
It is not possible to have a double saddle point outside of the real part of the Riemann surface.
The condition for coalescing saddle points leads to an algebraic equation of degree 8 for ξ 1 , ξ 2 and it precisely coincides with the equation listed in the appendix of [20], see also (6.10) below.
The real section of the degree 8 algebraic equation has two components in case α > 1, as shown in Figure 7. Both components are contained in the square −1 ≤ ξ 1 , ξ 2 ≤ 1. The outer component is a smooth closed curve that touches the square in the points (±1, 0), and (0, ±1). It is the boundary between the solid and liquid phases.

Gas phase
Our next result gives the limit of K N in the gas phase. Recall that K N is defined by (2.7). Theorem 2.10. Assume (ξ 1 , ξ 2 ) ∈ G. Suppose m, m , n, n are integers that vary with N in such a way that are fixed. Also assume that m + n and m + n are even. Then we have for N even, lim N →∞ K N (m, n; m , n ) = K gas (m, n; m , n ) (2.35) with K gas (m, n; m , n ) = (which we can also view as a closed loop on the first sheet of the Riemann surface).
The limit (2.36) does not depend on ξ 1 and ξ 2 as long as these are in the gas region. The kernel K gas still depends on the parameters α and β (cf. Remark 2.12) and therefore the 2-periodic structure is still present in this limit.
The proof of Theorem 2.10 is in Section 6.5.
Remark 2.11. In (2.36) we can see the exponential decay of correlations that is characteristic for the gas phase as follows. We combine the factor z (∆m)/2 with A −∆m (z), since by (2.15) We find from this and (2.10) and (2.15) that since λ 2 = λ −1 1 , see part (d) of Lemma 5.2 below. Thus (2.36) can be written as K gas (m, n; m , n ) = (2.37) By analyticity and Cauchy's theorem, we have the freedom to deform the contour γ as long as it goes around the interval [−β 2 , 0] and does not intersect (−∞, −α 2 ]. Since β < 1 < α we can deform it to a circle centered at zero of radius < 1 or to a circle of radius > 1. If ∆n ≥ 0 we deform to a circle |z| = r < 1 and if ∆n < 0 we deform to a circle |z| = r > 1. In both cases the factor z (∆n)/2 is exponentially small as |∆n| → +∞.
Since |λ 1 | > 1 on γ (as will follow from parts (d) and (f) of Lemma 5.2 below) the factor λ ±(∆m)/2 1 (z) is also exponentially small as |∆m| → +∞. It follows that the gas kernel (2.37) decays exponentially as |∆m| + |∆n| → ∞. Remark 2.12. Let's see what we have for the gas kernel (2.36) on the diagonal, i.e., for ∆m = ∆n = 0. then we obtain from (2.36) and (2.10) if m + n is even . (2.38) The diagonal entries of (2.38) are K gas (m, n; m, n) and K gas (m, n + 1; m, n + 1), see also (2.7). Thus for arbitary parity of m + n, K gas (m, n; m, n) = 1 2 where (2.39) follows from deforming the contour γ to the interval [−β 2 , 0] and making a change of variables. It is easy to check that and so (2.39) is between 0 and 1, as it is the particle density of the gas phase. Note also that the density depends on the parameters α and β.

Cusp points
On the boundary between the gas phase and the liquid phase, the gas kernel (2.36) is still the dominant contribution. This phenomenon was already observed by Chhita and Johansson [20] and further investigated by Beffara, Chhita and Johansson [9], who looked at the diagonal point ξ 1 = ξ 2 on the boundary and proved that after averaging there is an Airy like behavior in the first subleading term. We consider the cusp points, and show explicitly the appearance of a Pearcey like behavior in the subleading term of the kernel K N .
The four cusp points are located at positions (ξ 1 , ξ 2 ) = ± α−β α+β , 0 and (ξ 1 , ξ 2 ) = 0, ± α−β α+β in the phase diagram. We focus on the top cusp point with coordinates (ξ * 1 , ξ * 2 ) = 0, α−β α+β . At this cusp point the triple saddle point is located at the branch point −α 2 .   40) as N → ∞, with fixed u, u , v, v and with constants Then, in case m is even, with contours Σ and −Σ as shown in Figure 8. In case m is odd, we have  is, up to a gaussian, known as the Pearcey kernel. It is one of the canonical kernels from random matrix theory that arises typically as a scaling limit near a cusp point. It was first described by Brézin and Hikami [17] in the context of random matrices with an external source, see also [13]. The Pearcey process was given in [57,62]. More recent contributions are for example [1,10,40]. Note that the actual Pearcey kernel includes a gaussian in addition to the double integral in (2.44). Remarkably, this term is hidden in the gas kernel and can be retrieved by a steepest descent analysis of that kernel.
Theorem 2.14. Under the same assumptions as in Theorem 2.13 we have It is very curious that the double integral part of the Pearcey kernel appears in the scaling limit at the cusp point, but only in the subleading term. A similar phenomenon was already observed in [20] on the smooth parts of the liquid-gas boundary. The gas phase is dominant with a subleading Airy behavior. Also here the gaussian part of the Airy kernel is hiding in the gas kernel [20, §3.2]. With some effort we can also find this from our approach.

Non-intersecting paths
We discuss the non-intersecting paths on the Aztec diamond. What follows in this section is not new, and can be found in several places, see e.g. [42,43,44] and the recent works [7,45]. Note however, that we use a (random) double Aztec diamond to extend the paths, see Section 3.2 below, instead of a deterministic extension in [42,43].

Non-intersecting paths
The South, West and East dominos are marked by line segments as shown in Figure 9. The North dominos have no marking. There are also particles on the West and South dominos, but these will only play a role later on. We look at the line segment as part of paths that go from left to right and either go up (in a West domino), down (in an East domino), or go horizontal (in a South domino). Each segment enters a domino in the black square, and exits it from the white square within the domino.
We include the marking of the dominos into the tiling we obtain nonintersecting paths, starting at the lower left side of the Aztec diamond and ending at the lower right side. In the pictures that follow we forget about the black/white shading of the dominos.

West
North South East Figure 9: Line segments and particles on the dominos, that lead to nonintersecting paths in a domino tiling of the Aztec diamond.
Each path ends at the same height as where it started. Thus along each path the number of West dominos is the same as the number of East dominos. Also each path has in each row the same number of West dominos as the number of East dominos. Since we need this property, we state it in a separate lemma.

Double Aztec diamond
The lengths of the paths vary greatly. To obtain a more symmetric picture, which will be useful for what follows, we attach to the right bottom side  In other words: there are no dominos that are partly in the original Aztec diamond and partly in the newly attached Aztec diamond.
Proof. The smaller Aztec diamond is attached to the original Aztec diamond along its south-east boundary. In the checkerboard coloring there are only white squares adjacent to this boundary as is seen in Figure 1.
If some dominos are partly in the original Aztec diamond and partly in the new one, then these dominos would cover a number of white squares in the original Aztec diamond, and no black squares. That would leave us with more black squares than white squares that should be covered by dominos. This is impossible, since each domino covers exactly one black and one white square.
We thus cover the second Aztec diamond with dominos, independently from what we did in the original Aztec diamond. We also put in the markings with line segments. Note however that a horizontal domino at the very bottom of the second Aztec diamond is now a "North domino". A possible tiling is shown in Figure 10 together with the corresponding non-intersecting paths and the particle system that we discuss in Section 3.3.
The double Aztec diamond with partial overlap is considered in [2,3], where the phenomenon of a tacnode is studied. For us, the two Aztec diamonds do not overlap and there is no tacnode phenomenon.

Particle system
Next we put particles on the paths. On each West and South domino we put a particle at the left-most part of its marking as already shown in Figure 9. We do not put any particle on the East and North dominos, although there may be a particle on the right edge of an East (or West or South) domino if it is connected to a West or South domino. This a slight deviation from what we did before. We now put the particles on the boundary of the West and South dominos and not in the center of the black square. We also put a particle at the end of each path, see Figure 10.
Each path contains 2N + 1 particles. The particles are interlacing if we consider them along diagonal lines going from north-west to south-east. There are 2N +1 diagonal lines and each diagonal line contains 2N particles.
To find a clearer picture, we forget about the dominos and we rotate the figure over 45 degrees in clockwise direction. We further introduce a shear transformation so that the starting and ending points of each path remain at the same height. In a formula: (x, y) is mapped to (x + y, y).
Then we obtain a figure as in Figure 11 where each vertical line contains 2N particles. The paths now consist of diagonal (right-up) parts, horizontal parts and vertical parts. The diagonal parts come from the West dominos. The horizontal parts come from the South dominos and the vertical parts come from the East dominos. Each path contains 2N + 1 particles.

Modified paths on a graph
We are going to modify the paths, in such a way that each particle is preceeded by a horizontal step of a half unit (except for the initial particles).
The particles are on the integer lattice Z × Z. We put coordinates so that the initial vertices are at (0, j) for j = 0, . . . , 2N − 1 and the ending vertices are at (2N, j) for j = 0, . . . , 2N − 1. Each path consists of 2N parts, where the mth part goes from (m − 1, k) to (m, l) for some k, l ∈ Z with m ≤ k + 1. We modify this part by an affine transformation that results in a path from (m − 1, k) to (m − 1/2, l), followed by a horizontal step from (m − 1/2, l) to (m, l).
We also extend the particle system by putting particles at the new vertical lines. If there is no vertical part, then the path has a unique intersection point with the vertical line, and we put a particle there. If there is a vertical Figure 11: Particle system on vertical lines part of the path at that level, then we put the particle at the highest point, see Figure 12. Now there are 4N + 1 particles on each path.
The new paths have a two-step structure. Starting from an initial position, we either move horizontally to the right by half a unit and stay at the same height, or we go diagonally up by one unit in the vertical direction and horizontally by half a unit. We call this a Bernoulli step. In both cases we end at a particle on the line with horizontal coordinate 1/2.
Then we make a number of vertical down steps followed by a horizontal step by half a unit to the right. The number of down steps can be any non-negative integer, including zero. We call this a geometric step.
Then we repeat the pattern. We do a Bernoulli step, a geometric step, a Bernoulli step, etc. The final step in each path is a geometric step, which should take us to the same height as where the path started.
Another requirement is that the resulting paths are non-intersecting. Any such path structure is in one-to-one correspondence with a unique domino tiling of the double Aztec diamond.
The geometric steps cannot be too far down, since each path has to return to its initial height, and each up step is done by one unit only. In the example, the largest geometric down step is by two units, but in a larger size example, one could imagine that larger steps are possible.
The paths lie on an infinite directed graph that is also shown in Figure  12. We call it the Aztec diamond graph. Its set of vertices is ( 1 2 Z) × Z. The transition from m to m + 1 2 is a Bernoulli step. From (m + 1 2 , n) there are also two directed edges • From (m + 1 2 , n) to (m + 1 2 , n − 1) (vertical down step) • From (m + 1 2 , n) to (m + 1, n) (horizontal step) We go from level m + 1 2 to level m + 1 by making a number of vertical down steps (possibly zero) and then making the horizontal step to the right. The transition from m + 1 2 to m + 1 is a geometric step. For a given N the paths start at the vertices with coordinates (0, j), j = 0, . . . , 2N − 1, and end at the vertices with coordinates (2N, j), j = 0, . . . , 2N − 1. The paths lie on the graph and are not allowed to intersect, that is, the set of vertices for two different paths is disjoint. Since the graph is planar and paths cannot go to the left, the paths maintain their relative ordering. The path from (0, j) to (2N, j) stays below the path from (0, j +1) to (2N, j + 1) at all levels.

Weights
There is a one-to-one correspondence between tilings of the double Aztec diamond and non-intersecting paths on the Aztec diamond graph with prescribed initial and ending positions as described above.
In the two periodic Aztec diamond we assign a weight (2.2) to a tiling with the corresponding probability (2.3). To be able to transfer this to the paths, we recall that in a Bernoulli step a diagonal up-step corresponds to a West domino, and a horizontal step corresponds to a South domino. The vertical steps in a geometric step correspond to East dominos. The horizontal step that closes a geometric step was added artificially and does not correspond to a domino. We do not see the North dominos in the paths, and therefore we cannot transfer the weights on the dominos to weights on path segments directly. It is possible to assign weights to dominos in which North dominos have weight one and which is equivalent to (2.1) as it leads to the same probabilities (2.3) on tilings. This was also done in [20], but we present it in a different way here. Because of Lemma 3.1 each column has the same number of North and South dominos, and these dominos all have the same weight (2.1). We obtain the same weight (2.2) of a tiling if instead of assigning the same weight a or b to all the horizontal dominos in a column, we assign a 2 or b 2 to the South dominos and weight 1 to the North dominos in that column.
For symmetry reasons we apply the same operation to East and West dominos. Then instead of (2.1) we assign the following weight to a domino D, where we recall that α = a 2 and β = b 2 ,  Since North dominos have weight 1 we can transfer the weights (3.1) on dominos to weights on the edges of the Aztec diamond graph. The result is shown in Figure 14. The weights alternate per row.
Horizontal edges from (m, n) ∈ Z 2 to (m + 1/2, n) and diagonal edges from (m, n) to (m + 1/2, n + 1) have weight α if n is even and weight β if n is odd. All other edges have weight 1.

Transition matrices
We use the layered structure of the Aztec diamond graph to introduce transition matrices between levels. Here a level is just the horizontal coordinate. There are integer levels m and half-integer levels m + 1/2 with m ∈ Z.
The transition from level m to m + 1/2 is a Bernoulli step. Because of the weights, the transition matrix is, with x, y ∈ Z, Then T m,m+1/2 is 2-periodic, namely T m,m+1/2 (x+2, y +2) = T m,m+1/2 (x, y) for all x, y ∈ Z. As a matrix it is a block Laurent matrix (i.e., a block Toeplitz matrix that is infinite in both directions) with 2 × 2 blocks. The diagonal block is α α 0 β , the block on the first upper diagonal is 0 0 β 0 , and all other diagonals are zero. The associated symbol [14] is To go from level m + 1/2 to level m + 1, we make a number of vertical down steps (possibly zero) and then a horizontal step. All weights are 1 and so the transition matrix is This is a Laurent matrix, but we want to view it as a a block Laurent matrix with 2 × 2 blocks. The diagonal block is 1 0 1 1 , all blocks below the main diagonal are 1 1 1 1 and all blocks above the main diagonal are zero. The symbol is with |z| > 1. Then (the product is matrix multiplication) is the transition matrix from level m to level m + 1, and T is two periodic. The symbol for T is easily seen to be the product of (3.3) and (3.5) which agrees with (2.9). More generally, for any integers m < m we have a transition matrix T m −m to go from level m to level m with symbol A m −m . In particular T 2N is the transition matrix from level 0 to 2N with symbol A 2N . Now we want to invoke the Lindström-Gessel-Viennot lemma [37,52], see also [44,Theorem 3.1] for a proof, which gives an expression for the weighted number of non-intersecting paths on the graph with prescribed starting and ending positions. For us, the starting positions are (0, j), j = 0, . . . , 2N − 1 and the ending positions (2N, j), j = 0, . . . , 2N − 1. Since T 2N (j, k) is the sum of all weighted paths from (0, j) to (2N, k) and so by the Lindström-Gessel-Viennot lemma the partition function is a determinant Because of the layered structure in the graph, we can also look at the positions of the particles at intermediate levels m. We restrict to integer values m, but we could also include the half-integer values.
Given an admissible 2N -tuple of non-intersecting paths, we then find a point set configuration ( where x m 0 < x m 1 < · · · < x m N −1 are the vertical coordinates of the particles at level m. The probability measure on admissible tuples of non-intersecting paths yields a probability measure Then another application of the Lindström-Gessel-Viennot lemma yields that the joint probability for the particle configuration ( The point process (3.9) is determinantal. The correlation kernel is given by the Eynard-Mehta theorem [34], see also [15,16]. In gives the correlation kernel for the particles at level m. Note that G is a finite size submatrix of the two-sided infinite matrix T 2N and G −1 is the inverse of this matrix. To handle the correlation kernel in the large N limit, we need to find a suitable way to invert the matrix G. A fruitful idea is to biorthogonalize. This can be done with matrix valued orthogonal polynomials, and we will discuss this in greater generality in the next section.

The model
We analyze the following situation. We take an integer p ≥ 1, and we consider transition matrices T m : for every m and x, y ∈ Z.
The model also depends on integers N, L ∈ N and M ∈ Z. There will be L + 1 levels numbered as 0, 1, . . . , L. At each level m there are pN particles at integer positions denoted by The initial and ending positions (at levels 0 and L) are deterministic and are given by consecutive integers Our assumption for this section is the following.
is a multi-level particle system with joint probability where the transition matrices T m are p-periodic for every m. The constant Z N in (4.2) is a normalizing constant and δ j (x) = δ j,x is the Kronecker delta.
The determinantal factors with the delta-functions in (4.2) ensure the boundary conditions (4.1).
Assumption 4.1 is satisfied for the two periodic Aztec diamond by (3.9), provided we take p = 2, L = 2N , M = 0 and T m = T for each integer m, with T given by (3.6), using (3.2), (3.4). There are three crucial assumptions contained in Assumption 4.1.
• The transition matrices are p-periodic. It means that the T m are block Laurent matrices with p × p blocks.
• The initial and ending positions of the particles are at consecutive integers. This assumption allows to make a connection with matrix valued polynomials. Note that we allow for a shift pM in the positions at level L compared to the initial positions.
• The transition matrices are such that (4.2) is a probability. That is, (4.2) is always non-negative and we can find a normalization constant Z N such that all probabilities (4.2) add up to 1.
We made two other assumptions, namely • the number pN of particles at each level is a multiple of p, and • the shift pM in the positions of the particles at the final level is also a multiple of p, but these are less essential. They are made for convenience and ease of notation and could be relaxed if needed. For future analysis, we also assume Assumption 4.2. The symbols for the block Laurent matrices T m , m ∈ Z, are analytic in a common annular domain R 1 < |z| < R 2 in the complex plane.
For m < m we use for the transition matrix from level m to level m . The matrix multiplication is well defined because of Assumption 4.2. Every T m,m is a block Laurent matrix with the same period p. The Eynard-Mehta theorem [34] applies to (4.2). We present the Eynard-Mehta theorem as stated in [15]. We assume φ j , ψ j for j = 0, 1, . . . , N − 1 are given functions, and for functions φ, ψ : Z → R, and T : where φ j , ψ j for j = 0, 1, . . . , N − 1 are arbitrary given functions, is determinantal with correlation kernel.
where the Gram matrix G is defined by We apply Theorem 4.3 to (4.2) and we find the following correlation kernel (4.7).
Corollary 4.4. The multi-level particle system (4.2) is determinantal with correlation kernel Proof. This follows from Theorem 4.
The Gram matrix G in (4.6) is a finite section of the block Laurent matrix T 0,L . It has size pN × pN and we also view it as a block Toeplitz matrix of size N × N with blocks of size p × p. It is part of the conclusion of the Eynard Mehta Theorem 4.3 that G is invertible, and so it is in particular a consequence of Assumptions 4.1 and 4.2 that the matrix G is invertible.

Symbols and matrix biorthogonality
Associated with the block Laurent matrices T m and T m,m we have the symbols A m (z) and A m,m (z). According to Assumption 4.2 all symbols are analytic in an annular domain R 1 < |z| < R 2 . We have the identity for R 1 < |z| < R 2 . The series that define the symbol do not commute (in general) and thus the order of the factors in the product is important. We let γ be a circle of radius R ∈ (R 1 , R 2 ) with counterclockwise orientation. By Cauchy's theorem, we can recover the Laurent matrix entries from the symbols, and we have (4.10) In particular and if we restrict to 0 ≤ x, y ≤ N − 1 then the blocks (4.11) are the p × p blocks in the Gram matrix G, see (4.8).
We consider the matrix-valued weight on the contour γ. Clearly W 0,L also depends on M and N but we do not include this in the notation. It introduces a bilinear pairing between p × p matrix valued functions where Q t denotes the matrix transpose (no complex conjugation). The integral is taken entrywise, and so P, Q is again a p × p matrix. A matrix valued function is a polynomial of degree ≤ d if all its entries are polynomials of degree ≤ d.
For invertible matrices P and Q of size pN × pN we define and  Then P j and Q j for j = 0, 1, . . . , N − 1 are matrix valued polynomials of degrees ≤ N − 1.
Proposition 4.5. Let P and Q be invertible pN × pN matrices, and let P j , Q j be the matrix valued polynomials as in (4.14) and (4.15). Let G be the Gram matrix from (4.8). The following are equivalent.
Proof. Consider Then, by (4.12) and the definition (4.14)-(4.15) of the matrix valued polynomials, This is a a block Toeplitz matrix with p × p blocks. For 0 ≤ x, y ≤ N − 1, the xy-th block is 1 2πi γ A 0,L (z)z x−y−M dz z which by (4.8) and (4.11) is equal to the xy-th block of G. Thus which means that G −1 = Q t P if and only if X = I pN . This proves the proposition, since X = I pN is equivalent to the biorthogonality (4.16).
The property (4.16) is a matrix valued biorthogonality between the two sequences (P j ) N −1 j=0 and (Q j ) N −1 j=0 . The matrix valued biorthogonal polynomials are clearly not unique but depend on the particular factorization of G −1 .

Reproducing kernel
Let G −1 = Q t P be any factorization of G −1 and let P j , Q j be the matrix polynomials as in (4.14) and (4.15). We consider which is a bivariate polynomial of degree ≤ N − 1 in both w and z. (a) For every matrix valued polynomial P of degree ≤ N − 1 we have Proof. Parts (a) and (b) are immediate from the biorthogonality (4.16) and the fact that any matrix valued polynomials P of degree ≤ N − 1 can be written as P (z) = Let R N (w, z) be as in part (c), and suppose that for every matrix valued polynomial Q of degree ≤ N − 1. For a fixed w we note that z → R N (w, z) is a matrix valued polynomial of degree ≤ N − 1 and it can be written as a linear combination of P 0 , . . . , P N −1 with matrix coefficients. The matrix coefficients depend on w, and we get for some A j (w), j = 0, . . . , N − 1, Then taking (4.20) with Q = Q k , and using the biorthogonality (4.16), we get for every k = 0, . . . , N − 1. Thus R N (w, z) = R N (w, z) by (4.21) and (4.17). This proves part (c) in case the reproducing property of (b) is satisfied. The other case follows similarly.
Because of (4.18) and (4.19) we call R N (w, z) the reproducing kernel for the pairing (4.13).
From part (c) of Lemma 4.6 we find in particular that the sum in (4.17) does not depend on the particular choice of factorization G −1 = Q t P. It does depend on G as can be also seen from the calculation using (4.14), (4.15), and (4.17) which clearly only depends on G.

Main theorem
Now we are ready for the main theorem in this section.
x, y ∈ Z, 0 < m, m < L, (4.23) where R N (w, z) is the reproducing kernel (4.17) built out of matrix valued biorthogonal polynomials associated with the weight W 0,L (z) = Proof. We already know that the particle system is determinantal with kernel given by (4.7).
The first term in (4.7) gives rise to the first term on the right-hand side of (4.23) in view of (4.10) (note that x and y are interchanged in T m ,m (y, x) in (4.7)).
Let K 0 (m, x; m , y) be the second term in the right-hand side of (4.7). Instead of summation indices i and j in the double sum we use pν + k and pν + k , with 0 ≤ ν, ν ≤ N − 1, 0 ≤ k, k ≤ p − 1. Then from (4.7), T m ,L (py+i, pM +pν +k ) G −1 pν +k ,pν+k T 0,m (pν+k, px+j). (4.24) Using (4.10) we can write this in block form where the first factor in the right-hand side of (4.25) is a block row vector of length N with p × p blocks, and the last factor is a similar block column vector. We combine the integrals to obtain a double integral and then we use (4.22) to find which completes the proof.

Matrix valued orthogonal polynomials
The goal of this subsection and the next is to express the reproducing kernel (4.17) (or (4.22)) in terms of matrix valued orthogonal polynomials (MVOP) and use a Christoffel-Darboux formula for the sum (4.17). Such a formula is known for MVOP in various forms [25], [38], [4], [5]. We are however in a non-standard situation, with a non-Hermitian orthogonality, and the MVOP need not exist for every degree. Fortunately, the degree N MVOP will exist in the present situation, as we will explain in this subsection.
If we can find a factorization G −1 = Q t P leading to matrix valued polynomials (4.14)-(4.15) with P j (z) = Q j (z) and deg P j = j for every j = 0, 1, . . . , N − 1, then we would have a finite sequence of MVOP that are in fact orthonormal and then From (4.27) it would follow that R N (z, w) = R N (w, z) t and this is an identity that is not necessarily satisfied. Thus we cannot expect that the orthonormal MVOP exist. Instead we focus on monic MVOP. If the monic MVOP P j exists for every degree j, then we have with some matrix H j . If H j is invertible for every degree then the two sequences of matrix valued polynomials (H −1 j P j (z)) N −1 j=0 and (P j (z)) N −1 j=0 are biorthogonal, and thus The orthogonality (4.26) is non-Hermitian orthogonality, and it is not associated with a positive definite scalar product. Also existence and uniqueness of the monic MVOP is not guaranteed in general. However, the MVOP of degree N does exists, and this is a consequence of the fact that G is invertible.
Lemma 4.8. There is a unique monic matrix valued polynomial Proof. The conditions (4.30) give us p 2 N linear equations for the p 2 N unknown coefficients of a monic matrix valued polynomial of degree N . The linear system has matrix G, provided we number the coefficients and the conditions appropriately, and since G is invertible, the existence and uniqueness of P N follows.
More explicitly, write P N (z) = I p z N + N −1 j=0 C j z j with p × p matrices C j that are to be determined. The orthogonality conditions where G j,k denotes the j, kth block of the block Toeplitz matrix G, see also (4.11). Varying k = 0, 1, . . . , N − 1, we see that The matrix G is invertible, and thus the matrix coefficients C 0 , . . . , C N −1 are uniquely determined, and the monic MVOP of degree N exists uniquely.

Riemann-Hilbert problem and Christoffel-Darboux formula
The MVOP of degree N is characterized by a Riemann-Hilbert problem of size 2p × 2p. The RH problem asks for a 2p × 2p matrix valued function Y : C \ γ → C 2p×2p satisfying • Y is analytic, In the scalar valued case, i.e. p = 1, the RH problem is due to Fokas, Its, and Kitaev [35]. The matrix valued extension can be found in [19,29,38]. It is similar to the RH problem for multiple orthogonal polynomials [63].
The RH problem has a unique solution, since by Lemma 4.8 the monic MVOP of degree P N exists and is unique. The solution is where Q N −1 is a matrix valued polynomial of degree ≤ N − 1 such that (4.32) One can show that Q N −1 also uniquely exists, since the conditions (4.32) give a system of p 2 N linear equations for the p 2 N coefficients of Q N −1 , and the matrix of this system can be identified with G. Since G is invertible, there is a unique solution. If the leading coefficient of Q N −1 would be invertible (which is typically the case, but it is not guaranteed in general) then the monic MVOP P N −1 of degree N −1 would exist as well and In the following result we express the reproducing kernel R N (w, z) in terms of the solution of the RH problem. It can be viewed as a Christoffel-Darboux formula and in this form it is due to Delvaux [29]. It is similar to the Christoffel-Darboux formulas for multiple orthogonal polynomials [12,24] which also use the RH problem. Proposition 4.9. We have Proof. This is due to Delvaux [29, Proposition 1.10], see also [38]. Since the context and notation of [29] is somewhat different from the present setting, we give an outline of the proof. The right-hand side of (4.33) is a bivariate polynomial in z and w of degrees ≤ N − 1 in both variables, see Lemma 2.3 in [29] for details. We show that it satisfies the reproducing kernel property from Lemma 4.6 (b), and we follow the proof of [29, Proposition 2.4].
Let Q be a matrix valued polynomial of degree ≤ N − 1. We replace R N (w, z) by the right-hand side of (4.33) in the integral on the left of (4.19) and we obtain (4.31) we find is a matrix valued polynomial of degree ≤ N − 2. The first term in (4.34) then vanishes because of the orthogonality conditions (4.30) and (4.32) satisfied by P N and Q N −1 . The second term contains the second column of Y , see (4.31), and therefore (4.34) is equal to and this is indeed Q t (w), which proves that the right-hand side of (4.33) has the reproducing property of Lemma 4.6 (b) that characterizes R N (w, z) by Lemma 4.6 (c). The proposition follows.
We insert (4.33) into formula (4.23) and find a convenient formula for the correlation kernel in terms of the solution of the RH problem.
A possible asymptotic analysis of the kernel would consist of two parts. First we do an analysis of the RH problem that would give us the asymptotic behavior of the kernel (4.33). Then this is followed by an asymptotic analysis of the double contour integral in (4.23) by means of classical methods of steepest descent. We are able to this for the two periodic Aztec diamond.

Example 1: Aztec diamond
Let us put a weight 1 on a horizontal domino and a weight q on a vertical domino in a tiling of the Aztec diamond of size N . This corresponds to putting weights 1 on the horizontal edges in the Aztec diamond graph, see Figure 14, and weight q on the diagonal and vertical edges.
Assumption (4.1) holds with L = N , M = 0, and transition matrices that are independent of m, Then T is a Laurent matrix and the symbol is Hence, , and since L = N and M = 0, it follows that The (scalar) weight is rational with a pole at z = q and a zero at z = −1/q, both of order N . The orthogonal polynomials are properly rescaled Jacobi polynomials, but with parameters −N, N . The Jacobi polynomial of degree k with parameters α and β may be given by the Rodrigues formula If α and β are integers, and γ is a circle not going through ±1, then for j = 0, 1, . . ., we find by using (4.36) and after integrating by parts k times The weight (4.35) has a pole at z = q and a zero at z = −1/q. By an affine change of variables we map these to ±1, and then we use (4.37) with α = −N and β = N . It follows that Thus the orthogonal polynomials for the weight (4.35) are the rescaled Jacobi polynomials The uniform Aztec diamond has been analyzed in detail with Krawtchouk polynomials in [33,41,43], which are discrete orthogonal polynomials. We find it intriguing that an alternative approach with Jacobi polynomials seems also possible.

Example 2: Hexagon tiling
that is independent of m. It is a Laurent matrix with symbol A(z) = z + 1. Then A 0,L (z) = (z + 1) L and the weight is The scalar weight is again rational with one zero and one pole. The orthogonal polynomials are again rescaled Jacobi polynomials, but now with parameters −(M + N ) and L, namely for 0 ≤ k < N , An asymptotic analysis of the lozenge tilings of the hexagon based on the Jacobi polynomials (4.40) seems possible, but has not been pursued yet. See [51,53,54] for an asymptotic analysis of Jacobi polynomials with varying non-standard parameters.

Example 3: Aztec diamond with periodic weights
The third example is the main interest of this paper: the two periodic Aztec diamond of size 2N . We saw in Section 3 that the model gives rise to the multi-level particle system (3.9). This satisfies the Assumption 4.1 if we take p = 2 and L = 2N and M = 0. The transition matrices are independent of m, see (3.6) and the matrix symbol is given by (3.7). The weight matrix is W 0,L (z) = A 0,L (z) as in (2.15) and (2.27). Observe that W has no pole at the origin. Theorem 4.7 applies and it gives the form of the correlation kernel, in 2 × 2 matrix form, that will be stated in (5.4) below. It is equivalent to the form already announced in Section 2 in (2.29).
The correlation kernel contains the reproducing kernel R N (w, z) with respect to the varying weight W N , and that by Proposition 4.9 is expressed in terms of the RH problem for the MVOP of degree N . By Lemma 4.8 we conclude that the degree N monic MVOP with respect to the weight W N exists. What about degrees < N ? While it does not matter for the rest that follows, we can show that the MVOP of lower degrees indeed exist. Proof. Note that W N (z) = A 2N (z) z N , which we can also write as For k < N , we consider 2k non-intersecting paths on the same graph with L = 2N levels, starting at consecutive positions 0, 1, . . . , 2k − 1, and ending at shifted positions 2M, 2M + 1, . . . , 2M + 2k − 1. Provided that there are such non-intersecting paths, we have a determinantal point process as before, and we conclude by an application of Lemma 4.8 that the monic MVOP of degree k uniquely exists. It is readily seen that such paths indeed exist for the Aztec diamond graph. For example, by letting the 2k paths make a diagonal up-step 2M = 2N − 2k times followed by 2k horizontal steps, and there are no down steps. Then these paths are indeed non-intersecting.
The construction of non-intersecting paths does not work if k ≥ N + 1 and the MVOP does not exist for those degrees.
In the next section we continue with the analysis of the RH problem and we show that the correlation kernel (2.29) can be rewritten as (2.8).

Analysis of the RH problem
We consider an Aztec diamond of size 2N with two periodic weighting.

Correlation kernel
In the two periodic Aztec diamond we find the matrix symbol and the matrix valued weight is W N with W given by (4.41). Note that W N is a rational function with a pole at z = 1 only. The contour γ in the RH problem from Section 4.6 goes around 0 and lies in the domain |z| > 1. By analyticity, since W only has a pole at z = 1, we are free to deform the contour to a circle around 1. We use γ 1 to denote the circle of radius r < 1 around 1. We obtain the following RH problem for Y : • Y is analytic, • Y has asymptotic behavior Because of Theorem 4.7 and (4.33) we find the following correlation kernel for arbitrary integer levels m, m with 0 < m, m < 2N = L and M = 0, .
The contour γ 0,1 in (5.4) is a circle of radius > 1 + r around the origin, as before. The radius is large enough such that γ 1 lies inside γ 0,1 . The analysis of the correlation kernel (5.4) consists of two parts. First we apply a RH analysis to the RH problem for Y and then we use this for an asymptotic analysis of the double integral.
The RH analysis is remarkably simple. It is not an asymptotic analysis, since the outcome is an exact new formula for the correlation kernel.
Theorem 5.1. Assume 2y ≥ m and N is even. Then the correlation kernel .
where F is given by (2.10).
Passing from the non-intersecting path model back to the domino tilings of the Aztec diamond, we should make the change of variables m → m, 2x → m + n, m → m and 2y → m + n , that come from the shear transformation described in Section 3.3. Inserting these values in (5.5) we obtain the correlation kernel (2.8) and so Theorem 2.2 follows immediately from Theorem 5.1.
The rest of Section 5 is devoted to the proof of Theorem 5.1. We follow the general scheme of the analysis of RH problems, known as the Deift-Zhou steepest descent analysis [28], which was first applied to orthogonal polynomials in [26,27]. Extensions to larger size RH problems are for example in [13,30], see also the survey [50] and the references therein. However, the RH analysis in this section is not an asymptotic analysis, as it produces the exact formula (5.5).
We use z for a generic coordinate on R, and if we want to emphasize that z is on the jth sheet, we write z (j) , for j = 1, 2. We write λ for the function on R, is on the jth sheet, (5.8) see (2.17), and similarly for ρ. These are meromorphic functions on R, namely ρ = (α + β)z + y, λ = ρ 2 z(z − 1) 2 , see (2.11), (2.12) and (2.17).
Lemma 5.2. (a) ρ has a simple zero at z = 0, a double zero at z = 1 (2) (the point z = 1 on the second sheet), a triple pole at z = ∞, and no other zeros or poles, (b) λ has a double zero at z = 1 (2) , a double pole at z = 1 (1) , and no other zeros or poles, has a zero at z = 0, and a double zero at z = −1 (1) (if α > β).

First transformation Y → X of the RH problem
We use the matrix of eigenvalues (5.6) in the first transformation of the RH problem. We define which satisfies the following RH problem.
and from the definition (5.20) and the RH problem for X we obtain • U is analytic, • U has the jumps • U has asymptotic behavior To obtain the jump (5.22) on (−∞, −α 2 ] ∪ [−β 2 , 0] we also have to use the fact that λ 1,± = λ 2,∓ on these cuts. The asymptotic condition (5.23) requires some explanation. The first equality in (5.23) is clear from the definition (5.20) of U for |z − 1| > r, and the asymptotic behavior (5.16) of X. By (2.16) we have E(z)Λ ±N/2 (z) = W ±N/2 (z)E(z) so that by (4.41), (5.14), and (5.21) we get that and this leads to the second equality in (5.23).

Third transformation U → T
In the third transformation we turn the entries (z−1) ±2N in the jump matrix on γ 1 into an off-diagonal entry. It corresponds to the opening of lenses in a steepest descent analysis. We also remove the 24-entry in the jump matrix on γ 1 . We define Straightforward calculations, where we just use (5.24) and the RH problem for U , show that T satisfies • T has the jumps • T has asymptotic behavior
Since E is not invertible at z = 0, z = −α 2 , z = −β 2 , see also (5.18), we could have introduced singularities at these points. Therefore we look at the combined transformations Y → X → U → T → S in order to express S directly in terms of Y . For z outside of γ 1 we have by (5.15), (5.20), (5.24) and (5.27), Since EΛE −1 = W by (2.16), we simply have (recall N is even) This shows indeed that (5.27) does not introduce any singularities, since det W (z) = 1 for every z, and W (z) and W −1 (z) have poles at z = 1 only. Thus S has analytic continuation across (−∞, −α 2 ] and [−β 2 , 0) and satisfies the following RH problem that we obtain immediately from (5.27) and the RH problem for T .
• S has asymptotic behavior The RH problem is now normalized at infinity. Note also that the transformation (5.27) restores the property det S = 1, since det T = det X = (det E) 2 , see (5.19). The jump matrix in (5.29) is lower triangular, and the RH problem for S is normalized at infinity by (5.30). This means that we can solve the RH problem explicitly by a contour integral. We find

Proof of Theorem 5.1
We can now give the proof of Theorem 5.1.
Proof. We analyze the effect of the transformations on the correlation kernel (5.4). From (5.28) and (2.15) we have for z, w outside of γ 1 , which is part of the expression that appears in the double integral in (5.4). Because of (5.31) and Using (5.32) and (5.33) we see that the double integral in (5.4) is equal to We change the order of integration in (5.34) and evaluate the w-integral first. By a residue calculation Indeed the singularities at w = 0 and w = 1 in the integrand in the left-hand side of (5.35) are removable (we use (5.1), N is even, and 2y ≥ m ). The only singularity is at w = s and (5.35) indeed follows by Cauchy's formula since s ∈ γ 1 lies inside γ 0,1 . Using (5.35) in (5.34) and changing the integration variable s to w, we obtain the double integral in (5.5). The single integral in (5.5) is of course immediate from (5.4). This completes the proof of Theorem 5.1.

A consistency check
The RH analysis gives us explicit formulas, and in particular also for the 2 × 2 left upper block which by (4.31) should be the monic MVOP of degree N . Following the transformations Y → X → U → T → S and the expression (5.31) for S, we see that which is indeed a monic matrix valued polynomial of degree N since N is even. Note that W (z) has a double pole at z = 1, hence W −N/2 (z) has a pole of order N at z = 1, and the pole is compensated by the N th order zero of (z − 1) N . We then check that also has a removable singularity at z = 1. Hence it is also a matrix polynomial and the matrix orthogonality follows in a trivial way from Cauchy's theorem 1 2πi γ 1 P N (z)W N (z)Q(z)dz = 0 2 for every matrix valued polynomial Q (and not just for polynomials of degree ≤ N − 1). The degree N − 1 polynomial Q N −1 is in the left lower 2 × 2 block of Y , see (4.31). This is a polynomial of degree N − 1, but not necessarily a monic one. From the transformation in the RH analysis and (5.31), we find The analytic continuation to |z − 1| < r, is given by the Sokhotskii-Plemelj formula, Note that from (2.14) and (2.16) we have F (z)W −N/2 (z) = λ −N/2 1 (z)F (z) and this has a zero at z = 1 of order N because of Lemma 5.2 (b). The zero cancels the N th order pole in (z − 1) −N , and we see that the extra term is analytic for |z − 1| < r, which confirms that the expression defining Q N −1 is a polynomial of degree ≤ N − 1.
Let's verify the orthogonality (4.32) where we take a contour γ that lies in the exterior of γ 1 . For an integer k ≥ 0, we have by (5.37) We change the order of integration and use Cauchy's formula (the only pole is at z = s) to obtain Note that (I 2 − F (s))W N/2 (s) = (I 2 − F (s))λ N/2 2 (s) (again by (2.14) and (2.16)) and λ N/2 2 (s) has a zero of order N at s = 1 by Lemma 5.2 (b). Thus and combining this with (5.38) leads to Recall that W N/2 (s) is a rational matrix valued function whose only pole is at s = 1 and it is bounded at infinity. Then in (5.39) we move the contour γ 1 to infinity. There is no contribution from infinity if k ≤ N − 2, while for k = N − 1 there is a residue contribution at infinity and the expression (5.39) becomes −I 2 for k = N − 1. We conclude that (4.32) indeed holds.

Asymptotic analysis
In the final section of the paper we are analyzing the formula (2.8) in a scaling limit where N → ∞ and the coordinates (m, n) and (m , n ) scale linearly with N . We are going to distinguish the three phases of the model, and prove Theorems 2.10, 2.13, and 2.14.

Preliminaries
We first rewrite the formula (2.8) in a form that already contains the gas phase kernel (2.36) and double integrals with the phase functions Φ 1 and Φ 2 from (2.30), see Corollary 6.3. We may and do assume that the contour γ 0,1 is a contour in C\((−∞, −α 2 ]∪ [−β 2 , 0]) going around the interval [−β 2 , 0] once in positive direction. Note that n could go up to 2N − 1, and then O(w −N +n /2−1/2 ) = O(w −1 ). However, then we are close to the boundary of the Aztec diamond, and we do not consider this in what follows, since we focus on the gas phase. So we assume n ≤ 2N − 2, and then the integrand in double integral in (2.8) is O(w −3/2 ) as w → ∞.
Then for a fixed z ∈ γ 0,1 , we deform the contour γ 1 to a contour Γ going around the negative real axis, starting at −∞ in the upper half-plane and ending at −∞ in the lower half-plane, as in Figure 16. Since the integrand is O(w −3/2 ), by Lemma 6.1, there is no contribution from infinity, but there is a residue contribution from the pole at w = z. These residues combine to give the z-integral (we use that F (z) and A(z) commute) Together with the single integral in (2.8) this gives the limit (2.36) that we expect to get in the gas phase. We proved the following.
Thus to establish Theorem 2.10 we have to prove that in the gas phase the double integral in (6.1) tends to 0 as N → ∞ at an exponential rate.

2)
with contours as in Figure 16.

Saddle points
The large N behavior of the z-integrals in (6.2) is dominated by the factors e N Φ 1 (z) and e N Φ 2 (z)/2 that are exponential in N . Similarly the w part of the integrand is dominated by e −N Φ 1 (w)/2 . We study the saddle points, which in Definition 2.8 were already introduced as the zeros of the meromorphic differential Φ (z)dz from (2.31) defined on the Riemann surface R associated wth (2.11). Of course, Φ depends on ξ 1 , ξ 2 , and thus the saddle points depend on these parameters. Throughout we restrict to −1 < ξ 1 , ξ 2 < 1. The differential has simple poles at 1 (1) , 1 (2) , 0 and ∞ with residues given in the following table.
residue of residue residue residue pole of dz The residues of λ λ dz at z = 1 (1) and z = 1 (2) come from the double pole and double zero that λ has at these points, see Lemma 5.2 (b). The residues add up to zero, as it should be.
We assume α > 1 so that the genus of R is one. Then there are also four zeros of Φ dz counting multiplicities.
Recall that the real part of the Riemann surface consists of the cycles C 1 and C 2 as in (2.32). Proposition 6.4. For every ξ 1 , ξ 2 ∈ (−1, 1) there are at least two distinct saddle points on the cycle C 1 .
Proof. If C is a path from P to Q on the Riemann surface avoiding the poles, then by (2.30), for a choice of continuous branches of the logarithms along the path. Since ξ 1 and ξ 2 are real, it follows that the real part is well-defined, it depends on P and Q, but is otherwise independent of the path. Thus for a closed path C.
Observe that that there are no poles on the cycle C 1 , and Φ is real there. If there were no two distinct zeros on C 1 , then there would be no sign change, and the integral would be non-zero and real, which would contradict the condition (6.5).
The saddle points are explicit in case ξ 1 = 0, since then by (2.31) The equation 2 z−1 = 1+ξ 2 z has the unique solution This gives us two saddle points, namely the two points on R with (6.7) as z-coordinate. The other two saddles come from the branch points −α 2 , −β 2 , which are zeros of the differential dz. The branch point z = 0 is also a zero of dz, but this zero gets cancelled by the (double) pole of 1+ξ 2 z in (6.6). For special values of ξ 2 the saddles at z = z c (ξ 2 ) coincide with the saddle at −α 2 or −β 2 . This happens for the values ±ξ * 2 with Then depending on the value of ξ 2 , we are in the liquid or gas phase, or on the liquid-gas transition, as defined in Definition 2.9.
(a) (0, ξ 2 ) cannot be in the solid phase.
Proof. (a) It is clear from (6.7) that z c (ξ 2 ) < 0 and so there are no saddles on the positive real axis.
, and if ξ * 2 < ξ 2 < 1 then z c (ξ 2 ) ∈ (−∞, −α 2 ). Even though z c (ξ 2 ) is real, the two saddles with z coordinate equal to z c (ξ 2 ) are not on the real part of the Riemann surface, and thus in both cases we are in the liquid phase.
(c) If −ξ * 2 < ξ 2 < ξ * 2 then z c (ξ 2 ) ∈ (−α 2 , −β 2 ). Then the saddles with z coordinate equal to z c (ξ 2 ) are on the cycle C 1 . The branch points −α 2 and −β 2 are the other two saddles and they are also on the cycle. Thus all four saddles are on the cycle C 1 and they are distinct, and we are in the gas phase.

Algebraic equation
The condition of coalescing saddle points leads to an algebraic equation for ξ 1 and ξ 2 . We are able to calculate it with the help of Maple.
First of all, the saddle point equation Φ (z)dz = 0, see (2.31), leads us to consider 2 which after clearing denominators, and using (2.17) and (2.12), gives a polynomial equation in z and y = z(z + α 2 )(z + β 2 ). We eliminate the square root to obtain a polynomial equation in z of degree 4, which is By definition, the saddles are the four zeros of the polynomial (6.9). The discriminant with respect to z of (6.9) is a polynomial in ξ 1 and ξ 2 that has trivial factors ξ 2 1 and ξ 2 2 . The remaining factor is a degree 8 polynomial, which is symmetric in the two variables. Setting this to zero, we obtain the following equation for coalescing saddles: Up to a multiplicative constant, the equation (6.10) coincides with the one given by Chhita and Johansson [20,Appendix A]. See also [60,Section 8] for an equation that corresponds to (6.10) with α = 2 up to a change of variables. For α = 1, (6.10) reduces (up to a numerical factor) to and the real section is the unit circle.
Remark 6.6. The discriminant of (6.9) also vanishes for ξ 1 = 0 or ξ 2 = 0, and there is indeed a double root of (6.9) for these values of the parameters. However, they do not correspond to coalescing saddle points. For ξ 1 = 0, the double root is at z = z c (ξ 2 ) from (6.7). which corresponds to two different saddles on the Riemann surface R, unless ξ 2 = ± α−β α+β , see Lemma 6.5. For ξ 2 = 0, the double root turns out to be at z = −1, but the saddles are also on different sheets of R.
Proof. If ξ 1 = 0 then all statements of the proposition follow from Lemma 6.5. The proof in the general cases follows by a continuity argument, since the saddles depend continuously on the parameters ξ 1 , ξ 2 , and a saddle can only leave the real part of the Riemann surface if it coalesces with another saddle and then the pair can move away from the real part. This transition can thus only occur for ξ 1 , ξ 2 satisfying the algebraic equation (6.10). Note that this argument also applies to the point at infinity, since by definition (2.32) the point at infinity is included in the real part.
Combining with Lemma 6.5 we find that any point (ξ 1 , ξ 2 ) inside the inner component belongs to the gas phase, and any point in between the inner and outer component belongs to the liquid phase.
To treat the solid phase, we check that any (ξ 1 , ξ 2 ) close enough to a corner point is in the solid phase. We can see this from the equation (6.10). If ξ 1 = ±1 and ξ 2 = −1 then (6.9) has solutions z = 0, z = 1 (and two other solutions that are on the cycle C 1 ), and if ξ 1 = ±1 and ξ 2 = 1 then (6.9) has solutions z = ∞ and z = 1. Thus for each of the four corner points there are two distinct saddles on C 2 . Then by continuity this continues to be the case if (ξ 1 , ξ 2 ) is close enough to one of the corner points, and it continues to be so until the two saddles on C 2 coalesce, and this happens on the outer component.
This completes the proof of Proposition 6.7.

Gas phase: steepest descent paths
In the gas phase all four saddles are located on the cycle C 1 , and they are all simple. To prepare for the proof of Theorem 2.10, we need more precise information on the location of the saddles.
For ξ 1 = 0, we notice that both λ 1 and λ 2 are real and negative on the interval [−α 2 , −β 2 ], see (5.12), with a square root behavior at endpoints (which follows from (2.17) and the square roots in (2.12)). Thus λ 1,2 become infinite at the endpoints, and a closer inspection of (2.12), (2.17) shows that and 14) are infinite at the endpoints of the interval [−α 2 , −β 2 ] but with opposite signs. By continuity there is an odd number of zeros for each of them. There are exactly four simple saddles on the cycle C 1 as we are in the gas phase, and therefore one of Φ j , j = 1, 2, has three simple zeros and the other one has one simple zero. We already noted that for ξ 1 = 0 Re Φ j (z) = 2 log |z − 1| − (1 + ξ 2 ) log |z| + ξ 1 log |λ j (z)| (6.15) attains a local minimum at an interior saddle for j = 1, 2. Because of analytic dependence on parameters this continues to be the case for ξ 1 = 0, and in fact, since there is no coalescence of saddle points, it remains true for every (ξ 1 , ξ 2 ) ∈ G. Thus saddles z s,1 and z s,2 where Re Φ has a local minimum exist, and z s,j is on the jth sheet. Now, if ξ 1 > 0 then from (6.12) and (6.14) it follows that and so Re Φ 1 increases on an interval [−α 2 , −α 2 +δ] and decreases on [−β 2 − δ, −β 2 ] for some δ > 0. Since there is a local minimum at z s,1 on the first sheet, it should be that Re Φ 1 has two local maxima, say z s,3 < z s,4 , with −α 2 < z s,3 < z s,1 < z s,4 < −β 2 . In case ξ 1 < 0 we find in the same way that the local maxima are on the second sheet, with −α 2 < z s,3 < z s,2 < z s,4 < −β 2 .
The path of steepest descent from the saddle z s,j , j = 1, 2, is the curve γ sd,j through z s,j where the imaginary part of Φ j is constant and the real part decays if we move away from the saddle. Since Re Φ j on C 1 has a local minimum at z s,j , the path of steepest descent meets the real line at a right angle.
Emanating from z s,j , j = 1, 2 are also curves γ l,j and γ r,j where the real part is constant (l stands for left, and r stands for right). The curve γ l,j emanates from z s,j at angles ±3π/4. It consists of a part in the upper half plane and its mirror image with respect to the real line in the lower half plane. Similarly, γ r,j emanates from z s,j at angles ±π/4, and it is also symmetric in the real line. Near z s,j we have that γ l,j is to the left of the steepest descent path γ sd,j and γ r,j is to the right.
Then γ l,j and γ r,j are parts of the boundary of the domains: Lemma 6.9. Suppose (ξ 1 , ξ 2 ) ∈ G, and j = 1, 2. Then the following hold.
(b) The steepest descent path γ sd,j intersects the positive real line at 1, γ l,j intersects the positive real line at a value > 1 while γ r,j intersects the positive real line at a value < 1.
(c) Ω − j is a bounded open set with at most three connected components. One component (which we call the main component) contains γ sd,j \ {z s,j }. There are at most two other connected components, namely a component containing −β 2 (if Re Φ j (−β 2 ) < Re Φ j (z s,j )) and a component containing −α 2 (if Re Φ j (−α 2 ) < Re Φ j (z s,j )). The other components (if they exist) are at positive distance from the main component.

(d) Ω +
j is an open set with an unbounded component that contains a contour Γ 1,j that goes around (−∞, −α 2 ] and with a bounded component that contains a contour Γ 2,j going around the interval [−β 2 , 0]. Proof. The lemma is straightforward to verify if ξ 1 = 0, since in that case Re Φ 1 (z) = Re Φ 2 (z) = 2 log |z − 1| − (1 + ξ 2 ) log |z| (6.18) and z s,1 = z s,2 = z c (ξ 2 ) = − 1+ξ 2 1−ξ 2 which is in (−α 2 , −β 2 ) since we are in the gas phase. Then γ l,j , γ sd,j , and γ r,j are independent of j and we have a situation as in Figure 17. The curves γ l,j and γ r,j enclose the domain Ω − j that is shaded in the figure, and Ω − j has only one component in this case. The non-shaded domain Ω + j has an unbounded component with boundary γ l,j ∪(−∞, −α 2 ] and a bounded component with boundary γ r,j ∪[−β 2 , 0]. The contours Γ 1,j and Γ 2,j can be taken in Ω + j as specified in part (d) of the lemma.
For ξ 1 = 0 the behavior on the positive real axis is the same: for both j = 1 and j = 2, we have that Φ j (x) is real for x > 0 and it decreases from +∞ to −∞ on the interval [0, 1] and increases from −∞ to +∞ on [1, ∞). This is so because otherwise there would be a zero of the derivative, which would be a saddle point, but all saddle points are on the cycle C 1 since we are in the gas phase.
The behavior of Re Φ j on the cuts (−∞, −α 2 ] ∪ [−β 2 , 0] is exactly the Now let's follow the paths γ l,j and γ r,j , where the real part Re Φ j is constant, as they move away from z s,j into the upper half plane. These paths remain bounded, since Re Φ j (z) → +∞ as |z| → ∞, which follows from (6.15), Lemma 5.2 (c), and the fact that ξ 2 < 1. The two paths cannot come together in the upper half plane, since then they would enclose a domain on which Re Φ j is harmonic and constant on the boundary, which violates the maximum/minimum principle of harmonic functions. For the same reason they cannot meet at a point on the positive real axis.
Suppose now one of the paths comes to the cut (−∞, −α 2 ] at a point q. Then this path, together with its mirror image in the real line, encloses a bounded domain D and Re Φ j has the constant value Re Φ j (z s,j ) on its boundary. The value of Re Φ j is smaller on the interval [q, −α 2 ] because of (6.19). Also Re Φ j is harmonic on D \ [q, −α 2 ] and it follows from the maximum principle for harmonic functions that This is a contradiction, since z s,j is a saddle at which Re Φ j attains a local minimum when restricted to the cycle C 1 . We arrive at a similar contradiction if one of the paths γ l,j and γ r,j comes to the cut [−β 2 , 0].
Thus the two paths can only leave the upper half-plane via the positive real axis. Since Φ j is strictly decreasing on [0, 1] and strictly increasing on [1, ∞], we must conclude that γ l,j comes to the unique value in (1, ∞) and γ r,j to the unique value in (0, 1) where Φ j is equal to Re Φ j (z s,j ). Together with their mirror images in the real axis, they both form simple closed curves that go around the cut [−β 2 , 0]. They enclose a domain where Re Φ j is smaller than Re Φ j (z s,j ), i.e., it is contained in Ω − j , see (6.16). The path of steepest descent γ sd,j lies in Ω − j , and Re Φ j decreases along γ sd,j in the upper half plane. It will meet the positive real axis at 1, where Re Φ j is −∞, since if it would meet the positive real axis at some other point, this point would be a saddle, but there is no saddle on C 2 .
We now proved parts (a) and (b) of Lemma 6.9. We also proved that the component of Ω − j that is bounded by γ l,j and γ r,j contains the steepest descent curve γ sd,j \{z s,j }. We also see that Ω − j is bounded since Re Φ j (z) → +∞ as |z| → ∞.
Any other connected component of Ω − j has to intersect with the branch cuts (−∞, −α 2 ] ∪ [−β 2 , 0], since otherwise we have again a contradiction with the minimum principle for harmonic functions. If a component of Ω − j intersects (−∞, −α 2 ] then it will do so along an interval (q, −α 2 ] for some q < −α 2 , because of (6.19). Hence there can be at most one such component, and this exists if and only if if Re Φ j (−α 2 ) < Re Φ j (z s,j ). Similarly, there is at most one component that intersects [−β 2 , 0], and thus in total there are at most three components.
Finally, since z s,j is a simple saddle, and Re Φ j attains a minimum at z s,j when we restrict to the cycle C 1 , there is δ > 0 such that Thus the boundary of another component (if it exists) intersects the interval (−α 2 , −β 2 ) at a point different from z s,j , and thus the component stays at positive distance from the main component. This proves part (c) of the lemma. The component of Ω − j that contains −α 2 (if it exists) is bounded and belongs to the exterior of the closed contour γ l,j . However, by part (c), it remains at a positive distance from γ l,j , and therefor we can find a contour Γ 1,j as in the statement of part (d) that goes around the interval (−∞, −α 2 ] and avoids the closure of the domain Ω − j while staying to the left of γ l,j . Similarly, we can find Γ 2,j as in part (d). See Figures 18 and 19 for plots of γ sd,j and Γ 1,j and Γ 2,j in the situations where Ω − j has no component containing −α 2 or −β 2 ( Figure 18) and where Ω − j has a component containing −α 2 that Γ 1,j should avoid ( Figure 19).
This completes the proof of Lemma 6.9.
6.5 Gas phase: proof of Theorem 2.10 We are now ready to prove Theorem 2.10.
We are going to apply a similar argument to (6.21) and deform γ 0,1 to the steepest descent contour γ sd,2 passing through z s,2 . We again replace Γ by the union of two contours, one going around (−∞, −α 2 ] and the other around [−β 2 , 0], but now we are going to move these contours to the cuts, where we note that λ 1,± = λ 2,∓ and Φ 1,± = Φ 2,∓ plus a purely imaginary constant (depending on the precise branches of the logarithms that we choose in (2.30)). We also note that there is no pole at w = 0, and that I 2 − F (w) is the analytic continuation of F (w) to the second sheet. Then we deform the contours further to Γ 1,2 and Γ 2,2 as in Lemma 6.9 (d), and we see that (6.21) is equal to The minus sign comes since we use the orientation on Γ 1,2 and Γ 2,2 as in Figure 19. Analogous to (6.22) and (6.23) we now have and the integral (6.24) tends to 0 at an exponential rate as N → ∞. This completes the proof of Theorem 2.10.

Cusp points: proof of Theorem 2.13
For the proof of Theorem 2.13 we are going to use (6.2) once more. In the scaling (2.40) of the parameters we have that the saddle points coalesce to a triple saddle point at −α 2 . We are going to deform γ 0,1 so that it comes close to −α 2 . The main contribution to the integrals in (6.2) then comes from the triple saddle at −α 2 , as the integrands are exponentially small if w and/or z are outside of a small neighborhood of −α 2 . This follows as in the proof of the gas phase.
Hence, our task is to investigate how the entries in the integrals in (6.2) behave for z and w close to −α 2 . We do this in the following two lemmas and their corollaries. Besides the constants c 1 and c 2 from (2.41) we also use c 0 = α + β √ 2 .
We also need the behavior of F (z) near the saddle z = −α 2 . Lemma 6.12. We have the following.
(a) As z → −α 2 , Proof. Part (a) follows from (2.10) where we have to take into account that z(z + β 2 ) is the negative square root for z = −α 2 .
The contour Γ in (6.2) is split into two components, as in the proof of Theorem 2.10. We denote them by Γ 1 and Γ 2 as shown in Figure 20. We reverse the orientation on both Γ 1 ∪ Γ 2 and γ 0,1 , which does not change the values of the double integrals in (6.2).
As in the proof of Theorem 2.10 we then have Re Φ 1,2 (z) < Re Φ 1 (−α 2 ) < Re Φ 1 (w), whenever z ∈ γ 0,1 and w ∈ (Γ 1 \ {−α 2 }) ∪ Γ 2 . From these inequalities we find that it is only a neighborhood of −α 2 that contributes, and in particular the component Γ 2 that goes around [−β 2 , 0] will not contribute to the limit. Near −α 2 we introduce the change of variables z = −α 2 + α 2 c 0 sN −1/2 , w = −α 2 + α 2 c 0 tN −1/2 (6.43) with c 0 as in (6.25). Then where we used z ∼ −α 2 and the value (6.25) for c 0 . From (6.32), (6.41)-(6.42) and (6.44), we find that the leading order behavior for the sum of the two double integrals in (6.2) as N → ∞ is equal to where Σ 1 and Σ 2 are as in Figure 21. Σ 2 is the image of Γ 1 under the mapping (6.43). It is a contour going around the negative real axis. Σ 1 arises from applying (6.43) to the part of γ 0,1 in a small disk around −α 2 , and deforming and extending it so that Σ 1 is the imaginary axis with an indentation around 0. We make a further change of variables by mapping t = (t ) 2 , t ∈ iR in both integrals, and s = (s ) 2 , s ∈ Σ in the first integral, but s = (s ) 2 , s ∈ (−Σ) in the second double integral. Recall that Σ and −Σ are as in Figure 8. We find a sum of two double integrals, but the integrands are exactly the same. We pick up a Jacobian 4s t from the change of variables, and we use (s ) −1 + (t ) −1 s t (t ) 2 − (s ) 2 = 1 t − s to simplify the expression. Then the integrand on the right-hand side of (2.42) and (2.43) follows (provided we drop the accents from s and t ).
This completes the proof of Theorem 2.13.

Cusp points: proof of Theorem 2.14
Proof. The proof of Theorem 2.14 is based on the integral representation (2.37) for K gas . We already noted in Remark 2.11, that K gas (m, n; m , n ) tends to zero at an exponential rate whenever |∆m| + |∆n| → ∞. In the situation of Theorem 2.14 we have So we assume v > v . Then K gas (m, n; m , n ) still decays to 0 as N → ∞, but α −∆n increases since ∆n < 0. and the limit in (2.45) will exist, as we shown next.
We argue similar to the proof of Theorem 2.13. We start by deforming the contour γ to a contour going around the interval (−∞, −α 2 ] as before. That is, we deform γ to Γ 1 as in Figure 20 but with opposite direction. We can do this since for ∆n < 0 there is enough decay of the integrand in (2.37) as z → ∞.
After a change of variables (6.43), using the expansions (6.33) and (6.38), and finally changing s to s 2 , all as in the proof of Theorem 2.13, we find This concludes the proof of Theorem 2.14.

Acknowledgements
We thank Sunil Chhita and Kurt Johansson for many fruitful discussions. We are also grateful to Christophe Charlier for a careful reading of a preliminary manuscript. Most of this work was done when Arno Kuijlaars was a visiting professor at KTH in spring semester of 2017, funded by the Knut and Alice Wallenberg Foundation.