Biased \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2 \times 2$$\end{document}2×2 periodic Aztec diamond and an elliptic curve

We study random domino tilings of the Aztec diamond with a biased \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2 \times 2$$\end{document}2×2 periodic weight function and associate a linear flow on an elliptic curve to this model. Our main result is a double integral formula for the correlation kernel, in which the integrand is expressed in terms of this flow. For special choices of parameters the flow is periodic, and this allows us to perform a saddle point analysis for the correlation kernel. In these cases we compute the local correlations in the smooth disordered (or gaseous) region. The special example in which the flow has period six is worked out in more detail, and we show that in that case the boundary of the rough disordered region is an algebraic curve of degree eight.


Introduction
Domino tilings of the Aztec diamond, originally introduced in [12], form a popular arena for various interesting phenomena of integrable probability.A domino tiling of the Aztec diamond can be viewed as a perfect matching, also called dimer configuration, on the Aztec diamond graph.This is a particular bipartite subgraph of the square lattice (cf. Figure 2).By putting weights on the edges of the Aztec diamond graph, one defines a probability measure on the set of all perfect matchings, and hence all domino tilings, by saying that the probability of having a particular matching is proportional to the product of the weights of the edges in that matching.In recent years, several works have appeared on domino tilings of the Aztec diamond where the weights are doubly periodic.That is, they are periodic in two independent directions, and we will use the notation k × to indicate that they are kperiodic in one direction and -periodic in the other.In this paper, we will study a particular example of a 2 × 2 doubly periodic weighting that is a generalization of the model studied in [1,2,7,8,11,20].The difference is that we introduce an extra parameter that induces a bias towards horizontal dominos, and we refer to this model as the biased 2 × 2 periodic Aztec diamond.The model considered in [1,2,7,8,11,20] will be referred to as the unbiased 2 × 2 periodic Aztec diamond.
Doubly periodic weightings lead to rich behavior when the size of the Aztec diamond becomes large.The Aztec diamond can be partitioned into three regions: frozen, rough disordered (or liquid) and smooth disordered (or gaseous).They are characterized by the different local limiting Gibbs measures that one expects in these regions [22].The difference between the smooth and disordered regions is that the dimer-dimer correlations decay exponentially with their distance in the smooth disordered region and polynomially in the rough region.The three regions are clearly visible in Figure 1 where we have plotted a sample of our model for a large Aztec diamond.
From general arguments, that go back to [21], we know that the correlation functions in our model are determinantal.In order to perform a rigorous asymptotic study, one aims to find an expression for the correlation kernel that is amenable for an asymptotic analysis.For the unbiased 2 × 2 periodic Aztec diamond, a double integral representation was first found in [7] (more precisely, they were able to find the inverse Kasteleyn matrix [21]).Based on this expression, the boundary between the smooth and rough disordered region has been studied extensively in [1,2,20].Unfortunately, it is not obvious how the expression in [7] extends to the biased generalization that we consider in this paper.Instead, we follow the approach of [5].
In [5] the authors studied probability measures on particle configurations given by products of minors of block Toeplitz matrices.The biased 2 × 2 periodic Aztec diamond can be viewed as a special case of such a probability measure.The main result of [5] is an explicit double integral formula for the correlation kernel, provided one can find a Wiener-Hopf factorization for the product of the matrix-valued symbols for the block Toeplitz matrices.That Wiener-Hopf factorization can in principle be found by carrying out an iterative procedure, in which the total number of iterations is of the same order as the size of the Aztec diamond.In certain special cases, such as the unbiased 2 × 2 periodic Aztec diamond [5] and a family of 2 × k periodic weights [4], the procedure is periodic, and after a few iterations one ends up with the same parameters that one started with.This means that the Wiener-Hopf factorization has a rather simple form, and after inserting that expression in the double integral formula one obtains a suitable starting point for a saddle point analysis [5,11].However, generically, the iteration in [5] is too complicated to find simple expressions for the Wiener-Hopf factorization, and other ideas are needed.
The biased 2 × 2 periodic Aztec diamond is the simplest doubly periodic case in which it is difficult to trace the flow in [5].Our first main result is that the Wiener-Hopf factorization can alternatively be computed by following a linear flow on an explicit elliptic curve.This flow is rather simple and consists of repeatedly adding a particular point on the elliptic curve.For generic parameters, one expects the flow to be ergodic, but for special choices the flow will be periodic.We will identify a few explicit examples of these periodic cases, and perform an asymptotic study in the smooth disordered region for the general periodic situation.
The reason why the iterative procedure in our case is linearizable on an elliptic curve can be traced back to [26].In that work it was shown how an isospectral flow on certain quadratic matrix polynomials, obtained by repeatedly moving the right divisor of the polynomial with a given spectrum to the left side, is linearizable on the Jacobian (or the Prym variety) of the corresponding spectral curve.The main goal of [26] was to describe the dynamics of certain discrete analogs of classical integrable systems in terms of Abelian functions.Some of the key ideas used in that work had previously originated in constructing the so-called finite gap solutions of integrable PDEs, see their book-length exposition [3] with historic notes and references therein.The matrix case, which was most relevant for [26], had been originally developed in [9,10,16,23,24].
While our situation does not exactly fit into the formalism of [26], similar ideas do apply, and they led us to the linearization.We hope that they will also help with studying more general tiling models.
To conclude, let us mention that in [11] it was shown that the double periodicity leads to matrix-valued orthogonal polynomials.For the unbiased 2 × 2 periodic Aztec diamond, these matrix-valued orthogonal polynomials have a particularly simple structure.Somewhat surprisingly, they even have explicit integral expressions that lead to an explicit double integral representation for the correlation kernel.The expression in [11] was re-derived in [5].For the biased model it is interesting to see what the flow on the elliptic curve implies for the matrixvalued orthogonal polynomials, and if explicit expressions can be given in general and/or for the periodic case.Furthermore, it is interesting to compare our results with [6], in which matrix orthogonal polynomials were studied using an abelianization based on the spectral curve for the orthogonality weight.

Preliminaries
In this section we will introduce the dimer model that we are interested in, discuss several standard different representations from the literature and recall the determinantal structure of the correlation functions for a corresponding point processes.In our discussion we repeat necessary definitions from earlier works, in particular of [5,11,17,19], and we will make specific references to those works at several places to refer the reader for more details.We refer to [15] for a general introduction to random tilings.

