* If on the contrary, $w$ is “on the right” of $b$ (there are poles of type $\gamma$, but then none of type $\delta$), So the only poles of $g_{b,w}$ are either some $\alpha_j$ or some $\gamma_j$. The indicator function is equal to 1, so the final formula is $-\frac{1}{2i\pi}$ times an integral surrounding all the parameters from $\alpha$ and all the parameters from $\gamma$ of $g_{b,w}$ but since all the poles are encloses, the total is again equal to 0.

In this case again, $J_{b,w} = 0$ for all $w \in Q_N$.

* We can repeat the same arguments when $w$ is in the Aztec diamond, and $b$ is either in $Q_W$ or $Q_E$, by looking now first at the variable $u$ in the double integral. As a consequence we have the four cyan zero blocks in this infinite inverse.

* We now look at the formula when:

- $w$ is fixed in $Q_N$ (resp. $Q_S$) and $b$ is allowed to exit the Aztec region to be in $Q_W$, $Q_E$ and $Q_S$ (resp. and $Q_N$);

- $b$ is fixed in $Q_W$ (resp. $Q_E$) and $w$ is allowed to exit the Aztec region to be in $Q_N$, $Q_S$ and $Q_E$ (resp. and $Q_W$).

The previous arguments again show that in these situations, the value of $J_{b,w}$ is also 0. This is represented by the orange blocks.

* We then continue computing by taking $b$ in the Aztec region, and $w \in Q_W$: We can compute iteratively the value on every white vertex column by column, which leads at each step a linear equation with a single unknown. We do the same for $w \in Q_E$. This corresponds to the green block. We do the same for $w$ in the Aztec region, and $b \in Q_W \cup Q_N$. This corresponds to the purple region.