A doubly periodic dimer model
and white vertices and with edges E N between black and white vertices that are neighbors in the lattice graph (i.e., that have a difference of (±1, 0) or (0, ±1)).This gives the graph on the left of Figure 2. The picture on the right of Figure 2 is a perfect matching of this bipartite graph, also called a dimer configuration.A dimer model is a probability distribution on the space of all perfect matchings M of this graph G N such that the probability of a particular matching M is proportional to where w : E → (0, ∞) is a weight function.
In this paper, we will consider the weight functions defined as is shown in Figure 3.There are two parameters α, a ∈ (0, 1].The vertical and horizontal edges with a black vertex on top or on the right all have weights a and 1, respectively.For vertical edges with a black vertex on the bottom and horizontal edges with a black vertex on the left, the weight depends on the coordinates of that black vertex.These weights are given by aα and α, or by a/α and 1/α, depending on the coordinates of the black vertex in that edge.If the vertical coordinate of that vertex is 1  2 + k for an even k, then the weights are aα and α.If the vertical coordinate of that vertex is 1  2 + k for an odd k, then we have the same weight, but with α replaced by 1/α.The distribution of the weights is thus two periodic in two different directions; edges whose coordinates differ by a multiple of (2, 2) or (2, −2) have the same weight.
Note that only the parameter α is responsible for the double periodicity.Indeed, for α = 1 the weights no longer depend on the position of the black vertex in the center in Figure 3.We will be particularly interested in the doubly periodic situation and thus in the case 0 < α < 1.The effect of the extra parameter a is that all the vertical edges are given an extra factor a. If 0 < a < 1, this makes them less likely, and the model is biased towards horizontal edges.As we will see, adding this parameter has a profound effect on the integrable structure of this model.Moreover, we will see that the special case a = 1, studied by several authors [5,8,11], is a very particular point.
An alternative way of representing matchings is by drawing dominos.Indeed, each matching is equivalent to a domino tiling by drawing rectangles around the matched vertices as is shown in Figure 4.The dominos tile a planar region known as the Aztec diamond.We distinguish between four different types of dominos called the West, East, North and South dominos.The West dominos are the vertical dominos with a black vertex on the bottom, the East dominos are the vertical dominos with a black vertex on the top, the North dominos are the horizontal dominos with a black vertex on the right and, finally, the South dominos are the horizontal dominos with a black vertex on the left.In Figure 4 these four types of dominos are the furthermost ones in the corresponding corners.
Note that the weighting that we will consider is such that all North dominos have weight 1 and all East dominos have weight a.The weight of a West domino is either aα if the vertical coordinate of the lower left corner is even, or a/α if that coordinate is odd.Similarly, the weight of a South domino is either α if the vertical coordinate of the lower left corner is even, and 1/α if that coordinate is odd.For small a > 0 we expect to see more South and North dominos, as the West and East domino have small weight.

Non-intersecting paths
A useful alternative representation, that is easily obtained from the dominos, is the representation by DR-paths [19,18,29].By drawing an upright path across each West domino, a down-right path across each East domino, a horizontal across each South domino and nothing on a North domino, we obtain the picture given in Figure 5.There are four paths leaving from the lower left side of the Aztec diamond and ending at the lower right side.The paths also cannot intersect.Clearly, the paths determine the location of the East, West and South dominos, and therewith the entire tiling.One can therefore represent each dimer configuration with a collection of non-intersecting paths.
Instead of looking directly at the DR paths, however, we will consider closely related interpretation in terms of non-intersecting paths on a different graph.The reason for this is two-fold.First, the DR paths are rather uneven in length.The bottom path is much shorter than the top path.The second reason is that it turns out to be useful to add paths so that we have an infinite number of them.The auxiliary paths will have no effect on the model, but will give a very convenient integrable structure.
We start with a directed graph G p = ({0, 1, . . ., 2N }×Z, E p ) where we draw edges between the following vertices (we use the index p in G p and E p to distinguish this graph from the (0, 0) A part of the graph is shown in Figure 6.We then fix starting points (0, −j) for j = 0, . . ., M , and endpoints (2N, −j) for j = 0, . . ., M and consider collections of paths in the directed graph that connect the starting points with the endpoints, such that no paths have a vertex in common (i.e., they never intersect).
Note that if M ≥ 2N −2 only the N top paths and the N −1 bottom paths are non-trivial, but any path in between is, due to the non-intersecting condition, necessarily a straight line.In fact, even the top N and bottom N −1 paths have parts where they are necessarily straight lines.Indeed, in the region between the lines (m, −N + m/2) and (m, −M + N + m/2) for m = 0, . . ., 2N , all the paths are necessarily horizontal.
The connection with the dimer models is the following: If we remove all paths below the line (m, −N + m/2) then the configuration that remains is equivalent to the DR paths for the domino tilings of Aztec diamond.Indeed, by further removing all horizontal parts (m, u) → (m + 1, u) for odd m and concatenating the result, we obtain the picture in the middle of Figure 7.The coordinate transform (m, u) → (m, u − m) maps the middle picture to the DR-paths shown on the right of Figure 7.
The next step is to put a probability measure on the collection of non-intersecting paths that is consistent with the dimer model from Section 2.1.To make the correspondence, we note that each up-right diagonal edge in the graph G p corresponds to a West domino, each vertical edge to an East domino, and each horizontal edge (after removing the auxiliary horizontal edges at the odd steps) corresponds to a South domino.A careful comparison with the weights for the dimer models leads us to assigning weights to the underlying directed graph as follows: the horizontal edges (m, u) → (m + 1, u) for odd m are auxiliary and have weight 1, the vertical edges correspond to East dominos and have weight a, the horizontal edges (m, u) → (m + 1, u) for even m correspond to South dominos and have weight α if u is Figure 7: From the non-intersecting paths on the graph G p to the DR paths.The middle picture is obtained by removing the horizontal steps from the paths and the graph G p .In the second transformation (m, u) → (m, u − m) we obtain the rotated DR paths.even and weight 1/α if u is odd, and, finally, the up-right edges (m, u) → (m + 1, u + 1) for m even have weight aα if u is even and weight a/α if u is odd.This is also represented in the following finite weighted graph that is the building block for the rest of G p : Then the probability of having a particular configuration of non-intersecting paths is proportional to the product of the weights of all the edges in the corresponding dimer/domino configuration.

A determinantal point process
Let us now assign a point process to the above collections of paths.We place points on these paths by taking the lowest possible vertex on each vertical section (including those of length 0), as indicated in the right panel of Figure 6, (m, u j m ) for j = 1, . . ., M, m = 0, . . ., 2N, where u j 0 = u j 2N = −j + 1 and M ≥ N .Since the top N paths uniquely determine the dimer configuration, so do the points (m, u j m ).Further, our probability measure also turns the set of points with coordinates (m, u j m ) into a point process on {0, 1, . . ., 2N } × Z.We stress that we are only interested in the points (m, u j m ) with j ≤ N − m/2 + 1, as it is those that determine the tiling.The other points are auxiliary and only added for convenience.Indeed, by a theorem of Lindström-Gessel-Viennot (see, e.g., [25,14]) the probability of a given point configuration is proportional to where T m are the transition matrices defined by . (1) We will also use the notation By the Eynard-Mehta theorem (see, e.g., [13]), the point process is determinantal, meaning that there exists a kernel such that, for any (m k , u k ) ∈ {0, . . ., 2N } × Z and k = 1, . . ., n, Now we recall that we are only interested in the top N paths, and thus we will restrict u j to be in {−N + 1, . . ., 0}.Then the marginal densities are independent of M as long as M is sufficiently large and In [5] a double integral formula for the correlation kernel K N was given.That formula involves a solution to a Wiener-Hopf factorization.
Proposition 2.1.[5, Theorem 3.1] Suppose that we can find a factorization Then the kernel K N,M has the pointwise limit K N as M → ∞ given by where |a| 2 < ρ 1 < ρ 2 < 1/|a| 2 , 1 m <m = 1 if m < m and 0 otherwise, and the integration contours are positively oriented.
Remark 2.2.Proposition 2.1 is only part of Theorem 3.1 in [5].Indeed, the kernel in (3) is called K top in [5].We note that here we already shifted coordinates compared to [5].Also, in the formulation of Theorem 3.1 in [5] one needs a second factorization A(z) = Ã− (z) Ã+ (z).However, all that is needed for Proposition 2.1 is the existence of such a factorization, and that is guaranteed by Theorem 4.8 in [5].

The Wiener-Hopf factorization
The question remains how to find a Wiener-Hopf factorization that is explicit enough to be able to use (3) as a starting point for asymptotic analysis.The idea for finding a Wiener-Hopf factorization is simple (see also [5,Section 4.4]).Write where .
Next, we compute a factorization for P 1 (z) = P 1,− (z)P 1,+ (z) and set P 2 (z) = P 1,+ (z)P 1,− (z).At each step in the procedure we thus construct a new matrix valued function P k+1 (z) = P k,+ (z)P k,− (z) constructed by switching the order of the Wiener-Hopf factorization The result is that we find a Wiener-Hopf factorization for A(z) of the form An important point is that this procedure defines a flow P 0 (z) → P 1 (z) → P 2 (z) → . . .and to obtain explicit representations for the correlation kernel in (3) we need to have a sufficiently detailed description of this flow.As was pointed out in [5,Section 4], there is a general procedure to capture this flow.Generically, the description in [5] of the flow is rather difficult to control, but for specific values it can be written explicitly.Indeed, for a = 1, the double integral formula of [11] could be reproduced.See also [4] for other cases where it was tractable.It is important to note that in the cases of both [11] and [4], the flow was periodic, which is of great help, in particular for asymptotic analysis.For the model that we consider in this paper, however, it appears difficult to control this flow for a < 1, and the point of the present paper is to give an alternative more tangible description.We will show that the flow is equivalent to translations on an explicit elliptic curve.This will also help us to track other choices of parameters for which the flow is periodic.

Main results
We now present our main results.All proofs will be postponed to Sections 5.

An elliptic curve
Consider an elliptic curve E (over R) defined by the equation where α and a are the parameters from the dimer model in Section 2.1.One easily verifies that the curve crosses the x-axis precisely three times, once at the origin and at two further intersection points in (−∞, 0).The elliptic curve has therefore two connected components, and one of those, denoted by E − , lies entirely in the left half plane.Note also that (0, 0), (a 2 , a 2 ) and (a −2 , a −2 ) are the intersection points of the curve with the line y = x.The point (a −2 , a −2 ) will be of particular interest to us.It is well known that an elliptic curve carries an Abelian group structure, and we can add points on the curve.The point at infinity serves as the identity.We will be interested in a linear flow on the curve that is constructed by repeatedly adding the point (a −2 , a −2 ) starting from the initial parameters (x 0 , y . That is, we consider the flow and + represents addition on the elliptic curve.The flow can be nicely illustrated by the geometric description of the group addition on the curve.Starting from (x j , y j ) we compute (x j+1 , y j+1 ) as follows: the straight line passing though (x j , y j ) and (a −2 , a −2 ) passes through The flow on the elliptic curve.At each step we add the point (a −2 , a −2 ).This can be geometrically represented by drawing a straight line through (a −2 , a −2 ) and (x j , y j ).This line intersects the curve at a unique third point in E − .The point (x j+1 , y j+1 ) is then obtained from the intersection point by flipping the sign of the second coordinate.
a third point and (x j+1 , y j+1 ) is the reflection of that point with respect to the x axis (in other words, we flip the sign of the y-coordinate).See also Figure 8.It can happen that the line through (x j , y j ) and (a −2 , a −2 ) is tangent to E − at point (x j , y j ).In that case, (x j+1 , y j+1 ) is just the reflection of the (x j , y j ) with respect to the x axis.Note that the initial point (x 0 , y 0 ) lies on the oval E − , and from the geometric interpretation it is easy to see that every point (x j , y j ) is on the oval E − .Our first main result is that this flow uniquely determines the correlations for the biased Aztec diamond as described in Sections 2.1-2.4 above.But before we explain that, we first discuss properties of the flow that will be of interest to us.For generic choices of the parameters one can expect the flow to be ergodic on E − , but for certain special parameters (a −2 , a −2 ) will be a torsion point.In those cases the flow is periodic.This distinction has important implications for our asymptotic analysis of the tiling model.We will therefore discuss a few examples in which (a −2 , a −2 ) is a torsion point.
First, if we assume that α = 1, then our dimer model is an example of a Schur process [27], and we know that simpler double integral formulas for its correlation kernel can be given.This should mean that our flow has a particularly simple structure.Indeed, for α = 1, the oval E − reduces to a singleton E − = {(−1, 0)}, and the flow is constant.This can also be seen directly, from the fact that the two factors in the definition of P (z) commute.
The second case of interest is the unbiased case where a = 1.In that case, (a −2 , a −2 ) = (a 2 , a 2 ), and the elliptic curve is tangent to the line y = x at that point.For general a > 0 we have the relation (a 2 , a 2 )+(a −2 , a −2 ) = (0, 0) and thus, for a = 1, we have 2(a −2 , a −2 ) = (0, 0).It is also clear that (0, 0) is a point of order 2, and thus (a −2 , a −2 ) is of order 4.This implies that our flow is periodic and returns to its initial point after 4 steps.For completeness, we See the left panel of Figure 9 for an illustration.
The next example we would like to discuss is that of an order six torsion point.This happens when The flow on the elliptic curve is given by: Indeed, after six steps we have returned to our initial point.This case is illustrated on the right panel of Figure 9.
We found the relation ( 7) by computing the division polynomial of order 6 and requiring that (a −2 , a −2 ) is a zero of this polynomial.In fact, this provides a recipe for deriving relations between a and α such that (a −2 , a −2 ) is a torsion point of order m.We recall the notion of division polynomials in Appendix B and provide such relations for m = 4, 5, 6, 7, 8.

Correlation kernel
To explain the connection between the flow on the elliptic curve and the Wiener-Hopf factorization in Proposition 2.1 we define functions a, b, d : Since x < 0 for (x, y) ∈ E − , these functions are well-defined with no poles and take strictly positive values.Consider the maps and The first main result of this paper is that the factorization ( 4) is given by P k,± (z) = P ± (σ k (x, y)).
We will discuss this claim at length in Section 4 in a slightly more general setup and we refer to that section for more details.The claim is then a special case of Theorem 4.6.Of important to us now is that it, together with Proposition 2.1, implies the following.
Theorem 3.1.The correlation kernel K N from Proposition 2.1 can be written as where and and the contours of integration are counterclockwise oriented circles with radii ρ 1 and ρ 2 such that A proof of this theorem is given in Section 5.1.If (a −2 , a −2 ) is a torsion point of order d, the flow (x, y) → σ(x, y) is periodic, and the double integral formula can be rewritten in a useful way. and Then we can rewrite (9) as where Note that in (9) we have replaced the size of the Aztec diamond N by dN .This is not necessary and the upcoming analysis can also be performed for the general case.Since the difference will only involve non-essential cumbersome bookkeeping, we feel that working with dN instead of N makes for a cleaner presentation.

Asymptotics
The representation of the correlation kernel in ( 12) is a good starting point for an asymptotic study.We will compute the microscopic process in the limit N → ∞ near the point That is, we consider the limiting behavior of the correlation kernel as N → ∞, with m, m ∈ Z fixed.Note that the first coordinate of the point ( 13) is a multiple of 2d and the second coordinate is a multiple of 2. This restriction is made for clarity purposes and is not necessary.Note also that any finite shift from ( 13) can be absorbed into the variables 2m + ε, 2m + ε , 2x + j +2x + j in (14).

The spectral curve
To perform the asymptotic analysis it is convenient to diagonalize the matrices P (w), P (z), P + (w) and The spectral curve det(P (z) − λ) = 0 can be easily computed: The curve has branch points at z = 0, z = ∞, and at the zeros of the discriminant: These zeros are negative and will be denoted by x 1 and x 2 , ordered as x 1 < x 2 < 0. With these points, we define a Riemann surface R consisting of two sheets that we connect in the usual crosswise manner along the cuts (−∞, x 1 ) and (x 2 , 0).The sheets have 0 and ∞ as common points.See also Figure 10.We will write z (j) to indicate the point z on the sheet R (j) .Then we define the square root (R(z)) 1/2 on R such that (R(z (1) )) 1/2 > 0 for z (1) > 0. The spectral curve (15) then defines a meromorphic function on R given by with poles at 0 and ∞, and zeros at (a ±2 ) (2) .The restrictions of λ to R (j) will be denoted by λ j , i.e., λ j (z) = λ(z (j) ).
Next, consider the spectral curves for P det(P These spectral curves factorize (15) in the following way.
Lemma 3.3.The equations (18), (19) for µ and ν define meromorphic functions on R such that (λ(z) for z ∈ R. Then µ has a zero at (a 2 ) (2) and a pole at 0, both of the same order d, and ν has a zero at (a −2 ) (2) and a pole at ∞, both of the same order d.
With E(z) defined by we have and Here The proof of this lemma will be given in Section 5.2 One particular consequence of this lemma is that we can simultaneously diagonalize P (z) and P (d) ± (z).In the following theorem we use this to rewrite the correlation kernel in (12).Theorem 3.4.Assume (a −2 , a −2 ) is a torsion point of order d.Set, with E(z) as in (15), Then, where γ (1,2) 2 are the unit circles with counterclockwise orientation on the sheets R 1,2 , γ 1 is a counterclockwise oriented contour inside the contour γ (1) 2 on the sheet R 1 that goes around (a 2 ) (1) and the cut [x 2 , 0], and γ The proof of this Theorem will be given in 5.3.Note that A e (w) −1 is analytic at w = a 2 (even though A e (w) is not).Moreover, λ(w) −m µ(w) N −T has a zero at w = (a 2 ) (2) of order d(N − T ) − m , and this zero cancels the pole at w = (a 2 ) (2)  in the double integral in (26).The contour γ (2) Figure 11: The contours of integration in (26).The blue contour represents γ 1 and the orange contours are the unit circles on the two different sheets.
Remark 3.5.We note that the spectral curve det (P (z) − λ) = 0 and the elliptic curve E in (5) are related.Indeed, (5) can be written as In other words, the elliptic curve E equals the spectral curve after changing the spectral variable.

Saddle point equation and classification of different regions
The representation ( 26) is a very good starting point for asymptotic analysis.To illustrate this we will perform a partial asymptotic study, based on a saddle point analysis.We note that a similar analysis has been given in [4,11].An interesting feature is that our analysis will depend on the torsion d, but in such a way that we can treat all values of d simultaneously.
To perform a saddle point analysis of (26) we need to find the saddle points and the contours of steepest descent/ascent for the action defined by This is a multi-valued function, but the differential Φ (z)dz is single valued on R. Its zeros are the saddle points for Re Φ, and we will be especially interested in them.Let C 1 be the cycle on R defined by connecting the segments (x 1 , x 2 ) on R 1 and R 2 at the end points x 1 and x 2 .Similarly, let C 2 be the cycle that combines the copies of (0, ∞) on both sheets.Proposition 3.6.The differential Φ (z)dz has simple poles at 0, (a 2 ) (1) , (1/a 2 ) (2) and ∞.There are four saddle points (i.e., the critical points where Φ (z)dz = 0) counted according to multiplicity.There are at least two distinct saddle points on the cycle C 1 .
There are always two saddle points on the cycle C 1 , but it is the location of the two other saddle points that determines the phase at the point (τ, ξ).We say that (τ, ξ) is • in the frozen region, if we have two distinct saddle points on the cycle C 2 ; Smooth Rough Frozen Figure 12: Both pictures represent a partitioning of the region into the frozen region, the rough disordered region and the smooth disordered region.The picture on the left uses the natural coordinates (τ, ξ) corresponding to the point process associated with the nonintersecting paths.The picture on the right corresponds to the coordinates for the original dimer model.In both pictures we have a 2 = α/(1 + α + α 2 ) and α = 1 2 .
• in the smooth disordered region, if we have four distinct saddle points on the cycle C 1 ; • in the rough disordered region, if there is a saddle point in the upper half plane of R 1 or R 2 ; • on the boundary between the rough and smooth disorderd regions, when this saddle point from the upper half plane coalesces with its complex conjugate on the cycle C 1 ; • on the boundary between the rough and frozen regions, when the saddle point from the upper half plane coalesces with its complex conjugate on the cycle C 2 .
We note that the terminology rough, smooth and frozen goes back to at least [22].In that work also the alternatives gaseous for smooth disordered and liquid for rough disordered were mentioned.In the subsequent literature both these terms have been used.We chose to use terminology frozen, rough and smooth disordered.The difference between these regions is in the decay of the local correlations for the local Gibbs measure.In the frozen region, the randomness disappears.In the rough disordered region, the correlations between two points decay polynomially in their distance, whereas in the smooth disorder regions these correlations decay exponentially.Our list above suggests that these different behavior can be characterized in terms of the location of the two remaining saddle points.The following theorem justifies this characterization for the smooth disordered (or gaseous) region.
Theorem 3.7.Let (τ, ξ) be in the smooth disordered region.Then Note that from (28) we see that the limiting mean density in the smooth disordered region is given by lim and the right-hand side is independent of (T, X) (as long as it is in the smooth disordered region).
It is also not difficult to see that the right-hand side of (28) decays exponentially with the distance between (m, x) and (m , x ).Indeed, for m and m fixed, the right-hand side is the (x − x )-th Fourier coefficient of a function that is analytic in an annulus.Such coefficients decay exponentially with a rate that is determined by the width of the annulus.More generally, the exponential decay follows from a steepest descent analysis for the right-hand side of (28).
The proof of Theorem 3.7 will be given in Section 5.5 and it is based on a saddle point analysis of the integral representation (26).We are confident that such a saddle point analysis can be carried out similarly for the rough disordered and frozen regions.Since it requires nontrivial effort and since a full asymptotic study is not the main focus of this paper, we do not perform such an analysis here.

The boundary of the rough disordered region
We will now show that the boundary of the rough disordered region is an algebraic curve and discuss how this curve can be found explicitly in particular cases.
We start with the following proposition.
The proof of this proposition will be given in Section 5.6.By inserting the constants (29) into Φ (z), multiplying by (z −a 2 )(z −a −2 )R(z) 1/2 , and reorganizing the equation so that all terms with R(z) 1/2 are on the right, we see that Φ (z) = 0 can be written as Before we proceed, note that z = a −2 and z = a 2 are two solutions that we just introduced by multiplying by (z − a 2 )(z − a −2 ) and are not saddle points.
By squaring both sides of (31) and multiplying by z we find a polynomial equation of degree 6 in z with coefficients that are quadratic functions of τ and ξ.Since z = a ±2 are solutions that we are not interested in, we are left with an equation of degree four.There are four solutions to this equation, and each of them corresponds to exactly one point on the surface.This confirms that we indeed have four saddle points, which was part of the statement in Proposition 3.6.
This also allows to write an equation for the rough disordered boundary.Indeed, the coefficients of this fourth degree equation will be quadratic expressions in τ and ξ.We have a third order saddle in case the discriminant vanishes.The discriminant of a polynomial of degree four is a polynomial in its coefficients of degree six.Thus, the discriminant is a polynomial in τ and ξ of degree twelve.In the explicit cases that we tried, we found, with the help of computer software, that this degree twelve curve can be factorized into a curve of degree eight and remaining factors that are not relevant.This also matches with the findings of [7] and [5] for the special case a = 1.We have, however, only been able to verify that this holds numerically in special cases (one of them we will discuss in Appendix A) and do not have a proof that it holds generally.We leave this as an interesting open problem and post the following conjecture: Conjecture 3.9.The boundary of the rough disordered region is an algebraic curve in τ and ξ of degree eight.Remark 3.10.There is another way of parametrizing the boundary.Indeed, on the two components of the boundary of the rough region we have a coalescence of saddle points on the cycles C 1 or C 2 .This means that we have a double zero of the differential Φ (z)dz.This gives a way of parametrizing these curves.Indeed, Φ (z) = Φ (z) = 0 for z ∈ C 1 or C 2 gives a linear system of equations for µ and ξ that can be easily solved.
Another interesting consequence of ( 29) is that the saddle point equation Φ (z) = 0 only depends on the order d of the torsion via the constant γ 1 in (30).However, it is even possible to replace this with another expression that does not involve d: Lemma 3.11.The constants γ 1 and γ 2 from (30) are related via where R(x) > 0 for x ∈ (x 1 , x 2 ).
The proof of this lemma will ve given in Section 5.7.By replacing the equation for γ 1 in (30) by (32) we see that we have eliminated the dependence on d from the saddle point equation, and the saddle point equation makes sense for general parameters a and α.Although the arguments that we provide in this paper use the torsion at several places, it is natural to conjecture that the saddle point analysis and its consequences can be extended in this way.In particular, we conjecture the characterization of the different phases in Section 3.3.2and Theorem 3.7 to hold under this extension.We leave this as an open problem.

Overview of the rest of the paper and the proofs
In the remaining part of this paper we will prove the main results.In Section 4 we will show that the linear flow on the elliptic curve can be used to find a Wiener-Hopf factorization in Proposition 2.1.We will do this in a more general setup than only for the biased Aztec diamond.In Section 5 we will return to the biased Aztec diamond and prove Theorem 3.1 in Section 5.1, which is by then just an identification of the parameters in the discussion of Section 4. Then Lemma 3.3 and Theorem 3.4 are proved in Sections 5.2 and 5.3, respectively.The saddle point analysis starts with proving Lemma 3.6 in Section 5.4.After that, we perform a saddle point analysis in Section 5.5 and prove Theorem 3.7.Proposition 3.8 is proved in Section 5.6 and Lemma 3.11 in Section 5.7.In Appendix A we work out the example where (a −2 , a −2 ) is a torsion point of order six.We compute the boundary of the rough disordered region, and we provide numerical results supporting the saddle point analysis of Section 5.5.Finally, in Appendix B we will show how the notion of division polynomials can be used to find algebraic relations between a and α so that (a −2 , a −2 ) is a torsion point of order d.

The flow
In this section we introduce a flow on a space of matrices that will give a Wiener-Hopf factorization in the correlation kernel.We prove that this flow is equivalent to a linear flow on an elliptic curve using translations by a fixed point on that curve.

The space
First we have to define the space of matrices that we work on.To this end, we first introduce Clearly, the determinant det P (z) of any P ∈ S is a rational function in z with poles at z = 0 and z = ∞ and no other.Also, det P (z) will have two zeros z 1 and z 2 , and we will assume that 0 < z 1 < 1 < z 2 .
Then the winding number of det P (z) with respect to the unit circle equals zero.
The flow that we will define on S will be such that det P (z) and Tr P (z) will be invariant under it.We therefore introduce the sets Naturally, c 1 , c 2 and z 1 , z 2 be expressed in terms of a ij and b ij .Indeed, and z 1 , z 2 are the solutions to det P (z) = 0. Equivalently, z 1 and z 2 can be obtained from the following equations: which, combined with the condition 0 < z 1 < 1 < z 2 , determine z 1 and z 2 uniquely.We also note that the condition 0 < z 1 < 1 < z 2 is equivalent to requiring det P (1) > 0, because det P (z) → −∞ for z ↓ 0 and z → +∞.In terms of a ij and b ij this means that the condition is equivalent to a 11 a 22 > (a 12 + b 12 )(a 21 + b 21 ).
Note that this also shows that right-hand side of the second equation in ( 35) is positive, as it should be.It should also be noted that c 1 , c 2 , z 1 and z 2 cannot take arbitrary values.For instance, we have the following result.Lemma 4.1.We have Proof.The proof follows after inserting (34) and (35) into As we will see later, the inequality (36) is sufficient to ensure that S(z 1 , z 2 , c 1 , c 2 ) = ∅.We will give an explicit parametrization of S(z 1 , z 2 , c 1 , c 2 ) in terms of part of an elliptic curve.But first, let us define a flow on S(z 1 , z 2 , c 1 , c 2 ).

Definition of the flow
We will be interested in factorization of the matrices in S of a particular form.Start by introducing the sets and It is straightforward to verify that if Proposition 4.2.Let P ∈ S(z 1 , z 2 , c 1 , c 2 ).Then there exist unique Q ± ∈ S ± such that Proof.Note that To find Q ± we have to solve By comparing the coefficients on both sides we obtain six equations for the six unkowns a, b, c, d, z 1 and z 2 .The parameters z 1 , z 2 can be found from the condition det P (z 1 ) = det P (z 2 ) = 0. Then finding the remaining equation gives This determines the factorization Because of the special structure of S ± we have uniqueness of the factorization.However, for our purposes we need an additional degree of freedom by adding a diagonal factor.Indeed, if P = Q − Q + then P − = Q − D and P + = D −1 Q + for any diagonal matrix D also provides a factorization of P such that ).We will use this additional degree of freedom by requiring that In other words, we require that the leading term in the asymptotic behavior fo P − P + and P + P − match.In order to achieve this, we define where a, b are the parameters in Q − .The flow on S(z 1 , z 2 , c 1 , c 2 ) that we wish to study is then defined by iterating the map s, i.e., the flow is defined by the recurrence It turns out it is rather complicated to keep track of this dynamics, and our goal is to describe this dynamics in a way that it is easier to grasp.

Translations on an elliptic curve
Consider the elliptic curve E over R defined by (with c 1 , c 2 > 0 and 0 < z 1 < 1 < z 2 as before) We also assume, cf.Lemma 4.1, that This inequality implies that we have three points on the curve whose y coordinate is zero, (0, 0), (−t 1 , 0) and (−t 2 , 0), with t 1 , t 2 > 0.Moreover, the curve E has two connected components one in the left half plane and the other in the right half plane.It will also be important for us that the lines y = ±x lie above and below E − , meaning that y 2 − x 2 < 0 and thus |y|/|x| < 1.Indeed, the lines y = ±x intersect E at most at three points, and we already established that these points are on E + .This implies that E − has to lie fully below or above each of these lines and since (−t 1 , 0) and (−t 2 , 0) lie below the line y = −x and above the line y = x, so does E − .There is a classical construction of addition on an elliptic curve which we will use.We can add two points (x 1 , x 1 ), (x 2 , y 2 ) ∈ E as follows: generically, the line through (x 1 , y 1 ) and (x 2 , y 2 ) intersects the elliptic curve at exactly one point (x 3 , −y 3 ).Then we define (x 1 , y 1 ) + (x 2 , y 2 ) = (x 3 , y 3 ).One exception is when (x 2 , y 2 ) = (x 1 , y 1 ) (in which case the addition becomes a doubling of the point), but this can be defined by continuity.The other exception is (x 1 , y 1 )+(x 1 , −y 1 ) which we define to be the point at infinity.The addition turns E into group with the point at infinity as zero.
We will be mostly interested in translation by (z 2 , z 2 ) on E. Observe that if (x, y) ∈ E − then (x, y) + (z 2 , z 2 ) ∈ E − .We will define the translation operator It is not hard to put this into a concrete formula.Since it will be useful to have this formula at hand, and in order to simplify arguments later, we include it in the following lemma.
Proof.The line through the point (x, y) and (z 2 , z 2 ) is given by the formula Y = λ(X −z 2 )+z 2 where λ = y−z 2 x−z 2 .By substituting this into the equation for E, moving all terms to the righthand side and collecting the coefficient of X 2 we obtain and this equals −c 2 times the sum of the three zeros of the resulting cubic equation for X.
In other words, after setting (x * , −y * ) = (x, y) + (z 2 , z 2 ) we have Thus, Now use the fact that (x, y) ∈ E to find Inserting this back into y * = λ(x * − z 2 ) + z 2 we find and further simplification shows By flipping the sign of y * we thus obtain the statement.

Equivalence of the flows
Our main point is that the flows s and σ from Definition 4.3 and Lemma 4.4 are equivalent.
We start with the following.
To establish that π is a bijection we construct the inverse map as follows.It is not difficult to see that any matrix from the general space S can be written as in the right-hand side of (40) after choosing c 1 , c 2 , z 1 , z 2 as in (34) and (35) and u, x, y as By the assumptions a ij > 0 and b ij > 0 we see that u, c 1 , c 2 , z 1 z 2 > 0, hence x < 0 and |y| < |x|.We still need to verify that (x, y) lies on the elliptic curve.But this follows from the computation of the determinant (41).Indeed, since the determinant matches with det P (z) we must have that (x, y) ∈ E. Since we already know that x < 0 we find (x, y) ∈ E − , and we have thus proved the statement.
We now come to the key point of this section.

Now we can compute
To find (u , (x , y )) such that s(π(u, (x, y)) = π(u , (x , y )) we argue as follows.From (40) we see that u is the coefficient of z in the 12-entry.This gives u = u, so the parameter u is unchanged under the flow.Then x is the zero of the 12-entry viewed as a linear function in z and thus Next, by looking at the 22-entry of P + P − we find .

Wiener-Hopf factorizations
Let P (z) ∈ S with S as defined in (33) and n ∈ N. In this paragraph we will show how the flows above can be used to find an explicit Wiener-Hopf factorization First of all, as also discussed in Section 2.4, with P k (z) = s k (P (z)) and P k (z) = P k,− (z)P k,+ (z) as in Definition 4.3 we can take and Then, by Theorem 4.6 we can obtain an explicit representation in terms of the flow on the elliptic curve.To this end, we first define the functions (cf.( 42)) Using the parametrizaton for P (z) as in (40) we then have, by Theorem 4.6, Hence, For future reference, we note that the constant pre-factor is of no interest to us and will cancel out in the integrand for the double integral formula of Proposition 2.1 for the correlation kernel.It is thus the evolution of a(σ j (x, y)) that is of importance.Next, define the function Then we have Hence,

Proofs of the main results
We now return to the model of the biased doubly periodic Aztec diamond from Section 2.1 and prove our main results.

Proof of Theorem 3.1
Proof of Theorem 3.1.We recall from Proposition 2.1 that we are interested in finding a factorization for where .
Comparing this with the setting of Section 4.4 we see that we have the special case and thus the elliptic curve can be written as The flow starts with the initial parameters (x 0 , y 0 ) = (−1, −1−α 2 1+α 2 ).The theorem is a straightforward consequence of the factorization of Section 4.5.

Proof of Lemma 3.3
Proof of Lemma 3.3.It is readily verified that (22) holds.An important observation is that This implies that P (z) d , P  23) and (24).Furthermore, note that we can rewrite ( 23) and ( 24) as with E(z) as in (21).Now the entries of E(z) and E(z) −1 are meromorphic functions for ).From (48) we then see that µ 1,2 and ν 1,2 are also meromorphic

Proof of Proposition 3.6
Proof of Proposition 3.6.One can easily see that Φ (z)dz has simple poles at 0 and ∞.On the first sheet Φ (z) takes the form and we see that we have a simple pole at (a 2 ) (1) .On the second sheet we can use the relations and the pole at (a 2 ) (2) gets canceled at the cost of a new simple pole at (a −2 ) (2) .Thus, Φ (z)dz has four simple poles at said locations and thus also four zeros (since R(z) is of genus 1).We now show that there are at least two saddle points in C 1 , which can be done using the same argument as in [11,proof of Proposition 6.4].The point is that one can show that Indeed, since ν 1 (z) and µ 1 (z) are real-valued for z ∈ (x 1 , x 2 ), so is Φ (z (1,2) ) by ( 51) and ( 52), and thus Im As for the real part, note that x 1 Φ (z (2) )dz and, since Re Φ is single-valued on R, Therefore, also the real part of the left-hand side of (53) vanishes.By combining this with the fact that Φ (z)dz is real-valued and continuous on C 1 , we see that Φ (z)dz must change sign at least two times.This means that there are at least two zeros of Φ (z)dz.(Note that this argument does not work on C 2 since Φ (z)dz has two poles on C 2 .)

Asymptotic analysis in the smooth phase
We will work out the asymptotic analysis in the smooth phase.We prepare the proof of Theorem 3.7 by first performing two steps: 1. a preliminary deformation of paths.
2. a qualitative description of the paths of steepest descent and ascent leaving from the saddle points.
After these steps, the asymptotic analysis follows by standard arguments. γ 1 Figure 14: The first deformation of contours.The contours γ (2) 1 remain untouched.The blue contour γ(1) 1 is a deformation of the contour γ 1 .By deforming the contour like this, we pick up a residue at z = w.Note also that each blue contour can be deformed through the cuts and be entirely, or partly, on the second sheet.The orange contour is allowed to pass the cuts provided one does not pass through the origin while deforming.

A preliminary deformation
We will need the following lemma on the asymptotic behavior of the integrand in (26) near the poles at 0 and ∞.Lemma 5.1.We have that Proof.The behavior near w = ∞ follows readily after observing Similarly, the behavior near w = 0 follows after observing It is important to observe that we are considering (τ, ξ) in the parallellogram defined by τ = 0, τ = 1, ξ = τ /2 and ξ = (τ − 1)/2.By (13) this means that for any x , m ∈ Z we have that for N sufficiently large, and thus the left-hand side of (54) has no poles (and no residues) for either w = 0 or w = ∞.
We proceed with the first contour deformation.The contours γ (2) 1 will be untouched, but the contour γ (1) 1 is deformed to the contour γ(1) that goes around the cut (−∞, x 1 ) in clockwise direction, as indicated in Figure 14.While deforming we pick up possible residues at the pole at w = ∞ and at w = z for z ∈ γ (1) 2 .As mentioned above, with our choice of parameters there is no pole at w = ∞.The pole at w = z has a residue for z ∈ γ (1) 2 , and this gives us a contribution: This means that we can rewrite (26) as A e (w) This finishes the preliminary deformation.Before we continue to the steepest descent analysis, we mention that by (54) and (55), the integrand with respect to w has no pole at w = 0.This means that we can deform the contour γ (2) 1 to lie partly or even entirely on the first sheet.The integrand with respect to z does have poles at z = 0 and z = ∞ and thus, the contours γ 2 may be deformed over the surface R but cannot pass through the origin (without picking up a residue).

Description of the paths of steepest descent/ascent
By definition, we have four saddle points in the cycle C 1 , and in the interior of the smooth region these are distinct and simple.By viewing Re Φ as a function on the cycle C 1 , these saddle points will correspond to the locations of the two local minima and two local maxima of Re Φ.We will denote the saddles associated to the local minima by s 1 and s 3 , and the local maxima by s 2 and s 4 .We take the indexing such that when traversing the cycle C 1 starting from x 1 to x 2 on R 1 and then from x 2 to x 1 on R 2 , the first saddle point one encounters is s 1 , then s 2 and so on.Note also that both local minima are neighbors to both local maxima (on the cycle C 1 ) and therefore Re Φ(s 1,3 ) < Re Φ(s 2,4 ).
We proceed by giving a description of the contours of steepest descent and ascent for Re Φ leaving from these four saddles.Since each saddle point is simple, there will be two paths of steepest descent and two path of steepest ascent leaving from them.It is straightforward that the segment of C 1 between s 2j and s 2j±1 is a path of steepest descent for Re Φ leaving from s 2j and a path of steepest ascent leaving from s 2j±1 .What remains, is to identify the paths of steepest descent leaving from s 1 and s 3 and the paths of steepest ascent from s 2 and Figure 15: An illustration of the hypothetical case that the two saddle points s 2 and s 4 connect to the same saddle point (1/a 2 ) (2) .In this case, the four paths together form a contractible curve and enclose the (shaded) region that contains s 3 but not the cycle C 2 .This means that the steepest descent paths leaving s 3 will have to cross the paths from s 2 or s 4 , which is not possible.It is important to note that, even though Φ (z) is single-valued, Φ(z) is a multi-valued function, and we cannot replace the condition simply with Im Φ(z) = Im Φ(s j ).Indeed, because of the logarithmic terms the imaginary part Im Φ(z) jumps whenever we cross a cut (which we did not specify) for a logarithm.The real part Re Φ(z), however, is single-valued, and this will be important.We will also need the behavior near the logarithmic singularities of Re Φ(z) at z = 0, z = (a 2 ) (1) , z = (a −2 ) (2) and at z = ∞: from (27) (see also (51) and ( 52) and Re Φ((a 2 ) (1) By analyticity of Φ (z) the paths are a finite union of analytic arcs and ultimately have to end up at some special points.The only options for such special points are other saddle points or the poles of Φ .It takes only a short argument to exclude possibility that they will connect to another saddle point.Indeed, since Re Φ(s 1,3 ) < Re Φ(s 2,4 ) it is impossible to connect s 2j to s 2j±1 in this way.Moreover, it is obvious that a path of steepest descent (or ascent) from the global minimum (or maximum) cannot be connected to any other saddle point, hence s 1 cannot be connected to s 3 and s 2 cannot be connected to s 4 .We conclude so far that the path of steepest descent leaving from s 1,3 and the paths of steepest ascent from s 2,4 will have to end up at the four simple poles of Φ. From (57) we further deduce s 1,3 connect to ∞ and 0, and from (57) we find that s 2,4 connect to the simple poles at (a 2 ) (1)  and (1/a 2 ) (2) .
Observe that none of these paths can cross each other, since by analyticity of Φ such a crossing necessarily would be a saddle point (which we already excluded) or a pole.For the same reason, since Im z Φ (s)ds is constant on the cycles, the paths can never cross the cycles C 1 and C 2 as the point where it would cross would necessarily be a saddle point.The only point the paths have in common with the cycles are the saddle point at C 1 they started at and the pole at C 2 they end in.Hence, a path that starts at a saddle point on R 1 and continues in the upper half plane of R 1 always remain in the union of the upper half plane of R 1 and the lower half plane of R 2 glued together along the cuts (−∞, x 1 ) ∪ [x 1 , x 2 ].This important property shows, in particular, that steepest ascent/descent paths do not wind around the poles of Φ (z).
Next we argue that the paths of steepest ascent leaving from s 2 and s 4 cannot end in the same pole.Indeed, if they would, then all these four paths together would form a closed loop that is contractible and hence cuts the Riemann surface into two parts, as illustrated in Figure 15.The cycle C 2 lies fully in one of the parts.But s 1 and s 3 are in different parts, and hence there must be one of them that is in a part that is different from the part that contains the cycle C 2 .The steepest descent path that leaves that saddle point has to cross the closed loop in order to end up at a pole on C 2 , which is not possible, and we arrive at a contradiction.This means that the steepest ascent paths from s 2 and s 4 have to end up at different poles, one saddle connects to (a 2 ) (1) and the other to (1/a 2 ) (2) .A similar argument shows that one of the saddle points s 1 and s 3 connect, via steepest descent paths, to 0 and the other to ∞.
Let us summarize our findings above: Proposition 5.2.The steepest descent paths leaving from s 1 and s 3 and steepest ascent path from s 2 and s 4 form simple closed loops on R, such that no two loops intersect, each loop intersects each cycle C 1 and C 2 exactly once, and it does so at a saddle point for Re Φ in C 1 and a pole for Φ (z) in C 2 .
See Figure 16 and its caption for an illustration.

Proof of Theorem 3.7
Now we are ready for the Proof of Theorem 3.7.The starting point is the representation of the kernel after the preliminary deformation as given in (56).To prove the result, all that is needed is to show that the double integral tends to zero as N → ∞.This is rather straightforward after one has realized that the contours of the preliminary deformation strongly resemble the paths of steepest descent and ascent for the saddle point s j .Indeed, the two contours γ(1) 1 and γ (2) 1 can be deformed to go through the saddle points s 1 and s 3 and follow the paths of steepest descent, and the contours γ 2 can be deformed to the path of steepest ascent ending in z = (a 2 ) (1) and z = (1/a 2 ) (2) respectively.During this deformation, no additional residues are being picked up, and standard saddle point arguments show that there exists c > 0 such that as N → ∞.This finishes the proof.

Proof of Proposition 3.8
Before we come to the proof of Proposition 3.8 we need a few lemmas.We use the notation x for the largest integer smaller than x.
Figure 16: The seven pictures illustrate the possible locations (schematically) of the paths of steepest descent and ascent leaving from the four saddle points on the cycle C 1 in the smooth region.In (a) and (b) we have all four saddle points on the first sheet, in pictures (c)-(f) we have three saddle points on the first sheet and in picture (g) we have one point on the first sheet.It is also possible that all four saddle points are on the second sheet, and in that case the picture is similar to that of (a) and (b) with the two sheets switched (but keeping the poles a ±2 in place and slightly adjusting the contours accordingly).Similarly, for the case of three saddle point on the second sheet.All pictures can be reconstructed started from the picture in (a) by continuous deformations.For example, (b) can be obtained by moving the right most saddle point (and the orange contour) in (a) over the cycle C 1 , first passing the branch point x 1 to the second sheet and then passing the branch point x 2 back to the first sheet to become the left most saddle point at (b).The pictures (c) and (d) can be obtained from (a) by moving the right most and the left most points respectively to the second sheet, etc.We did not check whether all configurations indeed occur and perhaps some cases can be excluded, but our arguments hold for any of the above configurations.
Proof of Proposition 3.8.By (59) we have with R(z) = a 2 (z − x 1 )(z − x 2 )/z and the square root is taken so that R(z) > 0 for z > 0.
Here p(z) is a polynomial of degree at most d/2 and q is a polynomial of degree (d − 1)/2 with a zero at z = 0. Observe that ν 1 (z)ν 2 (z) can be written as where Since q(0) = 0 and R (z) has double pole at z = 0, r 1 and r 2 are polynomials and r 2 (0) = 0.The degree of r 1 and r 2 is at most d.By replacing R(z) by − R(z) in the derivation above we also find Therefore, we can write Since ν 2 (z) has a zero of order d at z = a −2 , this means that both r 1 and r 2 have a zero of order d − 1 at z = a −2 .This implies that ν 1 (z)ν 2 (z) can be written as where γ j ∈ R, for j = 1, 2, 3, are some real constants.By a similar reasoning, one can show that for some real parameters γj , for j = 1, 2, 3.
The next step is to compute the values of the constants γ j , γj for j = 1, 2, 3. To this end, add (65) and (66) to obtain On the other hand, we easily compute from ( 15)-( 17) that where we used z 2 R (z) = 4a 2 (z 2 − 1) in the last step.Comparing (67) and (68) leads to the following two equations: and From (70) we find and from (69) we find γ1 = γ 2 , γ2 = γ 2 (72) and Thus far, we have derived the first two identities in (30).

Proof of Lemma 3.11
Proof of Lemma 3.11.The cycle condition (53) implies that (using ( 29)) By a change of variable x → 1/x we find (using ( 16)) x 2 x 1 and after substituting this into the first integral, (75) reduces to after which the statement easily follows.

A Example: torsion point of order six
Let us now assume Then, as discussed in Section 3.1, (a −2 , a −2 ) is a torsion point of order six.Here we will compute the spectral curves for µ and ν, and derive a degree eight equation for the boundary of the rough disordered region.We will also show, numerically, how the steepest descent/ascent path can be chosen when (τ, ξ) are in the center of the diamond.

A.1 The flow on the matrices
The linear flow on the elliptic curve is given already in (8).The flow on the matrices and their decomposition can then be traced giving:

A.2 Contours of steepest descent/ascent
We will plot the contours of steepest descent/ascent for Re Φ, with Φ as in (27), for the special values a This is the midpoint of the Aztec diamond, where we indeed expect a smooth disordered region to appear.Indeed, with this choice of parameters, we find four saddle points on the cycle C 1 .Two of them are on the first sheet: and the other two on the second sheet: The branch points are at Observe that s 1 is close to x 1 and s 3 is close to x 2 , which we found to be typical for any choice of parameters.This has the unfortunate consequence that in numerical illustrations the saddle points s 1 and s 3 are hard to distinguish from the branch points x 1 and x 2 , respectively.The contours of steepest descent/ascent leaving from the saddle points locally coincide with the level lines of Im Φ(z).The problem is that Φ(z) has logarithmic terms making Im Φ(z) multi-valued and plotting the level sets Im Φ(z) = Im Φ(s j ) does not give the correct result.For this reason, we compute the vectorfield given by (Re Φ , − Im Φ ) and compute the streamlines using the function Streamplot in Mathematica.The results are given in Figure 17.
From Figure 17 one can see that the paths of ascent/descent are indeed as illustrated schematically in Figure 16(g).On the first sheet, the paths of steepest descent leaving from s 1 end up at ∞ and the paths of steepest ascent leaving from s 2 end up at a 2 .On the second sheet, the paths of steepest descent leaving from s 3 end up at 0 and the paths of steepest ascent leaving from s 4 end up at a −2 .The statement for s 3 is not immediately obvious from Figure 17(b) and this is why we zoom in around s 3 as in Figure 17(c).

A.3 Boundary of the rough disordered region
The last result for the case of order six that we will present here, is an explicit expression for the boundary of the rough disordered region.We follow the procedure indicated in Section 3.4 with a and α related by (76) and values for γ j , j = 1, 2, 3 as in (81).We start with (31), square both sides and remove the additional factors (z − a 2 )(z − a −2 )/z.Then we obtain an equation that is polynomial in z and of order four.The values of (τ, ξ) where the discriminant vanishes lead to a double saddle point and this will be the boundary of the liquid region.
The discriminant has degree twelve in τ and ξ, and for general parameters a (and α related by (76)) the expression is rather long.In order to obtain a shorter expression it will be convenient to perform the following change of variables (τ, ξ) = ((q + v + 1)/2, q/2).These coordinates change the parallellogram in the left panel into the tilted square in the right panel of Figure 12.
The first factor is only zero for (q, v) = (0, 0), and what is left is the boundary for the smooth disordered region.The second factor is an ellipse.Finally, for α ↓ 0 (and hence a ↓ 0 simultaneously) the curve reduces to 0 = (1 − 9q 2 ) 2 (4 − 9v 2 ) 2 , which gives a rectangular shape.In this case, there is no rough disordered region, but only a frozen region and a smooth disordered region.
In Figure 18 we have plotted the boundary of the rough disordered region for several particular values of a.

B Computation of torsion points
In Section 3.1 we gave a few examples of particular choices for the parameters α and a such that (1/a 2 , 1/a 2 ) is a torsion point.Here we will indicate how one can find such examples by recalling the notion of division polynomials.This is a standard construction for finding the torsion subgroups of the elliptic curve.We will base our discussion on [28, p. 105].
The torsion subgroup of order m consists of all zeros of ψ m , which are called division polynomials.Note that we are not looking for the entire subgroup, but for situations where 1/a 2 , (a + 1/a)(α + 1/α) 2a 2 is a torsion point.After substituting this point into ψ m we find a rational function in a and α.This rational function will have several zeros, but we are only interested in the zeros that satisfy 0 < α < 1 and 0 < a ≤ 1.For instance, for m = 3 we find the following equation This equation has no solutions such that 0 < α < 1 and 0 < a ≤ 1, and thus a third order torsion point cannot occur.
With the help of a computer code, we found the following equations that give proper solutions such that we have a torsion point of order m = 4, . . .Each term on the right-hand side is a factor of ψ m .

Figure 1 :
Figure 1: A sampling of the biased doubly periodic for a large Aztec diamond.The West and South dominos are colored yellow, and the North and East dominos are colored blue.The three different regions are clearly visible, with the smooth disordered region in the middle, surrounded by the rough disordered region and frozen regions in the corners.

Figure 2 :
Figure 2: The left picture is the bipartite graph G N , with N = 4, and the right picture is a perfect mathching of G N .

aα a 1 αα a 1 1 αFigure 3 :
Figure 3: The weights on the edges of G N .

Figure 4 :Figure 5 :
Figure 4: The right picture is the domino representation of the dimer configuration on the left

Figure 6 :
Figure 6: The left figure shows the underlying graph G p .The right figure shows the graph G p and a collection of non-intersecting paths starting in (0, −j) and ending in (2N, −j) for j = 0, . . .M , with N = 4 and M = 5.

Figure 9 :
Figure9: The picture on the left illustrates the flow in case (a −2 , a −2 ) is a torsion point of order four.The picture on the right shows the flow in case that point has order six.

Figure 10 :
Figure 10: The two sheeted Riemann surface R. The dashed lines represent the cycles C 1 and C 2 .
oriented contour on the sheet R 2 inside the contour γ (2) 2 that goes around the cut [x 2 , 0].See also Figure 11.

Figure 13 : 2 , z 2 = 2 and c 2 1 /c 2 = 7 .
Figure 13: An example with parameter z 1 = 1 2 , z 2 = 2 and c 2 1 /c 2 = 7.Under strict inequality in (39) we always have an oval in the left half plane.In case we have equality, the oval has shrunk to a point.

s 4 .
These paths will continue in the lower and upper half planes of the sheets R j and they are further characterized by the condition that Im z s j Φ (z)dz = 0.
(a) Level lines for Im Φ.The paths of steepest descent for Re Φ from s 1 connect to infinity on the first sheet.The paths of steepest ascent from s 2 end up in a 2 .(b) Level lines for Im Φ.The paths of steepest ascent from s 2 end up in a −2 .To see the paths of steepest descent from s 3 we need to zoom in.(c)Zooming in on the segment (x 2 , 0) shows that the path of steepest descent for Re Φ starting from the saddle point s 3 end in the origin.