* We then discuss the entries when $b$ and $w$ are in the same quadrant. Fix for example $w$ in $Q_N$. We know from computations above that the value of $J_{b,w}$ is 0 if $b$ is just outside of $Q_N$ (but connected to a white vertex $w'$ of $Q_N$). We can then solve iteratively for all the values of $J_{b,w}$ by solving $\sum_b K_{w',b} A_{b,w} = \delta_{w',w}$, as at each step there will be always a white vertex $w'$ such that we know the value of $J_{b,w}$ for two black neighbors of $w'$, but the value of the third value is not computed yet, but is completely determined by the linear equation above. In particular, we obtain that $J_{b,w} = 0$ if $b$ is not in the light cone of $w$, and is equal to $\frac{1}{K_{w,b}}$ if $(w, b)$ is the icy edge attached to $w$. We proceed in the same way for $Q_S$ and, after exchanging the roles of the colors, for $Q_E$ and $Q_W$. This determines the entries of the four diagonal blocks $D_N, D_S, D_W, D_E$.

* Finally, we proceed as before to compute the entries when $b \in Q_N \cup Q_S$ and $w \in Q_W \cup Q_E$, by propagating known values on the boundary to the bulk by the linear equation. The actual entries of that block will not matter for our purposes. $\square$

It turns out that this inverse of the Kasteleyn operator has a probabilistic meaning. More precisely, the determinant of minors of $J$ are related to local statistics of the Boltzmann measure on the Aztec diamond, as follows:
**Lemma 24.** Let $e_1 = (w_1, b_1), \dots, e_k = (w_k, b_k)$ be distinct edges of the Aztec diamond of size $n$, and $e_{k+1} = (w_{k+1}, b_{k+1}), \dots, e_l = (w_l, b_l)$ be distinct edges in the complement of $A_n$ in $G_n$. Then

$$ \left( \prod_{i=1}^{l} \tilde{K}_{w_i, b_i} \right) \det_{1 \le i,j \le l} J_{b_i, w_j} = \mathbb{P}_{A_n}[e_1, \dots, e_k \text{ are dimers}] \times \prod_{i=k+1}^{l} \mathbb{I}_{\{e_i \text{ is icy}\}} \quad (32) $$

where $\mathbb{P}_{A_n}$ is the Boltzmann probability measure on perfect matchings of the Aztec diamond of size $n$, computed via Theorem 5.

*Proof.* We can assume that the edges $e_{k+1}, \dots, e_l$ are ordered in such a way that they are grouped in four blocks corresponding to the four quadrants: $Q_N$, $Q_S$, $Q_W$, $Q_E$ (for edges joining a vertex of a quadrant to a vertex of a neighboring quadrant, consider them as belonging to either quadrant). In each block, order the edges in such a way that an edge $e'$ in the light-cone of a vertex of another edge $e$ of the same quadrant should come after $e$ in the list.

Due to the structure of the matrix $J$, and the order we put on the edges, the entries of the submatrix

$$ (J_{b_i, w_j})_{1 \le i, j \le l} $$

in columns indexed by white vertices in the quadrants $Q_N$ and $Q_S$ are zero, except maybe in the lower triangular part of the diagonal blocks. The diagonal entries for these columns are $\mathbb{I}_{\{(w_i, b_i) \text{ is icy}\}} \frac{1}{\tilde{K}_{w_i, b_i}}$. We can therefore expand the determinant

$$ \det_{1 \le i, j \le l} J_{b_i, w_j} $$

first along the columns corresponding to edges in the quadrants $Q_N$ and $Q_S$. Then, for similar reasons, we can expand the along the rows corresponding to edges in the quadrants $Q_W$ and $Q_E$. Therefore,

$$ \det_{1 \le i, j \le l} J_{b_i, w_j} = \det_{1 \le i, j \le k} J_{b_i, w_j} \times \prod_{i=k+1}^{l} \frac{\mathbb{I}_{\{(w_i, b_i) \text{ is icy}\}}}{\tilde{K}_{w_i, b_i}}, $$

which, once multiplied on both sides by $(\prod_{i=1}^l \tilde{K}_{w_i, b_i})$ is exactly Equation (32). $\square$

<figure><img src="image_14.png" /></figure>

The structure of the graph $G_n$ (as of the other infinite graphs from [Spe07]) is such that constraints at infinity propagate in all the quadrants up to the boundary of the Aztec region: suppose we are given a dimer configuration of $G_n$, which is such that in an annular region (large enough to contain the Aztec diamond part in its interior), all the dimers are icy edges. Then necessarily, all the icy edges in the four quadrants (and no other in these parts of the graph) are present in this dimer configuration. Such a configuration is called *ultimately frozen*. Using the determinantal formula from Theorem 5 with $\tilde{K}$ and $J$ turns out to define a Gibbs measure on dimer configurations on $G_n$, supported on ultimately frozen configurations.
**Lemma 24.** Let $e_1 = (w_1, b_1), \dots, e_k = (w_k, b_k)$ be distinct edges of the Aztec diamond of size $n$, and $e_{k+1} = (w_{k+1}, b_{k+1}), \dots, e_l = (w_l, b_l)$ be distinct edges in the complement of $A_n$ in $G_n$. Then

$$ \left( \prod_{i=1}^{l} \tilde{K}_{w_i, b_i} \right) \det_{1 \le i,j \le l} J_{b_i, w_j} = \mathbb{P}_{A_n}[e_1, \dots, e_k \text{ are dimers}] \times \prod_{i=k+1}^{l} \mathbb{I}_{\{e_i \text{ is icy}\}} \quad (32) $$

where $\mathbb{P}_{A_n}$ is the Boltzmann probability measure on perfect matchings of the Aztec diamond of size $n$, computed via Theorem 5.

*Proof.* We can assume that the edges $e_{k+1}, \dots, e_l$ are ordered in such a way that they are grouped in four blocks corresponding to the four quadrants: $Q_N$, $Q_S$, $Q_W$, $Q_E$ (for edges joining a vertex of a quadrant to a vertex of a neighboring quadrant, consider them as belonging to either quadrant). In each block, order the edges in such a way that an edge $e'$ in the light-cone of a vertex of another edge $e$ of the same quadrant should come after $e$ in the list.

Due to the structure of the matrix $J$, and the order we put on the edges, the entries of the submatrix

$$ (J_{b_i, w_j})_{1 \le i, j \le l} $$

in columns indexed by white vertices in the quadrants $Q_N$ and $Q_S$ are zero, except maybe in the lower triangular part of the diagonal blocks. The diagonal entries for these columns are $\mathbb{I}_{\{(w_i, b_i) \text{ is icy}\}} \frac{1}{\tilde{K}_{w_i, b_i}}$. We can therefore expand the determinant

$$ \det_{1 \le i, j \le l} J_{b_i, w_j} $$

first along the columns corresponding to edges in the quadrants $Q_N$ and $Q_S$. Then, for similar reasons, we can expand the along the rows corresponding to edges in the quadrants $Q_W$ and $Q_E$. Therefore,

$$ \det_{1 \le i, j \le l} J_{b_i, w_j} = \det_{1 \le i, j \le k} J_{b_i, w_j} \times \prod_{i=k+1}^{l} \frac{\mathbb{I}_{\{(w_i, b_i) \text{ is icy}\}}}{\tilde{K}_{w_i, b_i}}, $$

which, once multiplied on both sides by $(\prod_{i=1}^l \tilde{K}_{w_i, b_i})$ is exactly Equation (32). $\square$

<figure><img src="image_16.png" /></figure>

The structure of the graph $G_n$ (as of the other infinite graphs from [Spe07]) is such that constraints at infinity propagate in all the quadrants up to the boundary of the Aztec region: suppose we are given a dimer configuration of $G_n$, which is such that in an annular region (large enough to contain the Aztec diamond part in its interior), all the dimers are icy edges. Then necessarily, all the icy edges in the four quadrants (and no other in these parts of the graph) are present in this dimer configuration. Such a configuration is called *ultimately frozen*. Using the determinantal formula from Theorem 5 with $\tilde{K}$ and $J$ turns out to define a Gibbs measure on dimer configurations on $G_n$, supported on ultimately frozen configurations.
**Proposition 25.** Equation (32) from the previous lemma defines a determinantal point process on edges of $G_n$, which gives a Gibbs measure on the set of dimers configurations of $G_n$ for the specification given by Fock's weights on edges.

This probability measure $\mathbb{P}$ has the property to be frozen outside of the Aztec diamond: every icy edge appear with probability one, and the other edges of the four quadrants appear with probability zero, in a random configuration sampled according to that measure. In other words, $\mathbb{P}$ is the product measure of the Dirac mass supported on the set of icy edges in the complement of the Aztec diamond and the Boltzmann measure on dimer configurations of the Aztec diamond $\mathbb{P}_{A_n}$ given by Fock's weights.

*Proof.* The collection of relations (32) for all finite subsets of edges of $G_n$ is a consistent set of finite-dimensional probability distributions, by the multilinearity of the determinant.

By Kolmogorov's extension theorem, there is a unique probability measure on the subset of edges of $G_n$ with those finite dimensional marginals, which is automatically a determinantal process, supported on dimer configurations of $G_n$.

The fact that the configuration is frozen outside of the Aztec diamond is a consequence of Equation (32). This implies that it automatically satisfies the Dobrushin-Lanford-Ruelle [Dob68, LIR69] condition for large enough annuli, containing the Aztec diamond in their interior. Therefore, it is automatically Gibbs. $\square$

This inverse $J$ of the Kasteleyn matrix $\tilde{K}$ is of different nature than the inverses defined [BCdT23a]. Whereas the measures constructed from the latter would correspond to “almost linear” height profile, the one constructed in the proposition above has a different, extremal slope in each of the quadrants.

In both cases, the inverses are constructed as contour integrals involving the family of functions $g_{b,w}(u)$ in the kernel of $\tilde{K}$. This raises the question of constructing other inverses with that mechanism for other boundary conditions at infinity, for this family of graphs or generalizations, and further to classify inverses of $\tilde{K}$ with a probabilistic meaning.

# 6 Limit shapes

The formula for the inverse of the Kasteleyn matrix is very well suited for asymptotics analysis.

In this section, we discuss results about limit shapes which can be obtained directly from the analysis of this formula. Most of them are already present in the literature. We mostly focus on a setting where parameters $\alpha, \beta, \gamma, \delta$ are chosen in a periodic way (which is more general than having periodic weights on edges, see [BCdT23a, Section 4]).

From far away all typical random dimer configuration of a large Aztec look almost the same. This statement can be made rigorous by looking at the height function [Thu90],
**Proposition 25.** Equation (32) from the previous lemma defines a determinantal point process on edges of $G_n$, which gives a Gibbs measure on the set of dimers configurations of $G_n$ for the specification given by Fock's weights on edges.

This probability measure $\mathbb{P}$ has the property to be frozen outside of the Aztec diamond: every icy edge appear with probability one, and the other edges of the four quadrants appear with probability zero, in a random configuration sampled according to that measure. In other words, $\mathbb{P}$ is the product measure of the Dirac mass supported on the set of icy edges in the complement of the Aztec diamond and the Boltzmann measure on dimer configurations of the Aztec diamond $\mathbb{P}_{A_n}$ given by Fock's weights.

*Proof.* The collection of relations (32) for all finite subsets of edges of $G_n$ is a consistent set of finite-dimensional probability distributions, by the multilinearity of the determinant.

By Kolmogorov's extension theorem, there is a unique probability measure on the subset of edges of $G_n$ with those finite dimensional marginals, which is automatically a determinantal process, supported on dimer configurations of $G_n$.

The fact that the configuration is frozen outside of the Aztec diamond is a consequence of Equation (32). This implies that it automatically satisfies the Dobrushin-Lanford-Ruelle [Dob68, LIR69] condition for large enough annuli, containing the Aztec diamond in their interior. Therefore, it is automatically Gibbs. $\square$

This inverse $J$ of the Kasteleyn matrix $\tilde{K}$ is of different nature than the inverses defined [BCdT23a]. Whereas the measures constructed from the latter would correspond to “almost linear” height profile, the one constructed in the proposition above has a different, extremal slope in each of the quadrants.

In both cases, the inverses are constructed as contour integrals involving the family of functions $g_{b,w}(u)$ in the kernel of $\tilde{K}$. This raises the question of constructing other inverses with that mechanism for other boundary conditions at infinity, for this family of graphs or generalizations, and further to classify inverses of $\tilde{K}$ with a probabilistic meaning.

# 6 Limit shapes

The formula for the inverse of the Kasteleyn matrix is very well suited for asymptotics analysis.

In this section, we discuss results about limit shapes which can be obtained directly from the analysis of this formula. Most of them are already present in the literature. We mostly focus on a setting where parameters $\alpha, \beta, \gamma, \delta$ are chosen in a periodic way (which is more general than having periodic weights on edges, see [BCdT23a, Section 4]).

From far away all typical random dimer configuration of a large Aztec look almost the same. This statement can be made rigorous by looking at the height function [Thu90],
and saying that this random rescaled height function converges in probability to a con-
tinuous deterministic function, *the limit shape*, which maximizes a certain functional.

This has been established for the uniform measure on dimer configurations of simply
connected subgraphs of the square lattice [CKP01], generalized to more general periodic
weights [KOS06, Kuc17].

In the original setup of the uniform measure for the Aztec diamond, Jokush Propp
and Shor [JPS98] proved that the behaviour of the limit shape varies depending on the
position in the domain:

* Outside of the inscribed circle (the *frozen regions*): the limit shape is linear, and the corresponding dominos display a brickwall pattern, which has different orientations in each corner.

* inside the inscribed circle (the *liquid region*): the slope of the height function varies continuously, and all type of edges have a positive probability of appearance.

The inscribed circle separates *frozen regions* from a *temperate* (or *liquid*, using the terminology of [KOS06]). It is referred to as the *arctic circle*, and the main result of [JPS98] is known as the *arctic circle theorem*.

This has been extended to periodic weights for which the liquid region is an ellipse [Joh02], where the fundamental domain still contains a single pair of white and black vertices. Then it has been obtained for the 2-periodic case [CJ16] and biased variant [BD23], where the underlying spectral curve has genus 1, and more recently by Berggren and Borodin [BB23] for a generic arbitrary genus spectral curve, but with additional assumptions which seem purely technical.

We claim that Formula (17) allows us to recover and go beyond the results cited above: having the inverse Kasteleyn operator as an explicit contour integral is particularly well-suited for asymptotic analysis: it is then possible to extract from standard saddle point analysis the limiting behaviour of local probabilities, and from them reconstruct the rescaled expected height function which would give the limit shape.

We first give in Section 6.1 a short derivation of the arctic circle (ellipse) theorem, by the saddle point method, as presented in [BBC+17] for the uniform measure, inspired by its use in [OR03] for plane partitions. Then, in Section 6.2 and 6.3, we explore some limit shape results for Fock's weights given by genus 0 and 1 M-curves. Finally we discuss in Section 6.4 how the geometric arguments of [BB23] can be adapted to this more general context to give a less restrictive result, see also the forthcoming paper [BBS].

**6.1 A short derivation of the arctic circle theorem**

A short derivation of the arctic circle theorem for uniform weights is presented in [BBC+17].
We give briefly here a variant for Kenyon's *critical* (genus 0) weights, where all the pa-
rameters $\alpha_j$ (resp. $\beta_j, \gamma_j, \delta_j$) are equal to a single value $\alpha$ (resp. $\beta, \gamma, \delta$) satisfying the
and saying that this random rescaled height function converges in probability to a con-
tinuous deterministic function, *the limit shape*, which maximizes a certain functional.

This has been established for the uniform measure on dimer configurations of simply
connected subgraphs of the square lattice [CKP01], generalized to more general periodic
weights [KOS06, Kuc17].

In the original setup of the uniform measure for the Aztec diamond, Jokush Propp
and Shor [JPS98] proved that the behaviour of the limit shape varies depending on the
position in the domain:

* Outside of the inscribed circle (the *frozen regions*): the limit shape is linear, and the corresponding dominos display a brickwall pattern, which has different orientations in each corner.

* inside the inscribed circle (the *liquid region*): the slope of the height function varies continuously, and all type of edges have a positive probability of appearance.

The inscribed circle separates *frozen regions* from a *temperate* (or *liquid*, using the terminology of [KOS06]). It is referred to as the *arctic circle*, and the main result of [JPS98] is known as the *arctic circle theorem*.

This has been extended to periodic weights for which the liquid region is an ellipse [Joh02],
where the fundamental domain still contains a single pair of white and black vertices.
Then it has been obtained for the 2-periodic case [CJ16] and biased variant [BD23],
where the underlying spectral curve has genus 1, and more recently by Berggren and
Borodin [BB23] for a generic arbitrary genus spectral curve, but with additional assump-
tions which seem purely technical.

We claim that Formula (17) allows us to recover and go beyond the results cited above: having the inverse Kasteleyn operator as an explicit contour integral is particularly well-suited for asymptotic analysis: it is then possible to extract from standard saddle point analysis the limiting behaviour of local probabilities, and from them reconstruct the rescaled expected height function which would give the limit shape.

We first give in Section 6.1 a short derivation of the arctic circle (ellipse) theorem, by the saddle point method, as presented in [BBC+17] for the uniform measure, inspired by its use in [OR03] for plane partitions. Then, in Section 6.2 and 6.3, we explore some limit shape results for Fock's weights given by genus 0 and 1 M-curves. Finally we discuss in Section 6.4 how the geometric arguments of [BB23] can be adapted to this more general context to give a less restrictive result, see also the forthcoming paper [BBS].

**6.1 A short derivation of the arctic circle theorem**

A short derivation of the arctic circle theorem for uniform weights is presented in [BBC+17].
We give briefly here a variant for Kenyon's *critical* (genus 0) weights, where all the pa-
rameters $\alpha_j$ (resp. $\beta_j, \gamma_j, \delta_j$) are equal to a single value $\alpha$ (resp. $\beta, \gamma, \delta$) satisfying the
following cyclic order

$$ \alpha < \gamma < \beta < \delta. $$

By [CKP01], we know that the rescaled height function converges in probability, for the uniform topology, to a deterministic continuous function, with a gradient contained in some polygon. Our goal is to determine here the interface between the region where the gradient of the limiting height function is extremal and where it is not.

The expected gradient of the height function is directly related to the probability that an edge of given type/orientation at that position is a dimer. Let us look at a particular edge $e = (w, b)$ crossed by train-tracks with parameters $\delta$ and $\alpha$, where $w$ has coordinates $(w_x, w_y) = (2i, 2j+1)$ and $b$ has coordinates $(b_x, b_y) = (2i-1, 2j)$. By Equation (26), the probability to see this edge in a random dimer configuration is given by:

$$ K_{w,b} K_{b,w}^{-1} = (\alpha - \delta) \frac{1}{(2i\pi)^2} \int_{\mathcal{C}_2} \int_{\mathcal{C}_1} \frac{g_{b,0}(u) g_{0,w}(v)}{v-u} \left( \frac{(u-\beta)(v-\delta)}{(u-\delta)(v-\beta)} \right)^n du dv $$

where

$$ g_{0,b}(u) = \frac{(u-\delta)^{i-1}(u-\alpha)^j}{(u-\gamma)^i(u-\beta)^j}, $$

$$ g_{0,w}(v) = \frac{(v-\gamma)^i(v-\beta)^j}{(v-\delta)^i(v-\alpha)^{j+1}}. $$

We are interested in the behaviour of the integrand when $n$ is large, and $(i/n, j/n) \to (x, y)$, which can be rewritten in this asymptotic regime as:

$$ \exp(n(F(u; x, y) - F(v; x, y) + o(1))) \quad (33) $$

where

$$ F(u; x, y) = y \log(u - \alpha) + (1 - y) \log(u - \beta) - x \log(u - \gamma) - (1 - x) \log(u - \delta). $$

The critical points of $F$ are given by a quadratic equation in $u$: there are therefore two of them (counted with multiplicity). In order to apply the saddle point, one needs to move continuously the contour to make them pass through the critical points in the direction for which the critical point will be a maximum (resp. a minimum) for the real part of $F$ for the variable $u$ (resp. for $v$). In that configuration, the double integral tends to 0 exponentially fast as $N$ goes to infinity. See for example [OR03] where this technique has been introduced in the context of tilings. When the two critical points are not on the real line, one needs to make one of the contour cross the other, at least partially, to make the two contours pass through the two critical points (and cross orthogonally). By doing so, we get an additional contribution, given by the integral along the path between the two critical point of the residue of the integrand when $u = v$, which becomes the main contribution (the remaining double integral goes to zero by the saddle point method). See Figure 14.

The corresponding probability of the corresponding edge converges to a number strictly between 0 and 1: we are in the so-called *liquid region*, or *rough phase*. On the contrary, when the two critical points are both real, the integral becomes trivial, and the corresponding probability is either 0 or 1: we are then in a *frozen region*.
following cyclic order

$$ \alpha < \gamma < \beta < \delta. $$

By [CKP01], we know that the rescaled height function converges in probability, for the uniform topology, to a deterministic continuous function, with a gradient contained in some polygon. Our goal is to determine here the interface between the region where the gradient of the limiting height function is extremal and where it is not.

The expected gradient of the height function is directly related to the probability that an edge of given type/orientation at that position is a dimer. Let us look at a particular edge $e = (w, b)$ crossed by train-tracks with parameters $\delta$ and $\alpha$, where w has coordinates $(w_x, w_y) = (2i, 2j+1)$ and b has coordinates $(b_x, b_y) = (2i-1, 2j)$. By Equation (26), the probability to see this edge in a random dimer configuration is given by:

$$ K_{w,b} K_{b,w}^{-1} = (\alpha - \delta) \frac{1}{(2i\pi)^2} \int_{\mathcal{C}_2} \int_{\mathcal{C}_1} \frac{g_{b,0}(u) g_{0,w}(v)}{v-u} \left( \frac{(u-\beta)(v-\delta)}{(u-\delta)(v-\beta)} \right)^n du dv $$

where

$$ g_{0,b}(u) = \frac{(u-\delta)^{i-1} (u-\alpha)^j}{(u-\gamma)^i (u-\beta)^j}, $$

$$ g_{0,w}(v) = \frac{(v-\gamma)^i (v-\beta)^j}{(v-\delta)^i (v-\alpha)^{j+1}}. $$

We are interested in the behaviour of the integrand when $n$ is large, and $(i/n, j/n) \to (x, y)$, which can be rewritten in this asymptotic regime as:

$$ \exp(n(F(u; x, y) - F(v; x, y) + o(1))) \quad (33) $$

where

$$ F(u; x, y) = y \log(u - \alpha) + (1 - y) \log(u - \beta) - x \log(u - \gamma) - (1 - x) \log(u - \delta). $$

The critical points of $F$ are given by a quadratic equation in $u$: there are therefore two of them (counted with multiplicity). In order to apply the saddle point, one needs to move continuously the contour to make them pass through the critical points in the direction for which the critical point will be a maximum (resp. a minimum) for the real part of $F$ for the variable $u$ (resp. for $v$). In that configuration, the double integral tends to 0 exponentially fast as $N$ goes to infinity. See for example [OR03] where this technique has been introduced in the context of tilings. When the two critical points are not on the real line, one needs to make one of the contour cross the other, at least partially, to make the two contours pass through the two critical points (and cross orthogonally). By doing so, we get an additional contribution, given by the integral along the path between the two critical point of the residue of the integrand when $u = v$, which becomes the main contribution (the remaining double integral goes to zero by the saddle point method). See Figure 14.

The corresponding probability of the corresponding edge converges to a number strictly between 0 and 1: we are in the so-called *liquid region*, or *rough phase*. On the contrary, when the two critical points are both real, the integral becomes trivial, and the corresponding probability is either 0 or 1: we are then in a *frozen region*.
<figure><img src="image_15.png" /><figcaption>Figure 14: Left: representation of the two contours of integration when $\Sigma$ is the Riemann sphere for the inverse Kasteleyn matrix of Section 6.1. The curve $A_0$ is the unit circle in this case. Right: deformation of the contours in the case when the two critical points of $F$ are not on the unit circle (they are thus on different sides), to be in position to apply the saddle point method. Deforming the contours like this creates unwanted contribution to the integral of the residue at $u = v$ along the hatched segment joining the two critical points of $F$, which has to be subtracted.</figcaption></figure>

The *arctic curve* describing the transition from the liquid region to a frozen phase is thus given by the set of coordinates $(x, y)$ for which the two critical points merge on the real axis.

In our case, these two critical points merge when $x$ and $y$ satisfy the following equation:

$$r^2(x - \frac{1}{2})^2 + r^2(y - \frac{1}{2})^2 + 2r(2 - r)(x - \frac{1}{2})(y - \frac{1}{2}) = r - 1$$

where

$$r = \frac{(\beta - \alpha)(\delta - \gamma)}{(\beta - \gamma)(\delta - \alpha)} = (\alpha, \gamma; \beta, \delta) \in (1, +\infty)$$

is the cross ratio of the four points $\alpha, \gamma, \beta, \delta$.

We know by [KO06] that two sets of isoradial critical weights correspond are gauge equivalent if the corresponding train-track parameters are related by a Möbius transformation preserving the unit circle. The cross-ratio $r$ being Möbius invariant, it is expected that the limit shape if a function of this parameter.

When $r = 2$ (for which a representative is when $\alpha, \gamma, \beta, \delta$ cut the unit circle in four equal arcs), the corresponding Boltzmann measure is uniform, and the arctic curve is a circle. A cross-ratio $r \neq 2$ can be obtained by putting different weight on NE/SW and NW/SE dominos. Choosing the parameters $\alpha, \gamma, \beta, \delta$ to be $\pm e^{\pm i\theta/2}$, the ratio between the weights is given by $\tan \frac{\theta}{2}$ and $r = (\cos \frac{\theta}{2})^{-2}$. In this case, the arctic curve is an ellipse, as proved by Johansson [Joh02, Theorem 2.4].
<figure><img src="image_17.png" /><figcaption>Figure 14: Left: representation of the two contours of integration when $\Sigma$ is the Riemann sphere for the inverse Kasteleyn matrix of Section 6.1. The curve $A_0$ is the unit circle in this case. Right: deformation of the contours in the case when the two critical points of $F$ are not on the unit circle (they are thus on different sides), to be in position to apply the saddle point method. Deforming the contours like this creates unwanted contribution to the integral of the residue at $u = v$ along the hatched segment joining the two critical points of $F$, which has to be subtracted.</figcaption></figure>

The *arctic curve* describing the transition from the liquid region to a frozen phase is thus given by the set of coordinates $(x, y)$ for which the two critical points merge on the real axis.

In our case, these two critical points merge when $x$ and $y$ satisfy the following equation:

$$r^2(x - \frac{1}{2})^2 + r^2(y - \frac{1}{2})^2 + 2r(2 - r)(x - \frac{1}{2})(y - \frac{1}{2}) = r - 1$$

where

$$r = \frac{(\beta - \alpha)(\delta - \gamma)}{(\beta - \gamma)(\delta - \alpha)} = (\alpha, \gamma; \beta, \delta) \in (1, +\infty)$$

is the cross ratio of the four points $\alpha, \gamma, \beta, \delta$.

We know by [KO06] that two sets of isoradial critical weights correspond are gauge equivalent if the corresponding train-track parameters are related by a Möbius transformation preserving the unit circle. The cross-ratio $r$ being Möbius invariant, it is expected that the limit shape if a function of this parameter.

When $r = 2$ (for which a representative is when $\alpha, \gamma, \beta, \delta$ cut the unit circle in four equal arcs), the corresponding Boltzmann measure is uniform, and the arctic curve is a circle. A cross-ratio $r \neq 2$ can be obtained by putting different weight on NE/SW and NW/SE dominos. Choosing the parameters $\alpha, \gamma, \beta, \delta$ to be $\pm e^{\pm i\theta/2}$, the ratio between the weights is given by $\tan \frac{\theta}{2}$ and $r = (\cos \frac{\theta}{2})^{-2}$. In this case, the arctic curve is an ellipse, as proved by Johansson [Joh02, Theorem 2.4].
<figure><img src="image_16.png" /><figcaption>Figure 15: Arctic regions for (1-)periodic critical weights, for various values of the cross-ratio $r$ between $\alpha, \beta, \gamma$ and $\delta$. Here $r = (\cos \frac{\theta}{2})^{-2}$ and $\theta = k \frac{\pi}{20}$ for $k = 1, \dots, 10$. The green circle corresponding to the classical arctic circle for the uniform measure, for $k = 10$.</figcaption></figure>

## 6.2 The (periodic) critical Aztec diamond has no gaseous phase

We assume now that the parameters $\alpha, \beta$ (resp. $\gamma, \delta$) are periodic with period $l$ (resp. $k$), for some $k, l \ge 1$.

The integrand in the formula describing the probability of a single edge has the same form, as Equation (33), where now the function $F(u; x, y)$ is given by the formula:

$$F(u; x, y) = \frac{1}{l} \sum_{j=1}^{l} [y \log(u - \alpha_j) + (1 - y) \log(u - \beta_j)] \\ - \frac{1}{k} \sum_{j=1}^{k} [x \log(u - \gamma_j) + (1 - x) \log(u - \delta_j)].$$

The equation for critical points of $F(\cdot; x, y)$ can be rewritten as:

$$\frac{1}{l} \sum_{j=1}^{l} \left[ \frac{y}{u - \alpha_j} + \frac{1 - y}{u - \beta_j} \right] = \frac{1}{k} \sum_{j=1}^{k} \left[ \frac{x}{u - \gamma_j} + \frac{1 - x}{u - \delta_j} \right],$$

which by putting everything above the same denominator, becomes a polynomial equation in $u$ of degree $2k + 2l - 2$. It has thus $2k + 2l - 2$ complex roots, counted with multiplicity.

On the other hand, writing $u = e^{2i\bar{u}}, \alpha_j = e^{2i\bar{\alpha}_j}, \dots$, and using the fact that

$$\frac{2iu}{u - \alpha} = \frac{e^{i(\bar{u} - \bar{\alpha})}}{\sin(\bar{u} - \bar{\alpha})} = i + \cot(\bar{u} - \bar{\alpha}),$$
<figure><img src="image_18.png" /><figcaption>Figure 15: Arctic regions for (1-)periodic critical weights, for various values of the cross-ratio $r$ between $\alpha, \beta, \gamma$ and $\delta$. Here $r = (\cos \frac{\theta}{2})^{-2}$ and $\theta = k \frac{\pi}{20}$ for $k = 1, \dots, 10$. The green circle corresponding to the classical arctic circle for the uniform measure, for $k = 10$.</figcaption></figure>

## 6.2 The (periodic) critical Aztec diamond has no gaseous phase

We assume now that the parameters $\alpha, \beta$ (resp. $\gamma, \delta$) are periodic with period $l$ (resp. $k$), for some $k, l \ge 1$.

The integrand in the formula describing the probability of a single edge has the same form, as Equation (33), where now the function $F(u; x, y)$ is given by the formula:

$$F(u; x, y) = \frac{1}{l} \sum_{j=1}^{l} [y \log(u - \alpha_j) + (1 - y) \log(u - \beta_j)] \\ - \frac{1}{k} \sum_{j=1}^{k} [x \log(u - \gamma_j) + (1 - x) \log(u - \delta_j)].$$

The equation for critical points of $F(\cdot; x, y)$ can be rewritten as:

$$\frac{1}{l} \sum_{j=1}^{l} \left[ \frac{y}{u - \alpha_j} + \frac{1 - y}{u - \beta_j} \right] = \frac{1}{k} \sum_{j=1}^{k} \left[ \frac{x}{u - \gamma_j} + \frac{1 - x}{u - \delta_j} \right],$$

which by putting everything above the same denominator, becomes a polynomial equation in $u$ of degree $2k + 2l - 2$. It has thus $2k + 2l - 2$ complex roots, counted with multiplicity.

On the other hand, writing $u = e^{2i\bar{u}}, \alpha_j = e^{2i\bar{\alpha}_j}, \dots$, and using the fact that

$$\frac{2iu}{u - \alpha} = \frac{e^{i(\bar{u} - \bar{\alpha})}}{\sin(\bar{u} - \bar{\alpha})} = i + \cot(\bar{u} - \bar{\alpha}),$$
the same critical equation has the following form:

$$ \frac{1}{l} \sum_{j=1}^{l} [y \cot(\bar{u} - \bar{\alpha}_j) + (1-y) \cot(\bar{u} - \bar{\beta}_j)] = -\frac{1}{k} \sum_{j=1}^{k} [x \cot(\bar{u} - \bar{\gamma}_j) + (1-x) \cot(\bar{u} - \bar{\delta}_j)] $$

when $\bar{u}$ is real, both sides are real. The left-hand side takes all real values between $-\infty$ and $+\infty$ when $\bar{u}$ is in $(\bar{\alpha}_j, \bar{\alpha}_{j+1})$ or $(\bar{\beta}_j, \bar{\beta}_{j+1})$, for some $1 \le j \le l-1$, whereas the right-hand side stays bounded on those intervals. As a consequence, there is at least a solution of the critical equation on each of these intervals. Reasoning in the same way on the intervals $(\bar{\gamma}_j, \bar{\gamma}_{j+1})$, $(\bar{\delta}_j, \bar{\delta}_{j+1})$, $1 \le j \le k-1$, exchanging the roles of the left- and right-hand sides, this gives in total at least $2k+2l-4$ real roots. Therefore, there are at most two complex, non-real roots, for the equation in $\bar{u}$, which have to be complex conjugated. This means that for the original equation in $u = e^{2i\bar{u}}$, the two extra solutions are not on the unit circle and the reciprocal of the complex conjugate of one is equal to the other.

Repeating the same saddle-point analysis as above, one can see that the probability of the considered edge, will have a limit, which will be non trivial in $(0, 1)$ if and only the two extra critical points are not on the unit circle. When this is the case, let us call $u_0$ the one inside the unit disc. Then one sees that we can find a unique point $(x, y)$ which has $u_0$ (and the inverse of its complex conjugate) as critical points for $F$, by solving the system of two real linear equations obtained by separating the real and imaginary part of

$$ \frac{\partial}{\partial u} F(u_0; x, y) = 0. $$

whereas, the edge will be frozen if all the critical points are real. The argument that this linear system has a unique solution, which is furthermore in the square $[0, 1]$ is an slight adaptation of the proof of [BB23, Theorem 4.9].

The transition between the two regimes occur when the two extra complex critical points merge at some $u_0 = e^{2i\bar{u}_0}$ on the unit circle, where they become zeros of the first, and second derivative of $F(\cdot; x, y)$. The *arctic curve* separating the two regimes is obtained finding $x$ and $y$ such that

$$ \left\{ \begin{array}{l} \frac{\partial}{\partial u} F(u_0; x, y) = 0 \\ \\ \frac{\partial^2}{\partial u^2} F(u_0; x, y) = 0 \end{array} \right. \qquad (34) $$

for any point $u_0$ on the unit circle $S^1$. Since $F$ is an affine function of $x$ and $y$, this amounts to solving again a linear system in $x$ and $y$ with coefficients which are explicit functions of $u_0$, yielding an explicit parametrization of the arctic curve $u_0 \in S^1 \mapsto (x(u_0), y(u_0))$.

The result above can be summarized in the following proposition:

**Proposition 26.** For periodic “genus 0” weights, with arbitrary fundamental domain,

* The map from the unit disc to the liquid region is a homeomorphism.
the same critical equation has the following form:

$$ \frac{1}{l} \sum_{j=1}^{l} [y \cot(\bar{u} - \bar{\alpha}_j) + (1-y) \cot(\bar{u} - \bar{\beta}_j)] = -\frac{1}{k} \sum_{j=1}^{k} [x \cot(\bar{u} - \bar{\gamma}_j) + (1-x) \cot(\bar{u} - \bar{\delta}_j)] $$

when $\bar{u}$ is real, both sides are real. The left-hand side takes all real values between $-\infty$ and $+\infty$ when $\bar{u}$ is in $(\bar{\alpha}_j, \bar{\alpha}_{j+1})$ or $(\bar{\beta}_j, \bar{\beta}_{j+1})$, for some $1 \le j \le l-1$, whereas the right-hand side stays bounded on those intervals. As a consequence, there is at least a solution of the critical equation on each of these intervals. Reasoning in the same way on the intervals $(\bar{\gamma}_j, \bar{\gamma}_{j+1})$, $(\bar{\delta}_j, \bar{\delta}_{j+1})$, $1 \le j \le k-1$, exchanging the roles of the left- and right-hand sides, this gives in total at least $2k + 2l - 4$ real roots. Therefore, there are at most two complex, non-real roots, for the equation in $\bar{u}$, which have to be complex conjugated. This means that for the original equation in $u = e^{2i\bar{u}}$, the two extra solutions are not on the unit circle and the reciprocal of the complex conjugate of one is equal to the other.

Repeating the same saddle-point analysis as above, one can see that the probability of the considered edge, will have a limit, which will be non trivial in $(0, 1)$ if and only the two extra critical points are not on the unit circle. When this is the case, let us call $u_0$ the one inside the unit disc. Then one sees that we can find a unique point $(x, y)$ which has $u_0$ (and the inverse of its complex conjugate) as critical points for $F$, by solving the system of two real linear equations obtained by separating the real and imaginary part of

$$ \frac{\partial}{\partial u} F(u_0; x, y) = 0. $$

whereas, the edge will be frozen if all the critical points are real. The argument that this linear system has a unique solution, which is furthermore in the square $[0, 1]$ is an slight adaptation of the proof of [BB23, Theorem 4.9].

The transition between the two regimes occur when the two extra complex critical points merge at some $u_0 = e^{2i\bar{u}_0}$ on the unit circle, where they become zeros of the first, and second derivative of $F(\cdot; x, y)$. The *arctic curve* separating the two regimes is obtained finding $x$ and $y$ such that

$$ \left\{ \begin{array}{l} \frac{\partial}{\partial u} F(u_0; x, y) = 0 \\ \\ \frac{\partial^2}{\partial u^2} F(u_0; x, y) = 0 \end{array} \right. \qquad (34) $$

for any point $u_0$ on the unit circle $S^1$. Since $F$ is an affine function of $x$ and $y$, this amounts to solving again a linear system in $x$ and $y$ with coefficients which are explicit functions of $u_0$, yielding an explicit parametrization of the arctic curve $u_0 \in S^1 \mapsto (x(u_0), y(u_0))$.

The result above can be summarized in the following proposition:

**Proposition 26.** For periodic “genus 0” weights, with arbitrary fundamental domain,

* The map from the unit disc to the liquid region is a homeomorphism.
<figure><img src="image_17.png" /><figcaption>Figure 16: The arctic curve for critical periodic Aztec diamond (with period 3 vertically, and 2 horizontally).</figcaption></figure>

* the arctic curve has a single component and has an explicit parametrization by trigonometric (or rational) functions. It is thus a real algebraic curve of genus 0. In particular there is no gaseous phase.

* The tangency points with the sides of the of the square [0, 1] at the bottom (resp. top, left, right) to $u_0$ taking one of the values of $\alpha$ (resp. $\beta, \gamma, \delta$).

The genus 0 weights are not generic, and is not covered by the results of Berggren and Borodin [BB23]. The statement about the tangency points has been noticed by Dan Betea for Schur measures corresponding to the Aztec diamond, when parameters are periodic (which is essentially Stanley's weights in the periodic setting) [Bet15].

The discussion above is mainly about the arctic curve separating the liquid and the frozen region. Actually, the saddle point analysis allows to find the frequency of each type of domino near any point of the liquid region, from which we can reconstruct the expected slope, and then the limit shape.

## 6.3 The elliptic case

We briefly study here the case of elliptic weights, with periodic train-track parameters. The surface $\Sigma$ is now a torus $C/(Z + \tau Z)$ where $\tau \in i\mathbb{R}_^*$.

First, as in Section 6.1, let us consider the case where the parameters $\alpha_j$ (resp. $\beta_j, \gamma_j, \delta_j$) are equal to the same value $\alpha$ (resp. $\beta, \gamma, \delta$) belonging to $\mathbb{R}/Z$.

Rewriting the probability of a single edge using (25), with the same notation as for the critical case, we get:

$$K_{w,b}K_{b,w}^{-1} = \int_{C_1} \int_{C_2} \exp(n(F(u;x,y) - F(v;x,y)) + o(1)) \frac{du dv}{\theta_1(\pi(v-u))}$$

where now

$$F(u; x, y) = y \log \theta_1(\pi(u - \alpha)) + (1 - y) \log \theta_1(\pi(u - \beta)) \\ - x \log \theta_1(\pi(u - \gamma)) - (1 - x) \log \theta_1(\pi(u - \delta)).$$
<figure><img src="image_19.png" /><figcaption>Figure 16: The arctic curve for critical periodic Aztec diamond (with period 3 vertically, and 2 horizontally).</figcaption></figure>

* the arctic curve has a single component and has an explicit parametrization by trigonometric (or rational) functions. It is thus a real algebraic curve of genus 0. In particular there is no gaseous phase.

* The tangency points with the sides of the of the square [0, 1] at the bottom (resp. top, left, right) to $u_0$ taking one of the values of $\alpha$ (resp. $\beta, \gamma, \delta$).

The genus 0 weights are not generic, and is not covered by the results of Berggren and Borodin [BB23]. The statement about the tangency points has been noticed by Dan Betea for Schur measures corresponding to the Aztec diamond, when parameters are periodic (which is essentially Stanley's weights in the periodic setting) [Bet15].

The discussion above is mainly about the arctic curve separating the liquid and the frozen region. Actually, the saddle point analysis allows to find the frequency of each type of domino near any point of the liquid region, from which we can reconstruct the expected slope, and then the limit shape.

## 6.3 The elliptic case

We briefly study here the case of elliptic weights, with periodic train-track parameters. The surface $\Sigma$ is now a torus $C/(Z + \tau Z)$ where $\tau \in i\mathbb{R}_+^*$.

First, as in Section 6.1, let us consider the case where the parameters $\alpha_j$ (resp. $\beta_j, \gamma_j, \delta_j$) are equal to the same value $\alpha$ (resp. $\beta, \gamma, \delta$) belonging to $\mathbb{R}/Z$.

Rewriting the probability of a single edge using (25), with the same notation as for the critical case, we get:

$$K_{w,b}K_{b,w}^{-1} = \int_{C_1} \int_{C_2} \exp(n(F(u;x,y) - F(v;x,y)) + o(1)) \frac{du dv}{\theta_1(\pi(v-u))}$$

where now

$$F(u; x, y) = y \log \theta_1(\pi(u - \alpha)) + (1 - y) \log \theta_1(\pi(u - \beta)) \\ - x \log \theta_1(\pi(u - \gamma)) - (1 - x) \log \theta_1(\pi(u - \delta)).$$
<figure><img src="image_18.png" /><figcaption>Figure 17: Three arctic curves for a torus with $\tau = i\sqrt{3}$. Left: the limit shape for the unbiased 2 × 2-periodic Aztec diamond ($\alpha, \beta, \gamma, \delta$ cut the circle into four equal arcs). Center: the limit shape for the biased weights with $\rho = \frac{1}{6}$ ($\alpha = 0, \gamma = \frac{1}{6}, \beta = \frac{1}{2}, \delta = \frac{5}{6}$) and the same value of $\tau$. Right: $\alpha = 0, \beta = \frac{3}{5}, \gamma = \frac{3}{10}, \delta = \frac{17}{20}$</figcaption></figure>

The differential of $F$ is a meromorphic 1-form on the torus $\Sigma$. It has a divisor of degree $2g - 2 = 0$. Since $F$ is periodic in the horizontal direction of the torus, the integral of $dF$ along the cycle $\tau/2 + \mathbb{R}/\mathbb{Z}$ is equal to 0. Moreover, since $dF$ is real on this cycle, the intermediate value theorem implies that $dF$ has at least two zeros on this cycle. Since it has four simple poles, located at $\alpha, \beta, \gamma$ and $\delta$, it must have two additional zeros.

These extra zeros can be either both on the same real component of the torus (in a frozen or gaseous phase), or both non-real and symmetric, and complex conjugated (in the liquid phase).

The boundary of the liquid phase is thus obtained by looking at the transition between the two regimes, when the two extra critical points merge into a double critical point. As in the critical case, asking for $u$ to be a double critical point of $F$ yields a linear system in $x$ and $y$ with coefficients that are given in terms of elliptic functions of $u$. When $u$ runs along the two real components of the torus, we get two closed curves: the outer curve separates the liquid phase from the frozen ones in the corners; the inner one separates the liquid phase from the gas bubble near the center.

Note that having periodic (even constant here) train-track parameters does not imply that the edge weights are periodic [BCdT23a, Section 4]. In order for the edge weights to be periodic, for a fundamental domain containing two white and two black vertices, the parameters $\alpha, \beta, \gamma, \delta$ should satisfy:

$$ \alpha - \beta + \gamma - \delta = 0 \mod 1, \qquad \alpha - \beta - \gamma + \delta = 0 \mod 1. \qquad (35) $$

which means that $\beta - \alpha = \pm \frac{1}{2}$, and $\delta - \gamma = \pm \frac{1}{2}$. This case contains in particular the 2 × 2 periodic Aztec diamond, studied by

* Chhita and Johansson [CJ16] by taking $\alpha = 0, \beta = \frac{1}{2}, \gamma = \frac{1}{4}, \delta = \frac{3}{4}$ and letting $\tau$ vary on the imaginary axis,
<figure><img src="image_20.png" /><figcaption>Figure 17: Three arctic curves for a torus with $\tau = i\sqrt{3}$. Left: the limit shape for the unbiased $2 \times 2$-periodic Aztec diamond ($\alpha, \beta, \gamma, \delta$ cut the circle into four equal arcs). Center: the limit shape for the biased weights with $\rho = \frac{1}{6}$ ($\alpha = 0, \gamma = \frac{1}{6}, \beta = \frac{1}{2}, \delta = \frac{5}{6}$) and the same value of $\tau$. Right: $\alpha = 0, \beta = \frac{3}{5}, \gamma = \frac{3}{10}, \delta = \frac{17}{20}$</figcaption></figure>

The differential of $F$ is a meromorphic 1-form on the torus $\Sigma$. It has a divisor of degree $2g - 2 = 0$. Since $F$ is periodic in the horizontal direction of the torus, the integral of $dF$ along the cycle $\tau/2 + \mathbb{R}/\mathbb{Z}$ is equal to 0. Moreover, since $dF$ is real on this cycle, the intermediate value theorem implies that $dF$ has at least two zeros on this cycle. Since it has four simple poles, located at $\alpha, \beta, \gamma$ and $\delta$, it must have two additional zeros.

These extra zeros can be either both on the same real component of the torus (in a frozen or gaseous phase), or both non-real and symmetric, and complex conjugated (in the liquid phase).

The boundary of the liquid phase is thus obtained by looking at the transition between the two regimes, when the two extra critical points merge into a double critical point. As in the critical case, asking for $u$ to be a double critical point of $F$ yields a linear system in $x$ and $y$ with coefficients that are given in terms of elliptic functions of $u$. When $u$ runs along the two real components of the torus, we get two closed curves: the outer curve separates the liquid phase from the frozen ones in the corners; the inner one separates the liquid phase from the gas bubble near the center.

Note that having periodic (even constant here) train-track parameters does not imply that the edge weights are periodic [BCdT23a, Section 4]. In order for the edge weights to be periodic, for a fundamental domain containing two white and two black vertices, the parameters $\alpha, \beta, \gamma, \delta$ should satisfy:

$$ \alpha - \beta + \gamma - \delta = 0 \mod 1, \qquad \alpha - \beta - \gamma + \delta = 0 \mod 1. \qquad (35) $$

which means that $\beta - \alpha = \pm \frac{1}{2}$, and $\delta - \gamma = \pm \frac{1}{2}$. This case contains in particular the $2 \times 2$ periodic Aztec diamond, studied by

* Chhita and Johansson [CJ16] by taking $\alpha = 0, \beta = \frac{1}{2}, \gamma = \frac{1}{4}, \delta = \frac{3}{4}$ and letting $\tau$ vary on the imaginary axis,
<figure><img src="image_19.png" /><figcaption>Figure 18: Left: a 2x2 example with $\alpha = \{0.744904, 0.448131\}$, $\beta = \{2.35367, 1.58096\}$, $\gamma = \{1.21721, 0.983251\}$, $\delta = \{3.05117, 2.38214\}$. Right: a 3x2 example where $\alpha = \{0.599648, 0.0928766, 0.288755\}$, $\beta = \{2.29566, 1.7928, 2.20249\}$, $\gamma = \{1.18328\}$, $\delta = \{3.06226, 2.74699\}$. The fact that there is a single contact point on the left boundary (whereas there are two on the right side) is a consequence of the fact that $\gamma$ is 1-periodic (constant) whereas $\delta$ is 2-periodic.</figcaption></figure>

* Borodin and Duits [BD23], for a specific value of $t$, which in fact does not play a role in the computation of the limit shape.

These are the cases depicted on the left and middle of Figure 17. An example where the periodicity condition is *not* satisfied is shown on the right of that same figure. Note that in this case, the picture does not have axial symmetries anymore.

One can extend this study when the parameters $\alpha, \beta, \gamma$, and $\delta$ are $k \times l$ period, like in Section 6.2. With the additional constraint that the corresponding edge weights are periodic, this covers in particular the case $k = 2, l = 2$ in the study of Berggren and Borodin [BB23] (without the technical assumption they have on the weights to guarantee that the vertical tentacles of the amoeba of the corresponding spectral curve are separated), but also the extension to the case where the periodicity condition for edge weights is not met, like the one depicted on Figure 18. One can relate here, like in the critical case, the number of cusps along bottom, left, top, right side of the square with the number of distinct values for $\alpha, \gamma, \beta, \delta$ respectively, whereas there are four of them along the liquid/gas interface.

## 6.4 Higher genus

The study above can be done on surfaces with higher genus, with train-track parameters assigned in a periodic way, even though it can be technically challenging to work with prime forms to control the saddle point method. Imposing periodicity for weights
<figure><img src="image_21.png" /><figcaption>Figure 18: Left: a 2x2 example with $\alpha = \{0.744904, 0.448131\}$, $\beta = \{2.35367, 1.58096\}$, $\gamma = \{1.21721, 0.983251\}$, $\delta = \{3.05117, 2.38214\}$. Right: a 3x2 example where $\alpha = \{0.599648, 0.0928766, 0.288755\}$, $\beta = \{2.29566, 1.7928, 2.20249\}$, $\gamma = \{1.18328\}$, $\delta = \{3.06226, 2.74699\}$. The fact that there is a single contact point on the left boundary (whereas there are two on the right side) is a consequence of the fact that $\gamma$ is 1-periodic (constant) whereas $\delta$ is 2-periodic.</figcaption></figure>

* Borodin and Duits [BD23], for a specific value of $t$, which in fact does not play a role in the computation of the limit shape.

These are the cases depicted on the left and middle of Figure 17. An example where the periodicity condition is *not* satisfied is shown on the right of that same figure. Note that in this case, the picture does not have axial symmetries anymore.

One can extend this study when the parameters $\alpha, \beta, \gamma$, and $\delta$ are $k \times l$ period, like in Section 6.2. With the additional constraint that the corresponding edge weights are periodic, this covers in particular the case $k = 2, l = 2$ in the study of Berggren and Borodin [BB23] (without the technical assumption they have on the weights to guarantee that the vertical tentacles of the amoeba of the corresponding spectral curve are separated), but also the extension to the case where the periodicity condition for edge weights is not met, like the one depicted on Figure 18. One can relate here, like in the critical case, the number of cusps along bottom, left, top, right side of the square with the number of distinct values for $\alpha, \gamma, \beta, \delta$ respectively, whereas there are four of them along the liquid/gas interface.

## 6.4 Higher genus

The study above can be done on surfaces with higher genus, with train-track parameters assigned in a periodic way, even though it can be technically challenging to work with prime forms to control the saddle point method. Imposing periodicity for weights
from periodicity of train-track parameters implies that the fundamental domain is large enough (one needs to pick $g$ distinct integer points in the interior points of the corresponding Newton polygon), but we can relax this condition (and get just quasi-periodic weights), as in the genus 1 case above.

In [BB23], the authors give a similar looking formula for the inverse of the Kasteleyn matrix of the Aztec diamond with periodic weights. The main differences are the fact that the integrals are written with respect to the coordinates $z$ (and $w$) of the spectral curve instead of an intrinsic coordinate of the Riemann surface $\Sigma$, and the fact that the integrand is a block matrix<fnref n="1" />. Otherwise, unsurprisingly, the function $F$ appearing in the asymptotics of the integral is the same as here. We refer to the paper [BB23] for a detailed discussion of the saddle point method in the periodic case under additional technical assumptions. Working with formula (17) can be seen as a way to circumvent the technical assumptions they have to extend their results to non-generic cases, where the spectral curve is singular or the amoeba has a degenerate behaviour. In that case, the geometric correspondence is to be made between the limit shape of the Aztec diamond and the “upper part” $\Sigma^+$ of $\Sigma$, which is the connected component of $\Sigma \setminus \bigcup_{i=0}^g A_i$ whose oriented boundary is $A_0 - \sum_{i=1}^g A_i$ (instead of the amoeba of the spectral curve). Most of their geometric arguments can be carried or adapted in this slightly changed context to obtain the following:

**Proposition 27 ([BB23]).**

* The map $\psi$ from $(x, y)$ in the liquid region to the critical point of $F$ in the interior of $\Sigma^+$ is a homeomorphism.

* There is a bijective correspondence between gaseous bubbles and the real ovals $A_1, \dots, A_g$.

* There is a bijective correspondence between frozen regions and connected components of $A_0 \setminus \{\alpha_j, \beta_j, \gamma_j, \delta_j; 1 \le j \le n\}$.

The saddle point method allows one to not only compute the probability that a given edge occurs, but also to derive the asymptotics for the entry of inverse Kasteleyn matrix for vertices $b$ and $w$ at finite distance, near the same macroscopic point $(x, y)$ in the liquid region: it converges to $A_{b,w}^{u_0}$, where $u_0 = \psi(x, y)$ is the critical point of $F(\cdot; x, y)$ in $\Sigma^+$, and $A^{u_0}$ is the inverse from the family of inverses introduced in [BCdT23a] for the Fock Kasteleyn operator on any infinite minimal graph containing as a subgraph a neighborhood of the Aztec diamond containing $w$ and $b$ (for example the infinite square lattice). The inverse $A^{u_0}$ on the infinite square defines a Gibbs measure $\mathbb{P}^{u_0}$ on dimer configurations by the usual determinantal formula:

$$ \mathbb{P}^{u_0}((\mathbf{w}_1, \mathbf{b}_1), \dots, (\mathbf{w}_k, \mathbf{b}_k)) = \left( \prod_{j=1}^{k} \mathbf{K}_{\mathbf{w}_j, \mathbf{b}_j} \right) \det_{1 \le i,j \le k} (A_{\mathbf{b}_i, \mathbf{w}_j}^{u_0}) $$

<footnote n="1">¹In particular compare the term $E(u, v)$ in Equation (17) and the term $z_1 - z_2$ in their formula.</footnote>
from periodicity of train-track parameters implies that the fundamental domain is large enough (one needs to pick $g$ distinct integer points in the interior points of the corresponding Newton polygon), but we can relax this condition (and get just quasi-periodic weights), as in the genus 1 case above.

In [BB23], the authors give a similar looking formula for the inverse of the Kasteleyn matrix of the Aztec diamond with periodic weights. The main differences are the fact that the integrals are written with respect to the coordinates $z$ (and $w$) of the spectral curve instead of an intrinsic coordinate of the Riemann surface $\Sigma$, and the fact that the integrand is a block matrix<fnref n="1" />. Otherwise, unsurprisingly, the function $F$ appearing in the asymptotics of the integral is the same as here. We refer to the paper [BB23] for a detailed discussion of the saddle point method in the periodic case under additional technical assumptions. Working with formula (17) can be seen as a way to circumvent the technical assumptions they have to extend their results to non-generic cases, where the spectral curve is singular or the amoeba has a degenerate behaviour. In that case, the geometric correspondence is to be made between the limit shape of the Aztec diamond and the “upper part” $\Sigma^+$ of $\Sigma$, which is the connected component of $\Sigma \setminus \bigcup_{i=0}^g A_i$ whose oriented boundary is $A_0 - \sum_{i=1}^g A_i$ (instead of the amoeba of the spectral curve). Most of their geometric arguments can be carried or adapted in this slightly changed context to obtain the following:

**Proposition 27 ([BB23]).**

* The map $\psi$ from $(x, y)$ in the liquid region to the critical point of $F$ in the interior of $\Sigma^+$ is a homeomorphism.

* There is a bijective correspondence between gaseous bubbles and the real ovals $A_1, \dots, A_g$.

* There is a bijective correspondence between frozen regions and connected components of $A_0 \setminus \{\alpha_j, \beta_j, \gamma_j, \delta_j; 1 \le j \le n\}$.

The saddle point method allows one to not only compute the probability that a given edge occurs, but also to derive the asymptotics for the entry of inverse Kasteleyn matrix for vertices $b$ and $w$ at finite distance, near the same macroscopic point $(x, y)$ in the liquid region: it converges to $A_{b,w}^{u_0}$, where $u_0 = \psi(x, y)$ is the critical point of $F(\cdot; x, y)$ in $\Sigma^+$, and $A^{u_0}$ is the inverse from the family of inverses introduced in [BCdT23a] for the Fock Kasteleyn operator on any infinite minimal graph containing as a subgraph a neighborhood of the Aztec diamond containing $w$ and $b$ (for example the infinite square lattice). The inverse $A^{u_0}$ on the infinite square defines a Gibbs measure $\mathbb{P}^{u_0}$ on dimer configurations by the usual determinantal formula:

$$ \mathbb{P}^{u_0}((w_1, b_1), \dots, (w_k, b_k)) = \left( \prod_{j=1}^{k} K_{w_j, b_j} \right) \det_{1 \le i,j \le k} (A_{b_i, w_j}^{u_0}) $$

<footnote n="1">¹In particular compare the term $E(u, v)$ in Equation (17) and the term $z_1 - z_2$ in their formula.</footnote>
Therefore, one can obtain in this generalized setup the analogue of [BB23, Theorem 1.10]:

**Proposition 28.** The local statistics of dimers around a point $(x, y)$ of the liquid region converge to the infinite Gibbs measure $\mathbb{P}^{u_0}$ on the square lattice, with $u_0 = \psi(x, y)$.

We also note that the imaginary part of the equation used to determine the critical points of $F$, is the same as the master equation for the *tangent plane method* by Kenyon and Prause [KP22b, KP22a, KP24]: indeed, for a fixed $u$, the action $F$ is an affine function of $x$ and $y$, and the coefficients of $x$ and $y$ (respectively $\log z(u)$ and $\log w(u)$) which are analytic functions of $u$, and whose imaginary parts give the slope of the limit shape at the point $(x, y)$, see [BCdT23a, Section 4.5]. The remaining constant part is encoding the boundary conditions and should be directly related to the *intercept function*. We have thus two interpretations of the same equation: one coming from finding critical points of the action functional $F$, the other one by looking at *isothermal coordinates* in which the slope and the intercept of the asymptotic height function describing the limit shape are harmonic. Further investigation is needed to make a full connection between the two approaches. A starting point is the expression of the expression of the surface tension in terms of geometric quantities on $\Sigma$, obtained in [BCdT23a, Section 4.6], see also [BBS].

## References

[BB23] Tomas Berggren and Alexei Borodin. Geometry of the doubly periodic aztec dimer model, 2023.

[BBC$^{+}$17] Cédric Boutillier, Jérémie Bouttier, Guillaume Chapuy, Sylvie Corteel, and Sanjay Ramassamy. Dimers on rail yard graphs. *Ann. Inst. Henri Poincaré D*, 4(4):479–539, 2017.

[BBS] Alexander I. Bobenko, Nikolai Bobenko, and Yuri B. Suris. Limit shape for domino and lozenge tilings with Focks' weights. In preparation.

[BCdT22] Cédric Boutillier, David Cimasoni, and Béatrice de Tilière. Isoradial immersions. *Journal of Graph Theory*, 99(4):715–757, 2022.

[BCdT23a] Cédric Boutillier, David Cimasoni, and Béatrice de Tilière. Minimal bipartite dimers and higher genus harnack curves. *Probability and Mathematical Physics*, 4(1):151–208, March 2023.

[BCdT23b] Cédric Boutillier, David Cimasoni, and Béatrice de Tilière. Elliptic dimer models on minimal graphs and genus 1 harnack curves. *Communications in Mathematical Physics*, 400(2):1071–1136, February 2023.

[BD19] Tomas Berggren and Maurice Duits. Correlation functions for determinantal processes defined by infinite block toeplitz minors. *Advances in Mathematics*, 356:106766, 2019.
Therefore, one can obtain in this generalized setup the analogue of [BB23, Theorem 1.10]:

**Proposition 28.** The local statistics of dimers around a point $(x, y)$ of the liquid region converge to the infinite Gibbs measure $\mathbb{P}^{u_0}$ on the square lattice, with $u_0 = \psi(x, y)$.

We also note that the imaginary part of the equation used to determine the critical points of $F$, is the same as the master equation for the *tangent plane method* by Kenyon and Prause [KP22b, KP22a, KP24]: indeed, for a fixed $u$, the action $F$ is an affine function of $x$ and $y$, and the coefficients of $x$ and $y$ (respectively $\log z(u)$ and $\log w(u)$) which are analytic functions of $u$, and whose imaginary parts give the slope of the limit shape at the point $(x, y)$, see [BCdT23a, Section 4.5]. The remaining constant part is encoding the boundary conditions and should be directly related to the *intercept function*. We have thus two interpretations of the same equation: one coming from finding critical points of the action functional $F$, the other one by looking at *isothermal coordinates* in which the slope and the intercept of the asymptotic height function describing the limit shape are harmonic. Further investigation is needed to make a full connection between the two approaches. A starting point is the expression of the expression of the surface tension in terms of geometric quantities on $\Sigma$, obtained in [BCdT23a, Section 4.6], see also [BBS].

## References

[BB23] Tomas Berggren and Alexei Borodin. Geometry of the doubly periodic aztec dimer model, 2023.

[BBC$^{+}$17] Cédric Boutillier, Jérémie Bouttier, Guillaume Chapuy, Sylvie Corteel, and Sanjay Ramassamy. Dimers on rail yard graphs. *Ann. Inst. Henri Poincaré D*, 4(4):479–539, 2017.

[BBS] Alexander I. Bobenko, Nikolai Bobenko, and Yuri B. Suris. Limit shape for domino and lozenge tilings with Focks' weights. In preparation.

[BCdT22] Cédric Boutillier, David Cimasoni, and Béatrice de Tilière. Isoradial immersions. *Journal of Graph Theory*, 99(4):715–757, 2022.

[BCdT23a] Cédric Boutillier, David Cimasoni, and Béatrice de Tilière. Minimal bipartite dimers and higher genus harnack curves. *Probability and Mathematical Physics*, 4(1):151–208, March 2023.

[BCdT23b] Cédric Boutillier, David Cimasoni, and Béatrice de Tilière. Elliptic dimer models on minimal graphs and genus 1 harnack curves. *Communications in Mathematical Physics*, 400(2):1071–1136, February 2023.

[BD19] Tomas Berggren and Maurice Duits. Correlation functions for determinantal processes defined by infinite block toeplitz minors. *Advances in Mathematics*, 356:106766, 2019.
[BD23] Alexei Borodin and Maurice Duits. Biased 2× 2 periodic Aztec diamond and an elliptic curve. *Probability Theory and Related Fields*, pages 1–57, 2023.

[Ber21] Tomas Berggren. Domino tilings of the aztec diamond with doubly periodic weightings. *Ann. Probab.*, 49(4):1965–2011, 2021.

[Bet15] Dan Betea. private communication. https://icerm.brown.edu/materials/Slides/sp-s15-w3/Limit_shapes_in_the_Schur_process_]_Dan_Betea,_Universit%C3%A9_de_Paris_VI_(Pierre_et_Marie_Curie).pdf, 2015.

[CJ16] Sunil Chhita and Kurt Johansson. Domino statistics of the two-periodic Aztec diamond. *Adv. Math.*, 294:37–149, 2016.

[CJY15] Sunil Chhita, Kurt Johansson, and Benjamin Young. Asymptotic domino statistics in the Aztec diamond. *Ann. Appl. Probab.*, 25(3):1232–1278, 2015.

[CKP01] Henry Cohn, Richard Kenyon, and James Propp. A variational principle for domino tilings. *J. Amer. Math. Soc.*, 14(2):297–346 (electronic), 2001.

[CY14] Sunil Chhita and Benjamin Young. Coupling functions for domino tilings of Aztec diamonds. *Adv. Math.*, 259:173–251, 2014.

[DFSG14] Philippe Di Francesco and Rodrigo Soto-Garrido. Arctic curves of the octahedron equation. *J. Phys. A, Math. Theor.*, 47(28):34, 2014. Id/No 285204.

[DK21] Maurice Duits and Arno B. J. Kuijlaars. The two-periodic Aztec diamond and matrix valued orthogonal polynomials. *J. Eur. Math. Soc. (JEMS)*, 23(4):1075–1131, 2021.

[DLMF] NIST Digital Library of Mathematical Functions. https://dlmf.nist.gov/, Release 1.2.0 of 2024-03-15. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.

[Dob68] Roland L. Dobrušin. Description of a random field by means of conditional probabilities and conditions for its regularity. *Teor. Verojatnost. i Primenen*, 13:201–229, 1968.

[EKLP92] Noam Elkies, Greg Kuperberg, Michael Larsen, and James Propp. Alternating-sign matrices and domino tilings (part I). *Journal of Algebraic Combinatorics*, 1:111–132, 1992.

[Fay73] John D. Fay. *Theta functions on Riemann surfaces.*, volume 352. Springer, Cham, 1973.

[Foc15] Vladimir V. Fock. Inverse spectral problem for GK integrable system. *arXiv e-prints*, page arXiv:1503.00289, March 2015.
[BD23] Alexei Borodin and Maurice Duits. Biased 2× 2 periodic Aztec diamond and an elliptic curve. *Probability Theory and Related Fields*, pages 1–57, 2023.

[Ber21] Tomas Berggren. Domino tilings of the aztec diamond with doubly periodic weightings. *Ann. Probab.*, 49(4):1965–2011, 2021.

[Bet15] Dan Betea. private communication. https://icerm.brown.edu/materials/Slides/sp-s15-w3/Limit_shapes_in_the_Schur_process_]_Dan_Betea,_Universit%C3%A9_de_Paris_VI_(Pierre_et_Marie_Curie).pdf, 2015.

[CJ16] Sunil Chhita and Kurt Johansson. Domino statistics of the two-periodic Aztec diamond. *Adv. Math.*, 294:37–149, 2016.

[CJY15] Sunil Chhita, Kurt Johansson, and Benjamin Young. Asymptotic domino statistics in the Aztec diamond. *Ann. Appl. Probab.*, 25(3):1232–1278, 2015.

[CKP01] Henry Cohn, Richard Kenyon, and James Propp. A variational principle for domino tilings. *J. Amer. Math. Soc.*, 14(2):297–346 (electronic), 2001.

[CY14] Sunil Chhita and Benjamin Young. Coupling functions for domino tilings of Aztec diamonds. *Adv. Math.*, 259:173–251, 2014.

[DFSG14] Philippe Di Francesco and Rodrigo Soto-Garrido. Arctic curves of the octahedron equation. *J. Phys. A, Math. Theor.*, 47(28):34, 2014. Id/No 285204.

[DK21] Maurice Duits and Arno B. J. Kuijlaars. The two-periodic Aztec diamond and matrix valued orthogonal polynomials. *J. Eur. Math. Soc. (JEMS)*, 23(4):1075–1131, 2021.

[DLMF] NIST Digital Library of Mathematical Functions. https://dlmf.nist.gov/, Release 1.2.0 of 2024-03-15. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.

[Dob68] Roland L. Dobrušin. Description of a random field by means of conditional probabilities and conditions for its regularity. *Teor. Verojatnost. i Primenen*, 13:201–229, 1968.

[EKLP92] Noam Elkies, Greg Kuperberg, Michael Larsen, and James Propp. Alternating-sign matrices and domino tilings (part I). *Journal of Algebraic Combinatorics*, 1:111–132, 1992.

[Fay73] John D. Fay. *Theta functions on Riemann surfaces.*, volume 352. Springer, Cham, 1973.

[Foc15] Vladimir V. Fock. Inverse spectral problem for GK integrable system. *arXiv e-prints*, page arXiv:1503.00289, March 2015.
[GK13] Alexander B. Goncharov and Richard Kenyon. Dimers and cluster integrable systems. *Ann. Sci. Éc. Norm. Supér.*, 46(5):747-813, 2013.

[Hel00] Harald Helfgott. Edge effects on local statistics in lattice dimers: A study of the aztec diamond (finite case), 2000.

[Joh02] Kurt Johansson. Non-intersecting paths, random tilings and random matrices. *Probab. Theory Relat. Fields*, 123(2):225-280, 2002.

[JPS98] William Jokusch, James Propp, and Peter Shor. Random Domino Tilings and the Arctic Circle Theorem, 1998.

[Kas61] Pieter W. Kasteleyn. The statistics of dimers on a lattice: I. the number of dimer arrangements on a quadratic lattice. *Physica*, 27:1209-1225, December 1961.

[Ken97] Richard Kenyon. Local statistics of lattice dimers. *Ann. Inst. H. Poincaré Probab. Statist.*, 33(5):591-618, 1997.

[Ken02] Richard Kenyon. The Laplacian and Dirac operators on critical planar graphs. *Invent. Math.*, 150(2):409-439, 2002.

[KO06] Richard Kenyon and Andrei Okounkov. Planar dimers and Harnack curves. *Duke Math. J.*, 131(3):499-524, 2006.

[KOS06] Richard Kenyon, Andrei Okounkov, and Scott Sheffield. Dimers and amoebae. *Ann. of Math. (2)*, 163(3):1019-1056, 2006.

[KP22a] Richard Kenyon and István Prause. The genus-zero five-vertex model. *Probab. Math. Phys.*, 3(4):707-729, 2022.

[KP22b] Richard Kenyon and István Prause. Gradient variational problems in $\mathbb{R}^2$. *Duke Math. J.*, 171(14):3003-3022, 2022.

[KP24] Richard Kenyon and István Prause. Limit shapes from harmonicity: dominos and the five vertex model. *J. Phys. A, Math. Theor.*, 57(3):16, 2024. Id/No 035001.

[Kuc17] Nikolai Kuchumov. Limit shapes for the dimer model, 2017.

[Kup98] Greg Kuperberg. An exploration of the permanent-determinant method. *Electron. J. Combin.*, 5:Research Paper 46, 34, 1998.

[Law89] Derek F. Lawden. *Elliptic functions and applications*, volume 80 of *Applied Mathematical Sciences*. Springer-Verlag, New York, 1989.

[LIR69] Oscar E. Lanford III and David Ruelle. Observables at infinity and states with short range correlations in statistical mechanics. *Comm. Math. Phys.*, 13:194-215, 1969.
[GK13] Alexander B. Goncharov and Richard Kenyon. Dimers and cluster integrable systems. *Ann. Sci. Éc. Norm. Supér.*, 46(5):747-813, 2013.

[Hel00] Harald Helfgott. Edge effects on local statistics in lattice dimers: A study of the aztec diamond (finite case), 2000.

[Joh02] Kurt Johansson. Non-intersecting paths, random tilings and random matrices. *Probab. Theory Relat. Fields*, 123(2):225-280, 2002.

[JPS98] William Jokusch, James Propp, and Peter Shor. Random Domino Tilings and the Arctic Circle Theorem, 1998.

[Kas61] Pieter W. Kasteleyn. The statistics of dimers on a lattice: I. the number of dimer arrangements on a quadratic lattice. *Physica*, 27:1209-1225, December 1961.

[Ken97] Richard Kenyon. Local statistics of lattice dimers. *Ann. Inst. H. Poincaré Probab. Statist.*, 33(5):591-618, 1997.

[Ken02] Richard Kenyon. The Laplacian and Dirac operators on critical planar graphs. *Invent. Math.*, 150(2):409-439, 2002.

[KO06] Richard Kenyon and Andrei Okounkov. Planar dimers and Harnack curves. *Duke Math. J.*, 131(3):499-524, 2006.

[KOS06] Richard Kenyon, Andrei Okounkov, and Scott Sheffield. Dimers and amoebae. *Ann. of Math. (2)*, 163(3):1019-1056, 2006.

[KP22a] Richard Kenyon and István Prause. The genus-zero five-vertex model. *Probab. Math. Phys.*, 3(4):707-729, 2022.

[KP22b] Richard Kenyon and István Prause. Gradient variational problems in $\mathbb{R}^2$. *Duke Math. J.*, 171(14):3003-3022, 2022.

[KP24] Richard Kenyon and István Prause. Limit shapes from harmonicity: dominos and the five vertex model. *J. Phys. A, Math. Theor.*, 57(3):16, 2024. Id/No 035001.

[Kuc17] Nikolai Kuchumov. Limit shapes for the dimer model, 2017.

[Kup98] Greg Kuperberg. An exploration of the permanent-determinant method. *Electron. J. Combin.*, 5:Research Paper 46, 34, 1998.

[Law89] Derek F. Lawden. *Elliptic functions and applications*, volume 80 of *Applied Mathematical Sciences*. Springer-Verlag, New York, 1989.

[LIR69] Oscar E. Lanford III and David Ruelle. Observables at infinity and states with short range correlations in statistical mechanics. *Comm. Math. Phys.*, 13:194-215, 1969.
[Mum07] David Mumford. *Tata lectures on theta. II*. Modern Birkhäuser Classics. Birkhäuser Boston, Inc., Boston, MA, 2007. Jacobian theta functions and differential equations, With the collaboration of C. Musili, M. Nori, E. Previato, M. Stillman and H. Umemura, Reprint of the 1984 original.

[OR03] A. Okounkov and N. Reshetikhin. Correlation function of schur process with application to local geometry of a random 3-dimensional young diagram. *J. Amer. Math. Soc.*, 16(3):581–603, 2003.

[Per69] Jerome K. Percus. One more technique for the dimer problem. *J. Math. Phys.*, 10:1881, 1969.

[Pro97] James Propp. Talk given at the American Mathematical Society meeting in San Diego. *Slides available at http://jamespropp.org/san_diego.pdf*, 1997.

[Pro03] James Propp. Generalized domino-shuffling. *Theoret. Comput. Sci.*, 303(2-3):267–301, 2003. Tilings of the plane.

[Spe07] David E. Speyer. Perfect matchings and the octahedron recurrence. *J. Algebra. Comb.*, 25(3):309–348, 2007.

[TF61] Harold N. V. Temperley and Michael E. Fisher. Dimer problem in statistical mechanics-an exact result. *Philosophical Magazine*, 6(68):1061–1063, 1961.

[Thu90] William P. Thurston. Conway's tiling groups. *Amer. Math. Monthly*, 97(8):757–773, 1990.

[Thu17] Dylan P. Thurston. From dominoes to hexagons. In *Proceedings of the 2014 Maui and 2015 Qinhuangdao conferences in honour of Vaughan F. R. Jones' 60th birthday*, volume 46 of *Proc. Centre Math. Appl. Austral. Nat. Univ.*, pages 399–414. Austral. Nat. Univ., Canberra, 2017.

[Yan91] Bo-Yin Yang. *Two enumeration problems about the Aztec diamonds*. ProQuest LLC, Ann Arbor, MI, 1991. Thesis (Ph.D.)–Massachusetts Institute of Technology, available online at http://hdl.handle.net/1721.1/13937.
[Mum07] David Mumford. *Tata lectures on theta. II*. Modern Birkhäuser Classics. Birkhäuser Boston, Inc., Boston, MA, 2007. Jacobian theta functions and differential equations, With the collaboration of C. Musili, M. Nori, E. Previato, M. Stillman and H. Umemura, Reprint of the 1984 original.

[OR03] A. Okounkov and N. Reshetikhin. Correlation function of schur process with application to local geometry of a random 3-dimensional young diagram. *J. Amer. Math. Soc.*, 16(3):581–603, 2003.

[Per69] Jerome K. Percus. One more technique for the dimer problem. *J. Math. Phys.*, 10:1881, 1969.

[Pro97] James Propp. Talk given at the American Mathematical Society meeting in San Diego. *Slides available at http://jamespropp.org/san_diego.pdf*, 1997.

[Pro03] James Propp. Generalized domino-shuffling. *Theoret. Comput. Sci.*, 303(2-3):267–301, 2003. Tilings of the plane.

[Spe07] David E. Speyer. Perfect matchings and the octahedron recurrence. *J. Algebra. Comb.*, 25(3):309–348, 2007.

[TF61] Harold N. V. Temperley and Michael E. Fisher. Dimer problem in statistical mechanics-an exact result. *Philosophical Magazine*, 6(68):1061–1063, 1961.

[Thu90] William P. Thurston. Conway's tiling groups. *Amer. Math. Monthly*, 97(8):757–773, 1990.

[Thu17] Dylan P. Thurston. From dominoes to hexagons. In *Proceedings of the 2014 Maui and 2015 Qinhuangdao conferences in honour of Vaughan F. R. Jones' 60th birthday*, volume 46 of *Proc. Centre Math. Appl. Austral. Nat. Univ.*, pages 399–414. Austral. Nat. Univ., Canberra, 2017.

[Yan91] Bo-Yin Yang. *Two enumeration problems about the Aztec diamonds*. ProQuest LLC, Ann Arbor, MI, 1991. Thesis (Ph.D.)–Massachusetts Institute of Technology, available online at http://hdl.handle.net/1721.1/13937.