Research article

# Inhomogeneous field theory inside the arctic circle

Version 1 Released on 13 December 2015 under Creative Commons Attribution 4.0 International License

Nicolas Allegra 1,*, Jerome Dubail1, Jean-Marie Stéphan2, Jacopo Viti 3,*

## Authors' affiliations

1.  Institut Jean Lamour, Physique de la Matière et des Matériaux (IJL/P2M). CNRS : UMR7198 - Université de Lorraine
2.  Topology and Correlations in Condensed Matter Research Group, Department of Condensed Matter, Max Planck Institute for the Physics of Complex Systems
3.  Quantum Matter - Transport and Dynamics Research Group, Department of Condensed Matter, Max Planck Institute for the Physics of Complex Systems
4. *.   Unregistered author (unverified)

# Abstract

Motivated by quantum quenches in spin chains, a one-dimensional toy-model of fermionic particles evolving in imaginary-time from a domain-wall initial state is solved. The main interest of this toy-model is that it exhibits the arctic circle phenomenon, namely a spatial phase separation between a critically fluctuating region and a frozen region. Large-scale correlations inside the critical region are expressed in terms of correlators in a (euclidean) two-dimensional massless Dirac field theory. It is observed that this theory is inhomogenous: the metric is position-dependent, so it is in fact a Dirac field theory in curved space. The technique used to solve the toy-model is then extended to deal with the transfer matrices of other models: dimers on the honeycomb and square lattice, and the six-vertex model at the free fermion point ($\Delta=0$). In all cases, explicit expressions are given for the long-range correlations in the critical region, as well as for the underlying Dirac action. Although the setup developed here is heavily based on fermionic observables, the results can be translated into the language of height configurations and of the gaussian free field, via bosonization. Correlations close to the phase boundary and the generic appearance of Airy processes in all these models are also briefly revisited in the appendix.

## Introduction and overview

The arctic circle is the name given by Jockush, Propp & Shor to a phenomenon they discovered in 1995 [1] while studying dimer coverings on the so-called Aztec diamond [2]. The phenomenon consists in the appearance, in the thermodynamic limit, of a curve that separates two macroscopic domains, one in which the dimers appear to be frozen, meaning that one knows their exact position with a probability that goes quickly to one in the thermodynamic limit, and another domain where the dimers fluctuate. The same phenomenon was observed in other models [4], in particular in the six-vertex model with domain-wall boundary conditions [58], a system that already possessed a long history connected to integrability and to the evaluation of the norm of Bethe states [911]. In fact, the phenomenon itself appears to have been known as early as 1977, in a different context though: the statistics of large Young diagrams [12,13]. It was also well studied in the early eighties in the context of the physics of crystal growth, see for instance [14].

Since the late nineties, there has been an intense activity in this area, which lies at the frontier between physics, mathematics and computer science. Among the many results that were collected are the limiting distribution of dimers around the origin in the thermodynamic limit [16,17], the fluctuations of the boundary described by the Airy process [17], a general theory for dimer models [18] and the calculation of the corresponding arctic curves in connection with algebraic geometry [19], calculations of various correlation functions [2022], steps towards extension to interacting models [5,8,23] ( i.e. models that cannot be mapped to free fermions; for example: the six-vertex model at $\Delta \neq 0$, or dimers with aligning interactions [24]), and numerical investigations [2527]. There are also connections with the physics of glassy systems [28].

Our motivation for revisiting this ancient problem comes from the physics of quantum quenches in spin chains, and in particular their description by two-dimensional conformal field theory (CFT). This topic, launched by Calabrese & Cardy in 2005/2006 [29,30], received a lot of attention in recent years [3145]. The basic idea is the following: the usual quantum/classical correspondence tells us that, modulo Wick rotation, the real-time evolution of a one-dimensional quantum model is equivalent to a two-dimensional classical statistical model. In that setup, the initial state of the quantum system becomes a boundary condition in the classical model. Then the game is to chose the quantum system and the initial state such that (i) the bulk is critical, which is achieved when the Hamiltonian of the spin chain is gapless and (ii) the boundary condition flows (in the RG sense) towards a conformal boundary condition (see appendix 8. for some details). The second property is crucial and underlies most of the published works in that area. Fortunately, conformal boundary conditions are the rule rather than the exception: it is usually argued that, as soon as the correlations in the initial state are short range, the initial state corresponds to a conformal boundary condition in the thermodynamic limit (perturbed by irrelevant operators, to complete the picture). Good discussions of these ideas exist in the literature, see for instance [36].

What is the relation between the arctic circle and quantum quenches in spin chains ? Hopefully, this will become clear to the reader already in the next page. There, we introduce a quantum quench model which, we argue, may be regarded as a toy-model for arctic circle problems. The crucial ingredients are the ${\rm U}(1)$-symmetry, or conservation of the number of particles, together with the domain-wall initial state (DWIS), see below and Fig. 1. Because of particle conservation, the standard argument about the initial state renormalizing to a conformal boundary condition does not apply (see also the discussion in [46]). This is in very sharp contrast with what happens, for instance, in the quantum Ising chain, which does not possess a ${\rm U}(1)$ symmetry (the magnetization is not conserved), and where the argument can be safely exploited [47,48].

The question that motivates the rest of this introduction is: how does one apply the CFT framework of [30] to the quantum quench in the XX chain studied previously in [49]? We will soon see that the arctic circle shows up naturally in the process of answering that question. The purpose of the rest of the paper is then to revisit other arctic circle problems (dimer models and the six-vertex model) from that perspective.

Finally, let us emphasize that the model we deal with in the rest of this introduction can be mapped exactly on the PNG droplet model that appeared previously in the literature [5052]—this was pointed out to us by H. Spohn while we were completing this paper—. We shall not use this terminology here, because our own motivation really comes from quantum quenches in spin chains, so we prefer to view the model as the XX spin chain. Also, although there is some rather large overlap between this introduction and the results of [51], we think that our goal and our point of view are sufficiently different from the ones of [5052] so that they are worth explaining here independently.

### A toy-model for arctic-circle phenomena: the XX chain in imaginary-time

The model is a translation-invariant chain of free fermions living on lattice sites $x \in \mathbb{Z}+ \frac{1}{2}$. The Hamiltonian is diagonal in $k$-space, with a single band and a dispersion relation $\varepsilon(k)$, $$\label{eq:Ham} H \, = \, \int_{-\pi}^\pi \frac{dk}{2\pi} \, \varepsilon(k) \, c^\dagger(k) c(k) \, .$$ The fermion creation/annihilation modes obey the canonical anti-commutation relations $\{ c^\dagger(k) ,c(k')\} \, =\, 2\pi \, \delta(k-k')$ in $k$-space, or $\{ c^\dagger_x, c_{x'} \} = \delta_{x,x'}$ in real space. Our convention for Fourier transforms is $c^\dagger(k) \, = \, \sum_{x \in \mathbb{Z} + \frac{1}{2}} e^{i k x} c_x^\dagger$; notice that $c^\dagger(k+2\pi) = -c^\dagger (k)$. For simplicity, we specialize to the case of nearest-neighbor hopping only, such that the dispersion relation may be put in the form (possibly after a gauge transformation) \begin{equation*} \varepsilon(k) \, = \, - \cos k \,. \end{equation*} In other words, the Hamiltonian is the one of the XX chain, after a Jordan-Wigner transformation. We have normalized $\varepsilon(k)$ such that the maximal group velocity $v(k) \, = \, \frac{d}{dk} \varepsilon (k)$ is $1$ (at $k = \frac{\pi}{2}$).

Figure 1 Figure 1. The toy-model is a chain of free fermions hopping on the lattice $x \in \mathbb{Z} + \frac{1}{2}$, evolving from the DWIS, which is completely filled on the left, and completely empty on the right.
We focus on the evolution of the fermions from a DWIS, see Fig. 1. The DWIS $\mathinner{|{\Psi_0}\rangle}$ is completely filled on the left (density equal to one), and completely empty on the right (density equal to zero): $$\label{eq:definition_DWIS} \rho_x |\Psi_0\rangle \, = \, c^\dagger_x c_x \mathinner{|{\Psi_0}\rangle} \, = \, \Theta(-x) \mathinner{|{\Psi_0}\rangle} \, .$$ $\Theta$ is the Heaviside step function: $\Theta(-x) =1$ if $x<0$ and $\Theta(-x)=0$ if $x>0$. On top of this $\mathinner{|{\Psi_0}\rangle}$ is normalized: $\left< \Psi_0 \left| \Psi_0 \right> \right. = 1$. This completely determines $\mathinner{|{\Psi_0}\rangle}$, up to an irrelevant phase factor. The main originality of the model, compared to previous work on the DWIS in the XX chain, is that we focus on imaginary-time evolution. We want to compute correlators of local observables $\mathcal{O}(x,y)$, defined as the (imaginary-)time-ordered expectation values \begin{eqnarray*} \nonumber \left< \mathcal{O}_1(x_1,y_1) \dots \mathcal{O}_n(x_1,y_1) \right> & \equiv & \frac{ \mathinner{\langle{\Psi_0}|} e^{-2R \, H} \, \mathcal{T} \left[ \left( e^{(R+y_1) H} \mathcal{O}_1 (x_1) e^{-(R+y_1)H} \right) \dots \left( e^{(R+y_n) H} \mathcal{O}_n (x_n) e^{-(R+y_n)H} \right) \right]\mathinner{|{\Psi_0}\rangle}} { Z} , \end{eqnarray*} where $\mathcal{T}\left[ . \right]$ stands for time-ordering (the factors must appear from left to right with decreasing values of $y$). $Z$ is the partition function of the model, it ensures $\mathinner{\langle{1}\rangle}=1$. For our toy model it is given by the exact formula $$\label{eq:exactZ} Z \, = \, \mathinner{\langle{\Psi_0}|} e^{-2R \,H} \mathinner{|{\Psi_0}\rangle} \,=e^{R^2/2} .$$ A proof of (\ref{eq:exactZ}) is given in appendix. 11.. Aspects of the analytic continuation to real time problems for partition functions and correlations are discussed in appendix. 8.; in the rest of this introduction we stick however to the imaginary time point of view.

The constant $R>0$ is a parameter that sets the effective size of the system. Indeed, although the chain is infinite in the $x$-direction, it has a width $2R$ in the (imaginary-time) $y-$direction. We have the following cartoon in mind (see Fig. 2). One starts from the initial state $\mathinner{|{\Psi_0}\rangle}$ at imaginary time $y = -R$ and one lets the system evolve freely up to imaginary time $y = +R$, where it is conditioned to come back to the initial state $\mathinner{|{\Psi_0}\rangle}$.

Figure 2 Figure 2. (Left) Cartoon of the toy-model: the system is a strip of width $2R$, with boundary conditions corresponding to the DWIS. (Right) Density profile in the toy-model, measured numerically for $R = 128$.
At imaginary time $y = - R$ the degrees of freedom are frozen for any $x$, because all the sites are occupied on the left and empty on the right, so all correlations are trivial. Then the particles are released and can hop between left and right. So, as soon as $y >-R$, a small region around $x=0$ must appear, where the average density is between $0$ and $1$. In that region, there are density fluctuations, with non-trivial correlations. At $y = +R$, the system is conditioned to come back to its initial state, so we expect the fluctuating region to shrink back to a point as $y$ increases and approaches $+R$. In fact, it is easy to see that the model is symmetric under $y \rightarrow -y$, so the fluctuating region is reflection-invariant across the horizontal axis.

Our model possesses the same phenomenology as the arctic circle problem in its original formulation for the Aztec diamond dimer model [1]. This should not be a surprise, when looking back at (\ref{eq:exactZ}). Indeed, from general statistical mechanical arguments, the partition function is expected to scale as the exponential of the fluctuating area, which is proportional to $R^2$. The connection with dimers, and with other models, will be developed further in the main parts of the paper. In this introduction, our purpose is to use this toy-model as an illustration of the approach that will be followed later, in the main parts. As can be seen in Fig. 2, in the scaling limit $R \rightarrow \infty$, the exact shape of the fluctuating region is a disc of radius $R$ centered at the origin. Let us be more precise about what we mean by scaling limit here and in the rest of the paper. We are interested in expectation values and correlation of observables of the form $\left< \mathcal{O}_1(x_1,y_1) \dots \mathcal{O}_n(x_n,y_n) \right>$ in the scaling limit $$\label{eq:scaling_regime} x_j/R \quad {\rm and} \quad y_j/R \quad {\rm fixed \; for \; any \;} j \, , \qquad \quad{\rm and} \quad R \rightarrow \infty \, .$$ In this limit, the density inside the arctic circle takes a very simple form, interpolating from $\left< \rho \right> = 1$ on the left to $\left< \rho \right> = 0$ on the right (this formula is derived below): \begin{equation*} (x^2+y^2 < R^2)\qquad \quad \left< \rho_{x,y} \right> \, = \, \frac{1}{\pi} {\rm arccos} \frac{x}{\sqrt{R^2 - y^2}} \, . \end{equation*} The main focus of this paper, however, is neither on the shape of the critically fluctuating region, nor on the density profile inside it. Instead, what we are really eager to understand is the field theory inside the disc. What do we mean by that ? Simply the fact that the critically fluctuating region is, by definition, the region where one finds observables with long-range (connected) correlations. Such correlations must be describable by a euclidean quantum field theory with a local action; moreover this theory should be massless. In fact, there presumably is a theory that describes both the interior and the exterior of the disc, with a non-zero mass term only outside the disc. But, in the scaling limit (\ref{eq:scaling_regime}), any such mass term with $m>0$ is automatically renormalized to $m = +\infty$, so the theory outside the disc is trivial in that regime. The theory inside the disc, on the contrary, is non-trivial in the scaling limit; our goal is to find out what it is.

In fact, it is rather easy to guess what the field theory should be. Since we are dealing with a quadratic Hamiltonian for fermionic particles, we are looking for a free fermion theory. Moreover, a key aspect of the model is the particle conservation, meaning that it is invariant under ${\rm U}(1)$ transformations $c^\dagger \rightarrow e^{i \varphi} c^\dagger$, $c \rightarrow e^{-i \varphi} c$. The field theory must inherit this symmetry from the lattice model. So, there is not much choice, it has to be a massless Dirac theory in two dimensions, or possibly a direct product of several of those (meaning a theory with more than one species of massless fermions). The basic object in any such theory is the two-component Dirac spinor, $\Psi = \left( \begin{array}{c} \psi \\ \overline{\psi} \end{array} \right)$, $\Psi^\dagger = \left( \begin{array}{cc} \psi^\dagger & \overline{\psi}^\dagger \end{array} \right)$. We pick our conventions for the (euclidean) gamma matrices as $\gamma^1 \, = \, -\sigma^2$ and $\gamma^2 \, = \, \sigma^1$ (the $\sigma^j$'s are the Pauli matrices), such that $\gamma^5 \, = \, i \gamma^2 \gamma^1 \, = \, \sigma^3$ is diagonal and refer to $\psi$ and $\overline{\psi}$ as the chiral and anti-chiral components. The latter are the local fields in the continuous theory. At a given point, one must be able to express the lattice degrees of freedom in terms of the continuous fields. Such relations hold in the sense that lattice correlators are equal, in the scaling limit, to field theory correlators of carefully chosen fields at the same point. There is always some freedom in the choice of the relation between the lattice and the continuous degrees of freedom, because local perturbations may be absorbed either in this relation, or in the action, in the form of perturbations by irrelevant operators. Here, what we want to do is to fix this relation once and for all, and then work out the action of the theory. The lattice fermion $c^\dagger_{x,y}$ must be a linear combination of $\psi^\dagger$ and $\overline{\psi}^\dagger$, so let us fix $$\label{eq:lattice_continuous_XX} \left\{ \begin{array}{rcl} \displaystyle c^\dagger_{x,y} &=& \displaystyle \frac{1}{\sqrt{2\pi}}\, \psi^\dagger(x,y) \, + \, \frac{1}{\sqrt{2\pi}} \, \overline{\psi}^\dagger(x,y) \\ \displaystyle c_{x,y} &=& \displaystyle \frac{1}{\sqrt{2\pi}} \, \psi(x,y) \, + \, \frac{1}{\sqrt{2\pi}} \, \overline{\psi}(x,y) \, . \end{array} \right.$$ The normalization factors $1/\sqrt{2\pi}$ are introduced for later convenience; global phase factors may always be absorbed in the definition of the continuous fields. The point is that we chose the expression of the lattice degree of freedom in terms of the continuous fields such that it is the same everywhere, there are no position-dependent phase or scale factors. The question is now: what is the action ? The system is clearly not homogeneous—because the density is position-dependent—so there is no reason for the field theory to be homogeneous. Therefore, we are looking for some massless Dirac action in which the parameters vary with position. But what parameter is there to vary in the Dirac action, if not the mass ?

The simplest possible answer is: the metric. After all, it makes sense that the metric of a field theory describing a non-homogeneous system could vary. Also, there might be some background gauge potentials. Thus, we are led to contemplate the Dirac action in fully covariant form, allowing the possibility of both vector (v) and axial (a) background gauge potentials $$\label{eq:Dirac} S \, = \, \frac{1}{2\pi} \int d^2 {\rm x} \sqrt{ g} \, e^{\mu}_{\phantom{a}a} \left[ \overline{\Psi} \gamma^a \left( \frac{1}{2} \overset{\leftrightarrow}{\partial_\mu} \, + \, i A^{({\rm v})}_\mu \, + \, i A^{({\rm a})}_\mu \gamma^5 \right) \Psi \right] \, .$$ Here, the 'Dirac adjoint' $\overline{\Psi}$ is $\Psi^\dagger \gamma^2$, $e^\mu_{\phantom{a}a}({\rm x})$ is the tetrad, and $(d^2 {\rm x} \sqrt{ g}) = (d^2{\rm x}\, |e|^{-1})$ is the volume element. The spin connection drops out of the two-dimensional Dirac action ( cf. [53] for details). This action (\ref{eq:Dirac}) is a priori the most general possibility for the field theory describing the long-range correlations inside the critically fluctuating region (apart from the possibility of having more species of massless fermions, as already mentioned). The goal is now to properly identify the tetrads and the background gauge potentials in our toy-model. Before we proceed to this identification, let us stress the important difference between the vector and the axial potentials: the ${\rm U}(1)$-symmetry of the lattice model ( i.e. particle number conservation) is a true symmetry. These ${\rm U}(1)$ gauge transformations act on the vector part, not the axial one. Physical observables must be gauge-invariant. In fact, as we will see, both the vector and axial parts are pure gauge: $A_\mu^{({\rm v})} \, = \, \partial_\mu f$ for some function $f$ (and the same with $A_\mu^{({\rm a})}$ for a different function), so the vector part can simply be gauged away, and $A_\mu^{({\rm v})}$ can really be dropped. It can never appear in physical observables. The situation is different for the axial part, as it does not correspond to a symmetry of the model. There is no axial symmetry on the lattice; it happens to be a symmetry of the continuous Dirac Lagrangian, but it is well-known that the 'symmetry' is not preserved by any UV regularization procedure required to define the path integral: the 'symmetry' is anomalous. Stated in a more pedestrian way: the two gauge potentials $A_\mu^{({\rm v})}$ and $A_\mu^{({\rm a})}$ are not on the same footing. Physical observables must not depend on $A_\mu^{({\rm v})}$, but they can (and they do) depend on $A_\mu^{({\rm a})}$.

### Stationary phase and long-range correlations in the toy-model

The parameters in the action would be easily obtained if we knew the long-range behavior of fermionic correlators $\left< c^\dagger_{x_1,y_1} \dots c^\dagger_{x_n,y_n} c_{x'_1,y'_1} \dots c_{x'_n,y'_n} \right>$ in the strip of width $2R$. As we shall see in part 2., these correlators can be computed exactly. They turn out to be directly related to the correlators in the DWIS, through a rather simple linear transformation, \begin{eqnarray} \label{eq:correl_XX} \nonumber && \left< c^\dagger_{x_1,y_1} \dots c^\dagger_{x_n,y_n} c_{x'_1,y'_1} \dots c_{x'_n,y'_n} \right> \\ && = \, \int \prod_j \frac{dk_j}{2\pi} e^{-i k_j x_j + y_j \varepsilon(k_j) - i R \tilde{\varepsilon}(k_j)} \prod_p \frac{dk'_p}{2\pi} e^{i k'_p x'_p - y'_p \varepsilon(k'_p) + i R \tilde{\varepsilon}(k'_p)} \mathinner{\langle{\Psi_0}|} c^\dagger (k_1) \dots c^\dagger(k_n) c(k'_1) \dots c(k'_n) \mathinner{|{\Psi_0}\rangle} \, . \end{eqnarray} Here $\tilde{\varepsilon}(k)$ is the Hilbert transform of the dispersion relation $\varepsilon (k)$. For $\varepsilon(k) = -\cos k$, this is simply $\tilde{\varepsilon}(k) = - \sin k$. We will come back to the Hilbert transform at length in the main parts of the paper, so we do not insist on it in this introduction. Of course, since we are dealing with a free fermion model, Eq. (\ref{eq:correl_XX}) requires a proof only for $n=1$; the case $n > 1$ then follows immediately from applying Wick's theorem on both sides, \begin{eqnarray*} \det_{1 \leq i,j\leq n} \left( \, \left< c^\dagger_{x_i,y_i} c_{x'_j,y'_j} \right> \, \right) &=& \det_{1\leq i,j \leq n} \left( \, \int \frac{dk}{2\pi} \frac{dk'}{2\pi} e^{-i k x_i + y_i \varepsilon(k) - i R \tilde{\varepsilon}(k)} e^{i k' x'_j - y'_j \varepsilon(k') + i R \tilde{\varepsilon}(k')} \mathinner{\langle{\Psi_0}|} c^\dagger (k) c(k') \mathinner{|{\Psi_0}\rangle} \, \right) . \end{eqnarray*} For later purposes, notice that the propagator in the DWIS appearing in the r.h.s is $$\label{eq:prop_DWIS} \mathinner{\langle{\Psi_0}|} c^\dagger(k) c(k') \mathinner{|{\Psi_0}\rangle} \, = \, \sum_{x \in \mathbb{Z} + \frac{1}{2}} e^{ i (k-k') x} \Theta(-x) \, = \, \frac{1}{ 2 i \sin \left( \frac{k-k'}{2} - i 0 \right)} \, .$$ The formula (\ref{eq:correl_XX}) and the role played by $\tilde{\varepsilon}(k)$ are a key aspect of this paper, analyzed in great detail in part 2.. For now, let us take this formula for granted, and proceed to the analysis of the scaling limit of fermionic correlators. First we introduce some useful piece of notation. We write $$\label{eq:notation_doteq} \left\{ \begin{array}{rcl} c^\dagger_{x,y} & \doteq & \displaystyle \int_{-\pi}^\pi \frac{dk}{2\pi} e^{-i k x + y \varepsilon(k) - i R \tilde{\varepsilon}(k)} c^\dagger (k) \\ \\ c_{x,y} & \doteq & \displaystyle \int_{-\pi}^\pi \frac{dk}{2\pi} e^{i k x - y \varepsilon(k) + i R \tilde{\varepsilon}(k)} c (k) , \end{array} \right.$$ where we use the symbol $\doteq$ to emphasize that the relation must be understood in the sense of Eq. (\ref{eq:correl_XX}). It is not an equality between operators. Instead, it holds inside correlators, the l.h.s being inserted inside the bracket $\left< .\right>$, and the r.h.s between $\mathinner{\langle{\Psi_0}|}$ and $\mathinner{|{\Psi_0}\rangle}$. Eq. (\ref{eq:notation_doteq}) is really nothing more than the equality (\ref{eq:correl_XX}), but it is written in a more compact way.

Figure 3 Figure 3. The critically fluctuating region $x^2+y^2 < R^2$ is mapped to an infinite strip of width $\pi$ by $z(x,y) \, = \, {\rm arccos}\frac{x}{\sqrt{R^2-y^2}} - i\, {\rm arcth} \frac{y}{R}$.
Next, we observe that the scaling limit (\ref{eq:scaling_regime}) of the correlators may be obtained with the steepest descent method. In the relations (\ref{eq:notation_doteq}), the contour of integration over $k$ may be deformed such that it passes through the points where the phase is stationary; these are the solutions of \begin{equation*} \frac{x}{R} \, + \, i \frac{y}{R} \, \frac{d\varepsilon(k)}{dk} \, + \, \frac{d \tilde{\varepsilon}(k)}{dk} \, = \, \frac{1}{R} \frac{d}{dk} \left[ k x + i y \varepsilon (k) + R \tilde{\varepsilon}(k) \right] \, = \, 0 \, . \end{equation*} For $\varepsilon(k) = -\cos k$ and $\tilde{\varepsilon}(k) = - \sin k$, this is a quadratic equation in $e^{i k}$, which is easily solved. As long as $x^2+y^2 < R^2$, there are two simple roots: $k \, = \, z(x,y)$ and $k \, =\, -z^*(x,y)$, where $z^*$ is the complex conjugate of $z$, with \begin{eqnarray*} z\, = \, z(x,y) \, \equiv \, {\rm arccos} \frac{x}{\sqrt{R^2-y^2}} \, - \, i \, {\rm arcth}\frac{y}{R} \, . \end{eqnarray*} The integral over $k$ is then evaluated by the stationary phase approximation around these points: we approximate the argument of the exponential by \begin{equation*} \left\{ \begin{array}{rcll} e^{-i k x \, + \, y \varepsilon(k) \, - \, i R \tilde{\varepsilon}(k)} & \simeq & e^{ - i \varphi \, - \, \frac{i}{2}(k-z)^2\, e^{\sigma}} & \qquad {\rm around} \quad k = z \\ e^{-i k x \, + \, y \varepsilon(k) \, - \, i R \tilde{\varepsilon}(k)} & \simeq & e^{ i \varphi^* \, + \, \frac{i}{2}(k+z^*)^2\, e^{\sigma}} & \qquad {\rm around} \quad k = -z^* \, , \end{array} \right. \end{equation*} where we have defined \begin{eqnarray*} \varphi \, = \, \varphi(x,y) & \equiv & z(x,y)\, x \, + \, i y \, \varepsilon(z(x,y)) \, + \, R \, \tilde{\varepsilon}(z(x,y)) \, = \,z(x,y)\, x \, - \, \sqrt{R^2-x^2-y^2} \, ,\\ e^{\sigma} \, = \, e^{\sigma(x,y)} & \equiv & \, i y \, \frac{d^2 \varepsilon}{dk^2} (z(x,y)) \, + \, R \, \frac{d^2 \tilde{\varepsilon}}{dk^2}(z(x,y)) \, = \, \sqrt{R^2 - x^2 -y^2} \, . \end{eqnarray*} Finally, one performs the integral in (\ref{eq:notation_doteq}), which is now gaussian; it yields the asymptotic formulas $$\label{eq:notation_doteq_asympt} \left\{ \begin{array}{rcl} c^\dagger_{x,y} & \underset{ R \rightarrow \infty} {\doteq} & \displaystyle \frac{ e^{-i \varphi(x,y) - i \frac{\pi}{4}}} {\sqrt{2\pi}} e^{-\frac{1}{2} \sigma(x,y)} c^\dagger (z) \, + \, \frac{ e^{i \varphi^*(x,y) + i \frac{\pi}{4}}} {\sqrt{2\pi}} e^{-\frac{1}{2} \sigma(x,y)} c^\dagger (-z^*) \\ \\ c_{x,y} & \underset{ R \rightarrow \infty} {\doteq} & \displaystyle \frac{ e^{i \varphi(x,y) + i \frac{\pi}{4}}} {\sqrt{2\pi}} e^{-\frac{1}{2} \sigma(x,y)} c (z) \, + \, \frac{ e^{-i \varphi^*(x,y) - i \frac{\pi}{4}}} {\sqrt{2\pi}} e^{-\frac{1}{2} \sigma(x,y)} c (-z^*) \, . \end{array} \right.$$ What exactly do we mean by $c^\dagger(z)$, $c(z)$ here, when $z \notin \mathbb{R}$? Again, these relations must be understood in the sense of Eqs. (\ref{eq:correl_XX})-(\ref{eq:notation_doteq}): the correlator of fermions in the l.h.s of (\ref{eq:notation_doteq_asympt}) are equal to expressions that involve $\mathinner{\langle{\Psi_0}|} c^\dagger(z) c(z')\mathinner{|{\Psi_0}\rangle} = 1/(2i \sin \frac{z-z'}{2} )$, the analytic continuation of (\ref{eq:prop_DWIS}).

### Identifying the parameters in the Dirac action

Now let us come back to Eqs. (\ref{eq:lattice_continuous_XX})-(\ref{eq:Dirac}) and to the discussion around them. The resemblance with Eq. (\ref{eq:notation_doteq_asympt}) is striking, and clearly suggests that the continuous field $\psi^\dagger$ is, roughly speaking, nothing but $e^{- i \varphi - i \frac{\pi}{4}} e^{-\frac{1}{2} \sigma} c^\dagger$ (and similar relations for $\psi$, $\overline{\psi}^\dagger$ and $\overline{\psi}$). It is then very easy to extract the parameters of the Dirac action in the toy-model. One simple way of doing this is to match the local behaviors of the lattice and continuous propagators. The behavior of the lattice propagator follows from Eqs. (\ref{eq:notation_doteq_asympt}), which give its scaling limit in the form of a sum of four terms. When $(x',y')$ approaches $(x,y)$, it is dominated by only two of the terms, that diverge as $$\label{eq:prop_short} \left< c_{x,y}^\dagger c_{x',y'} \right> \, \simeq \, \frac{e^{-\frac{\sigma + \sigma'} {2}}} {2\pi i} \left( \frac{e^{-i (\varphi - \varphi' )}} {z-z'} - \frac{e^{i (\varphi^* - \varphi'^* )}} {z^*-z'^*} \right) \, .$$ Looking back at Eq. (\ref{eq:lattice_continuous_XX}), one sees that the propagator of the continuous fields must have the local behavior \begin{equation*} \left< \psi^\dagger(x,y) \psi(x',y') \right> \, \simeq \, e^{-\frac{\sigma + \sigma'} {2}} \frac{e^{-i (\varphi - \varphi' )}} {i(z-z')} \qquad \qquad \qquad \left< \overline{\psi}^\dagger(x,y) \overline{\psi}(x',y') \right> \, \sim \, e^{-\frac{\sigma + \sigma'} {2}} \frac{e^{i (\varphi^* - \varphi'^* )}} {-i(z^*-z'^*)} \, , \end{equation*} implying that they are Green's functions for the following operators \begin{equation*} \left\{ \begin{array}{cll} e^\sigma \left( i \partial_{\bar{z}} + \frac{i}{2} (\partial_{\bar{z}} \sigma) - (\partial_{\bar{z}} \varphi) \right) \big\langle \psi^\dagger(z,\bar{z}) \psi(z',\bar{z}') \,\big\rangle & = & \pi \, \delta^{(2)} (z-z') \\ e^\sigma \left( -i \partial_z - \frac{i}{2} (\partial_z \sigma) - (\partial_z \varphi^*) \right) \big\langle \,\overline{\psi}^\dagger(z,\bar{z}) \overline{\psi}(z,\bar{z}') \,\big\rangle & = & \pi \, \delta^{(2)} (z-z') \, . \end{array} \right. \end{equation*} We view the latter as arising from the action \begin{equation*} \begin{array}{c} \displaystyle S \, = \, \frac{1}{\pi} \int d^2 z \, e^{\sigma} \left[ \left( \begin{array}{cc} \psi^\dagger & \overline{\psi}^\dagger \end{array} \right) \left( \begin{array}{cc} - \frac{i}{2} \overset{\leftrightarrow}{\partial}_{\bar{z}} - (\partial_{\bar{z}} \varphi) & 0 \\ 0 & \frac{i}{2} \overset{\leftrightarrow}{\partial}_{z} - (\partial_{z} \varphi^*) \end{array}\right) \left( \begin{array}{c} \psi \\ \overline{\psi} \end{array} \right) \right] \\ \\ \left( = \, \frac{1}{\pi} \int d^2 z \, e^{\sigma} \, \left[ \left( [ i \partial_{\bar{z}} + \frac{i}{2} (\partial_{\bar{z}} \sigma) - (\partial_{\bar{z}} \varphi )] \psi^\dagger \right) \psi + \left( [ -i \partial_{z} - \frac{i}{2} (\partial_{z} \sigma) - (\partial_{z} \varphi^*) ] \overline{\psi}^\dagger \right) \overline{\psi} \right] \quad + \quad {\rm boundary \; terms} \right) \, . \end{array} \end{equation*} This, as anticipated in the discussion above, is a particular case of the generic Dirac action (\ref{eq:Dirac}). To match the two expressions, one needs to write $(\ref{eq:Dirac})$ in the coordinate system provided by the map $(x,y) \mapsto z(x,y)$, \begin{equation*} \left\{ \begin{array}{rcl} {\rm x}^1 \, + \, i \, {\rm x}^2 & = & z(x,y) \\ {\rm x}^1 \, - \, i \, {\rm x}^2 & = & z^*(x,y) \, , \end{array} \right. \end{equation*} where $0 <{\rm x}^1 < \pi$ and ${\rm x}^2 \in \mathbb{R}$ (Fig. 3). The tetrad in (\ref{eq:Dirac}) must be taken as the diagonal matrix $$e^\mu_{\phantom{b}a} \, = \, e^{-\sigma} \delta^\mu_{\phantom{b}a} \, ,$$ so the corresponding metric is $d s^2 = e^{2 \sigma} [(d {\rm x}^1)^2 + (d {\rm x}^2)^2]$. In other words, what the map $(x,y) \mapsto z(x,y)$ does for us is that it provides isothermal coordinates, for free. The background gauge potentials are readily identified as $$A^{({\rm v})}_\mu \, = \, -i \partial_\mu {\rm Im} \,\varphi \qquad \qquad \qquad A^{({\rm a})}_\mu \, = \, -\partial_\mu {\rm Re}\, \varphi \, .$$ They are both pure gauge. It may look problematic that $A^{({\rm v})}_\mu$ is imaginary, but again, it can never appear in physical observables, which must be gauge invariant. [Notice that the propagator itself, as defined so far, is not a physical observable, as it is not gauge invariant. To make it gauge invariant, one would have to define it with the insertion of a Wilson line.] The axial potential, on the other hand, can (and does) appear in physical quantities. One example is the density: it may be obtained from Eq. (\ref{eq:prop_short}) by taking $y'=y$ first, and then $x'\rightarrow x$: $\left< \rho_{x,y} \right> \, = \, -\frac{1}{\pi} \frac{\partial_x {\rm Re} \varphi} {e^{\sigma} \partial_x z}$. The formula for the density can be further simplified using $e^{\sigma} \partial_x z=-1$. This identity can be checked directly using the formulas for $z$ and $e^\sigma$ or by observing that the stationary phase equation is of the form $F(x,y,z)=0$, where $F(x,y,z)=x+iyv(z)+ R\tilde v(z)$ and $v(z)=\frac{d\varepsilon}{dz}$, $\tilde{v}(z)=\frac{d\tilde{\varepsilon}}{dz}$. We also have $e^{\sigma}=\partial_z F$. By the chain rule and the fact that $dF=0$ when $z=z(x,y)$ we obtain $\partial_x z(x,y)=-\partial_x F/\partial_z F$ and the result follows from $\partial_x F=1$. Coming back to the density we obtain \begin{eqnarray*} \left< \rho_{x,y} \right> \, = \,-\frac{1}{\pi} A^{({\rm a})}_{x} \, , \end{eqnarray*} and plugging the explicit formula for $A^{({\rm a})}_{x}$, one recovers $\left< \rho_{x,y} \right> = \frac{1}{\pi} {\rm arccos} \frac{x}{\sqrt{R^2-y^2}}$ as claimed above.

### Discussion

Let us pause and take a look back at the logic of the previous paragraphs. We defined a two-dimensional toy-model that displays an arctic circle; it is discrete in one direction (space coordinate $x$) and continuous in the other direction (imaginary-time coordinate $y$). It was argued, on general grounds (fermion statistics, ${\rm U}(1)$ symmetry, locality), that the long-range correlations inside the fluctuating region should be described by a massless Dirac theory. Then:

• independently of these field theory considerations, we exhibited an exact formula for the lattice fermion propagator $\left< c_{x,y}^\dagger c_{x',y'} \right>$, to be proved later (part 2.)
• the scaling limit of this propagator was analyzed with the stationary phase approximation. We found that the stationary points are the solution of a quadratic equation. If one solution is called $z$, the other is $-z^*$. The mapping $(x,y) \mapsto z$ is a diffeomorphism that sends the interior of the arctic circle to a certain subset of $\mathbb{C}$ (here, an infinite strip), so it can be used as a coordinate system. Moreover, in the analysis of the stationary points, two other functions $\sigma$ and $\varphi$ show up naturally
• trying to match the local asymptotic behavior of the lattice fermion propagator with the one coming from the generic Dirac action (\ref{eq:Dirac}), we were able to fix the free parameters in this action, namely the tetrad and the background gauge potentials. It turns out that the coordinate system provided by $(x,y) \mapsto z$ is isothermal, with a conformal factor given by $e^\sigma$, and that the gauge potentials are given by the real- and imaginary-part of $\partial_\mu \varphi$

Some important aspects were neglected in the present discussion, for instance the boundary conditions that should be imposed on the continuous fields (we will focus on the boundary in part 3.4.). But, once again, the point we wish to emphasize in this introduction is that, if one is willing to describe the fluctuating region by a field theory, then the theory must be inhomogeneous, in the sense that, once a relation between the lattice and the continuous degrees of freedom is fixed (such as Eq. (\ref{eq:lattice_continuous_XX}) above), then the parameters in the action vary with position. Since the components of the metric themselves are such parameters, it is natural that one ends up with a theory in curved euclidean space. [In the isothermal coordinate system above, $\sigma \, = \, \log ( R \,\sin ({\rm x}^1) / \cosh ({\rm x}^2) )$, which gives the non-zero scalar curvature $S \, = \, - e^{2 \sigma} \Delta \sigma \, = \, \left( \sin^2 ({\rm x}^1) + \cosh^2({\rm x}^2) \right) / (R^2 \sin^4 ({\rm x}^1))$ , so it is really a curved space.]

The reader might be puzzled by our logic: usually, when one discusses field theory in relation with lattice models, the goal is to use arguments or calculation tools that are simpler in the continuous setting to make interesting predictions about lattice observables. Here, even though we initially promised a field theory investigation, technically, we have proceeded backwards: we first provided the solution of the lattice model with formula (\ref{eq:correl_XX}), and then imposed some field theory interpretation on top of it. It is our opinion that the observation that 'arctic circle problems' are related to field theories in curved space deserves to be highlighted, and that the precise derivation of the action deserves a proper exposition. We also believe that some of the ideas exposed here will be useful to tackle other problems involving inhomogeneous systems, including inhomogeneous quantum quenches, particles trapped by inhomogeneous potentials as well as closely related random matrix questions; this will be discussed elsewhere.

### Organization of the paper

The goal of the rest of the paper is to extend the approach we just presented to other free-fermion models that have appeared previously in the arctic circle context: dimer models on the honeycomb and square lattice, and the six-vertex model at $\Delta = 0$. We also need to complete the derivation of the key formula (\ref{eq:correl_XX}), which was left aside in this introduction. This is achieved in part 2.: we obtain the propagator in $k$-space for the toy-model, in the form of a closed, simple, exact formula. We also extend the calculation to the case of two-band models, needed later to solve the six-vertex model and dimers on the square lattice. The most important aspect of part 2. is the appearance of the Hilbert transform of the dispersion relation; it is instrumental throughout the paper. Part 3. contains the results about the dimer model on the honeycomb lattice. Because this is a one-band model, the honeycomb model is the simplest possible extension of the toy-model; the results are quite similar to the ones we sketched in this introduction. In parts 4. and 5., we again repeat the analysis, for the dimer model on the square lattice, and for the six-vertex model. Essentially, these are nothing but the toy-model on steroids. Intermediate calculations become increasingly cumbersome, while all main conclusions remain qualitatively identical: the large-scale correlations inside the fluctuating region are those of a Dirac theory, and the parameters in the action (\ref{eq:Dirac}) are identified from a pair of stationary points. In part 6. we translate our results into the language of discrete height configurations and their continuum limit, the gaussian free field. While the approach taken in this paper is really tied to free fermions, we think that presenting the results in the alternative (bosonic) way should be useful to facilitate connections with other references, for example [18].

Some additional information is gathered in five appendices. In appendix. 8., we explain how the results of the toy-model may be analytically continued to real time. Appendix. 9. deals with the behavior of the correlations close to the arctic curve, which were neglected in the main text. We provide a simple heuristic explanation for the genereric appearance of the Airy kernel as well as the Tracy-Widom distribution in that case. In appendix. 10. we compute the various Hilbert transforms needed throughout the paper, while all the technical steps needed to obtain exact finite-size formulae for the two-band models are gathered appendix. 12.. Finally, it is explained in appendix. 11. how the well known partition functions for all the models studied here may be recovered from bosonization arguments.

## Propagator in imaginary-time and appearance of the Hilbert transform

In this part, we prove the formula (\ref{eq:correl_XX}) of the introduction. More precisely, we prove the $n=1$ case: we compute the propagator inside the strip. The result involves the Hilbert transform of the dispersion relation $\varepsilon (k)$, which comes from the DWIS. Indeed, the DWIS acts as a projector onto the negative real axis in real-space, which, when Fourier transformed, gives the Hilbert transform. We derive this key formula in full details in the case of a single band, and then sketch the derivation for the two-band case, which differs from the former case only by relatively minor details. Not surprisingly, the main technical tool in this section is Wick's theorem; we use it in a relatively non-standard bosonized form though, which considerably simplifies the calculations. This bosonization trick will not be needed in the rest of the paper; only the results are, and the latter will be recalled when needed. Therefore, the impatient reader may want to skip this technical part and go directly to section 3..

### Imaginary time propagator for one-band models

This subsection is devoted to the calculation of the propagator in $k$-space for the one-band Hamiltonian (\ref{eq:Ham}). Let us recall that this propagator is defined as \begin{equation*} \left< c^\dagger(k,y) c(k',y') \right> \, \equiv \, \left\{ \begin{array}{c} \displaystyle \frac{\mathinner{\langle{\Psi_0}|} e^{ -(R-y) H} c^\dagger(k) e^{-(y -y') H} c(k') e^{- (R+ y') H} \mathinner{|{\Psi_0}\rangle}}{\mathinner{\langle{\Psi_0}|} e^{-2 R\, H} \mathinner{|{\Psi_0}\rangle}} \qquad {\rm if} \;(y>y') \\ \\ \displaystyle - \,\frac{\mathinner{\langle{\Psi_0}|} e^{ -(R-y') \, H} c(k) e^{-(y' -y) H} c^\dagger(k') e^{- (R+ y) \, H} \mathinner{|{\Psi_0}\rangle}}{\mathinner{\langle{\Psi_0}|} e^{- 2 R \, H} \mathinner{|{\Psi_0}\rangle}} \qquad {\rm if} \;(y<y') \,.\end{array} \right. \end{equation*} We are going to show that $$\label{eq:ktr} \left< c^\dagger(k,y) c(k',y') \right> \, = \, \frac{e^{y \varepsilon(k) -i R \tilde{\varepsilon}(k)} e^{-y' \varepsilon(k') + i R \tilde{\varepsilon}(k')}} { 2i \sin \left( \frac{k-k'}{2} - i 0 \right)} \, ,$$ where $\tilde{\varepsilon}(k)$ is the Hilbert transform of the dispersion relation $\varepsilon(k)$: $$\label{eq:Hilbert} \tilde{\varepsilon}(k) \, = \, {\rm p.v.} \int_{-\pi}^\pi \frac{dk'}{2\pi} \varepsilon(k') ~ {\rm cot} \left( \frac{k-k'}{2}\right).$$ The formula (\ref{eq:ktr}) is quite remarkable, and it completely solves the toy-model, as we explained in the introduction. It also solves the dimer model on the honeycomb lattice, as we explain in part 3.. To prove this formula, we proceed as follows. First, we note that, since \begin{equation*} \left\{ \begin{array}{ccccc} c^\dagger(k,y) & = & e^{y H} c^\dagger (k) e^{-y H} & = & e^{ y \varepsilon(k)} c^\dagger (k)\\ c(k',y') & = & e^{y' H} c (k') e^{-y' H} & = & e^{- y' \varepsilon(k)} c (k')\, , \end{array} \right. \end{equation*} the $y$- and $y'$-dependence of (\ref{eq:ktr}) is trivial. Instead, what we really need to show is $$\label{eq:ktr0} \left< c^\dagger(k,0) c(k',0) \right> \, = \, \frac{e^{-i R \left[ \tilde{\varepsilon}(k) - \tilde{\varepsilon}(k') \right]}} { 2i \sin \left( \frac{k-k'}{2} - i 0 \right)} \, .$$ We introduce the normal order $:.:$ adapted to our initial state $\mathinner{|{\Psi_0}\rangle}$, namely \begin{equation*} : c^\dagger_x c_x : \;\equiv \; \left\{ \begin{array}{rcl} c^\dagger_x c_x && {\rm if} \quad x>0 \\ - c_x c_x^\dagger && {\rm if} \quad x<0 \, , \end{array} \right. \end{equation*} such that \begin{eqnarray*} \nonumber :c^\dagger (k) c(k'): & = & c^\dagger(k) c(k') \, - \, \mathinner{\langle{\Psi_0}|} c^\dagger (k) c(k') \mathinner{|{\Psi_0}\rangle} \\ & = & c^\dagger(k) c(k') \, - \, \frac{1}{ 2i \sin \left( \frac{k-k'}{2} -i 0\right)} \, . \end{eqnarray*} It is probably possible to get to formula (\ref{eq:ktr0}) using Wick's theorem for fermions directly, but we haven't been able to do so; the reader may try and convince themselves that the computation is rather cumbersome. So, instead, we chose to use standard bosonization formulas that make the various computational steps much lighter. For this purpose, we introduce a (real) free bosonic field $\varphi(k)$, and we bosonize the fermion creation/annihilation operators, as well as the Hamiltonian, according to the rules $$\label{eq:bosonization} ({\rm bosonization}) \qquad \left\{ \begin{array}{rcl} c^\dagger(k) & \rightarrow & : e^{i \varphi(k)}: \\ c(k) & \rightarrow & : e^{-i \varphi(k)}: \\ \int \frac{dk}{2\pi} \varepsilon(k) \, c^\dagger(k) c(k) & \rightarrow & \int \frac{dk}{2\pi} \varepsilon(k) \, \partial \varphi(k) \, . \end{array} \right.$$ The propagator of the bosonic field is chosen as $$\label{eq:bosprop} \mathinner{\langle{\Psi_0}|} \varphi(k - i \epsilon ) \varphi(k') \mathinner{|{\Psi_0}\rangle} \, = \, - \log \left[ 2 i \sin \left( \frac{k- i \epsilon -k'}{2} \right) \right]$$ such that $\mathinner{\langle{\Psi_0}|} : e^{i \varphi(k- i \epsilon)} :\, :e^{-i \varphi(k')} : \mathinner{|{\Psi_0}\rangle} \, = \, \mathinner{\langle{\Psi_0}|} c^\dagger(k) c(k') \mathinner{|{\Psi_0}\rangle}$ when $\epsilon \rightarrow 0^+$. When computing correlators or commutators, one needs to carefully regularize the bosonic expressions. This is done by adding some small imaginary part to the argument $k$, which is larger, or smaller, depending on the order of appearance of the operators in the different expressions. An instructive example consists in recovering the commutator $\left[ H , c^\dagger(k) \right] \, = \, \varepsilon(k) \, c^\dagger(k)$ in bosonized form: \begin{eqnarray} \nonumber \left[ H , c^\dagger(k) \right] & \rightarrow & \lim_{\epsilon \rightarrow 0^+} \int \frac{dq}{2\pi} \varepsilon(q) \left[ \partial \varphi(q- i \epsilon ) \,: e^{i \varphi(k)} :\; -\; :e^{i \varphi(k)}: \, \partial \varphi(q+ i \epsilon ) \right] \\ \nonumber & = & \lim_{\epsilon \rightarrow 0^+} \int \frac{dq}{2\pi} \varepsilon(q) \left[ i \partial_{q} \mathinner{\langle{\Psi_0}|} \varphi(q- i \epsilon ) \varphi(k) \mathinner{|{\Psi_0}\rangle} \; -\; i\partial_{q} \mathinner{\langle{\Psi_0}|} \varphi(k ) \varphi(q+i \epsilon) \mathinner{|{\Psi_0}\rangle} \, \right] \, : e^{i \varphi(k)} : \\ \nonumber &=& \lim_{\epsilon \rightarrow 0^+} \int \frac{dq}{2\pi} \varepsilon(q) \left[ -\frac{i}{2}\, {\rm cot} \left(\frac{q - i \epsilon - k}{2} \right) \; +\; \frac{i}{2}\, {\rm cot} \left(\frac{q + i \epsilon - k}{2} \right) \, \right] \, : e^{i \varphi(k)} : \\ \nonumber &=& \int \frac{dq}{2\pi} \varepsilon(q) \, 2\pi\delta(q-k) \, : e^{i \varphi(k)} : \;\, = \; \varepsilon(k) \, :e^{i \varphi(k)}: \, , \end{eqnarray} which is the expected result. In the first line, we used the bosonization formulas (\ref{eq:bosonization}), introducing the regulator $\varepsilon > 0$ according to the respective position of the two operators; in the second line we applied Wick's theorem; in the third line we used the bosonic propagator (\ref{eq:bosprop}); in the last line we used $$\label{eq:delta_pv} \lim_{\epsilon \rightarrow 0^+} {\rm cot} \left(\frac{k - k' \pm i \epsilon} {2} \right) \; = \; \mp \, 2\pi i\, \delta(k-k') \, + \, {\rm p.v.} \left[ {\rm cot} \left(\frac{k - k'}{2} \right) \right] \, .$$ With these bosonization tricks at hand, it is a relatively pleasant exercise to prove formula (\ref{eq:ktr0}). Indeed, the quantity that we need to evaluate becomes \begin{eqnarray} \label{eq:interm0} \nonumber && \frac{\mathinner{\langle{\Psi_0}|} e^{-R H} c^\dagger(k) c(k') e^{-R H} \mathinner{|{\Psi_0}\rangle}}{ \mathinner{\langle{\Psi_0}|} e^{-2 R H} \mathinner{|{\Psi_0}\rangle}} \\ && \rightarrow \; \frac{ \mathinner{\langle{\Psi_0}|}: e^{ -R \int \frac{dq}{2\pi} \varepsilon(q) \, \partial \varphi(q - 2i \epsilon)} : \;:e^{i \varphi(k- i \epsilon)}: \; :e^{-i \varphi(k')}: \; : e^{ -R \int \frac{dq'}{2\pi} \varepsilon(q') \,\partial \varphi(q'+ i \epsilon)} :\mathinner{|{\Psi_0}\rangle} } { \mathinner{\langle{\Psi_0}|} :e^{ - R \int \frac{dq}{2\pi} \varepsilon(q) \, \partial \varphi(q)} : \; :e^{ - R \int \frac{dq'}{2\pi} \varepsilon(q') \, \partial \varphi(q')} : \mathinner{|{\Psi_0}\rangle}} \; , \end{eqnarray} and we may compute this by applying the celebrated (exact) fusion formula for vertex operators \begin{equation*} : e^{\alpha \varphi(k)} : \, : e^{ \beta \varphi(k')}: \;\, = \; e^{\alpha \beta \mathinner{\langle{\Psi_0}|} \varphi(k) \varphi(k') \mathinner{|{\Psi_0}\rangle}} \, :e^{ \alpha \varphi(k) + \beta \varphi(k')} : \end{equation*} which follows from Wick's theorem. In (\ref{eq:interm0}), we start by fusing the two central vertex operators that come from $c^\dagger(k)$ and $c(k')$. This gives \begin{equation*} \mathinner{\langle{\Psi_0}|} c^\dagger(k) c(k') \mathinner{|{\Psi_0}\rangle} \, \times \, \frac{ \mathinner{\langle{\Psi_0}|}: e^{ -R \int \frac{dq}{2\pi} \varepsilon(q) \, \partial \varphi(q - 2i \epsilon)} : \;:e^{i \varphi(k- i \epsilon) -i \varphi(k')}: \; : e^{ -R \int \frac{dq'}{2\pi} \varepsilon(q') \,\partial \varphi(q'+ i \epsilon)} :\mathinner{|{\Psi_0}\rangle} } { \mathinner{\langle{\Psi_0}|} :e^{ - R \int \frac{dq}{2\pi} \varepsilon(q) \, \partial \varphi(q)} : \; :e^{ - R \int \frac{dq'}{2\pi} \varepsilon(q') \,\partial \varphi(q')} : \mathinner{|{\Psi_0}\rangle}} . \end{equation*} Next, we fuse the vertex operator involving the integral over $q$ with the one involving the integral over $q'$, both in the numerator and in the denominator. Clearly, the scalar factor that comes out of the fusion is the same in both cases, so it cancels when we take the ratio, and we do not need to evaluate it. The denominator is then the expectation value of a single normal-ordered exponential, which is one by definition. Thus, we are left with \begin{equation*} \mathinner{\langle{\Psi_0}|} c^\dagger(k) c(k') \mathinner{|{\Psi_0}\rangle} \, \times \, \mathinner{\langle{\Psi_0}|}: e^{ -R \int \frac{dq}{2\pi} \varepsilon(q) \, \partial \varphi(q - 2i \epsilon) \, -\, R \int \frac{dq'}{2\pi} \varepsilon(q') \,\partial \varphi(q'+ i \epsilon)} : \;:e^{i \varphi(k- i \epsilon) -i \varphi(k')}: \mathinner{|{\Psi_0}\rangle} . \end{equation*} We use the fusion formula one last time to get \begin{eqnarray*} \nonumber \mathinner{\langle{\Psi_0}|} c^\dagger(k) c(k') \mathinner{|{\Psi_0}\rangle} \, \times \, e^{ -R \int \frac{dq}{2\pi} \varepsilon(q) \, i \partial_{q} \mathinner{\langle{\Psi_0}|} \varphi(q - 2i \epsilon) \varphi(k- i \epsilon) \mathinner{|{\Psi_0}\rangle} \, -\,R \int \frac{dq'}{2\pi} \varepsilon(q') \,i \partial_{q'} \mathinner{\langle{\Psi_0}|} \varphi(q'+ i \epsilon) \varphi(k- i \epsilon) \mathinner{|{\Psi_0}\rangle}} \\ \times \; e^{ R \int \frac{dq}{2\pi} \varepsilon(q) \, i \partial_{q} \mathinner{\langle{\Psi_0}|} \varphi(q - 2i \epsilon) \varphi(k') \mathinner{|{\Psi_0}\rangle} \, +\,R \int \frac{dq'}{2\pi} \varepsilon(q') \,i \partial_{q'} \mathinner{\langle{\Psi_0}|} \varphi(q'+ i \epsilon) \varphi(k') \mathinner{|{\Psi_0}\rangle}} \, . \end{eqnarray*} Finally, we plug in the explicit form of the bosonic propagator (\ref{eq:bosprop}), \begin{eqnarray*} \nonumber \mathinner{\langle{\Psi_0}|} c^\dagger(k) c(k') \mathinner{|{\Psi_0}\rangle} \, \times \, e^{ i R \int \frac{dq}{2\pi} \varepsilon(q) \left[ \frac{1}{2} {\rm cot} \left( \frac{q - k - i \epsilon}{2} \right) \, +\, \frac{1}{2}{\rm cot} \left( \frac{q - k + 2 i \epsilon} {2} \right) \right]} \; e^{ - i R \int \frac{dq}{2\pi} \varepsilon(q) \left[ \frac{1}{2} {\rm cot} \left( \frac{q - k' - 2 i \epsilon}{2} \right) \, +\, \frac{1}{2} {\rm cot} \left( \frac{q - k' + i \epsilon} {2} \right) \right]} \, , \end{eqnarray*} and we take the limit $\epsilon \rightarrow 0^+$, using (\ref{eq:delta_pv}). The integrals in the exponentials become the principal value that appears in the definition of the Hilbert transform (\ref{eq:Hilbert}). We arrive at \begin{equation*} \frac{\mathinner{\langle{\Psi_0}|} e^{-R H} c^\dagger(k) c(k') e^{-R H} \mathinner{|{\Psi_0}\rangle}}{ \mathinner{\langle{\Psi_0}|} e^{-2 R H} \mathinner{|{\Psi_0}\rangle}} \; = \; e^{- i R\left[ \tilde{\varepsilon}(k) - \tilde{\varepsilon}(k') \right]} \; \mathinner{\langle{\Psi_0}|} c^\dagger(k) c(k') \mathinner{|{\Psi_0}\rangle}, \end{equation*} as claimed above.

### Imaginary time propagator for two-band models

In parts 4. and 5., we will see that, in our setup, the dimer model on the square lattice and the six-vertex model are both models with two fermionic degrees of freedom per unit cell. Therefore, their transfer matrices $T$—as well as their 'Hamiltonians' $H \, =\, -\log T$—have two bands, instead of one. So, in order to prepare for these models, we need to generalize the previous calculation to the two-band case. We have not presented the six-vertex and square-lattice dimer models yet, and actually, for the purposes of this section, we do not need them. Instead, we introduce only the few features that are strictly necessary for now; the details that are more specific to either model will be given later, in parts 4. and 5..

First, we are still dealing with fermionic particles on the lattice $\mathbb{Z} + \frac{1}{2}$, and with the DWIS $\mathinner{|{\Psi_0}\rangle}$. The transfer matrix $T$, however, is not invariant under translations of one site, but only of two sites. Therefore, in $k$-space, the transfer matrix can be put in the form (possibly up to an unimportant global numerical factor) \begin{equation*} T \, = \, \exp \left[ - \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \frac{dk}{2\pi} \left( \begin{array}{cc} c^\dagger(k) & c^\dagger(k+\pi) \end{array} \right) h(k) \left( \begin{array}{c} c(k) \\ c(k+\pi) \end{array} \right) \right] \, , \end{equation*} where $h(k)$ is some $2 \times 2$ matrix whose elements depend on the details of the model. Now let us make the following assumptions. They are all satisfied by the models tackled later; these assumptions allow to express the results in a more compact way, which will be most useful in parts 4. and 5.. We assume that:

1. the transfer matrix $T$ is a hermitian operator; equivalently $h(k)$ is a hermitian $2 \times 2$ matrix
2. the two (real) eigenvalues of $h(k)$ are opposite to each other, namely there is a unitary $2\times 2$ matrix $U(k)$ and a real-valued function $\varepsilon (k)$ such that \begin{equation*} h(k) \, = \, U^t(k) \left( \begin{array}{cc} \varepsilon(k) & 0 \\ 0 & - \varepsilon (k) \end{array} \right) U^*(k) \end{equation*}
3. the function $\varepsilon (k)$, which is initially defined on $[ - \frac{\pi}{2} , \frac{\pi}{2} ]$, may be analytically continued to the real axis $\mathbb{R}$. Moreover, this function $\varepsilon : \mathbb{R} \rightarrow \mathbb{R}$ is anti-periodic with period $\pi$: $\varepsilon(k + \pi ) = - \varepsilon (k)$. Similarly, the matrix-valued function $U(k)$ can be chosen such that it can be analytically continued to $\mathbb{R}$, with the property $$\label{eq:periodU} U(k + \pi) \, = \, \left( \begin{array}{cc} 0 & 1 \\ -1 & 0 \end{array} \right) U(k) \left( \begin{array}{cc} 0 & -1 \\ 1 & 0 \end{array} \right) \, .$$

The second and third assumptions are not really crucial to the expansion of the formalism of the previous section to two-band models, but they hold for the two models of interest in this paper, and they will considerably simplify the calculations in parts 4. and 5.. Therefore, we chose to rely quite heavily on these properties already in the present section, in order to set up the proper notations for later parts; but let us stress that, if we were interested in more general two-band models in which these two properties do not hold, it would still be possible to generalize the calculation of the previous section. In fact, these additional assumptions reflect the very special structure of both the dimer model on the square lattice and the six-vertex model, compared to more generic two-band models. These properties come from the fact these two models are somehow one-band models in disguise. This is rather clear for the six-vertex model, because it is well-known that, when formulated in the standard language of the row-to-row transfer matrix—as opposed to the diagonal-to-diagonal transfer matrix which we use in this paper, see part 5.—the six-vertex model at $\Delta = 0$ is a genuine one-band model. Our setup, however, based on the DWIS, makes it more natural for us to work with the diagonal-to-diagonal transfer matrix, and therefore, with a two-band model. Presumably, using the row-to-row transfer matrix instead would lead us to a framework that should resemble the one of Okounkov and Reshetikhin [21]. But let us not pursue this direction here, and come back to our purpose, which is to get the imaginary-time propagator of a model that satisfies the three assumptions above.

The canonical transformation \begin{equation*} \left( \begin{array}{c} c^\dagger_+(k) \\ c^\dagger_-(k) \end{array} \right) \, = \, U(k) \left( \begin{array}{c} c^\dagger(k) \\ c^\dagger (k+\pi) \end{array} \right) \end{equation*} diagonalizes the transfer matrix, \begin{equation*} T \, = \, \exp \left[ - \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \frac{dk}{2\pi} \varepsilon(k) \left( c^\dagger_+(k) c_+(k) - c^\dagger_-(k) c_-(k) \right) \right] \, . \end{equation*} Moreover, thanks to property (\ref{eq:periodU}), we have \begin{equation*} \left\{ \begin{array}{rcl} c^\dagger_+(k+\pi) & = & c^\dagger_-(k) \\ c^\dagger_-(k+\pi) &=& -c^\dagger_+(k) \, , \end{array} \right. \end{equation*} and this, together with the anti-periodicity of $\varepsilon(k)$, allows us to rewrite the transfer matrix as \begin{equation*} T \, = \, \exp \left[ - \int_{-\pi}^{\pi} \frac{dk}{2\pi} \varepsilon(k) c^\dagger_+(k) c_+(k) \right] \, . \end{equation*} It is then obvious how to adapt the bosonization trick learned above. Essentially, we just need to replace the modes $c^\dagger(k)$ and $c(k)$ by $c^\dagger_+(k)$ and $c_+(k)$. We introduce a free real boson $\varphi_+$, with the propagator \begin{equation*} \left< \varphi_+ (k-i \varepsilon ) \varphi_+ (k') \right> \, = \, \log \left[ \mathinner{\langle{\Psi_0}|} c^\dagger_+(k) c_+(k') \mathinner{|{\Psi_0}\rangle} \right] \, . \end{equation*} The r.h.s needs to be computed. This, in fact, is the only important difference with the previous section: instead of being simply $1/(2 i \sin ( \frac{k-k'}{2} - i0 ) )$, the correlator in the r.h.s depends on the change of basis $U(k)$. More precisely, we have \begin{equation*} \left( \begin{array}{cc} \mathinner{\langle{\Psi_0}|} c_+^\dagger(k) c_+(k') \mathinner{|{\Psi_0}\rangle} & \mathinner{\langle{\Psi_0}|} c_+^\dagger(k) c_-(k') \mathinner{|{\Psi_0}\rangle} \\ \mathinner{\langle{\Psi_0}|} c_-^\dagger(k) c_+(k') \mathinner{|{\Psi_0}\rangle} & \mathinner{\langle{\Psi_0}|} c_-^\dagger(k) c_-(k') \mathinner{|{\Psi_0}\rangle} \end{array} \right) \, = \, U(k) \left( \begin{array}{cc} \mathinner{\langle{\Psi_0}|} c^\dagger(k) c(k') \mathinner{|{\Psi_0}\rangle} & \mathinner{\langle{\Psi_0}|} c^\dagger(k) c(k'+\pi) \mathinner{|{\Psi_0}\rangle} \\ \mathinner{\langle{\Psi_0}|} c^\dagger(k+\pi) c(k') \mathinner{|{\Psi_0}\rangle} & \mathinner{\langle{\Psi_0}|} c^\dagger(k+\pi) c(k'+\pi) \mathinner{|{\Psi_0}\rangle} \end{array} \right) U^\dagger(k') \, , \end{equation*} where the matrix in the middle in the r.h.s is $\left( \begin{array}{cc} 1/( 2 i \sin \frac{k-k'} {2} ) & -1 /( 2 i \cos \frac{k-k'} {2} ) \\ 1/( 2 i \cos \frac{k-k'} {2} ) & 1 /(2 i \sin \frac{k-k'} {2} ) \end{array} \right)$. So we find that the propagator in the DWIS is given by the $(1,1)$ matrix element $$\label{eq:fU} \mathinner{\langle{\Psi_0}|} c^\dagger_+(k) c_+(k') \mathinner{|{\Psi_0}\rangle} \, = \, \left[ U(k) \left( \begin{array}{cc} \frac{1}{ 2 i \sin \left( \frac{k-k'} {2} \right)} & \frac{-1}{ 2 i \cos \left( \frac{k-k'} {2} \right)} \\ \frac{1}{ 2 i \cos \left( \frac{k-k'} {2} \right)} & \frac{1}{ 2 i \sin \left( \frac{k-k'} {2} \right)} \end{array} \right) U^\dagger(k') \right]_{11} \, .$$ The rest of the calculation is straightforward. We apply the bosonization rules \begin{equation*} \left\{ \begin{array}{ccc} c^\dagger_+ (k) & \rightarrow & e^{ i \varphi_+(k)} \\ c_+ (k) & \rightarrow & e^{ -i \varphi_+(k)} \\ T & \rightarrow & \exp \left( - \int_{-\pi}^\pi \frac{dk}{2\pi} \varepsilon (k) \, \partial \varphi_+(k) \right) \, , \end{array} \right. \end{equation*} and we compute the imaginary-time propagator as in the previous section; the computational steps are exactly the same, so we do not repeat them here. In the end, we find that the propagator at $y=y'=0$ is: \begin{eqnarray*} \nonumber \left< c^\dagger_+ (k,0) c_{+} (k',0) \right> & \equiv & \frac{ \mathinner{\langle{\Psi_0}|} T^N c^\dagger_+ (k) c_{+}(k') T^N \mathinner{|{\Psi_0}\rangle}} {\mathinner{\langle{\Psi_0}|} T^{2N} \mathinner{|{\Psi_0}\rangle}} \\ \nonumber &=& e^{- i N \int \frac{dq}{2\pi} \varepsilon (q) \, \left[ \partial_{q} \log \mathinner{\langle{\Psi_0}|} c^\dagger_+(k-i \epsilon) c_+(q-2i \epsilon) \mathinner{|{\Psi_0}\rangle} + \partial_{q} \log \mathinner{\langle{\Psi_0}|} c^\dagger_+(k-i \epsilon) c_+(q+i \epsilon) \mathinner{|{\Psi_0}\rangle} \right]} \\ \nonumber && \quad \times \; e^{i N \int \frac{dq}{2\pi} \varepsilon (q) \, \left[ \partial_{q} \mathinner{\langle{\Psi_0}|} c^\dagger_+(q-2 i \epsilon) c_+(k') \mathinner{|{\Psi_0}\rangle} + \partial_{q} \log \mathinner{\langle{\Psi_0}|} c^\dagger_+(q+ i \epsilon) c_+(k') \mathinner{|{\Psi_0}\rangle} \right]} \\ && \quad \quad \times \; \mathinner{\langle{\Psi_0}|} c^\dagger_+(k) c_+(k') \mathinner{|{\Psi_0}\rangle} \, . \end{eqnarray*} Finally, since the matrix elements of $U(k)$ are analytic, the only possible singularities in $\mathinner{\langle{\Psi_0}|} c^\dagger_+(k) c_+(k') \mathinner{|{\Psi_0}\rangle}$ come from the entries $1/(2i \sin ( \frac{k-k'}{2} ))$ in (\ref{eq:fU}); therefore the kernel $\partial_{q} \log \mathinner{\langle{\Psi_0}|} c^\dagger_+(q) c_+(k) \mathinner{|{\Psi_0}\rangle}$ must have a single poles at $k = k' + 2\pi \mathbb{Z}$. Hence, the expressions in the two exponentials are again principal values. Bringing back the (trivial) dependence on $y$ and $y'$, we end up with the formula $$\label{eq:prop_interm_2} \left< c^\dagger_+ (k,y) c_{+} (k',y') \right> \, =\, e^{ y \varepsilon(k) - i N \tilde{\varepsilon}_{U} (k)} \, e^{ -y' \varepsilon(k') + i N \tilde{\varepsilon}_{U} (k')} \, \mathinner{\langle{\Psi_0}|} c^\dagger_+(k) c_+(k') \mathinner{|{\Psi_0}\rangle} \, ,$$ where the function $\tilde{\varepsilon}_{U}(k)$ is defined as $$\tilde{\varepsilon}_{U}(k) \, \equiv \, {\rm p.v.} \int_{-\pi}^\pi \frac{dq}{2\pi} \varepsilon (q) \; 2 \partial_{q} \log \mathinner{\langle{\Psi_0}|} c^\dagger_+(q) c_+(k) \mathinner{|{\Psi_0}\rangle} \, .$$ It does not seem possible to go further without being more specific about the matrix elements of $U(k)$. In the concrete models of parts 4. and 5., we will compute these matrix elements, and the associated kernel $2 \partial_{q}\log \mathinner{\langle{\Psi_0}|} c^\dagger_+(q) c_+(k) \mathinner{|{\Psi_0}\rangle}$. It will turn out that, both for the dimer model on the square lattice and for the six-vertex model, this kernel takes a particularly simple form. We will see that, in both cases, $\tilde{\varepsilon}_{U}(k)$ is the Hilbert transform of $\varepsilon (k)$, modulo a conjugation by a change of variables. Let us mention that a quite similar bosonization procedure can be used to compute the partition function of these models (see appendix 11. for details).

## Dimers on the honeycomb lattice

In the introduction, we focused on a simple toy model, the XX chain in imaginary time. A particularly neat feature of that model is the simple dispersion relation $\varepsilon(k)=-\cos k$, which implies that the Hilbert transform is also very simple: $\tilde{\varepsilon}(k)=-\sin k$. The correlations are then analyzed using the steepest descent method. Crucial in the whole analysis are of course the stationary points, which satisfy the equation $x+i y v(k)+R\tilde{v}(k)=0$, where $v(k)=\frac{d\varepsilon}{dk}$ is the group velocity, and $\tilde{v}(k)=\frac{d\tilde{\varepsilon}}{dk}$ is its Hilbert transform. Solving this equation is elementary, as it is quadratic in $e^{ik}$. In particular the arctic curve is the set $(x,y)$ where the two solutions are degenerate; it turns out to be the circle $x^2+y^2=R^2$. Of course, this circle can be turned into an ellipse by a trivial rescaling of the dispersion relation.

There is, in fact, a whole (one-parameter) family of dispersion relations that lead to a quadratic equation for the stationary phase. It is obtained by considering a group velocity of the form $$\label{eq:vk} v(k)=\frac{i/2}{1+u e^{ik}}-\frac{i/2}{1+u e^{-ik}}=\frac{u \sin k}{1+u^2+2u \cos k},$$ which is a sum of two simple poles in $e^{ik}$ and $e^{-ik}$. Here $u$ is some real parameter in the interval $(0,1)$; the XX case is recovered by considering $v(k)/u$ and by taking the limit $u\to 0$, or, equivalently, by rescaling the vertical coordinates as $y\to y/u$, $R\to R/u$ and then by taking the limit $u\to 0$. It turns out that there are a few models that realize the more general dispersion relation corresponding to (\ref{eq:vk}), and that are therefore essentially the same model, as they can be mapped onto each other. The first representative is the multilayer PNG droplet [51]. Another representative, on which we focus in this part, is the honeycomb dimer model on a strip, with boundary conditions at the top and bottom corresponding to the DWIS. In this model, the parameter $u$ is the fugacity of the dimers on a subset of the lattice edges (Fig. 4).

Figure 4 Figure 4. Mapping onto a system of non-intersecting lattice paths (particle trajectories). (i): Reference configuration with staggered dimers (shown in red). (ii): A configuration of dimers (shown in thick blue) is superimposed on the reference configuration. This generates four lattice paths which propagate from bottom to top. (iii): Rules for the trajectories. A particle is defined as a vertical link not occupied by a real (thick blue) dimer. Note the distinction between even and odd rows.
The dimer model can be mapped onto a system of free fermions as follows. We consider the ordered configuration of staggered dimers that is displayed in Fig. 4 (i); we refer to it as the reference configuration. Then, any other dimer configuration may be superimposed on the reference configuration, as shown in Fig. 4 (ii). This procedure generates a set of lattice paths, all going upwards. It is convenient to think of these as trajectories of particles, where a particle is defined as a vertical link not occupied by a real dimer. The particle trajectories have the following important properties:

• There is a one to one correspondence between the particle trajectories and the dimer coverings.
• The number of particles is conserved from one row to the next.
• The trajectories do not cross.

The last property implies that the statistics of the particles does not matter. We are thus free to consider them as fermions, and look at the transfer matrix as an operator acting on the fermion Fock space. The rules for the fermion trajectories are summarized in Fig. 4 (iii). Each fermion can go either to the nearest left or to the nearest right site from one row to the next. However, the two jumps are not symmetric: for even rows we assign a weight $u$ to a right jump, and $1$ to the left jump. For odd rows the rule is reversed, a weight $u$ is assigned to the left jump and a weight $1$ to the right jump.

Before we turn to the construction of the transfer matrix, let us explain the geometry induced by the DWIS boundary conditions. In dimer language, the DWIS $\mathinner{|{\Psi_0}\rangle}=\prod_{x<0} c_x^\dagger \mathinner{|{0}\rangle}$ corresponds to all vertical links on the right occupied by dimers and all links on the left empty. As a result, many other dimer occupancies are automatically set to one, see Fig. 5(i). An example of corresponding fermions trajectories is pictured in Fig. 5 (ii).

Figure 5 Figure 5. (i) Geometry obtained by enforcing the DWIS as bottom and top boundary conditions. The corresponding dimers cross either of the dotted lines, and are shown in green. Other dimer occupancies are also set to one as a result; these are also shown in green. (ii) Example of a possible set of particle trajectories, from bottom to top.

### The transfer matrix

The transfer matrix acts on the fermion Fock space that represents the configurations of particles along a given row, namely on linear combinations of states of the form $\prod_{i} c_{x_i}^\dagger \mathinner{|{0}\rangle}$. It is however more elegant to look at the action of the transfer matrix on a single fermion creation operator. Given the rules shown in Fig. 4 (iii), a transfer matrix satisfying the rules \begin{eqnarray*} T_{\rm e} c_x^\dagger T_{\rm e}^{-1}&=& c_x^\dagger+u c_{x+1}^\dagger,\\ T_{\rm o} c_x^\dagger T_{\rm o}^{-1}&=&u c_{x-1}^\dagger + c_{x}^\dagger, \end{eqnarray*} generates all possible dimer configurations, with the correct weights. The action on any many-fermion configuration may be reconstructed using $T\prod_{i} c_{x_i}^\dagger \mathinner{|{0}\rangle} =\left(\prod_{i} T c_{x_i}^\dagger T^{-1} \right)T\mathinner{|{0}\rangle}$ and the fact that $T\mathinner{|{0}\rangle}=\mathinner{|{0}\rangle}$ for either $T=T_{\rm e}$ or $T=T_{\rm o}$. In particular, the partition function is simply $Z=\mathinner{\langle{\Psi_0|\left(T_{\rm o} T_{\rm e}\right)^{N}|\Psi_0}\rangle}$. Notice that $T_{\rm e}$ and $T_{\rm o}$ are not Hermitian; however $T^2=T_{\rm o} T_{\rm e}$ is. Therefore, $T^2$ is the natural hermitian transfer matrix we focus on, and to which one can readily apply the formalism of Part 2.. The action of $T_{\rm e}$, $T_{\rm o}$, $T^2$ on the Fourier modes of the fermions (recall that they are defined as $c^\dagger(k)=\sum_{x\in \mathbb{Z}+1/2} e^{ikx} c_x^\dagger$) is diagonal: \begin{eqnarray*} T_{\rm e} \,c^\dagger(k)\,T_{\rm e}^{-1}&=&\left(1+u e^{ik}\right) c^\dagger(k),\\ T_{\rm o} \,c^\dagger(k)\,T_{\rm o}^{-1}&=& \left(1+ue^{-ik}\right) c^\dagger(k),\\ T^2 \,c^\dagger(k)\, T^{-2}&=&\left(1+u^2 +2u \cos k\right) c^\dagger(k). \end{eqnarray*} More explicitely, $T^2$ reads \begin{equation*} T^2=\exp\left(-2 \int \frac{dk}{2\pi} \varepsilon(k) c^\dagger(k)c(k)\right) \;, \qquad \quad {\rm with}\qquad \varepsilon(k)=-\frac{1}{2}\log \left(1+u^2+2u \cos k\right). \end{equation*} Notice that the derivative of the dispersion relation $\varepsilon(k)$ is exactly Eq. (\ref{eq:vk}).

### Fermionic correlators: exact finite-size formulas

In the rest of this part, we assume for simplicity that $N$ is even. The vertical axis is labelled by integers $-N\leq y\leq N$, while the horizontal axis is still labelled by half integers $x\in \mathbb{Z}+1/2$. The expectation value of an observable $\mathcal{O}$ at position $x$ and $y$ is \begin{eqnarray} \label{eq:hex_eveny} \mathinner{\langle{\mathcal{O}_{x,y}}\rangle}=\left\{ \begin{array}{ccccc} \displaystyle{\frac{ \mathinner{\langle{\Psi_0}|} (T^2)^{\frac{N-y}{2}} \mathcal{O}_x (T^2)^{\frac{N+y}{2}} \mathinner{|{\Psi_0}\rangle}} { \mathinner{\langle{\Psi_0}|} T^{2N} \mathinner{|{\Psi_0}\rangle} } } &&,&& y \,\textrm{ even}\\ \\ \displaystyle{\frac{ \mathinner{\langle{\Psi_0}|} (T^2)^{\frac{N-y-1}{2}} T_{\rm o}\mathcal{O}_x T_{\rm o}^{-1} (T^2)^{\frac{N+y+1}{2}} \mathinner{|{\Psi_0}\rangle}} { \mathinner{\langle{\Psi_0}|} T^{2N} \mathinner{|{\Psi_0}\rangle}}} &&,&& y \textrm{ odd} \, . \end{array} \right. \end{eqnarray} Similar formulae hold for more than one local observable. In particular, the fermion propagator is, using the result of part 2., $$\label{eq:prop_honey} \left< c_{x,y}^\dagger c_{x',y'} \right> \, = \, \int \frac{dk}{2\pi} e^{- i k x + y \varepsilon(k) - i N \tilde{\varepsilon}(k)} r_y(k) \int \frac{dk'}{2\pi} e^{ i k' x' - y' \varepsilon(k) + i N \tilde{\varepsilon}(k')} s_{y'}(k') \, \mathinner{\langle{\Psi_0}|} c^\dagger(k) c(k') \mathinner{|{\Psi_0}\rangle} \, ,$$ where the propagator in the DWIS is still $\mathinner{\langle{\Psi_0}|} c^\dagger(k) c(k') \mathinner{|{\Psi_0}\rangle} = 1/(2 i \sin \frac{k-k'}{2})$. The factors $r_y(k)$ and $s_y(k)$ are \begin{equation*} r_y(k) = s_y(k) = 1 \qquad {\rm if} \quad y \quad {\rm is \;even} \, , \qquad \qquad \qquad \quad r_y(k) \,=\, \frac{1}{s_y(k)} \,=\, \left( \frac{1+ u e^{- ik}}{ 1+ u e^{i k}} \right)^{\frac{1}{2}} \quad {\rm if} \quad y \quad {\rm is \;odd} \, . \end{equation*} These coefficients follow from the exponentials in (\ref{eq:hex_eveny}): for $y$ odd, one uses the relations $T_{\rm o} c^\dagger (k)T_{\rm o}^{-1}=(1+u e^{-ik})c^\dagger (k)$ and $T_{\rm o} c(k)T_{\rm o}^{-1} =(1+u e^{-ik})^{-1}c(k)$, and one takes into account the small shift $y\rightarrow y+1$ in (\ref{eq:hex_eveny}) compared to the formulae in part 2.. This gives $r_y(k) = 1/s_y(k) = e^{\varepsilon(k)} ( 1+u e^{-i k})$, as claimed. In what follows, it is convenient to use again the notation from the introduction, see Eq. (\ref{eq:notation_doteq}). We rewrite (\ref{eq:prop_honey}) as $$\label{eq:ident_hex} \left\{ \begin{array}{rcl} c^\dagger_{x,y} & \doteq & \displaystyle \int_{-\pi}^\pi \frac{dk}{2\pi} e^{-i k x + y \varepsilon (\kappa) - i N \tilde{\varepsilon}(\kappa)} \, r_{y} (k) \, c^\dagger(k) \\ \\ c_{x,y} & \doteq & \displaystyle \int_{-\pi}^\pi \frac{dk}{2\pi} e^{i k x-y \varepsilon (\kappa) + i N \tilde{\varepsilon}(\kappa)} \, s_{y} (k) \, c(k) \, , \end{array} \right.$$ where the dot means that this is an equality that holds only when the l.h.s is evaluated within the bracket $\left< . \right>$, and the r.h.s is the expectation value in the DWIS $\mathinner{\langle{\Psi_0}|} . \mathinner{|{\Psi_0}\rangle}$. Finally, we note that, here, the Hilbert transform of the dispersion relation is (see the Appendix 10.) \begin{equation*} \varepsilon(k)=-\frac{1}{2}\log \left(1+u^2+2u \cos k\right) \qquad\Rightarrow\qquad \tilde{\varepsilon}(k)=\frac{i}{2}\log \left(\frac{1+u e^{ik}}{1+u e^{-ik}}\right) \, . \end{equation*}

Figure 6 Figure 6. The critically fluctuating region $X^2+y^2 < N^2$ is mapped to an infinite strip of width $\pi$ by $( x,y) \mapsto z(x,y)$, see Eq. (\ref{eq:zhoney}). The red curve is $x = \frac{N-|y|}{2}$; it is tangent to the ellipse at the two red points. Those two points are mapped to ${\rm x}^2 = \pm \infty$ on the strip. The two points $(x,y) = (-\frac{u^2}{1-u^2}N,\pm N)$ in orange are mapped to $({\rm x}^1,{\rm x}^2) = (\pi, \mp {\rm arcth}\left( \frac{1-u^2}{1+u^2}\right))$. The strip itself can be mapped conformally to the upper half-plane, $\zeta = e^{i z}$. Because this last mapping is conformal, it makes no difference to use the coordinate $z$ or $\zeta$ as an isothermal coordinate system. Here we chose to use $z$.

### Fermionic correlators: scaling limit

Exactly like in the introduction, we apply the steepest descent method to obtain the large scale correlations in our honeycomb dimer model. The stationary points in (\ref{eq:ident_hex}) are the solutions to ($v(k) = \frac{d}{dk} \varepsilon(k)$ and $\tilde{v}(k) = \frac{d}{dk} \tilde{\varepsilon}(k)$): \begin{equation*} x+iy v(k)+N\tilde{v}(k)=0 \, . \end{equation*} This equation possesses two roots, $k=z(x,y)$ and $k=-z^*(x,y)$, with $$\label{eq:zhoney} z(x,y)=\arccos \left(\frac{(1+u^2)x-N u^2}{u\sqrt{(N-2x)^2-y^2}}\right)-i\, {\rm arctanh}\left( \frac{y}{N-2x}\right) \, ,$$ or equivalently $\zeta (x,y) \, = \, e^{i z(x,y)} \, = \, \frac{ \frac{1+u^2}{u} x - N u + i\sqrt{R^2 - X^2 - y^2}} {N-2 x- y}$, where $X=\frac{1-u^2}{u}x+N u$. The two solutions are usually distinct; they are degenerate only if \begin{equation*} X^2+y^2 \, = \, N^2 \, . \end{equation*} This curve, an ellipse, is the arctic curve of the model. Inside this ellipse the long-range correlations are critical. Around the points where the phase is stationary, the argument in the exponential in (\ref{eq:ident_hex}) is approximated by \begin{equation*} \left\{ \begin{array}{ccl} e^{-ikx+y\varepsilon(k)-i N \tilde{\varepsilon}(k)} & \simeq & e^{-i \varphi-\frac{i}{2}(k-z)^2 e^{\sigma+i\theta}} \qquad \,\textrm{around}\quad k=z \\ e^{-ikx+y\varepsilon(k)-i N\tilde{\varepsilon}(k)}& \simeq & e^{i \varphi^*+\frac{i}{2}(k+z^*)^2 e^{\sigma-i\theta}}\qquad \textrm{around}\quad k=-z^* \, , \end{array} \right. \end{equation*} where $$\label{eq:phi_honey} \left\{ \begin{array}{rcl} \varphi & \equiv &\displaystyle z(x,y)x+iy \varepsilon(z(x,y))+N\tilde{\varepsilon}(z(x,y)) \\ \\ e^{\sigma +i \theta} &\equiv & \displaystyle i y\, \frac{d^2\varepsilon}{dk^2}(z(x,y)) \, + \, N\,\frac{d^2\tilde{\varepsilon}}{dk^2}(z(x,y)) \\ &=&\frac{u\sqrt{(N-2x)^2-y^2}\sqrt{N^2-X^2-y^2}}{(1-u^2)\sqrt{N^2-y^2}}\times\exp\left[-i\arccos\left(\frac{N(N-2x)-\frac{1+u^2}{1-u^2}y^2}{\sqrt{N^2-y^2}\sqrt{(N-2x)^2-y^2}}\right)\right] \, . \end{array} \right.$$ Performing the gaussian integration, we find that the scaling limit of (\ref{eq:ident_hex}) is $$\label{eq:scaling_ident_hex} \left\{\begin{array}{rcl} c^\dagger_{x,y} &\doteq& \displaystyle \frac{e^{-i \frac{\pi}{4}}}{\sqrt{2\pi}} e^{- \frac{\sigma + i \theta}{2}} e^{- i \varphi} r_y(z) c^\dagger(z) \, + \, \frac{e^{i \frac{\pi}{4}}}{\sqrt{2\pi}} e^{- \frac{\sigma - i \theta}{2}} e^{i \varphi^*} r_y(-z^*) c^\dagger(-z^*) \\ \\ c_{x,y} &\doteq& \displaystyle \frac{e^{i \frac{\pi}{4}}}{\sqrt{2\pi}} e^{- \frac{\sigma + i \theta}{2}} e^{ i \varphi} s_y(z) c(z) \, + \, \frac{e^{-i \frac{\pi}{4}}}{\sqrt{2\pi}} e^{- \frac{\sigma - i \theta}{2}} e^{-i \varphi^*} s_y(-z^*) c(-z^*) \, . \end{array} \right.$$

### Identification of the Dirac theory inside the arctic ellipse

Again, we proceed by repeating the procedure of the introduction. There are two important differences though. The first is in the relation between the lattice degrees of freedom and the continuous fields. Here we chose $$\begin{array}{lcr} \left\{ \begin{array}{rcl} c^\dagger_{x,y} &=& \displaystyle \frac{1}{\sqrt{2\pi}} \left[ \psi^\dagger(x,y) \, + \, \overline{\psi}^\dagger(x,y) \right] \\ c_{x,y} &=& \displaystyle \frac{1}{\sqrt{2\pi}} \left[ \psi(x,y) \, + \, \overline{\psi}(x,y) \right] \end{array}\right. &\qquad& (y \quad{\rm even}) \\ \\ \left\{ \begin{array}{rcl} c^\dagger_{x,y} &=& \displaystyle \frac{1}{\sqrt{2\pi}} \left[ \psi^\dagger(x,y+1) + u \, \psi^\dagger(x+1,y+1) \, + \, \overline{\psi}^\dagger(x,y+1) + u \, \overline{\psi}^\dagger(x+1,y+1)\right] \\ c_{x,y} &=& \displaystyle \frac{1}{\sqrt{2\pi}} \left[ \psi(x,y-1)+ u \, \psi(x+1,y-1) \, + \, \overline{\psi}(x,y-1) + u \, \overline{\psi}(x+1,y-1) \right] \end{array}\right. &\qquad& (y \quad{\rm odd}) \, . \end{array}$$ We make this choice so as to naturally incorporate the coefficients $r_y(k)$ and $s_y(k)$ above; the field $\psi^\dagger$ is then essentially $e^{-i\frac{\pi}{4}} e^{-\frac{\sigma + i\theta}{2}} e^{-i \varphi} c^\dagger(z)$, and similar relations for $\overline{\psi}^\dagger$, $\psi$ and $\overline{\psi}$, see Eq. (\ref{eq:scaling_ident_hex}). It follows that the short-distance behavior of the propagator of the continuous fields is \begin{equation*} \left< \psi^\dagger(x,y) \psi(x',y') \right> \, \sim \, e^{-\frac{\sigma + i \theta}{2} -\frac{\sigma' + i \theta'}{2}} \frac{e^{-i (\varphi - \varphi')}}{i (z-z')} \quad , \qquad \qquad \left< \overline{\psi}^\dagger(x,y) \overline{\psi}(x',y') \right> \, \sim \, e^{-\frac{\sigma - i \theta}{2} -\frac{\sigma' - i \theta'}{2}} \frac{e^{i (\varphi^* - \varphi'^*)}}{i (z^*-z'^*)} \, . \end{equation*} Here we see the second main difference with the case treated in the introduction: the appearance of the angle $\theta = \theta(x,y)$, which was absent from the propagator in the XX chain, see Eq. (\ref{eq:prop_short}). This angle, however, does not change the fact that $\left< \psi^\dagger(z) \psi (z') \right>$ is the propagator of a euclidean Dirac theory with action (\ref{eq:Dirac}). The only difference with the introduction is that, now, the action written in $z$-coordinates must be \begin{equation*} \displaystyle S \, = \, \frac{1}{\pi} \int d^2 z \, e^{\sigma+i\theta} \left[ \left( \begin{array}{cc} \psi^\dagger & \overline{\psi}^\dagger \end{array} \right) \left( \begin{array}{cc} - \frac{i}{2} \overset{\leftrightarrow}{\partial}_{\bar{z}} - (\partial_{\bar{z}} \varphi) & 0 \\ 0 & \frac{i}{2} \overset{\leftrightarrow}{\partial}_{z} - (\partial_{z} \varphi^*) \end{array}\right) \left( \begin{array}{c} \psi \\ \overline{\psi} \end{array} \right) \right] \, , \end{equation*} while in the introduction there was no phase $e^{i \theta}$. But, comparing this to the covariant expression (\ref{eq:Dirac}), we see that this phase factor can be absorbed in the choice of the tetrad, $$\label{eq:tetrad_honey} e^{\mu}_{\phantom{a}a} \, = \, e^{- \sigma} \left( \cos \theta \, \delta^\mu_{\phantom{a}a} \, - \, \sin \theta\, \varepsilon^\mu_{\phantom{a}a} \right)$$ where $\varepsilon^{1}_{\phantom{1}2} = -\varepsilon^{2}_{\phantom{1}1} = 1$ and $\varepsilon^{1}_{\phantom{1}1} = \varepsilon^{2}_{\phantom{1}2} =0$. In complex coordinates (using complex coordinates both in the tangent space and in the moving frame) the components of the tetrad are $e^z_{\phantom{a}z} = e^{-\sigma + i \theta}$, $e^{\bar{z}}_{\phantom{a}\bar{z}} = e^{-\sigma - i \theta}$ and $e^z_{\phantom{a}\bar{z}} = e^{\bar{z}}_{\phantom{a}z}=0$. Multiplied by the Jacobian from the volume element $\sqrt{g} d^2 {\rm x} = |e|^{-1} d^2 {\rm x}$, this gives $|e|^{-1} e^{z}_{\phantom{z}z} = e^{\sigma + i\theta}$ in the action above. Thus, we end up with exactly the same conclusion as in the introduction: the long-range correlations inside the arctic ellipse follow from the Dirac action (\ref{eq:Dirac}), with the tetrad (\ref{eq:tetrad_honey}) and with background gauge potentials that are again derivatives of the function $\varphi(x,y)$, \begin{equation*} A^{({\rm v})}_\mu \, =\, - i \partial_{\mu} {\rm Im}\, \varphi \qquad {\rm and}\qquad A^{({\rm a})}_\mu \, =\, - \partial_{\mu} {\rm Re} \,\varphi \, . \end{equation*} Finally, there is one point that we left aside in the introduction, that we wish to discuss now: the boundary conditions obeyed by the continuous fields. We use the coordinate system $({\rm x}^1, \,{\rm x}^2) \, = \, ({\rm Re} z , \, {\rm Im} z)$; the boundaries are then at ${\rm x}^1 = 0$ and ${\rm x}^1=\pi$. We claim that the boundary conditions are \begin{equation*} \begin{array}{lll} {\rm x}^1 \, = \, 0 \,: &\quad & \left\{ \begin{array}{rcl} e^{i\frac{\pi}{4} + i {\rm Re}\, \varphi} \psi^\dagger & = & e^{-i \frac{\pi}{4}} e^{-i {\rm Re}\, \varphi} \overline{\psi}^\dagger \\ e^{-i\frac{\pi}{4} - i {\rm Re}\, \varphi} \psi & = & e^{i \frac{\pi}{4}} e^{i {\rm Re} \,\varphi} \overline{\psi} \\ \end{array} \right. \\ \\ {\rm x}^1 \, = \, \pi \quad {\rm and} \quad |{\rm x}^2|< {\rm arcth}\left( \frac{1-u^2}{1+u^2} \right) : &\quad & \left\{ \begin{array}{rcl} e^{i\frac{\pi}{4} + i {\rm Re} \,\varphi} \psi^\dagger & = & -\, e^{-i \frac{\pi}{4}} e^{-i {\rm Re} \, \varphi} \overline{\psi}^\dagger \\ e^{-i\frac{\pi}{4} - i {\rm Re}\, \varphi} \psi & = & \, -\, e^{i \frac{\pi}{4}} e^{i {\rm Re} \, \varphi} \overline{\psi} \end{array} \right. \\ \\ {\rm x}^1 \, = \, \pi \quad {\rm and} \quad |{\rm x}^2|> {\rm arcth}\left( \frac{1-u^2}{1+u^2} \right) : &\quad & \left\{ \begin{array}{rcl} e^{i\frac{\pi}{4} + i {\rm Re}\, \varphi} \psi^\dagger & = & e^{-i \frac{\pi}{4}} e^{-i {\rm Re}\, \varphi} \overline{\psi}^\dagger \\ e^{-i\frac{\pi}{4} - i {\rm Re}\, \varphi} \psi & = & e^{i \frac{\pi}{4}} e^{i {\rm Re} \, \varphi} \overline{\psi} \, . \end{array} \right. \end{array} \end{equation*} To see that, one can focus on the propagator in the DWIS, $\mathinner{\langle{\Psi_0}|} c^\dagger(z) c(z') \mathinner{|{\Psi_0}\rangle}$, which is equal to $1/(2i \sin \frac{z-z'}{2})$. It is clear from that formula that $c^\dagger (z) = c^\dagger (-z^*)$ if ${\rm x}^1 = {\rm Re} \, z = 0$. Then it follows from (\ref{eq:scaling_ident_hex}) that $e^{i \frac{\pi}{4}} e^{\frac{\sigma + i \theta}{2}} e^{i \varphi} \psi^\dagger \, = \, e^{-i \frac{\pi}{4}} e^{\frac{\sigma-i \theta}{2}} e^{-i \varphi^*} \overline{\psi}^\dagger$. Moreover, one can check that $\theta = 0$ if ${\rm x}^1 = 0$; this gives the boundary condition relating $\psi^\dagger$ and $\overline{\psi}^\dagger$ above. The other cases are similar; we use the fact that $c^\dagger (z) = -c^\dagger (-z^*)$ if ${\rm x}^1 = \pi$, and that $\theta =0$ or $\pi$ if ${\rm x}^1 = \pi$ and $|{\rm x}^2 | < {\rm arcth} \left( \frac{1-u^2}{1+u^2} \right)$ or $|{\rm x}^2 | > {\rm arcth} \left( \frac{1-u^2}{1+u^2}\right)$ respectively (see Fig. 6). Notice that these boundary conditions look complicated only because we want to keep the dependence on ${\rm Re} \varphi$; clearly, if we took the liberty of making an axial gauge transformation $\psi \rightarrow e^{ i {\rm Re}\, \varphi} \psi$, $\overline{\psi} \rightarrow e^{- i {\rm Re}\, \varphi} \overline{\psi}$, then the phases above would drop. The remaining factors $e^{\pm i \frac{\pi}{4}}$ can be traced back to the fact that, in the $z$-coordinate system, the boundary is vertical. Using the $\zeta$-coordinates instead (see Fig. 6), the boundary would simply be the real axis, and these phase factors would drop as well. Therefore, the only meaningful information in the boundary conditions above is the sign.

### Density profile

What is the density profile inside the ellipse ? Like in the introduction, one way of obtaining the density is to take the limit $x' \rightarrow x$ of the propagator $\left< c^\dagger_{x,y} c_{x',y} \right>$. From the previous paragraphs, we know that the short-distance behavior of the propagator is \begin{equation*} \left< c^\dagger_{x,y} c_{x',y} \right> \, \simeq \, \frac{e^{-\frac{\sigma+\sigma'}{2}}}{2\pi i} \left[ \frac{e^{-i \frac{\theta+\theta'}{2}} e^{-i (\varphi - \varphi')} } {z-z'} - \frac{e^{i \frac{\theta+\theta'}{2}} e^{i (\varphi^* - \varphi'^*)}} {z^*-z'^*} \right] \, . \end{equation*} Expanding the $\varphi$-exponentials to the first order, one finds \begin{equation*} \left< c^\dagger_{x,y} c_{x',y} \right> \, \underset{x' \rightarrow x}{ \simeq} \, \frac{e^{-\sigma}}{\pi} \left[ {\rm Im} \left( \frac{e^{-i \theta}}{z-z'} \right) - {\rm Re} \left( \frac{e^{-i \theta} (\varphi - \varphi')}{z-z'} \right) \, + \, \dots \right] \, . \end{equation*} It looks like this expression diverges; however it does not: in fact, one can check from the above explicit formulas that $e^{\sigma + i \theta} \partial_x z = -1$ (see also the Introduction). Using this, we see that ${\rm Im} \left( \frac{e^{-\sigma -i \theta}}{z-z'} \right) \, \simeq \, {\rm Im} \left( \frac{e^{-\sigma-i \theta}}{\partial_x z} \frac{1}{x-x'} \right) \, = \, 0$. Then, exactly like in the introduction, we are left with \begin{eqnarray*} \left< \rho_{x,y} \right> &=& \underset{x' \rightarrow x}{\rm lim} \left< c^\dagger_{x,y} c_{x',y} \right> \, = \, - \frac{1}{\pi} {\rm Re} \left( \frac{ \partial_x \varphi} {e^{\sigma+i \theta} \partial_x z} \right) \\ &=& -\frac{1}{\pi} A_x^{({\rm a})} \, . \end{eqnarray*} We then evaluate $\varphi(x,y)$ explicitly and find that, inside the ellipse, \begin{eqnarray} \label{eq:Rephi_honeycomb} \nonumber {\rm Re} \, \varphi(x,y) &=& x \, {\rm arccos}\left( \frac{(1+u^2)x - N u^2}{u \sqrt{(N-2x)^2-y^2}} \right) \, + \, \frac{|y|}{2}\, {\rm arccos} \left( \frac{ N(N-2x) - \frac{1+u^2}{1-u^2} y^2} { \sqrt{ (N^2-y^2) ((N-2x)^2 -y^2)}} \right) \\ && \qquad \qquad - \, \frac{N}{2} \,{\rm arccos} \left( \frac{ (1-2u^2) (N-x)^2 + x^2-y^2 } { \sqrt{ (N^2-y^2) ((N-2x)^2 -y^2)}} \right) \, . \end{eqnarray} Differentiating with respect to $x$, one gets $A_x^{({\rm a})} = - \partial_x {\rm Re} \varphi$. The formula (\ref{eq:Rephi_honeycomb}) holds inside the ellipse only, but similar tricks work also outside the ellipse. The final result is $$\label{eq:densityhexa} \left< \rho_{x,y} \right> \;=\; \left\{\begin{array}{ccl} \frac{1}{\pi}\arccos \left(\frac{(1+u^2)x-N u^2}{u\sqrt{(N-2x)^2-y^2}}\right)&&{\rm if} \quad X^2+y^2<N^2 \, , \\ \\ 1&&{\rm if} \quad X\leq -\sqrt{N^2-y^2} \, , \\ \\ 1&& {\rm if} \quad X\geq \sqrt{N^2-y^2} \quad \textrm{ and} \quad 2x<N-|y| \quad \textrm{ and} \quad \frac{y}{N} \geq \frac{1-u^2}{1+u^2} \, , \\ \\ 0 && \textrm{otherwise} \,. \end{array} \right.$$ Numerical checks of the formula are shown in Fig. 7, and show excellent agreement with (\ref{eq:densityhexa}).

Figure 7 Figure 7. Left: Numerical density profile in a system of height $N=128$, for an anisotropy $u=1/2$. The simulation were performed on a finite system of size $L=512$. The color convention is the following: Green is $\rho=1$, Blue is $\rho=0$, and intermediate values are shown with a mixture of green, yellow and red. The fluctuating region is in perfect qualitative agreement with Eq. (\ref{eq:densityhexa}), which gives the inside of an ellipse $(N/2+3x/2)^2+y^2<N^2$. Right: Quantitative comparison between a numerical simulation for the same system size and the exact asymptotic solution (\ref{eq:densityhexa}) for $y=0$ (blue stars), $y=-N+N/4=-3N/4$ (green triangles), $y=-N+N/16=-15N/16$ (red circles). The black lines are the analytical solution. The agreement is excellent.
Alternative derivation of the density profile. The reader may be unsatisfied with the way we obtained the density profile, by a limiting process which may look suspicious. Of course, other derivations of the result (\ref{eq:densityhexa}) exist; let us quickly sketch another one, that perhaps looks more mathematically sound. From the previous paragraphs, one gets that the average density $\left< \rho(x,y) \right> = \left< c_{x,y}^\dagger c_{x,y}\right>$ is given by the integral \begin{equation*} \left< \rho_{x,y} \right> \, = \, \int \frac{dk}{2\pi} \int \frac{dk'}{2\pi} e^{- i (k-k') x + y (\varepsilon(k)-\varepsilon(k')) - i N (\tilde{\varepsilon}(k)-\tilde{\varepsilon}(k'))} \frac{1}{2 i \sin\left( \frac{k-k'}{2} \right)} \, . \end{equation*} Trying to analyze this integral by stationary phase, one faces the issue that the two stationary points (for $k$ and $k'$) are identical, so the denominator is singular. To circumvent this issue, one can differentiate with respect to $y$ under the integral (notice that $y$, although discrete in our model, may be treated as a continuous parameter in the integral in the r.h.s). This differentiation removes the singularity in the integrand: for a point $(x,y)$ inside the ellipse, we find \begin{eqnarray*} \frac{\partial \left< \rho_{x,y} \right>} {\partial y} & = & \int \frac{dk}{2\pi} \int \frac{dk'}{2\pi} e^{- i (k-k') x + y (\varepsilon(k)-\varepsilon(k')) - i N (\tilde{\varepsilon}(k)-\tilde{\varepsilon}(k'))} \frac{\varepsilon(k)-\varepsilon(k')}{2 i \sin\left( \frac{k-k'}{2} \right)} \\ & \underset{N \rightarrow \infty}{\simeq} & \frac{e^{-\sigma-i \theta}}{2\pi i} \,\frac{d \varepsilon}{dk} (z)\, + \, \frac{e^{-\sigma+i\theta}}{2\pi i} \, \frac{d \varepsilon}{dk} (-z^*) \\ & = & \frac{1}{\pi} \, {\rm Im} \left[ e^{-\sigma-i\theta}\,\frac{d \varepsilon}{dk} (z) \right] \\ &=& \frac{1}{\sqrt{R^2-X^2-y^2}}\times\frac{y(x-u^2(N-x))}{\pi u((N-2x)^2-y^2)}\, . \end{eqnarray*} Integrating back, we get the result (\ref{eq:densityhexa}) inside the ellipse, up to an undetermined additive constant, that may depend on $x$. This constant can be fixed by setting $y=0$ and using a continuity argument. Outside the ellipse, one may repeat similar arguments; the derivative vanishes everywhere, except on the segments $N-2x=\pm y$, see Fig. 6.

## Dimers on the square lattice

In this section, the method introduced in section 2. is applied to the dimer model on the Aztec diamond, where the arctic circle phenomenon has been originally discovered [1]. As we shall see, it is possible to recover this exact geometry using our formalism. This case is slightly more intricate than the dimers on honeycomb, as the corresponding fermionic problem has two bands. We will consider a slight generalization where all dimers do not have the same weights. This is done by assigning a fugacity $u$ to all horizontal dimers, and a fugacity $v$ to all vertical dimers, so that the partition function reads $$Z=\sum_{\{c\}} u^{\# {\rm\, horizontal \;dimers}} v^{\# {\rm \,vertical\; dimers}} .$$ The sum runs over all dimer coverings of the square lattice. Mapping the dimer model onto a system of free fermions is standard, see e.g. Refs. [5456], and can be done as follows. First, we consider a particular configuration, where all dimers are arranged in a staggered fashion. This is shown in Fig. 8(i); in the following we call it Reference configuration. Then, any dimer configuration will be compared to this reference configuration, by superimposition of the two. For each dimer configuration, this procedure generates a set of lattice paths, an example being shown in Fig. 8(ii). It is convenient to think of these lattice paths as trajectories of particles propagating (say) upwards, where a particle is defined as a vertical link occupied by a reference dimer only, or by a “real” dimer only. While moving from one row to the next, the particles have to follow the rules shown in Fig. 8(iii) (Each particle is represented as a black zigzag line in the Figure).

Figure 8 Figure 8. Mapping onto a system of non intersecting lattice paths (or particle trajectories). (i): Reference configuration with staggered dimers (shown in red). (ii): A configuration of dimers (shown in thick blue) is superimposed on the reference configuration. This generates four lattice paths which propagate from bottom to top. (iii): Rules for the trajectories. Recall a particle is defined as a vertical link occupied by a reference dimer only, or a vertical link occupied by a real dimer only. In the latter case the particle, represented by a black zigzag line, has to go straight ahead (bottom right). In the former there are three possibilities, the particle can either go straight ahead (bottom left), to the left (upper left) or to the right (upper right). A weight $u$ (resp. $v$) is assigned to horizontal (resp. vertical) dimers.
The particle trajectories have the following important properties:

• There is a one to one correspondence between the valid particle trajectories and the valid dimer coverings.
• The number of particle is conserved.
• The trajectories do not cross.

The later property means that we need not be careful about the statistics of the particles. As for dimers on the honeycomb lattice, we are allowed to represent them using fermionic operators, and use a transfer matrix acting on those operators. From now one we will call the particle fermions. We explicitly construct the transfer matrix in the following subsection. Before doing so, let us explain how the celebrated Aztec diamond may be recovered using our formalism, and the DWIS $\mathinner{|{\Psi_0}\rangle}=\prod_{x<0}c_x^\dagger \mathinner{|{0}\rangle}$. In terms of dimer occupancies, the DWIS corresponds to an alternation $\ldots 101010$ for sites with index $x<0$ and $010101\ldots$ for sites with index $x>0$. This is shown in Fig. 9(i). The (imaginary-)time evolution is generated by several application of the transfer matrix, and the system is projected back onto the initial state $\mathinner{|{\Psi_0}\rangle}$. Due to the staggered nature of the DWIS, many other dimer occupancies are automatically set to one, by construction of the dimer model. They are shown in thick green in Fig. 9(i). As can be seen, the remaining edges that are left free define precisely an Aztec diamond. In the following, we use this simple observation to compute all the large scale correlations, and derive the effective low energy theory. The phenomenology here is strikingly different from the homogenous case (see e. g. [46,57]) where the standard machinery of boundary CFT can be safely used. Let us mention also that a slightly different transfer matrix approach has been used in [58] to deal with more generic lattices.

### The transfer matrix

Figure 9 Figure 9. (i): Aztec diamond obtained when the DWIS $\mathinner{|{\Psi_0}\rangle}$ is set as bottom and top boundary condition. The corresponding dimers cross either of the two dotted line, and are shown in thick green. Other dimers occupancies are set to one as a result; they are also shown in green. The remaining edges define the inside of an Aztec diamond, which is highlighted in thick black. (ii): Example of fermion trajectories generated by successive application of the transfer matrix on the DWIS.
Let us now write the transfer matrix in explicit form. First, we are still dealing with fermionic particles on the lattice $\mathbb{Z} + \frac{1}{2}$, and with the DWIS $\mathinner{|{\Psi_0}\rangle}$. The transfer matrix, however, is not invariant under translations of one site, but only of two sites. To have a consistent labeling of the edges, it is necessary to distinguish between even and rows, so that there are in fact two transfer matrices $T_{\rm e}$ and $T_{\rm o}$. Their respective action on the many-fermion states follows from the rules shown in Fig. 8(iii). We have \begin{eqnarray}\label{eq:even_action} &&T_{\rm e}\ c^\dagger_{2j+1/2}\ T_{\rm e}^{-1}=v \,c^\dagger_{2j+1/2},\nonumber \\ &&T_{\rm e}\ c^\dagger_{2j-1/2}\ T_{\rm e}^{-1}=u \,c^\dagger_{2j-3/2}+v \,c^\dagger_{2j-1/2}+u\, c^\dagger_{2j+1/2}, \end{eqnarray} and \begin{eqnarray}\label{eq:odd_action} &&T_{\rm o}\ c^\dagger_{2j+1/2}\ T_{\rm o}^{-1}=u \,c^\dagger_{2j+3/2}+v\, c^\dagger_{2j+1/2}+u \,c^\dagger_{2j-1/2},\nonumber \\ &&T_{\rm o}\ c^\dagger_{2j-1/2}\ T_{\rm o}^{-1}=v \,c^\dagger_{2j-1/2}. \end{eqnarray} The action of the transfer matrix on the DWIS is illustrated in Fig. 9(ii), where an example of a possible trajectory for the fermions is shown. To diagonalize the transfer matrix, it is convenient to go to momentum space. Recall our convention for the Fourier transform is $$\label{eq:fouriertransform} c^\dagger(k)=\sum_{x\in \mathbb{Z}+1/2}e^{ikx} c_x^\dagger \qquad,\qquad c^\dagger(k+\pi)=\sum_{x\in \mathbb{Z}+1/2} e^{i\pi x}e^{ikx}c_x^\dagger.$$ The action of the transfer matrix on those modes reads $$T_{\rm e} \left(\begin{array}{c}c^\dagger(k)\\ c^\dagger(k+\pi)\end{array}\right) T_{\rm e}^{-1}=\left(\begin{array}{cc} v+u\cos k & -iu \cos k \\ -iu \cos k & v-u\cos k \end{array}\right) \left(\begin{array}{c}c^\dagger(k)\\ c^\dagger(k+\pi)\end{array}\right)$$ and $$T_{\rm o} \left(\begin{array}{c}c^\dagger(k)\\ c^\dagger(k+\pi)\end{array}\right) T_{\rm o}^{-1}=\left(\begin{array}{cc} v+u\cos k & iu \cos k \\ iu \cos k & v-u\cos k \end{array}\right) \left(\begin{array}{c}c^\dagger(k)\\ c^\dagger(k+\pi)\end{array}\right).$$ $T_{\rm e}$ and $T_{\rm o}$ do not satisfy the assumption 2.2. of section 2., because they are not hermitian, however $T^2=T_{\rm o}T_{\rm e}$ is. Its action on the Fourier modes is given by $$T^2 \left(\begin{array}{c}c^\dagger(k)\\ c^\dagger(k+\pi)\end{array}\right) T^{-2}=M(k)\left(\begin{array}{c}c^\dagger(k)\\ c^\dagger(k+\pi)\end{array}\right)\quad,\quad M(k)=\left(\begin{array}{cc} u^2\cos^2 k+(v+u \cos k)^2 & 2i u^2 \cos ^2 k \\ -2iu^2 \cos^2 k & u^2\cos^2 k+(v-u \cos k)^2 \end{array}\right).$$ $M(k)$ is hermitian, so may be diagonalized by unitary matrix: $$\label{eq:umat} U(k)=\left(\begin{array}{cc} \cos \alpha(k) & i\sin \alpha(k) \\ i\sin\alpha(k) & \cos \alpha(k) \end{array}\right).$$ We parametrize the angle $\alpha$ as $$\alpha(k)=\frac{1}{2}\arcsin \left(\frac{u}{w}\cos \kappa\right)\qquad,\qquad w=\sqrt{u^2+v^2}.$$ This choice of parametrization might appear somewhat artificial. However, as we shall see, it makes many subsequent computations much lighter. It is easy to check that $U(k)M(k)U^\dagger(k)$ is diagonal provided the condition $$\label{eq:param} \tan \kappa =\tan \kappa(k)=\frac{v}{w}\tan k$$ is fulfilled. The previous equation defines a map $k\mapsto \kappa(k)$ which is a bijection of $[-\pi/2,\pi/2]$ onto itself. $\kappa$ will be the natural variable to perform the steepest descent method. Using this change of basis we have $$U(k)M(k)U^\dagger(k)=\left(\begin{array}{cc} \lambda(k) & 0 \\ 0 & 1/\lambda(k) \end{array}\right)\qquad,\qquad \lambda(k)=\left(u \cos k+\sqrt{v^2+u^2 \cos^2 k}\right)^2= v^2\frac{w+u \cos \kappa}{w-u \cos \kappa}.$$ Introducing a new set of fermions through $$\left(\begin{array}{c}c_{+}^\dagger(k)\\ c_{-}^\dagger(k)\end{array}\right)=U(k) \left(\begin{array}{c}c^\dagger(k)\\ c^\dagger(k+\pi)\end{array}\right)$$ allows to rewrite the transfer matrix as \begin{eqnarray}\label{eq:tmsquare} T^2&=&\exp\left[-\int_{-\pi/2}^{\pi/2}\frac{dk}{2\pi} 2\varepsilon(k)\left(c_{+}^\dagger(k)c_+(k)-c_{-}^\dagger(k)c_-(k)\right)\right],\\ \label{eq:square_energy} \varepsilon(k)&=&-\frac{1}{2}\log \lambda(k)= -\frac{1}{2}\log \left(v^2\frac{w+u \cos \kappa}{w-u \cos \kappa}\right). \end{eqnarray} The form (\ref{eq:tmsquare}) obviously satisfies the condition 2.2.. The last requirement, 2.2., may be fulfilled by extending the map $k\mapsto \kappa(k)$ to $\mathbb{R}$ in the following way: \begin{eqnarray} \kappa(k+ p\pi)&=&\kappa(k)+p\pi\qquad\,,\qquad \forall p\in \mathbb{Z} \end{eqnarray} In particular, this implies $c^\dagger_+(k+\pi)=c_-^\dagger(k)$ and $\varepsilon(k+\pi)=-\varepsilon(k)$, so that the transfer matrix may be rewritten as $$T^2=\exp\left[-\int_{-\pi}^\pi \frac{dk}{2\pi}2\varepsilon(k)c_+^\dagger(k)c_+(k)\right].$$

### The Hilbert transform

Now that all the assumptions 2.2.,2.2.,2.2. are satisfied, we are ready to apply the formalism developed in section 2.. First, we need the propagator in the DWIS. We find, after a long calculation, \begin{eqnarray}\nonumber \left\langle \Psi_0|c_+^\dagger(k)c_{+}(k')|\Psi_0\right\rangle&=&e^{i\frac{\beta(\kappa)-\beta(\kappa')}{2}} \frac{\cos \frac{\kappa-\kappa'}{2}}{i \sin (k-k')}\\ \label{eq:correl_cp_square} &=&\left(\frac{d\kappa}{dk}\right)^{1/2}\left(\frac{d\kappa'}{dk'}\right)^{1/2} \frac{e^{i\frac{\beta(\kappa)-\beta(\kappa')}{2}}}{2i \sin \frac{\kappa-\kappa'}{2}}. \end{eqnarray} where $\kappa=\kappa(k)$, $\kappa'=\kappa(k')$, $\beta=\beta(k)$, and $\beta'=\beta(k')$. The ancillary angle $\beta$ is given by $$\beta(k)=\arctan\left(\frac{u}{v} \sin \kappa(k)\right).$$ The transform needed to obtain the full propagator $\mathinner{\langle{c_+^\dagger(k,y)c_{+}(k',y')}\rangle}$ is $$\tilde{\varepsilon}(k)={\rm p.v.}\int \frac{dk'}{2\pi}\varepsilon(k')2\partial_{k'} \log \mathinner{\langle{\Psi_0|c_+^\dagger(k')c_{+}(k)|\Psi_0}\rangle}.$$ We first remark that the phase in (ref.) only gives a constant – independent on $k$ – contribution to $\tilde{\varepsilon}(k)$. Since only the difference $\tilde{\varepsilon}(k)-\tilde{\varepsilon}(k')$ enters any correlation function, this contribution may simply be ignored. Therefore, we set $\beta,\beta'=0$ in the remainder of the present subsection. Using Eq. (\ref{eq:correl_cp_square}), we find that the appropriate transform is a Hilbert transform, modulo conjugation by the change of variables $k\mapsto\kappa(k)$: \begin{eqnarray} \nonumber \tilde{\varepsilon}(k) & = & {\rm p.v.} \int_{-\pi}^\pi \frac{dk'}{2\pi} \varepsilon(k') 2\partial_{k'} \log \mathinner{\langle{\Psi_0}|} c^\dagger_+(k) c_+(k') \mathinner{|{\Psi_0}\rangle} \\ \nonumber &=& {\rm p.v.} \int_{-\pi}^\pi \frac{d\kappa'}{2\pi} \varepsilon(k') 2\partial_{\kappa'} \log \left[ \frac{1}{2 i \sin \left( \frac{\kappa - \kappa'}{2} \right)} \right] \, + \, {\rm p.v.} \int_{-\pi}^\pi \frac{d k'}{2\pi} \varepsilon(k') \partial_{k'} \log \left[ \frac{d\kappa'}{d k'} \right] \\ \nonumber &=& {\rm p.v.} \int_{-\pi}^\pi \frac{d\kappa'}{2\pi} \varepsilon(\kappa') \cot \left( \frac{\kappa - \kappa'}{2} \right) \quad + \quad 0\, . \end{eqnarray} In the second line, we have split the function $\log \mathinner{\langle{\Psi_0}|} c^\dagger_+(k) c_+(k') \mathinner{|{\Psi_0}\rangle}$ into a sum of three terms, one of which trivially vanishes; then in the first term we changed variables $k' \rightarrow \kappa'$, using the fact that the Jacobian $\frac{dk'}{d\kappa'}$ in the integration measure cancels the one coming from the chain rule for the derivative $\partial_{k'} = \frac{d\kappa'}{dk'} \partial_{\kappa'}$. In this term, we also use the notation $\varepsilon (\kappa')$ for $\varepsilon(k')$. In the third line, the second term is zero because of the anti-periodicity of $\varepsilon(k)$, $\varepsilon(k+\pi) = -\varepsilon (k)$. Indeed, the derivative $\frac{d\kappa}{dk}$ is periodic (with period $\pi$), therefore the integrand $\varepsilon(\kappa') \partial_{k'} \log \left[ \frac{d\kappa'}{d k'} \right]$ is anti-periodic, and its integral from $-\pi$ to $0$ cancels the one from $0$ to $+\pi$.

We call $\tilde{\varepsilon}(\kappa)$ the Hilbert transform of the function (of $\kappa$) $\varepsilon(\kappa)$. For the function of interest here, we have the result $$\label{eq:hilbertsquare} \varepsilon(\kappa)=-\frac{1}{2}\log \left(v^2\frac{w+u \cos \kappa}{w-u \cos \kappa}\right) \qquad \Leftrightarrow \qquad \tilde{\varepsilon}(\kappa)=\frac{i}{2}\log\left(\frac{v+i u \sin \kappa}{v-i u \sin \kappa}\right)-\log v=-\arctan \left(\frac{u}{v} \sin \kappa\right)-\log v,$$ which is derived in the appendix 10.. We note that the function $\tilde{\varepsilon}(\kappa)$ is automatically anti-periodic with period $\pi$, a property it inherits from $\varepsilon(\kappa)$.

### Fermionic correlators: exact finite-size formulas

Here $N$ is assumed to be even. The vertical axis is labelled by integers $-N\leq y\leq N$, while the horizontal axis is still labelled by half integers $x\in \mathbb{Z}+1/2$. This way the left/right symmetry is $x\to -x$ and the up/down symmetry is $y\to -y$. For an observable $\mathcal{O}_{x,y}$ at position $x$ and $y$, we have \begin{eqnarray} \mathinner{\langle{\mathcal{O}_{x,y}}\rangle}=\left\{ \begin{array}{ccccc} \displaystyle{\frac{\mathinner{\langle{\Psi_0}|}(T^2)^{\frac{N-y}{2}} \mathcal{O}_x (T^2)^{\frac{N+y}{2}}\mathinner{|{\Psi_0}\rangle}} {\mathinner{\langle{\Psi_0|T^{2N}|\Psi_0}\rangle}}}&&,&& y \,\textrm{ even}\\ \\ \displaystyle{\frac{\mathinner{\langle{\Psi_0}|}(T^2)^{\frac{N-y-1}{2}} T_{\rm o}\mathcal{O}_x T_{\rm o}^{-1} (T^2)^{\frac{N+y+1}{2}}\mathinner{|{\Psi_0}\rangle}}{\mathinner{\langle{\Psi_0|T^{2N}|\Psi_0}\rangle}}} &&,&& y \textrm{ odd} \end{array} \right. \end{eqnarray} We will now use the notation introduced in the introduction (see Eq. \ref{eq:notation_doteq}), and also used in Sec. 3.. We write $$\label{eq:ident_square} \left\{ \begin{array}{rcl} c^\dagger_{x,y} & \doteq & \displaystyle \int_{-\pi}^\pi \frac{dk}{2\pi} e^{-i k x + y \varepsilon (\kappa) - i N \tilde{\varepsilon}(\kappa)} \, r_{x,y} (k) \, c^\dagger_+(k) \\ \\ c_{x,y} & \doteq & \displaystyle \int_{-\pi}^\pi \frac{dk}{2\pi} e^{i k x-y \varepsilon (\kappa) + i N \tilde{\varepsilon}(\kappa)} \, s_{x,y} (k) \, c_+(k) \, , \end{array} \right.$$ where the dot means that this is an equality that holds only when taking expectation values in the DWIS. The main difference with the toy model is the emergence of the extra factors $r_{x,y}(k)$, and $s_{x,y}(k)$. These coefficients are here because we are dealing with a two band problem; they only depend on the parity of $x$ and $y$. The aim of this subsection is to compute them.

We start with the case $y\in 2\mathbb{Z}$, which is easiest. Using Eqs. (\ref{eq:fouriertransform}), (\ref{eq:umat}) and extending the Brillouin zone to $[-\pi,\pi]$ we obtain $$c^\dagger_{2j\pm 1/2}=\int_{-\pi}^{\pi}\frac{dk}{2\pi} e^{-i k (2j\pm 1/2)} \sqrt{1\mp (u/w)\cos \kappa}\; c_{+}^\dagger(k).$$ Now due to (\ref{eq:prop_interm_2}), it is easy to see that the coefficients we are looking for are given by $$r_{x,y}(k)\;=\;s_{x,y}(k)\;=\; \sqrt{1\mp \frac{u}{w}\cos \kappa},$$ where the upper sign corresponds to $x\in 2\mathbb{Z}+1/2$ and the lower sign corresponds to $x\in 2\mathbb{Z}-1/2$. Note that the equality between $r_{x,y}(k)$ and $s_{x,y}(k)$ simply follows from hermiticity.

The case $y\in 2\mathbb{Z}+1$ is more tricky. Indeed, one has to shift $y$ by one, and conjugate with one application of $T_{\rm o}$, which is not Hermitian. For $x\in 2\mathbb{Z}+1/2$ we have – see Eq. (\ref{eq:odd_action}) – $T_{\rm o} c_{x}^\dagger T_{\rm o}^{-1}=v c_{x}^\dagger +u(c_{x-1}^\dagger +c_{x+1}^\dagger)$, so that $$c_{x,y}^\dagger \;\doteq \; v \,c_{x,y+1}^\dagger +u\left[c_{x-1,y+1}^\dagger+c_{x+1,y+1}^\dagger\right],$$ which yields \begin{eqnarray} r_{x,y}(k)&=&e^{\varepsilon(\kappa)} \left[v \,r_{x,y+1}+2 u \cos k \,r_{x+1,y+1}\right]\\ &=& \sqrt{1+ \frac{u}{w}\cos \kappa}. \end{eqnarray} For $x\in 2\mathbb{Z}-1/2$ we have $T_{\rm o} c_{x}^\dagger T_{\rm o}^{-1}=v c_{x}^\dagger$, and we find $$r_{x,y}(k)=\sqrt{1- \frac{u}{w}\cos \kappa}.$$ Finally, it remains to perform a similar identification for $s_{x,y}(k)$. To do so, we first need to determine the action of the odd transfer matrix on the annihilation operators. This is explained in Appendix. 12.1.. We find \begin{eqnarray}\label{eq:odd_action_2} &&T_{\rm o}\ c_{2j+1/2}\ T_{\rm o}^{-1}=\frac{1}{v}\, c_{2j+1/2},\nonumber \\ &&T_{\rm o}\ c_{2j-1/2}\ T_{\rm o}^{-1}=\frac{1}{v^2}\left(-u c_{2j-3/2}+ v c_{2j-1/2}-u c_{2j+1/2}\right), \end{eqnarray} and a similar calculation gives $r_{x,y}(k)=s_{x,y}(k)$ for $y$ odd. In total, all possible cases may be summarized by the following simple formula $$\label{eq:rsk_result} r_{x,y}(k)\,=\,s_{x,y}(k)\,=\, \sqrt{1-\sigma \frac{u}{w}\cos \kappa(k)}\qquad,\qquad \sigma=(-1)^{x-y-1/2},$$ which depends only on the choice of sublattice. In the following, we will refer to the case $x-y\in 2\mathbb{Z}+1/2$, for which $\sigma=+1$ as the even sublattice. The odd sublattice, $x-y\in 2\mathbb{Z}-1/2$, has $\sigma=-1$.

### Fermionic correlators: scaling limit

Once again, we are interested in the scaling limit of the correlators. These are obtained by the steepest descent method, exactly like in the introduction. We focus on the argument of the exponentials in Eq. (\ref{eq:ident_square}). The stationary points are the solutions of \begin{equation*} \frac{x}{u} \, +\, i \frac{y}{v} \sin \kappa - \frac{N}{w} \cos \kappa \, = \, \frac{1}{u}\,\frac{d}{dk} \left[ k x \,+\, i y \varepsilon(\kappa) + N \tilde{\varepsilon}(\kappa) \right] \, = \, 0 \, . \end{equation*} Here we have used the chain rule $\frac{d}{d k} = \frac{d\kappa}{dk} \frac{d}{d\kappa}$ and the formulas (see Eq. (\ref{eq:hilbertsquare})) \begin{equation*} \frac{d \varepsilon(\kappa)}{d \kappa} \, = \, \frac{ u w \, \sin \kappa} {w^2 - u^2 \cos^2 \kappa} \, , \qquad \frac{d \tilde{\varepsilon}(\kappa)}{d \kappa} \, = \, \frac{- u v \, \cos \kappa} {w^2 - u^2 \cos^2 \kappa} \, ,\qquad \frac{d k}{d \kappa} \, =\, \frac{v w}{ w^2 - u^2 \cos^2 \kappa} \, . \end{equation*} Generically, there are two non-degenerate stationary points \begin{equation*} \kappa(x,y) \, =\, z(x,y) \qquad {\rm and} \qquad \kappa(x,y) \, = \, -z^*(x,y) \, , \end{equation*} where $$\label{eq:z_square} z(x,y) \, = \, \arccos \frac{x/u}{\sqrt{ (\frac{N}{w} )^2 - (\frac{y}{v} )^2}} \,-\, i \,{\rm arctanh} \left( \frac{y/v}{N/w} \right) \, .$$ By 'generically', we mean 'unless $\frac{x^2}{u^2} + \frac{y^2}{v^2} = \frac{N^2}{w^2}$'. This curve is, of course, the arctic curve for the dimer model on the Aztec diamond. Now let us assume that $(x,y)$ is not on the arctic curve; then around the two stationary points, the exponential is approximated by \begin{equation*} \left\{ \begin{array}{rcl} e^{ -i N \left[ k \frac{x}{N} + i \frac{y}{N} \varepsilon(k) + \tilde{\varepsilon}(k) \right]} & \simeq & e^{ -i \varphi(x,y ) - i \frac{1}{2} (k - k(z) )^2 \, e^{\sigma (x,y)} \left(\frac{dk}{d\kappa}\right)^{-1} } \qquad {\rm around} \quad k = k(z) \\ \\ e^{ -i N \left[ k \frac{x}{2N} + i \frac{y}{N} \varepsilon(k) + \tilde{\varepsilon}(k) \right]} & \simeq & e^{ i \varphi^*(x,y ) - i \frac{1}{2} (k - k(-z^*) )^2 \, e^{\sigma (x,y)} \left(\frac{dk^*}{d\kappa^*}\right)^{-1} } \qquad {\rm around} \quad k = k(-z^*) , \end{array} \right. \end{equation*} where $$\label{eq:phi_sigma_square} \left\{ \begin{array}{rcl} \varphi \,= \, \varphi(x,y) & \equiv & k(z(x,y)) x + i y \varepsilon(z(x,y)) + N \tilde{\varepsilon}(z(x,y)) \\ \\ e^{\sigma} \, = \, e^{\sigma(x,y)} &\equiv& \left(\frac{dk}{d\kappa}\right) \frac{d^2}{dk^2} \left[ i y \varepsilon(k) + N \tilde{\varepsilon}(k) \right] \, = \, u \, \frac{d}{d\kappa} \left[ -i \frac{y}{v} \, \sin \kappa + \frac{N}{w} \, \cos \kappa \right] \, =\, u \, \sqrt{ \left(\frac{N}{w}\right)^2 - \left(\frac{x}{u}\right)^2 - \left(\frac{y}{v}\right)^2} \, . \end{array} \right.$$ Performing the gaussian integration, we find that the identification (\ref{eq:ident_square}) becomes $$\label{eq:ident_square_2} \left\{ \begin{array}{rcl} c^\dagger_{x,y} & \underset{N \rightarrow \infty}{\doteq} & \displaystyle \frac{ e^{-i \frac{\pi}{4}}} {\sqrt{2\pi}} e^{ -i \varphi} e^{-\frac{\sigma} {2}} \, r_{x,y} (z) \, \left( \frac{dk}{d \kappa} \right)^{\frac{1}{2}} c^\dagger_+(z) \; + \; \frac{ e^{i \frac{\pi}{4}}} {\sqrt{2\pi}} e^{ i \varphi^*} e^{-\frac{\sigma}{2}} \, r_{x,y} (-z^*) \, \left( \frac{dk^*}{d \kappa^*} \right)^{\frac{1}{2}} c^\dagger_+(-z^*) \\ \\ c_{x,y} & \underset{N \rightarrow \infty}{\doteq} & \displaystyle \frac{ e^{i \frac{\pi}{4}}} {\sqrt{2\pi}} e^{ i \varphi} e^{-\frac{\sigma} {2}} \, s_{x,y} (z) \, \left( \frac{dk}{d \kappa} \right)^{\frac{1}{2}} c_+(z) \; + \; \frac{ e^{-i \frac{\pi}{4}}} {\sqrt{2\pi}} e^{ -i \varphi^*} e^{-\frac{\sigma} {2}} \, s_{x,y} (-z^*) \, \left( \frac{dk^*}{d \kappa^*} \right)^{\frac{1}{2}} c_+(-z^*) \, . \end{array} \right.$$ Notice that the appearance of $\left( \frac{dk}{d \kappa} \right)^{\frac{1}{2}} c^\dagger_+(z)$ and $\left( \frac{dk}{d \kappa} \right)^{\frac{1}{2}} c_+(z)$ is particularly convenient, since the roots of the Jacobians are exactly what is needed to cancel the ones in (\ref{eq:correl_cp_square}). Their propagator is $$\label{eq:someprop} \mathinner{\langle{\Psi_0}|} \left( \frac{dk}{d \kappa} \right)^{\frac{1}{2}} c^\dagger_+(z) \left( \frac{dk'}{d \kappa'} \right)^{\frac{1}{2}} c_+(z') \mathinner{|{\Psi_0}\rangle} \, =\, \frac{e^{i\frac{\beta(z)-\beta(z')}{2}}}{2 i \sin \left( \frac{z-z'}{2}\right)} \, .$$ We note that it is also possible to get rid of the the phase factor $e^{i\frac{\beta-\beta'}{2}}=e^{-i\frac{\tilde{\varepsilon}-\tilde{\varepsilon}'}{2}}$ in (\ref{eq:someprop}), and absorb it in the phase $\varphi(x,y)$. We will make use of this observation in the next subsection.

### Identification of the Dirac theory inside the arctic circle

At this point, it is quite clear that the Dirac action in the dimer model is going to be similar to the one discussed in the introduction. There are however a couple differences. First, the relation between the lattice and continuous degrees of freedom is more complicated than in the case of the XX chain or honeycomb dimers, see Eq. (\ref{eq:lattice_continuous_XX}). $$\label{eq:square_c_psi} \left\{ \begin{array}{rcl} c_{x,y}^\dagger &=& \displaystyle \sum_{n \in \mathbb{Z}} \gamma_{x,y,n} \left[ \frac{1}{\sqrt{2\pi}} \psi^\dagger(x+n,y) + \frac{1}{\sqrt{2\pi}} \overline{\psi}^\dagger(x+n,y) \right] \\ c_{x,y} &=& \displaystyle \sum_{n \in \mathbb{Z}} \gamma_{x,y,n} \left[ \frac{1}{\sqrt{2\pi}} \psi(x+n,y) + \frac{1}{\sqrt{2\pi}} \overline{\psi}(x+n,y) \right] \end{array} \right.$$ where the $\gamma_{x,y,n}$ are the $n$-th Fourier coefficients of $r_{x,y}$, i.e. $\gamma_{x,y,n}=\int_{-\pi}^{\pi} \frac{dk}{2\pi}e^{-ikn}r_{x,y}$. Since the $r_{x,y}$ are analytic in $k$, the coefficients decay exponentially fast with $n$. In that sense the identification (\ref{eq:square_c_psi}), despite the infinite number of terms, is still local. The propagators of $\psi^\dagger$, $\psi$ and $\overline{\psi}^\dagger$, $\overline{\psi}$ are \begin{eqnarray*} \left< \psi^\dagger(x,y) \psi(x',y') \right> & = & e^{-i (\varphi - \varphi')} e^{-\frac{\sigma+\sigma'}{2}} \mathinner{\langle{\Psi_0}|} \left( \frac{dk}{d \kappa} \right)^{\frac{1}{2}} c^\dagger_+(z) \left( \frac{dk'}{d \kappa'} \right)^{\frac{1}{2}} c_+(z') \mathinner{|{\Psi_0}\rangle} \, = \, \frac{e^{-i (\varphi+\frac{\tilde{\varepsilon}}{2} - \varphi'-\frac{\tilde{\varepsilon}'}{2})} e^{-\frac{\sigma+\sigma'}{2}}}{2i \sin \left( \frac{z-z'}{2} \right)} \\ \left< \overline{\psi}^\dagger(x,y) \overline{\psi}(x',y') \right> & = & e^{-i (\varphi - \varphi')} e^{-\frac{\sigma+\sigma'}{2}} \mathinner{\langle{\Psi_0}|} \left( \frac{dk^*}{d \kappa^*} \right)^{\frac{1}{2}} c^\dagger_+(-z^*) \left( \frac{dk'^*}{d \kappa'^*} \right)^{\frac{1}{2}} c_+(-z'^*) \mathinner{|{\Psi_0}\rangle} \, = \, \frac{ e^{i (\varphi^*+\frac{\tilde{\varepsilon}^*}{2} - \varphi'^*-\frac{\varepsilon'^*}{2})} e^{-\frac{\sigma+\sigma'}{2}}}{2i \sin \left( \frac{z^*-z'^*}{2} \right)} \, . \end{eqnarray*} where $\tilde{\varepsilon}=\tilde{\varepsilon}(z)$. These are exactly the same propagators as in the introduction, with functions $\varphi$ and $\sigma$ that are now given by (\ref{eq:phi_sigma_6v}), and $\varphi$ replaced by $\varphi+\tilde{\varepsilon}/2$. We thus find that the field theory inside the arctic circle is again the generic Dirac action (\ref{eq:Dirac}), with a diagonal tetrad $e^\mu_a \, = \, e^{-\sigma} \delta_{a}^\mu$, and with gauge fields $A^{({\rm a})}_\mu \, = \, -\partial_\mu {\rm Re} [\varphi+\tilde{\varepsilon}/2]$ and $A^{({\rm v})}_\mu \, = \,- i\partial_\mu {\rm Im} [\varphi+\tilde{\varepsilon}/2]$. Let us finally comment on the boundary conditions for the fields $\psi,\psi^\dagger, \overline{\psi}, \overline{\psi}^\dagger$. As is explained in detail in Sec. 3.4., they follows from the relation $1/\sin \left(\frac{z-z'}{2}\right)=\pm 1/\sin\left(\frac{z^*-z'}{2}\right)$, where the upper sign corresponds to the boundary ${\rm Re}\, z=0$, and the lower sign corresponds to the boundary ${\rm Re}\, z=\pi$.

### Density profile

For the sake of completeness, we also show the density profile $\mathinner{\langle{\rho_{x,y}}\rangle}=\mathinner{\langle{c_{x,y}^\dagger c_{x,y}}\rangle}$ inside the arctic ellipse. A possible way to derive it is to take the derivative with respect to $y$: \begin{eqnarray*} \frac{\partial \mathinner{\langle{\rho_{x,y}}\rangle}}{\partial y} & = & \int \frac{dkdk'}{(2\pi)^2} e^{- i (k-k') x + y (\varepsilon(\kappa)-\varepsilon(\kappa')) - i N (\tilde{\varepsilon}(\kappa)-\tilde{\varepsilon}(\kappa'))} r_{x,y} (k) s_{x,y} (k') \left( \frac{d\kappa}{dk} \right)^{\frac{1}{2}} \left( \frac{d\kappa'}{dk'} \right)^{\frac{1}{2}} \frac{(\varepsilon(\kappa)-\varepsilon(\kappa'))e^{i \frac{\beta(\kappa)-\beta(\kappa')}{2}}}{2 i \sin\left( \frac{\kappa-\kappa'}{2} \right)} \\ & \underset{N \rightarrow \infty}{\simeq} & \frac{e^{-\sigma}}{2\pi i} r_{x,y}(z) s_{x,y} (z) \,\frac{d \varepsilon}{d \kappa} \, + \, \frac{e^{-\sigma}}{2\pi i} r_{x,y}(-z^*) s_{x,y} (-z^*) \, \frac{d \varepsilon}{d \kappa} (-z^*) \\ & = & \frac{e^{-\sigma}}{\pi} \, {\rm Im} \left[ r_{x,y}(z) s_{x,y} (z) \, \frac{d \varepsilon}{d \kappa} (z) \right] \, . \end{eqnarray*} Then, we use Eq. (\ref{eq:rsk_result}): $$r_{x,y}(k) s_{x,y}(k)= 1\mp \frac{u}{w}\cos \kappa(k),$$ where the upper (resp. lower) sign corresponds to $x-y\in 2\mathbb{Z}+1/2$ (resp. $x\in 2\mathbb{Z}-1/2$), and recall $\frac{d\varepsilon}{d\kappa}=\frac{u w \sin \kappa}{w^2-u^2 \cos^2 \kappa}$. We obtain \begin{eqnarray*} \frac{\partial \mathinner{\langle{\rho_{x,y}}\rangle}}{\partial y} & = &\frac{1}{\pi\sqrt{(N/w)^2-(x/u)^2-(y/v)^2}}\times {\rm Im}\,\left[\frac{\sin z(w\mp u\cos z)}{w^2-u^2 \cos^2 z}\right]\\ & = &\frac{1}{\pi\sqrt{(N/w)^2-(x/u)^2-(y/v)^2}}\times {\rm Im}\,\left[\frac{\sin z}{w\pm u \cos z}\right]\\ &=&\frac{1}{\pi\sqrt{(N/w)^2-(x/u)^2-(y/v)^2}}\times \frac{-y(w^2 x\pm u^2 N)}{uvw((N\pm x)^2-y^2)}. \end{eqnarray*} This integrates to $$\mathinner{\langle{\rho_{x,y}}\rangle}=\frac{1}{\pi}\arccos \left(\frac{w x/u\pm u N/w}{\sqrt{(N\pm x)^2-y^2}}\right).$$ Now recalling that a fermion is a dimer on the even sublattice, but the absence of a dimer on the odd sublattice, the density of dimers on vertical links is given by $\rho_d(x,y)=\rho_{x,y}$ on the even sublattice, and $\rho_d(x,y)=1-\rho_{x,y}$ on the odd. Hence $$\rho_d(x,y)=\frac{1}{\pi}\arccos \left(\frac{\pm w x/u+u N/w}{\sqrt{(N\pm x)^2-y^2}}\right).$$ The results are summarized in Table. I.

Figure 10 Figure 10. Numerical fermionic density plots corresponding to an Aztec diamond of size $N=256$. Left: Isotropic case $u/v=1$, where the arctic curve is a circle. Right: $u/v=1/2$. The arctic curve is an ellipse with excentricity $e=\sqrt{3}/2$.
We finally note that the density of vertical dimers may also be expressed as $$\rho_d(x,y)=\frac{1}{2}-\frac{1}{\pi} \arctan\left(\frac{\frac{N u}{v w}\pm\frac{w x}{u v}}{\sqrt{N^2/w^2-x^2/u^2-y^2/v^2}}\right),$$ which makes the arctic ellipse slightly more apparent. Setting $u=v=1$, $w=\sqrt{2}$ in this formula, we recover the result of Ref. [16].

 $\quad x - y\in 2\mathbb{Z}-1/2\quad$ $\quad x - y \in 2\mathbb{Z}+1/2 \quad$ (odd sublattice) (even sublattice) $\rho_{x,y}$ $\quad\frac{1}{\pi}\arccos\left(\frac{w x/u- u N/w}{\sqrt{(N-x)^2-y^2}}\right)\quad$ $\quad\frac{1}{\pi}\arccos\left(\frac{w x/u+ u N/w}{\sqrt{(N+x)^2-y^2}}\right)\quad$ $\rho_{d}(x,y)$ $\quad\frac{1}{\pi}\arccos\left(\frac{-w x/u+ u N/w}{\sqrt{(N-x)^2-y^2}}\right)\quad$ $\quad\frac{1}{\pi}\arccos\left(\frac{w x/u+ u N/w}{\sqrt{(N+x)^2-y^2}}\right)\quad$
Table I. Average density of fermions $\rho(x,y)$ in the scaling limit, and corresponding density of dimers $\rho_{d}(x,y)$ on vertical edges.

## Six-vertex model at $\Delta = 0$ with domain-wall boundary conditions

It is conventional to draw the vertices of the six-vertex model as

with weights $a,b,c \in \mathbb{R}$. The configurations of the six-vertex model are in one-to-one correspondence with (directed) self-avoiding trajectories on the square lattice as illustrated in Fig. 11.(c). The trajectories can touch—a situation that corresponds to the first vertex with weight $a$—, but they cannot overlap along an edge. By convention, we decide that the trajectories never cross. This is convenient, because it means we do not have to be too careful about the fermionic statistics of the particles whose trajectories we are describing. [We would have to be more careful if we took periodic boundary conditions.] We view the system as a set of fermionic particles on the lattice $\mathbb{Z} + \frac{1}{2}$ that evolve in the $y$-direction, which is imaginary time. The six vertices drawn above are then operators acting on two neighboring sites $x$ and $x+1$, from bottom to top. They can be written respectively as \begin{equation*} \begin{array}{rcl} {\rm weight} \; a \, : && \left\{ \begin{array}{c} c^\dagger_x c_x \, c^\dagger_{x+1} c_{x+1} \\ \\ (1-c^\dagger_x c_x ) \, ( 1-c^\dagger_{x+1} c_{x+1} ) \end{array} \right. \\ \\ {\rm weight} \; b \, : && \left\{ \begin{array}{c} c^\dagger_x c_{x+1} \\ \\ c^\dagger_{x+1} c_x \end{array} \right. \\ \\ {\rm weight} \; c \, : && \left\{ \begin{array}{c} ( 1-c^\dagger_x c_x ) c^\dagger_{x+1} c_{x+1} \\ \\ c^\dagger_x c_x ( 1-c^\dagger_{x+1} c_{x+1} ) \; . \end{array} \right. \end{array} \end{equation*} Summing these six terms, we build the $R$-matrix of the six-vertex model, \begin{eqnarray} \label{eq:R_matrix} \nonumber R_{x,x+1} & = & a \left[ c^\dagger_x c_x c^\dagger_{x+1} c_{x+1} + (1-c^\dagger_x c_x ) \, ( 1-c^\dagger_{x+1} c_{x+1} )\right] + b \left[ c^\dagger_x c_{x+1} + c^\dagger_{x+1} c_x \right] \\ \nonumber && \qquad \qquad \quad + \, c \left[ c^\dagger_x c_x (1-c^\dagger_{x+1} c_{x+1} ) + ( 1-c^\dagger_x c_x ) c^\dagger_{x+1} c_{x+1} \right] \\ &=& \, a \, \exp \left[ \log \left(\frac{b+c}{a}\right) \left( c^\dagger_x c_{x+1} \, + \, c^\dagger_{x+1} c_x \right) \right] \,. \end{eqnarray} This last equality holds provided that $$\Delta \, = \, \frac{a^{2}+b^{2}-c^{2}}{2ab} \, = \, 0 \,,$$ and this is what we assume from now on. Also, to keep formulas as simple as possible, we assume that $a>0$, $b>0$ and $c>0$; this will save us the trouble of having to keep track of the signs. Now, since the $R$-matrix is gaussian at $\Delta = 0$, the model is free: the transfer matrix is gaussian as well, and all observables can be obtained from the propagator using Wick's theorem. The techniques developed in previous parts apply, as we explain below.

However, before we turn to the technicalities about the transfer matrix and Hilbert transforms, let us briefly explain why, in our setting, we recover the six-vertex model with domain-wall boundary conditions—a model that has been studied extensively in the literature both analytically and numerically [5,59,8,23,26]—. We focus once again on the DWIS $\mathinner{|{\Psi_0}\rangle}$, which is completely filled with particles on the left, and completely empty on the right. The evolution of the system is generated by iterations of the transfer matrix. After a number of iterations, we project the one-dimensional system of particles back to the DWIS (Fig. 11.(a)). Because of the initial and final states, most vertices on the left and right of the system are forced to be in the configurations with weight $a$, as illustrated in Fig. 11.(b). Only a central region, which turns out to be a square, is not determined a priori by those initial and final states. The square in the middle hosts configurations of the six-vertex model with boundary conditions (Fig. 11.(b)) that are known traditionally as 'domain-wall boundary conditions'. Therefore, our setup, which combines the DWIS with the transfer matrix, is well suited to revisit this popular problem.

Figure 11 Figure 11. Six-vertex model on the square lattice. We regard the horizontal axis as space, and the vertical axis as imaginary time. (Left) One starts form the domain-wall initial state (DWIS) at the bottom, then we apply the transfer matrix several times, and finally, we project back to the DWIS at the top. (Middle) Because of the DWIS, a large number of edges are constrained to be either empty or filled; the only vertices that are not fixed by the choice of boundary conditions are the ones inside the square region in the middle of the figure. This central square is what is known in the literature as the 'six-vertex model with domain-wall boundary conditions'. (Right) A configuration of vertices compatible with the boundary conditions: note that one can interpret the filled edges as trajectories in space-time. This is the idea of the mapping to a fermion model.

### The transfer matrix

We need to define two transfer matrices, one for the even rows and another for the odd rows, $$\left\{ \begin{array}{rcl} \displaystyle T_{{\rm e}} & \equiv & \displaystyle \prod_{ x \in 2\mathbb{Z} + \frac{1}{2}} R_{x,x+1} \\ \displaystyle T_{{\rm o}} & \equiv & \displaystyle \prod_{x \in 2 \mathbb{Z} - \frac{1}{2}} R_{x,x+1} \, . \end{array} \right.$$ Notice that the terms commute within each product; but, of course, $T_{{\rm e}} T_{{\rm o}} \neq T_{{\rm o}} T_{{\rm e}}$. The weighted configurations of the model are obtained by applying iteratively these two operators, one after the other. One could focus on one of the double-row transfer matrices $T_{\rm e}T_{\rm o}$ or $T_{\rm o}T_{\rm e}$, which generate the imaginary time evolution. However the drawback of both those operators is that they are not hermitian. In order to apply the formalism developed in part 2., it is more convenient to deal with a hermitian generator of translations in imaginary time. Thus, we define our (double-row) transfer matrix as the operator $$T^2 \, = \, T_{{\rm e}}^{\frac{1}{2}} \,T_{{\rm o}} \, T_{{\rm e}}^{\frac{1}{2}}$$ whose hermiticity follows from $T_{{\rm e}}^\dagger \, = \, T_{{\rm e}}$ and $T_{{\rm o}}^\dagger \, = \, T_{\rm o}$. Again, the upperscript/exponent '$2$' in $T^2$ is there to remind us that it is an operator that adds two rows to the system, instead of one. We will always deal with $T^2$, there will be no object '$T$'. Now let's diagonalize $T^2$. In $k$-space, we have \begin{eqnarray*} \left\{ \begin{array}{rcl} \displaystyle T_{{\rm e}} &=& \displaystyle \exp \left[ \log \left(\frac{b+c}{a} \right) \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \frac{dk}{2\pi} \left( \begin{array}{cc} c^\dagger(k) & c^\dagger (k+\pi) \end{array} \right) \left( \begin{array}{cc} \cos k & \sin k \\ \sin k & - \cos k \end{array} \right) \left( \begin{array}{c} c(k) \\ c(k+\pi) \end{array} \right) \right] \\ \\ \displaystyle T_{{\rm o}} &=& \displaystyle \exp \left[ \log \left(\frac{b+c}{a} \right) \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \frac{dk}{2\pi} \left( \begin{array}{cc} c^\dagger(k) & c^\dagger (k+\pi) \end{array} \right) \left( \begin{array}{cc} \cos k & -\sin k \\ -\sin k & - \cos k \end{array} \right) \left( \begin{array}{c} c(k) \\ c(k+\pi) \end{array} \right) \right] \, . \end{array} \right. \end{eqnarray*} [Here and in what follows, we drop the unimportant normalization constant $a^{{\rm number \; of \; sites}}$ that comes from the factor $a$ in the r.h.s of (\ref{eq:R_matrix}).] Then the problem boils down to the diagonalization of a $2 \times 2$ matrix; after some easy algebra, we find \begin{equation*} T^2 \, = \, \exp \left[ -\int_{-\frac{\pi}{2}}^\frac{\pi}{2} \frac{dk}{2\pi} \,2\varepsilon(k)\, \left( c_{+}^\dagger(k) c_{+}(k) \, - \, c_{-}^\dagger(k) c_{-}(k) \right) \right] \end{equation*} where $$\left\{ \begin{array}{rcl} \label{eq:diag_6v} 2\varepsilon(k) &=& \displaystyle - \log\left[ \frac{c + b \cos \kappa(k)} {c - b \cos \kappa(k)} \right] \\ \\ \displaystyle \left( \begin{array}{c} c_{+}^\dagger(k) \\ c_{-}^\dagger(k) \end{array} \right) &=& \displaystyle U(k) \left( \begin{array}{c} c^\dagger(k) \\ c^\dagger(k+\pi) \end{array}\right) \, = \, \displaystyle \left( \begin{array}{cc} \cos \left(\frac{k-\kappa(k)}{2}\right) & \sin \left(\frac{k-\kappa(k)}{2}\right) \\ \\ - \sin \left(\frac{k-\kappa(k)}{2}\right) & \cos \left(\frac{k-\kappa(k)}{2}\right) \end{array} \right) \left( \begin{array}{c} c^\dagger(k) \\ c^\dagger(k+\pi) \end{array} \right) \, , \end{array} \right.$$ and $\kappa(k)$ is the real-valued function $\kappa : [-\frac{\pi}{2},\frac{\pi}{2}] \rightarrow [-\frac{\pi}{2}, \frac{\pi}{2}]$ defined by $$\label{eq:pseudo_mom} \tan \kappa(k) \, =\, \frac{a}{c} \, \tan k \, .$$ The mapping $k\mapsto \kappa(k)$ is always bijective (since we assume $a > 0$, $c>0$), so it is possible to switch from the variable $k$ to the new variable $\kappa$ and to reparametrize all functions of $k$, like $\varepsilon(k) \rightarrow \varepsilon(\kappa)$ or $\tilde{\varepsilon}(k) \rightarrow \tilde{\varepsilon}(k)$ (introduced in the next paragraph). This change of variables $k \rightarrow \kappa$ is instrumental in the appearance of the Hilbert transform in the six-vertex model, as we explain in the next section. To keep notations as light as possible, it will be implicit in most formulas, unless specified otherwise. Let us make one last important remark about $\kappa$: the mapping $k \mapsto \kappa(k)$ can be analytically continued to $\mathbb{R}$, and it has the property $\kappa ( k + p \pi) \,= \, \kappa (k) + p \pi$ for any $p \in \mathbb{Z}$. In the formulas below, we use the function $\kappa : \mathbb{R} \rightarrow \mathbb{R}$ and this property without recalling it.

Finally, notice that the three assumptions of part 2. are satisfied: in particular, $\varepsilon(k)$ is an anti-periodic function of $k$ with period $\pi$ (namely, $\varepsilon(k+ \pi) \,=\, - \varepsilon(k)$), and the relation (\ref{eq:periodU}) holds, such that $c^\dagger_+(k+\pi) = c^\dagger_-(k)$.

### $k \mapsto \kappa$ and the Hilbert transform

We are now in position to apply the formalism of part 2. to the transfer matrix $T^2$. First, we need the propagator in the DWIS, $\mathinner{\langle{\Psi_0}|} c^\dagger_+(k) c_+(k') \mathinner{|{\Psi_0}\rangle}$; it is computed very easily from Eqs. (\ref{eq:diag_6v})-(\ref{eq:pseudo_mom}), leading to the explicit result $\cos (\frac{\kappa-\kappa'}{2} ) / [i \,\sin (k - k')]$, where $\kappa = \kappa(k)$ and $\kappa' = \kappa(k')$. For our purposes though, it is more convenient to write this equation in the form $$\label{eq:correl_cp_6v} \mathinner{\langle{\Psi_0}|} c^\dagger_+(k) c_+(k') \mathinner{|{\Psi_0}\rangle} \, = \, \left( \frac{d\kappa}{d k} \right)^{\frac{1}{2}} \left( \frac{d\kappa'}{d k'} \right)^{\frac{1}{2}} \frac{ 1} {2 i \sin (\frac{\kappa - \kappa'}{2})} \, .$$ As anticipated in part 2., the transform $\tilde{\varepsilon}_{U}(k)$ that appears is nothing but the Hilbert transform, modulo conjugation by the change of variable $k \rightarrow \kappa$: \begin{eqnarray} \nonumber \tilde{\varepsilon}_{U}(k) & = & {\rm p.v.} \int_{-\pi}^\pi \frac{dk'}{2\pi} \varepsilon(k') 2\partial_{k'} \log \mathinner{\langle{\Psi_0}|} c^\dagger_+(k) c_+(k') \mathinner{|{\Psi_0}\rangle} \\ \nonumber &=& {\rm p.v.} \int_{-\pi}^\pi \frac{d\kappa'}{2\pi} \varepsilon(k') 2\partial_{\kappa'} \log \left[ \frac{1}{2 i \sin \left( \frac{\kappa - \kappa'}{2} \right)} \right] \, + \, {\rm p.v.} \int_{-\pi}^\pi \frac{d k'}{2\pi} \varepsilon(k') \partial_{k'} \log \left[ \frac{d\kappa'}{d k'} \right] \\ \nonumber &=& {\rm p.v.} \int_{-\pi}^\pi \frac{d\kappa'}{2\pi} \varepsilon(\kappa') \cot \left( \frac{\kappa - \kappa'}{2} \right) \quad + \quad 0\, . \end{eqnarray} In the second line, we have split the function $\log \mathinner{\langle{\Psi_0}|} c^\dagger_+(k) c_+(k') \mathinner{|{\Psi_0}\rangle}$ into a sum of three terms, one of which trivially vanishes; then in the first term we changed variables $k' \rightarrow \kappa'$, using the fact that the Jacobian $\frac{dk'}{d\kappa'}$ in the integration measure cancels the one coming from the chain rule for the derivative $\partial_{k'} = \frac{d\kappa'}{dk'} \partial_{\kappa'}$. In this term, we also use the notation $\varepsilon (\kappa')$ for $\varepsilon(k')$. In the third line, the second term is zero because of the anti-periodicity of $\varepsilon(k)$, $\varepsilon(k+\pi) = -\varepsilon (k)$. Indeed, the derivative $\frac{d\kappa}{dk}$ is periodic (with period $\pi$), therefore the integrand $\varepsilon(\kappa') \partial_{k'} \log \left[ \frac{d\kappa'}{d k'} \right]$ is anti-periodic, and its integral from $-\pi$ to $0$ cancels the one from $0$ to $+\pi$.

We call $\tilde{\varepsilon}(\kappa) = \tilde{\varepsilon}_U(k(\kappa))$ the Hilbert transform of the function $\varepsilon(\kappa)$—viewed as a function of $\kappa$, not $k$—. For the function $\varepsilon (\kappa)$ of interest here, we have (see the appendix 10.): $$\label{eq:eps_epst_6v} 2\varepsilon(\kappa) \, = \, - \log\left[ \frac{c + b \cos \kappa}{c - b \cos \kappa} \right] \quad \Leftrightarrow \quad 2\tilde{\varepsilon}_{U}(k) \,=\, 2 \tilde{\varepsilon}(\kappa) \, = \, i \log\left[ \frac{a \, + \, i \, b\, \sin \kappa} { a \, - \, i \, b\, \sin \kappa } \right] \, .$$ [Here we are using $a^2+b^2=c^2$, and also the assumption that $a>0$, $b>0$, $c>0$, which was announced earlier; there could be a minus sign in the r.h.s otherwise.] We note that the function $\tilde{\varepsilon}(\kappa)$ is automatically anti-periodic with period $\pi$, a property it inherits from $\varepsilon(\kappa)$.

### Fermionic correlators: exact finite-size formulas

In this section we give the explicit formulas for the fermionic correlators. These formulas are exact in finite size. The most convenient way to express these results is, like in the introduction, to exhibit the relation between correlators $\left< c^\dagger_{x_1,y_1} \dots c^\dagger_{x_n,y_n} c_{x'_1,y'_1} \dots c_{x'_n,y'_n} \right>$ in the two-dimensional geometry, and (known) correlators in the DWIS $\mathinner{\langle{\Psi_0}|} c^\dagger(k_1) \dots c^\dagger(k_n) c(k'_1) \dots c(k'_n) \mathinner{|{\Psi_0}\rangle}$, based on the results of part 2.. Using the same notation as in the introduction—see Eq. (\ref{eq:notation_doteq})—this relation can be expressed through the identification $$\label{eq:ident_6v} \left\{ \begin{array}{rcl} c^\dagger_{x,y} & \doteq & \displaystyle \int_{-\pi}^\pi \frac{dk}{2\pi} e^{-i k x + y \varepsilon (\kappa) - i N \tilde{\varepsilon}(\kappa)} \, r_{x,y} (k) \, c^\dagger_+(k) \\ \\ c_{x,y} & \doteq & \displaystyle \int_{-\pi}^\pi \frac{dk}{2\pi} e^{i k x-y \varepsilon (\kappa) + i N \tilde{\varepsilon}(\kappa)} \, s_{x,y} (k) \, c_+(k) \, . \end{array} \right.$$ The functions $r_{x,y}(k)$ and $s_{x,y}(k)$ can be computed by combining the results of the previous sections. The calculation itself is a little bit cumbersome, with various cases that need to be distinguished, so we defer it to the appendix 12.. The final result, however, takes a simple form. We find $$\label{eq:rsk_vertex} r_{x,y}(k)=s_{x,y}(k)=\sqrt{1-i \,(-1)^{x-y}\, \frac{b}{a}\sin \kappa} \, .$$ The reason why various cases need to be distinguished is the alternation of $T_{\rm e}$ and $T_{\rm o}$ in our transfer matrix setup. For an observable $\mathcal{O}_{x,y}$ at point $(x,y) \in (\mathbb{Z} + \frac{1}{2})\times (\mathbb{Z} + \frac{1}{2})$, the formula for the expectation value depends on the parity of $y-\frac{1}{2}$. More explicitly, assuming that $N$ is odd as in the figure 11 [the case $N$ even would require small shifts in the powers of $T^2$], the formula reads $$(N \;{\rm odd}) \qquad \quad \left< \mathcal{O}_{x,y} \right> \, = \, \left\{ \begin{array}{lll} \displaystyle \frac{ \mathinner{\langle{\Psi_0}|} T_{{\rm e}}^{1/2} T^{2(\frac{N-y-1/2}{2})} \, T_{{\rm e}}^{1/2} \mathcal{O}_x T_{{\rm e}}^{-1/2} \, T^{2(\frac{N+y+1/2}{2})} T_{{\rm e}}^{1/2} \mathinner{|{\Psi_0}\rangle} } {\mathinner{\langle{\Psi_0}|} T_{{\rm e}}^{1/2} T^{2N} T_{{\rm e}}^{1/2} \mathinner{|{\Psi_0}\rangle}} & \quad & {\rm if} \quad y \in 2 \mathbb{Z} + \frac{1}{2} \\ \\ \displaystyle \frac{ \mathinner{\langle{\Psi_0}|} T_{{\rm e}}^{1/2} T^{2(\frac{N-y+1/2}{2})} \, T_{{\rm e}}^{-1/2} \mathcal{O}_x T_{{\rm e}}^{1/2} \, T^{2(\frac{N+y-1/2}{2})} T_{{\rm e}}^{1/2} \mathinner{|{\Psi_0}\rangle} } {\mathinner{\langle{\Psi_0}|} T_{{\rm e}}^{1/2} T^{2N} T_{{\rm e}}^{1/2} \mathinner{|{\Psi_0}\rangle}} & \quad & {\rm if} \quad y \in 2 \mathbb{Z} - \frac{1}{2} \, . \end{array} \right.$$ So it is clear that, in the derivation of the functions $r_{x,y}(k)$ and $s_{x,y}(k)$, one needs to treat the two cases $y \in 2 \mathbb{Z} \pm \frac{1}{2}$ separately, because they involve different conjugations by $T_{{\rm e}}^{\pm 1/2}$. The cases $x \in 2 \mathbb{Z} \pm \frac{1}{2}$ are distinguished for convenience; the formulas are more compact if one treats them separately.

### Fermionic correlators: scaling limit

Once again, we are interested in the scaling limit of the corelators. These are obtained by the steepest descent method, exactly like in the introduction. We focus on the argument of the exponentials in Eq. (\ref{eq:ident_6v}). The stationary points are the solutions of \begin{equation*} \frac{x}{b} \, +\, \left( -i \frac{y}{a} \sin \kappa + \frac{N}{c} \cos \kappa \right) \, = \, \frac{1}{b}\,\frac{d}{dk} \left[ k x \,+\, i y \varepsilon(\kappa) + N \tilde{\varepsilon}(\kappa) \right] \, = \, 0 \, . \end{equation*} Here we have used the chain rule $\frac{d}{d k} = \frac{d\kappa}{dk} \frac{d}{d\kappa}$ and the formulas \begin{equation*} \frac{d \varepsilon(\kappa)}{d \kappa} \, = \, \frac{- b c \, \sin \kappa} {c^2 - b^2 \cos^2 \kappa} \, , \qquad \frac{d \tilde{\varepsilon}(\kappa)}{d \kappa} \, = \, \frac{ a b \, \cos \kappa} {c^2 - b^2 \cos^2 \kappa} \, ,\qquad \frac{d k}{d \kappa} \, =\, \frac{a c}{ c^2 - b^2 \cos^2 \kappa} \, . \end{equation*} Generically, there are two non-degenerate stationary points \begin{equation*} \kappa(x,y) \, =\, z(x,y) \qquad {\rm and} \qquad \kappa(x,y) \, = \, -z^*(x,y) \, , \end{equation*} where $$\label{eq:z_6v} z(x,y) \, = \, - \frac{\pi}{2} \,- \, {\rm arcsin} \frac{x/b}{\sqrt{ (\frac{N}{c} )^2 - (\frac{y}{a} )^2}} \,-\, i \,{\rm arcth} \left( \frac{y/a}{N/c} \right) \, .$$ By 'generically', we mean 'unless $\frac{x^2}{b^2} + \frac{y^2}{a^2} = \frac{N^2}{c^2}$'. This ellipse is, of course, the arctic curve for the six-vertex model with domain-wall boundary conditions. Now let us assume that $(x,y)$ is inside the ellipse; then around the two stationary points, the exponential is approximated by \begin{equation*} \left\{ \begin{array}{rcl} e^{ -i N \left[ k \frac{x}{N} + i \frac{y}{N} \varepsilon(k) + \tilde{\varepsilon}(k) \right]} & \simeq & e^{ -i \varphi(x,y ) - i \frac{1}{2} (k - k(z) )^2 \, e^{\sigma (x,y)} \left(\frac{dk}{d\kappa}\right)^{-1} } \qquad {\rm around} \quad k = k(z) \\ \\ e^{ -i N \left[ k \frac{x}{2N} + i \frac{y}{N} \varepsilon(k) + \tilde{\varepsilon}(k) \right]} & \simeq & e^{ i \varphi^*(x,y ) - i \frac{1}{2} (k - k(-z^*) )^2 \, e^{\sigma (x,y)} \left(\frac{dk^*}{d\kappa^*}\right)^{-1} } \qquad {\rm around} \quad k = k(-z^*) \, , \end{array} \right. \end{equation*} where $$\label{eq:phi_sigma_6v} \left\{ \begin{array}{rcl} \varphi \,= \, \varphi(x,y) & \equiv & k(z(x,y)) x + i y \varepsilon(z(x,y)) + N \tilde{\varepsilon}(z(x,y)) \\ \\ e^{\sigma} \, = \, e^{\sigma(x,y)} &\equiv& \left(\frac{dk}{d\kappa}\right) \frac{d^2}{dk^2} \left[ i y \varepsilon(k) + N \tilde{\varepsilon}(k) \right] \, = \, b \, \frac{d}{d\kappa} \left[ -i \frac{y}{a} \, \sin \kappa + \frac{N}{c} \, \cos \kappa \right] \, =\, b \, \sqrt{ \left(\frac{N}{c}\right)^2 - \left(\frac{x}{b}\right)^2 - \left(\frac{y}{a}\right)^2} \, . \end{array} \right.$$ Performing the gaussian integration, we find that the identification (\ref{eq:ident_6v}) becomes $$\label{eq:ident_6v_2} \left\{ \begin{array}{rcl} c^\dagger_{x,y} & \underset{N \rightarrow \infty}{\doteq} & \displaystyle \frac{ e^{-i \frac{\pi}{4}}} {\sqrt{2\pi}} e^{ -i \varphi} e^{-\frac{\sigma} {2}} \, r_{x,y} (z) \, \left( \frac{dk}{d \kappa} \right)^{\frac{1}{2}} c^\dagger_+(z) \; + \; \frac{ e^{i \frac{\pi}{4}}} {\sqrt{2\pi}} e^{ i \varphi^*} e^{-\frac{\sigma} {2}} \, r_{x,y} (-z^*) \, \left( \frac{dk^*}{d \kappa^*} \right)^{\frac{1}{2}} c^\dagger_+(-z^*) \\ \\ c_{x,y} & \underset{N \rightarrow \infty}{\doteq} & \displaystyle \frac{ e^{i \frac{\pi}{4}}} {\sqrt{2\pi}} e^{ i \varphi} e^{-\frac{\sigma} {2}} \, s_{x,y} (z) \, \left( \frac{dk}{d \kappa} \right)^{\frac{1}{2}} c_+(z) \; + \; \frac{ e^{-i \frac{\pi}{4}}} {\sqrt{2\pi}} e^{ -i \varphi^*} e^{-\frac{\sigma} {2}} \, s_{x,y} (-z^*) \, \left( \frac{dk^*}{d \kappa^*} \right)^{\frac{1}{2}} c_+(-z^*) \, . \end{array} \right.$$ Notice that the appearance of $\left( \frac{dk}{d \kappa} \right)^{\frac{1}{2}} c^\dagger_+(z)$ and $\left( \frac{dk}{d \kappa} \right)^{\frac{1}{2}} c_+(z)$ is particularly convenient, since the roots of the Jacobians are exactly what is needed to cancel the ones in (\ref{eq:correl_cp_6v}). Their propagator is $$\label{eq:prop_6v_wJac} \mathinner{\langle{\Psi_0}|} \left( \frac{dk}{d \kappa} \right)^{\frac{1}{2}} c^\dagger_+(z) \left( \frac{dk'}{d \kappa'} \right)^{\frac{1}{2}} c_+(z') \mathinner{|{\Psi_0}\rangle} \, =\, \frac{1}{2 i \sin \left( \frac{z-z'}{2}\right)} \, .$$

### Identification of the Dirac theory inside the arctic circle

At this point, it is quite clear that the Dirac action in the six-vertex model is going to be the same as in the introduction, with the new $\varphi$ and $\sigma$ given in Eq. (\ref{eq:phi_sigma_6v}). The difference is rather in the relation between the lattice and continuous degrees of freedom, which is more complicated than in the case of the XX chain, see Eq. (\ref{eq:lattice_continuous_XX}). Indeed, we want a relation between the lattice fermions $c^\dagger_{x,y}$ and the two components of the Dirac spinor $\Psi^\dagger = \left( \begin{array}{cc} \psi & \overline{\psi} \end{array} \right)$ that is local, and that does not depend on position, apart from possible sub-lattice distinctions. Looking back at Eqs. (\ref{eq:ident_6v})-(\ref{eq:rsk_vertex}), we see that this is obtained by expanding the factors $r_{x,y}(k) = s_{x,y}(k)$ in Fourier series: \begin{equation*} r_{x,y}(k) \, = \,s_{x,y}(k) \, = \, \sum_{n \in \mathbb{Z}} \alpha_{x,y,n} \, e^{i n k} \, . \end{equation*} The coefficients $\alpha_{x,y,n}$ decay exponentially with $|n|$ because $r_{x,y}(k)$ is analytic in $k$ (and, of course, periodic with period $2 \pi$). Plugging this expansion in Eq. (\ref{eq:ident_6v}) the coefficients given by Eq. (\ref{eq:rsk_vertex}), and using the analysis of the previous section, we arrive at $$\label{eq:6v_c_psi} \left\{ \begin{array}{rcl} c_{x,y}^\dagger &=& \displaystyle \sum_{n \in \mathbb{Z}} \alpha_{x,y,n} \left[ \frac{1}{\sqrt{2\pi}} \psi^\dagger(x-n,y) + \frac{1}{\sqrt{2\pi}} \overline{\psi}^\dagger(x-n,y) \right] \\ c_{x,y} &=& \displaystyle \sum_{n \in \mathbb{Z}} \alpha_{x,y,n} \left[ \frac{1}{\sqrt{2\pi}} \psi(x+n,y) + \frac{1}{\sqrt{2\pi}} \overline{\psi}(x+n,y) \right] \, . \end{array} \right. \\ \\$$ The propagators of $\psi^\dagger$, $\psi$ and $\overline{\psi}^\dagger$, $\overline{\psi}$ are \begin{eqnarray*} \left< \psi^\dagger(x,y) \psi(x',y') \right> & = & e^{-i (\varphi - \varphi')} e^{-\frac{\sigma+\sigma'}{2}} \mathinner{\langle{\Psi_0}|} \left( \frac{dk}{d \kappa} \right)^{\frac{1}{2}} c^\dagger_+(z) \left( \frac{dk'}{d \kappa'} \right)^{\frac{1}{2}} c_+(z') \mathinner{|{\Psi_0}\rangle} \, = \, e^{-i (\varphi - \varphi')} e^{-\frac{\sigma+\sigma'}{2}} \frac{1}{2i \sin \left( \frac{z-z'}{2} \right)} \\ \left< \overline{\psi}^\dagger(x,y) \overline{\psi}(x',y') \right> & = & e^{-i (\varphi - \varphi')} e^{-\frac{\sigma+\sigma'}{2}} \mathinner{\langle{\Psi_0}|} \left( \frac{dk}{d \kappa} \right)^{\frac{1}{2}} c^\dagger_+(z) \left( \frac{dk'}{d \kappa'} \right)^{\frac{1}{2}} c_+(z') \mathinner{|{\Psi_0}\rangle} \, = \, e^{-i (\varphi - \varphi')} e^{-\frac{\sigma+\sigma'}{2}} \frac{1}{2i \sin \left( \frac{z-z'}{2} \right)} \, . \end{eqnarray*} These are exactly the same propagators as in the introduction, with functions $\varphi$ and $\sigma$ that are now given by (\ref{eq:phi_sigma_6v}). We thus find one again that the field theory inside the arctic circle is the generic Dirac theory with action (\ref{eq:Dirac}), with a diagonal tetrad $e^\mu_a \, = \, e^{-\sigma} \delta_{a}^\mu$, and with gauge fields $A^{({\rm a})}_\mu \, = \,- \partial_\mu {\rm Re}\, \varphi$ and $A^{({\rm v})}_\mu \, = \,- i\partial_\mu {\rm Im} \,\varphi$. The boundary conditions for the fields $\psi$, $\psi^\dagger$, $\overline{\psi}$, $\overline{\psi}^\dagger$ follow again from $1/\sin\left( \frac{z-z'}{2} \right) = \pm 1/\sin\left( \frac{-z^*-z'}{2} \right)$ at the two boundaries ${\rm Re} \,z = 0$ and ${\rm Re} \, z = \pi$ repectively; see the more detailed discussion in section 3.4..

### Density profile

Let us conclude this part on the six-vertex model with the density profiles, which we can now derive easily, using the above formulas. We start from the definition of the density operator, and from the above formula with the propagator, \begin{eqnarray*} \left< \rho_{x,y} \right> &= & \, \left< c_{x,y}^\dagger c_{x,y} \right> \\ & = & \int_{-\pi}^\pi \frac{dk}{2\pi} \int_{-\pi}^\pi \frac{dk'}{2\pi} e^{- i k x + y \varepsilon(\kappa) - i N \tilde{\varepsilon}(\kappa)} e^{ i k' x - y \varepsilon(\kappa') +i N \tilde{\varepsilon}(\kappa')} r_{x,y} (k) s_{x,y} (k') \left( \frac{d\kappa}{dk} \right)^{\frac{1}{2}} \left( \frac{d\kappa'}{dk'} \right)^{\frac{1}{2}} \frac{1}{2 i \sin\left( \frac{\kappa-\kappa'}{2} \right)} \, . \end{eqnarray*} Again, we want to evaluate the double-integral with the steepest descent method. There is a problem with the denominator, due to the fact that the stationary points are going to be exactly the same for $k$ and $k'$ (or equivalently $\kappa$ and $\kappa'$). Therefore, we apply the same trick as in part 3.: we treat $y$ as a continuous variable and differentiate with respect to $y$ inside the integral. That way, we obtain \begin{eqnarray*} \partial_y \left< \rho_{x,y} \right> & = & \int_{-\pi}^\pi \frac{dk}{2\pi} \int_{-\pi}^\pi \frac{dk'}{2\pi} e^{- i k x + y \varepsilon(\kappa) - i N \tilde{\varepsilon}(\kappa)} e^{ i k' x - y \varepsilon(\kappa') +i N \tilde{\varepsilon}(\kappa')} r_{x,y} (k) s_{x',y} (k') \left( \frac{d\kappa}{dk} \right)^{\frac{1}{2}} \left( \frac{d\kappa'}{dk'} \right)^{\frac{1}{2}} \frac{\varepsilon(\kappa)-\varepsilon(\kappa')}{2 i \sin\left( \frac{\kappa-\kappa'}{2} \right)} \\ & \underset{N \rightarrow \infty}{\simeq} & \frac{e^{-\sigma}}{2\pi i} r_{x,y}(z) s_{x,y} (z) \,\frac{d \varepsilon}{d \kappa} \, + \, \frac{e^{-\sigma}}{2\pi i} r_{x,y}(-z^*) s_{x,y} (-z^*) \, \frac{d \varepsilon}{d \kappa} (-z^*) \\ & = & \frac{e^{-\sigma}}{\pi} \, {\rm Im} \left[ r_{x,y}(z) s_{x,y} (z) \, \frac{d \varepsilon}{d \kappa} (z) \right] \, . \end{eqnarray*} Now let's focus first on the case $x - y \in 2\mathbb{Z}$, for which $r_{x,y} = s_{x,y} = \sqrt{1 - i \frac{b}{a} \sin \kappa}$. Then, using the above formulas for $e^\sigma$ and for $\frac{d\varepsilon}{d\kappa}$, we get \begin{eqnarray*} \partial_y \left< \rho_{x,y} \right> & = & \frac{1/\pi}{b \sqrt{\left( \frac{N}{c} \right)^2 - \left( \frac{x}{b} \right)^2 - \left( \frac{y}{a} \right)^2}} \, {\rm Im} \left[ \left( 1-i \frac{b}{a} \sin z \right) \frac{- b c \sin z}{(a-i b \sin z) (a+ i b \sin z)} \right] \\ &=& \frac{-1/\pi}{\frac{a}{c} \sqrt{\left( \frac{N}{c} \right)^2 - \left( \frac{x}{b} \right)^2 - \left( \frac{y}{a} \right)^2}} \, {\rm Im} \left[ \frac{ \sin z}{a+ i b \sin z} \right] \\ &=& \frac{-1/\pi}{\frac{a}{bc} \sqrt{\left( \frac{N}{c} \right)^2 - \left( \frac{x}{b} \right)^2 - \left( \frac{y}{a} \right)^2}} \frac{ \frac{x (x+y)}{b^2} - \frac{N^2}{c^2}} { N^2 - (x+y)^2} \, . \end{eqnarray*}

Figure 12 Figure 12. Left: Numerical density profile in a system of height $N=256$ at the free fermion point $\Delta \propto a^2+b^2- c^2 = 0$. Left: $b/a=1$, the arctic curve is a circle. Right: $b/a=1/2$, the arctic curve is an ellipse with excentricity $\sqrt{3}/2$.
To go from the second to the last line, it is useful to notice that $\sin z \, = \, \frac{N}{|c|} \frac{ i \frac{c x y}{a b N} - \sqrt{ \left( \frac{N}{c} \right)^2 - \left( \frac{x}{b} \right)^2 - \left( \frac{y}{a} \right)^2} } { \left(\frac{N}{c}\right)^2 - \left(\frac{y}{a}\right)^2 }$ and $\left| a + i b \sin z \right|^2 \, = \, \frac{ N^2 - (x+y)^2} { \left(\frac{N}{c} \right)^2 - \left( \frac{y}{a} \right)^2}$, and then use those two relations. This can be integrated; it gives $$\left< \rho_{x,y} \right> \, = \, \frac{1}{\pi} {\rm arccos} \left( \frac{ a x/b - b y/a} {\sqrt{ N^2 - (x+y)^2}} \right) \, ,$$ up to some additive function of $x$, which is argued to be zero by looking at the density along the arctic curve $\frac{x^2}{b^2}+\frac{y^2}{a^2}=\frac{N^2}{c^2}$, where it must be either one or vanish. The other case, namely $x-y \in 2\mathbb{Z}+1$, is almost identical, up to a few sign changes. The results are summarized in table II; they agree perfectly with our numerical checks displayed in Fig. 12.

 $x - y \in 2 \mathbb{Z}$ $x - y \in 2 \mathbb{Z}+1$ (SW - NE edges) (NW - SE edges) $\left< \rho_{x,y} \right>$ $\displaystyle \quad \frac{1}{\pi} {\rm arccos} \left( \frac{ a x/b - b y/a} {\sqrt{ N^2 - (x+y)^2}} \right) \quad$ $\displaystyle \quad \frac{1}{\pi} {\rm arccos} \left( \frac{ a x/b + b y/a} {\sqrt{ N^2 - (x-y)^2}} \right) \quad$
Table II. Average occupation of the edges of the six-vertex model in the scaling limit.

## Height mapping and the gaussian free field

In previous parts, we have made extensive use of fermionic degrees of freedom. The fact that all the above models can be mapped to free fermions is instrumental; it is the key that allows to compute correlators exactly on the lattice. Starting from the solutions of these models formulated in fermion language, it is not surprising that we end up with a rather generic free fermion field theory: the Dirac theory in curved two-dimensional euclidean space and in background gauge potentials. In this part we adopt a different point of view. As is well-known in the literature, it is possible to map dimer models and the six-vertex models onto height models, and to make the connection with bosonic free field theory, also called gaussian free field in the mathematical literature ( cf. [60] for a rigorous overview). In the physics literature, such methods are often called Coulomb gas methods (see Ref. [61] for a review); those have played a very important role in the historical development of exact results and field theoretic methods in two-dimensional statistical mechanics. Therefore, we believe it is of some interest to reformulate some of the above results in the bosonic language, which may be more familiar to some readers. When formulated in bosonic language, we expect the continuous field theory describing the correlations inside the arctic circle to be of the form $$\label{eq:boson_curved} S \, = \, \frac{1}{8\pi} \int \sqrt{g} d^2 {\rm x} \, g^{\mu \nu} \, \partial_\mu h \, \partial_\nu h \, ,$$ Going from fermionic degrees of freedom to bosonic ones is a standard procedure in two-dimensional physics, called bosonization (the same procedure we used as a mathematical trick in part 2.). There are a few subtle points, however, that we find sufficiently non-standard to deserve further discussion. Among those, the fate of the background gauge potentials $A^{({\rm v})}$ and $A^{({\rm a})}$ and of the tetrad; those appear in the fermionic action, but are absent from the bosonic one. Thus, for completeness, this part focuses on the mapping to the height field, with some level of detail.

We shall illustrate the procedure with the dimer model on the honeycomb lattice. It is not difficult to adapt this discussion to deal with dimers on the square lattice or with the six-vertex model.

### From dimers to the height field: a quick review

Dimer configurations can be mapped on configurations of a discrete height field. The construction is well-known, see e.g. Refs. [15,62], so we shall go through it quickly.

Figure 13 Figure 13. (Left) A dimer configuration satisfying our boundary conditions, and the corresponding height configuration. ${\bf x}_0$ is the marked hexagon, where the height is set to $h_{{\bf x}_0} =0$. (Right) Same configuration drawn with rhombi, which help visualize the height configuration.
Dimers live on the edges of the honeycomb lattice, while the height field lives on the faces. One needs to mark one of the faces, say at a chosen position ${\bf x}_0$, and decide that the height there is $h_{{\bf x}_0} =0$. Then the height is determined on all other faces by applying the following local rules. If two adjacent hexagons touch along a vertical edge, then, going through the edge in the East$\rightarrow$West direction, the height increases by $+\frac{4\pi}{3}$ if it is occupied by a dimer, and drops by $\frac{2\pi}{3}$ otherwise. The rule is the same for diagonal edges, going Northeast$\rightarrow$Southwest, and Southeast$\rightarrow$Northwest. See Fig. 13.(a) for an illustration. The height field $h$ is related to the fermions as follows. Recall from Part 3. that a fermion on a vertical edge corresponds to the absence of a dimer. Thus the height operator is related to the fermion density through \begin{equation*} h_{x+\frac{1}{2},y} - h_{x-\frac{1}{2},y} \, = \, \frac{4\pi}{3}(1- c^\dagger_{x,y} c_{x,y}) - \frac{2\pi}{3} c^\dagger_{x,y} c_{x,y} \, . \end{equation*} The dimer configurations are frozen at infinity, so the height operator is a scalar as soon as $x$ is large enough: $h_{x,y} = \frac{4\pi}{3} x$ for any $x$ larger than some fixed $M$ (for instance, $M = N/2$ works, see Fig. 5). It follows that \begin{eqnarray*} h_{x,y} & = & h_{M,y} - \sum_{n=0}^{M-x-1} \left( h_{x+n+1,y} - h_{x+n,y} \right) \, = \, \frac{4\pi}{3} x + 2\pi \sum_{n=0}^{M-x-1} c^\dagger_{x+\frac{1}{2}+n,y} c_{x+\frac{1}{2}+n,y} \, . \end{eqnarray*} The average height $\left< h \right>$ can then be evaluated using the results of Part 3., \begin{eqnarray} \label{eq:average_height} \nonumber \left< h_{x,y} \right> & = & \frac{4\pi}{3} x + 2\pi \sum_{n=0}^\infty \left< \rho_{x+\frac{1}{2}+n,y} \right> \\ \nonumber & \simeq & \frac{4\pi}{3} x + 2\pi \int^\infty_x dx' \left< \rho_{x',y} \right> \\ & = & \left\{ \begin{array}{lll} \frac{4\pi} {3} x - 2 \,{\rm Re} \,\varphi (x,y) && {\rm if} \quad X^2 + y^2 \leq N^2 \, , \\ -\frac{2\pi} {3} x && {\rm if} \quad X \leq - \sqrt{ N^2- y^2} \, , \\ -\frac{2\pi} {3} x + \pi (N - |y|) && {\rm if} \quad X \geq \sqrt{ N^2- y^2} \quad {\rm and}\quad 2x \leq N - |y| \quad {\rm and} \quad \frac{|y|}{N} \geq \frac{1-u^2}{1+u^2} \, , \\ \frac{4\pi} {3} x && {\rm otherwise} \, . \end{array} \right. \end{eqnarray} where the function ${\rm Re}\, \varphi(x,y)$ is the one given by Eq. (\ref{eq:Rephi_honeycomb}). This height profile is shown in Fig. 14(b), together with a typical dimer configuration for a large system (Fig. 14(a)).

Figure 14 Figure 14. (Left) A typical configuration for a large width $2N=100$; the fluctuating region is an ellipse. (Right) The average height, as given by Eq. (\ref{eq:average_height}), for $u=1/2$.

### Large-scale fluctuations of the height field

It is customary to decompose the height field $h$ into its classical part $\left< h \right>$, given by (\ref{eq:average_height}), and its zero-average fluctuating part $\tilde{h}$: \begin{eqnarray*} h_{x,y} & = & \left< h_{x,y} \right> \, + \, \tilde{h}_{x,y} \, . \end{eqnarray*} Inside the fluctuating region, we know that we can make the replacements $c^\dagger_{x,y} = \frac{1}{\sqrt{2\pi}} (\psi^\dagger(x,y) + \overline{\psi}^\dagger(x,y))$ and $c_{x,y} = \frac{1}{\sqrt{2\pi}} (\psi(x,y) + \overline{\psi}(x,y))$; thus $\tilde{h}_{x,y}$ is a primitive integral of the form \begin{eqnarray} \label{eq:tildeh} \nonumber \tilde{h}_{x,y} &\simeq & -\int^x dx' : [\psi^\dagger(x',y) + \overline{\psi}^\dagger(x',y)] [\psi(x',y) + \overline{\psi}(x',y)] : \\ &\simeq & -\int^x dx' \left( : \psi^\dagger(x',y) \psi(x',y) : \, + \, : \overline{\psi}^\dagger(x',y) \overline{\psi}(x',y) : \right) \, . \end{eqnarray} Here, normal ordering is defined such that $\left< : \psi^\dagger_{x,y} \psi_{x',y'}: \right> \, = \, 0$. Going from the first to the second line, we have neglected the cross-terms, which are oscillating fast because of the chiral background gauge potential $A^{({\rm a})}$. Now, a crucial observation, that follows from the relation between $\tilde{h}$ and $:\psi^\dagger \psi:$ and $:\overline{\psi}^\dagger\overline{\psi}:$, is that the field $\tilde{h}_{x,y}$ is gaussian. Namely, all its connected correlators vanish: $$\left< \tilde{h}_{x_1,y_1}\tilde{h}_{x_2,y_2} \dots \tilde{h}_{x_n,y_n} \right>_{\rm conn.} \, = \, 0 \, , \quad \qquad {\rm if} \quad n > 2\, .$$ This follows immediately from the fact that the connected correlators $\left< :\psi^\dagger_{x_1,y_1} \psi_{x_1,y_1}: \dots :\psi^\dagger_{x_n,y_n} \psi_{x_n,y_n}: \right>_{\rm conn.}$ themselves vanish if $n>2$. [One way to prove this is to observe that the only singular term in the operator product expansion of $:\psi^\dagger_{x_i,y_i} \psi_{x_i,y_i}:$ with $:\psi^\dagger_{x_j,y_j} \psi_{x_j,y_j}:$ is a scalar, so its connected correlation with any other terms (present if $n>2$) is zero.] The fluctuations of the height field around its mean value $\left< h_{x,y} \right>$ are therefore the ones of a gaussian field theory with the action (\ref{eq:boson_curved}); the only thing we need to do is to compute its propagator. In the process we expect to learn why the background gauge potential $A^{({\rm v})}$ and the tetrads are absent from the action. The derivative of the propagator is, locally, \begin{eqnarray*} \partial_x \partial_{x'} \left< \tilde{h}_{x,y} \, \tilde{h}_{x',y'} \right> & \simeq & \left< \left( : \psi^\dagger_{x,y} \psi_{x,y} : + : \overline{\psi}^\dagger_{x,y} \overline{\psi}_{x,y} :\right) \left( : \psi^\dagger_{x',y'} \psi_{x',y'} : + : \overline{\psi}^\dagger_{x',y'} \overline{\psi}_{x',y'} :\right) \right> \\ &= & - \left<:\psi^\dagger_{x,y}\psi_{x',y'}:\right>\left<:\psi^\dagger_{x',y'}\psi_{x,y}:\right> \, + \, {\rm c.c.} \\ & =& - \,\frac{e^{-\sigma - i \theta - \sigma' - i \theta'}}{(z-z')^2} \, + \, {\rm c.c,} \end{eqnarray*} where recall that $z=z(x,y)$ is given by Eq. (\ref{eq:zhoney}). By locally we mean distances large compared to the lattice spacing ($|x-x'|$, $|y-y'|\gg 1$) but small compared to the full system ($|x-x'|$, $|y-y'|\ll N$). Using the fact that $e^{\sigma + i \theta} \partial_x z(x,y) \, = \, -1$, the previous equation can be easily integrated, leading to \begin{eqnarray}\label{eq:proplog} \left< \tilde{h}_{x,y} \, \tilde{h}_{x',y'} \right> \, = \, - \log |z-z'|^2 \qquad,\qquad |z-z'|\ll N \quad,\quad |z-z'|\gg 1 \end{eqnarray} up to some unimportant additive constant. This result is fully consitent with the action (\ref{eq:boson_curved}), with $h$ replaced by $\tilde{h}$. Indeed, written in the coordinate system provided by $z(x,y)$, the propagator for the action (\ref{eq:boson_curved}) is the Green function for the following operator $$\partial_z \partial_{\bar{z}} \left\langle\tilde{h}(z,\bar{z})\tilde{h}(z',\bar{z}')\right\rangle\;=\;4\pi\, \delta^{(2)}(z-z'),$$ whose solution is exactly (\ref{eq:proplog}). Note also that $\tilde{h}_{x,y}$ is clearly gauge invariant, so the background gauge potential $A^{({\rm v})}$ cannot play any role. What about the tetrad ? Since $\tilde{h}_{x,y}$ is a scalar ( i.e. spinless) field, it should not play any role either; we see that it has dropped from the propagator (\ref{eq:proplog}). In summary, the long distance correlations in the fluctuation region are fully described by a gaussian action in curved space, Eq. (\ref{eq:boson_curved}). The tetrad as well as the vector gauge potential $A_\mu^{(\textrm{v})}$ essentially disappear from the bosonic formalism, as they should. However the axial gauge potential still plays a crucial role, as it fixes the classical part of the field $\mathinner{\langle{h}\rangle}$, the average height (\ref{eq:average_height}).

### Probable connection with other works

We should mentioned here that the gaussian free field has been rigorously derived in a large variety of mathematical models of random surfaces (see [6366] for some recent references) and we hope that a connection with our work would lead to a better understanding of the universality of these models. In other works (see e.g. Refs. [19,67,18]), the arctic curve is obtained by minimizing a functional of the form $$\mathcal{A} \, = \, \int_{{\rm strip}} dx dy \, \sigma \left( \partial_x h , \partial_y h \right) \, ,$$ where $\sigma$ is the so-called surface tension. It is the free energy per unit area of the homogeneous dimer model with a fixed average slope $(\partial_x h , \partial_y h)$. In that formalism, the average height configuration is the one that minimizes the functional $\mathcal{A}$. The fluctuations around that minimal configuration are given, at the leading order, by the hessian of $\sigma$. Thus, we expect that applying this formalism to our setting (the particular dimer model above on the strip) should lead to the formulas (\ref{eq:average_height}) and (\ref{eq:proplog}) above. Unfortunately, as far as we know this result is not directly available in the literature. We thus lead this particular exercise for future work.

## Conclusion

In summary: motivated by the quantum quench studied by Antal, Rácz, Rácoz and Schütz in [49] and by the imaginary-time setup of Calabrese and Cardy [30], we studied the XX chain evolving from a domain-wall initial state (DWIS) in imaginary time. As it turns out that this problem is naturally related to arctic circle problems, we revisited a number of those from the perspective of 1d quantum systems evolving with a transfer matrix. Our main message is that the long-range correlations inside the fluctuating region are described by a massless field theory with a local action; since the system is inhomogeneous, the few parameters that enter this action must vary with position. In particular, the components of the metric themselves vary, so the field theory is, generically, a field theory in curved space-time.

All the models that we investigated in this paper could be mapped on free fermions. Obviously, this is a very strong restriction, and it is natural to ask whether our message is supposed to generalize to interacting problems. For instance, it would be very interesting to see what the action of the field theory becomes when one replaces the XX chain in the introduction by the XXZ chain (say in the critical regime, with $-1 < \Delta \leq 1$) [6870]. There, the natural conjecture is that the continuous theory is a gaussian free field (like the one in part 6.), but this time with a stiffness that depends on position as well. The reason is that, as is well-known, the density fluctuations of one-dimensional quantum systems with conserved particle number can always be described by a scalar field theory which is gaussian up to irrelevant perturbations. These systems are known as 'Luttinger liquids' (see [71] for a review), and they are all described by a free boson action, with (at least) two free parameters: the velocity $v$ of the low-energy excitations, and the Luttinger parameter $K$, related to the stiffness of the gaussian free field. The variation of the velocity with position in space-time is, in the language of this paper, equivalent to the variation of the metric: $ds^2 \sim v(x,t)^2 dt^2 - dx^2$. What is new compared to free fermion systems is the variation of $K$ with position in space-time. Importantly, since the Luttinger liquid action is only an effective theory in which the continuous free fields correspond to some renormalized version of the lattice operators, we should expect that the relation between the lattice degrees of freedom and the continuous fields should vary with positions as well, in contrast with (\ref{eq:lattice_continuous_XX}) and other similar relations in this paper. We hope to come back to these questions in future papers. Of course, there already are some known results on arctic circle problems in interacting systems, such as the six-vertex model at $\Delta \neq 0$. For instance, the partition function is known thanks to the work of Korepin and Zinn-Justin [5] and the equation that determines the arctic curve itself was conjectured by Colomo and Pronko in [8]. To our knowledge, however, not much is known about the long-range correlations inside the fluctuating region, and we hope that the present paper, together with the perspective of possible applications to quantum quenches, motivates further work in that direction.

Figure 15 Figure 15. Typical density profile in imaginary time for a free fermion chain with increased lattice periodicity ($N=512$ sites, and $R = 128$; the conventions are the same as in the introduction). Here we added a staggered chemical potential term $\sum_x (-1)^x c^\dagger_x c_x$ to the XX Hamiltonian. This model is gapped at half-filling. The density profile now contains more phases: the central region is half-filled, and is gapped. The colored region around it is critical, and the green and blue regions on the sides are again completely frozen.
Another possible extension of the present paper, not quite as ambitious as trading the XX chain for the XXZ one, is to look at free fermion models with more bands and with energy gaps. This is achieved by increasing the size of the unit cell in the model, namely by breaking the symmetry group of translations $\mathbb{Z}$ to a subgroup $n \mathbb{Z}$. For instance, one could replace the XX chain (which has a single band) by the so-called Su-Schrieffer-Heeger model (with two bands separated by a gap) [72], or by a model with a staggered chemical potential term. The ground state of either model, at half-filling, is gapped, so it has short-range correlations. On the other hand, at other filling fractions, it is gapless. So we expect that, evolving in imaginary time from the DWIS, we should see a more complex density profile in space-time, with a central region that is massive, some intermediate region that is critical, and frozen regions outside. Some of our preliminary numerical results are displayed in Fig. 15. This model is a particularly simple realization of a situation with three phases, a phenomenon already studied by [15], and also revisited recently in [73]. This phenomenon appears to be quite generic when one starts from a quantum 1d system that is gapped; this is the case for free fermion cases, as the previous references illustrate, but also when the gap is due to interactions [26] rather than to broken translation symmetry. We plan to come back to this in full detail elsewhere.

\acknowledgments

We thank Pasquale Calabrese for discussions in 2012 that motivated this paper, and Leonid Glazman for discussions about Ref. [75] which inspired the setup of part 2.. We are grateful to J. Bouttier, F. Colomo, V. Korepin, P. Krapivsky, N. Reshetikhin, G. Schütz, H. Spohn and A. Sportiello for their encouragements to publish these results, and for providing guidance through the literature on the arctic circle. We also thank M. Haque and G. Misguich for stimulating discussions about real time dynamics. We acknowledge hospitality and financial support from the GGI in Florence in 2012 (JD, JV) and in 2015 (NA, JMS, JV), and from the Visitors Program of the MPIPKS Dresden (JD, 2014). Part of this work was supported by the Mission Interdisciplinaire du CNRS (JD) and by the Conseil Régional de Lorraine (JD).

## Analytic continuation and quantum quenches

We discuss here the analytic continuation to real time of our main results for the XX chain. As emphasized in the introduction, the connection between imaginary and real-time dynamics is one of the motivations behind the present paper. We explain how this can be done, and check that one gets results consistent with previous studies [49,76,77]. This means our results can be used to determine the real time dynamics as a biproduct. We mainly focus on partition functions and correlators. It is our hope that the above connection can be used to make non trivial predictions about more complicated quantities such as particle fluctuations or entanglement entropies; this will be discussed elsewhere, however.

The present appendix also serves as an opportunity to discuss a well known [78,76,79,77] heuristic method to obtain results about the real time dynamics, called semi-classics. Here we will go a bit further that that, and try to reverse the logic. Is it possible to interpret imaginary time correlations using semiclassical arguments ? We will see that this can be done in the XX chain.

The section is organized as follows. In 8.1. we check the analytic continuation of partition functions, asymptotic forms of simple correlations such as the density profile, as well as the metric of the underlying Dirac action. We then explain the logic behind (real-time) semiclassical methods in 8.2., and extend those to determine the analogous quantities in imaginary time. Finally, we check in 8.3. that the analytic continuation may also be consistently performed at the level of exact finite-size correlators.

### Partition functions, asymptotics of the correlation functions and the metric

Let us start with the simplest quantity one can compute, the partition function $Z(R)=\mathinner{\langle{\Psi_0|e^{-2RH}|\Psi_0}\rangle}=e^{R^2/2}$ (the proof of this formula may be found in appendix. 11.). Setting $R=it/2$ in the previous formula, one recovers the correct return probability $R(t)=\left|\mathinner{\langle{\Psi(0)|\Psi(t)}\rangle}\right|=e^{-t^2/8}$ after the quench [77]. Another simple physical quantity is the density profile in the slab. Recall that it is given by the formula $$\label{eq:xxdensity} \mathinner{\langle{\rho_{x,y}}\rangle}=\frac{1}{\pi}\arccos \left(\frac{x}{\sqrt{R^2-y^2}}\right)\qquad,\qquad x^2+y^2\leq R^2,$$ inside the fluctuating region in the scaling limit. To obtain the real time density profile, it is necessary to first replace $y$ by $it$, and then take the limit $R\to 0^+$. We obtain $$\label{eq:exactdensityt} \rho_x(t)=\frac{1}{\pi}\arccos \left(\frac{x}{t}\right)\qquad,\qquad x\leq t,$$ which is the result of Ref. [49]. Here by scaling limit we mean large $t$ and $x$, but with $x/t$ finite. More generally, one can easily check that all equal (real) time correlators can be recovered from their imaginary time counterpart. Said differently, $$\left\langle c_x^\dagger(t) c_{x'}(t)\right\rangle_{\rm real}\;=\;\lim_{R\to 0^+}\,\left.\left\langle c_{x,y}^\dagger c_{x',y}\right\rangle\right|_{y=it}\,.$$ Let us finally comment on the metric. We have seen that it is given in imaginary time by \begin{eqnarray} ds^2&=&e^{2\sigma} dz d\bar{z}\\ &=&\left(dx+\frac{xy}{R^2-y^2}dy\right)^2+\frac{R^2\left(R^2-x^2-y^2\right)} {\left(R^2-y^2\right)^2}dy^2. \end{eqnarray} Now performing the analytic continuation as explained above, we obtain $$ds^2\,=\,\left(dx-v(x,t)dt\right)^2,$$ where $v(x,t)=x/t$ is the speed of propagation.

### Semi-classical interpretation

The XX chain has dispersion relation $\varepsilon(k)=-\cos k$, and group velocity $v(k)=\frac{d\varepsilon(k)}{dk}=\sin k$. However in this subsection, we will try to keep the notations $\varepsilon(k)$ and $v(k)$ for as long as possible, as it makes the physical explanations much more transparent. One particularly nice feature of the real time results is that they can be interpreted using a semiclassical picture, as pointed out in [76], and studied in detail in [77]. The key idea behind such methods is to assign to each “particle” both a well-defined position $x$ and a well-defined momentum $k\in [-\pi,\pi]$. Doing that violates a basic principle of quantum mechanics, the Heisenberg's uncertainty relations. Therefore, such an assumption can only be approximate; it is often called semiclassical (or hydrodynamic) approximation in the physics literature. To see how this can be done, let us write the occupation number of a single $k$ mode in real space. We have \begin{eqnarray} \mathinner{\langle{c^\dagger (k) c(k)}\rangle}=\sum_{j \in \mathbb{Z}} \sum_{j'\in \mathbb{Z}} e^{ik(j-j')}\mathinner{\langle{c_j^\dagger c_{j'}}\rangle}. \end{eqnarray} [In the following we use different conventions than in the main text, and label the sites by integers. The DWIS becomes $\mathinner{|{\Psi_0}\rangle}=\prod_{x\leq 0}c_x^\dagger \mathinner{|{0}\rangle}$.] The idea is to rewrite this as sums of terms centered at positions $x=\frac{j+j'}{2}$ (or $x=\frac{j+j'+1}{2}$, depending on the parity of $j+j'$). Writing $j=x+y$ and $j'=x-y$, we obtain \begin{eqnarray} \mathinner{\langle{c^\dagger (k) c(k)}\rangle}&=&\sum_{x\in \mathbb{Z}}\left(\sum_{y \in \mathbb{Z}} e^{ik 2y}\mathinner{\langle{c^\dagger_{x+y}c_{x-y}}\rangle}+e^{ik (2y+1)}\mathinner{\langle{c^\dagger_{x+y+1}c_{x-y}}\rangle}\right)\\ \label{eq:wigner} &\equiv& \sum_{x\in \mathbb{Z}} W(x,k). \end{eqnarray} The function $W(x,k)$ defined by Eq. (\ref{eq:wigner}) is a discrete analog of the Wigner function; in the following we will also use that name for the discrete case. $W(x,k)$ is, roughly speaking, the probability of finding a fermion at position $x$ and momentum $k$ in the phase space. All correlation functions can be, in principle, expressed in terms of $W$. For example, the mean density is $\mathinner{\langle{\rho_x}\rangle}=\int_{-\pi}^{\pi} \frac{dk}{2\pi}W(x,k)$. For our quench protocol we want to know how $W(x,k,t)$ evolves in time. In the domain wall initial state we have $$W(x,k,t=0)=\Theta(-x),$$ where we use the slightly unusual convention $\Theta(x\geq 0)=1$ and $\Theta(x<0)=0$, for the Heaviside theta function. This yields $\mathinner{\langle{\rho_x}\rangle}=\Theta(-x)$, as expected. We then assume that each particle at momentum $k$ propagates ballistically (classically) with group velocity $v(k)=\frac{d\varepsilon(k)}{dk}$. This means the Wigner function at time $t$ is simply given by $$W(x,k,t)=W(x-v(k)t,k,0)=\Theta(-x+v(k)t).$$ The density profile is then $\mathinner{\langle{\rho_x(t)}\rangle}=\int_{-\pi}^{\pi} \frac{dk}{2\pi}\Theta(-x+v(k)t)$, which reproduces Eq. (\ref{eq:exactdensityt}) exactly for dispersion relation $\varepsilon(k)=-\cos k$. Said differently, semiclassics is exact in the scaling limit $x,t\to \infty$ with $x/t$ fixed. For this particular quench, such an observation can be proved by expressing the correlators in terms of Bessel functions [76]. We refer also to Ref. [77] for a general discussion of these methods in free fermions systems. All this essentially follows from the stationary phase equation, which is here given by $$\label{eq:statphasereal} x-t v(k)=0,$$ and can be shown to dominate the correlation functions. The ballistic spreading of the particles essentially follows from (\ref{eq:statphasereal}) (see Fig. 16 for an illustration).

Figure 16 Figure 16. Semiclassical picture for the real time dynamics. Left: Domain wall initial state. The “probability density” (Wigner function $W(x,k)$) to find a fermion at position $x$ and momentum $k$ is uniform equal to one for $x<0$ and zero for $x>0$. Right: same picture at time $t>0$, assuming that each “fermion” at momentum $k$ moves freely with speed $v(k)=\sin k$.
Now can such (heuristic) pictures be used to determine the density profile in imaginary time, at least for cosine dispersion relation? As explained in the introduction the stationary phase equation now reads $$x+i y v(k)+R\tilde{v}(k)=0,$$ where $\tilde{v}(k)=\frac{d\tilde{\varepsilon(k)}}{dk}$ is the derivative of the Hilbert transform of the dispersion relation. We now use the following trick. We suppose for a moment that $y$ may be replaced by $it$ with $t$ real. We then assume that the particle move at speed $v(k)$, with an extra contribution proportional to $\tilde{v}(k)$ in the imaginary-time direction. Therefore, the new “Wigner function” $W_R(x,k,t)$ corresponding to the imaginary time problem satisfies $$W_R(x,k,t)=W_R(x-v(k)t-R\tilde{v}(k),k,0),$$ and the density is given by \begin{eqnarray} \mathinner{\langle{ \rho_{x,t}}\rangle} &=& \int_{-\pi}^{\pi} \frac{dk}{2\pi} W_R(x,k,t)\\ &=& \int_{-\pi}^{\pi} \frac{dk}{2\pi} \Theta \left(-x+t\sin k +R \cos k \right)\\ &=&\frac{1}{\pi}\arccos \left(\frac{x}{\sqrt{R^2+t^2}}\right). \end{eqnarray} Now plugging $t=-iy$ in the above formula, we recover the exact result (\ref{eq:xxdensity}).

### Finite-size correlation functions

We now present more general arguments about the analytic continuation. Let $\mathcal{O}_k (x_k)$ be local operators at positions $x_k\in\mathbb R$, with expectation value in time given by $$\label{QC} \langle\mathcal{O}_1(x_1,t_1)\dots\mathcal{O}_n(x_n,t_n) \rangle\equiv\lim_{R\rightarrow 0} \frac{\langle\psi_0|\mathcal{T}\bigl[(e^{iH(t_1+iR)}\mathcal{O}(x_1)e^{-iH(t_1-iR)})\dots (e^{iH(t_n+iR)} \mathcal{O}(x_n)e^{-iH(t_n-iR)})\bigr]|\psi_0\rangle} {\langle\psi_0|e^{-2RH}|\psi_0\rangle},$$ where $R>0$ is a dumping factor introduced to ensure convergence and $\mathcal{T}$ the time-ordering prescription. Explicit manipulations of \eqref{QC} are performed in imaginary time, i.e. by setting $t_k=-i y_k$, and eventually taking the limit $R\rightarrow 0$ after analytic continuation back in real time. We should stress that along the procedure analyticity is assumed a priori. Notice also that the effective geometry of the time-evolution in \eqref{QC} – before taking the limit $R\rightarrow 0$ – is the slab in shown Fig. 2. Computations in the slab geometry and analytic continuation back to real time were performed in [30] to extract the large-time behavior of correlation functions in a protocol satisfying the conditions

• the bulk of the slab is described by a gapless quantum field theory (CFT),
• the boundary state $|\psi_0\rangle$ renormalizes to a conformal invariant boundary condition.

As emphasized in the introduction, the dynamics from the domain wall initial state in the XX model provides a particularly neat example of violation of the second condition above (see also [44] for a discussion) that is also exactly solvable in the slab geometry. We start by establishing under which conditions the two-point correlation functions can be brought into a form suitable for analytic continuation to real time. We have $$\label{corryy1} \langle c^{\dagger}(x,y)c(x',y')\rangle= \int_{-\pi}^{\pi}\frac{dk}{2\pi}\int_{-\pi}^{\pi}\frac{dk'}{2\pi}\frac{e^{-ixk-y\cos k+iR\sin k+ix'k'+ y'\cos k'-iR\sin k'}}{2i\sin\left(\frac{k-k'}{2}-i0\right)}.$$ Introducing the notations $\alpha=\operatorname{arcth}(y/R)$ (resp. $\alpha'$ ) and $\tilde{R}=\sqrt{R^2-y^2}$ (resp. $\tilde{R'}$), we can rewrite \eqref{corryy1} as $$\label{Cf_tilde} \langle c^{\dagger}(x,y)c(x',y')\rangle =\int_{-\pi}^{\pi}\frac{dk}{2\pi}\int_{-\pi}^{\pi}\frac{dk'}{2\pi}\frac{e^{-ixk+i\tilde{R}\sin(k+i\alpha)+ix'k' -i\tilde{R}'\sin(k'+i\alpha')}}{2i\sin\left(\frac{k-k'}{2}-i0\right)}.$$ To show the main idea of the calculation we assume $\alpha, \alpha'\geq 0$, the other cases can be worked out in a similar fashion. First we look at the correlation functions as contour integrals in the complex plane of $z$ and $z'$ defined as follows $$\label{Cf_an} \langle c^{\dagger}(x,y)c(x',y')\rangle =e^{-\alpha x+\alpha'x'}\int_{C_{\alpha}}\frac{dz}{2\pi}\int_{C_{\alpha'}}\frac{dz'}{2\pi} \frac{e^{-ixz+i\tilde{R}\sin z+ix'z' -i\tilde{R}'\sin z'}}{2i\sin\left(\frac{z-z'-i\alpha+i\alpha'}{2}-i0\right)}.$$ The contours $C_{\alpha}$, $C_{\alpha'}$ are shifted along the imaginary axis by a constant amount $i\alpha$ and $i\alpha'$ with respect to the segment $[-\pi,\pi]$, see Fig. 17.

Figure 17 Figure 17. Left. The contour $C_{\alpha}$ is denoted in red. The contour can be lowered back to the real axis provided the gray area is free of singularity. The pole $z=k'+i\alpha+i0$ for $k'\in[-\pi,\pi]$ is also denoted. Right. Contour of integration $C_{\alpha'}$ and pole in the variable $z'=k+i(\alpha'-\alpha)-i0$.
The contributions of the vertical edges in Fig. 17 cancel due to periodicity of the integrand under $z\rightarrow z+2\pi$ and $z'\rightarrow z'+2\pi$. Therefore we can lower the contours back to the real axis provided no singularities are encountered in the process. Let us first look at the integral over $z$: The only pole with $|\textrm{Re}\, z|\leq \pi$ is located at $z_*=z'+i\alpha-i\alpha'+i0$, but remembering that $z'$ belongs to the contour $C_{\alpha'}$ we have $z_*=k'+i\alpha+i0$, with $k'\in[-\pi,\pi]$. Therefore for $\alpha\geq 0$ we can deform the contour $C_{\alpha}$ into the real interval $[-\pi,\pi]$. Consider now the integral over $z'$. The only singularity we can encounter is located at $z'_*=z-i(\alpha-\alpha')-i0$, but now the variable $z$ lies on the real axis and $z'_*=k-i(\alpha-\alpha')-i0$, for $k\in[-\pi,\pi]$. The point $z'_*$ will be outside the gray region of Fig. 17 if and only if $\alpha-\alpha'\geq 0$ since we have assumed $\alpha$ and $\alpha'$ non negative. Summarizing, in the case $y$ and $y'$ non-negative we can analytically continue the integral \eqref{Cf_an} to $$\label{Cf_fin} \langle c^{\dagger}(x,y)c(x',y')\rangle =e^{-\alpha x+\alpha'x'}\int_{-\pi}^{\pi}\frac{dk}{2\pi}\int_{-\pi}^{\pi}\frac{dk'}{2\pi} \frac{e^{-ixk+i\tilde{R}\sin k+ix'k' -i\tilde{R}'\sin k'}}{2i\sin\left(\frac{k-k'-i\alpha+i\alpha'}{2}-i0\right)}.$$ if and only if $y'\leq y$. A case where \eqref{Cf_fin} can be obtained in closed form is $y=y'$, one gets $$\label{samey} \langle c^{\dagger}(x,y)c(x',y)\rangle=\frac{\tilde{R} e^{-\alpha(r-s)}}{2(r-s)} [J_{r-1}(\tilde{R})J_{s}(\tilde{R})-J_{r}(\tilde{R})J_{s-1}(\tilde{R})]$$ where $x=r-1/2$, $x'=s-1/2$ for $r,s\in\mathbb Z$ and $J_{\nu}(t)=\int_{-\pi}^{-\pi}\frac{dk}{2\pi} e^{-i\nu x+i\sin k t}$ is the Bessel function, in agreement with [51]. One can prove \eqref{samey} by taking the derivative with respect to $\tilde{R}$ in \eqref{Cf_fin}. Now we are allowed to analytically continue directly \eqref{samey} for $y=it$ and set $\tilde{R}$ to zero; in this limit $\alpha=i\pi/2$ and we are back with exactly the same correlation functions of the real time quench problem first found in [49].

## Airy kernel, Tracy-Widom distribution and Airy process

The goal of this appendix is to provide a concise, non-rigorous, derivation of the boundary scaling behavior of the propagator, and its link with the celebrated Tracy-Widom distribution [80]. This behavior at the edge is obtained in a scaling regime that is different from the one in which one sees the bulk correlations—recall (\ref{eq:scaling_regime})—therefore it is a topic that is not really the focus of this paper. However, once the fermion propagator is calculated, it is really a minor effort to get to the Tracy-Widom distribution. We therefore discuss these aspects now, for completeness. We also briefly review the appearance of the Airy process introduced by Prähofer and Spohn in [50].

### Fermion propagator close to the boundary: generic appearance of the Airy kernel

Once again, we want to analyze the behavior of some integral using the stationary phase approximation. To illustrate the general features, consider the following warm-up exercise: let us determine the behavior of $$I(x,y)=\int_{-\pi}^{\pi}\frac{dk}{2\pi} e^{i(kx+R\tilde{\varepsilon}(k)+iy\varepsilon(k))}$$ for a point that is near the arctic curve. [To fix the ideas, one can focus on the XX chain of the introduction, for which $\varepsilon(k) = -\cos k$ and $\tilde{\varepsilon}(k) = -\sin k$; then the arctic curve is the circle $x^2 + y^2 = R^2$. We will however try to keep the notations more general, in order to emphasize that the discussion of this appendix applies more generally, and in particular it can be adapted straightforwardly to all the models treated in the main text, by chosing the appropriate dispersion $\varepsilon(k)$, or $\varepsilon(\kappa)$.] Let us parametrize the arctic curve locally as $(x,y) = (x_{\rm a.c}(y), y)$. Then when we say 'near the arctic curve', we mean, more precisely, that we are interested in the following limit: $$\label{eq:edge_scaling} \frac{y}{R} \quad {\rm fixed} , \qquad \quad \frac{x-x_{{\rm a.c.}}(y)}{R^{1/3}} \quad {\rm fixed} , \qquad \qquad {\rm and} \quad R \rightarrow \infty .$$ The reason why the power $R^{1/3}$ shows up will become clear below. First, notice that such points are 'near the boundary' in the sense that the distance from the point $(\frac{x}{R},\frac{y}{R})$ to the arctic curve goes to zero, so in the bulk scaling regime (\ref{eq:scaling_regime}), these are really points on the boundary. Now, the crucial difference with bulk scaling is that, in the limit (\ref{eq:edge_scaling}), the two solutions to the stationary phase equation $x + R \frac{d}{dk} \tilde{\varepsilon}+iy \frac{d}{dk}\varepsilon$, which we called $k=z$ and $k=-z^*$ throughout the paper, are not distinct; instead, the two stationary points coalesce at a degenerate point $z_d(y) = z = -z^*$. Around this critical point, the argument in the exponential vanishes not only at first order, but also at second order. Therefore, one needs to expand it to third order [81], leading to \begin{eqnarray} \label{Airy_analytic} \nonumber I(x,y) & \stackrel{\text{near the b.}}{\longrightarrow} & e^{i z_d(y) (x-x_{{\rm a.c.}})+ i\varphi(z_d(y) )} \int_{-\infty}^{\infty}\frac{d\tilde{k}}{2\pi} e^{i\left[\tilde{k}\left(x-x_{\rm a.c.}(y) \right)+ \frac{\tilde{k}^3 C(y)} {3}\right]} \\ && = \; e^{i z_d(y) (x-x_{{\rm a.c.}})+ i\varphi(z_d(y) )}~ C(y)^{-1/3} \text{Ai}(X) \, , \end{eqnarray} where $\tilde{k} \equiv k-z_d(y)$ and $C(y) \equiv \frac{1}{2} \frac{d^3}{dk^3} \left( R\tilde{\varepsilon}+iy \varepsilon \right)_{|_{k = z_d(y)}}$. Notice that $C(y)$ must be of order $O(R)$. [One can check that $C(y)$ is real and non-zero, for all the models we looked at in the main text. Moreover, we can assume without loss of generality that $C(y) > 0$ (if $C(y) <0$, all formulas are identical up to a change $C(y) \rightarrow -C(y)$ and $X \rightarrow -X$). For instance, for the XX chain, $C(y)=\sqrt{R^2-y^2}$.] In the first line, we extended the integral to the real line, since the integrand is oscillating fast at large $\tilde{k}$ and therefore gives only subleading contributions. In going from the first to the second line, we recognized the Airy function, $\operatorname{Ai}(X)=\int_{\mathbb R}\frac{dK}{2\pi} e^{iKX+iK^3/3}$, and used the rescaled variables $X= C(y)^{-1/3}(x-x_{\rm a.c.}(y))$ and $K = C(y)^{1/3} \tilde{k}$.

After this instructive exercise, let us turn to the task of determining the behavior of the propagator for two points near the boundary. Recall that the propagator takes the form $$\label{cf_exact} \langle c^{\dagger}_{x,y} c_{x',y'}\rangle=\int_{-\pi}^{\pi} \int_{-\pi}^{\pi}\frac{dk~dq}{(2\pi)^2}~ \frac{e^{-i(kx+R\tilde{\varepsilon}(k)+iy\varepsilon(k))+i(q x'+R\tilde{\varepsilon}(q) +iy'\varepsilon(q))}}{2i\sin\left(\frac{k-q}{2}-i0\right)}.$$ Since we want the two points $(x,y)$ and $(x',y')$ to be near the arctic curve, the stationary points in $k$ and $q$ are the degenerate critical points above $z_d(y)$ and $z_d(y')$. When $y\not=y'$, $z_d(y) \neq z_d(y')$, so the denominator in \eqref{cf_exact} does not vanish, and the asymptotics of the two-point functions immediately follows from the result of the previous exercise. It is given by the product of Airy functions $$\label{Airy_1} \langle c^{\dagger}_{x,y}c_{x',y'}\rangle_{|_{y\not = y'}} \stackrel{\text{near the b.}}{\longrightarrow} \frac{[I(x,y)]^* I(x',y')}{2i\sin\left(\frac{z_c(y)-z_c(y')}{2}\right)}.$$ The case $y\rightarrow y'$ is more interesting. Although it could be derived directly from \eqref{Airy_1}, taking the limit $y' \rightarrow y$, we prefer to follow a more direct path. When $y=y'$ we can expand the denominator in \eqref{cf_exact} to the first order both in $k$ and $q$, close to $z_c(y)$. This gives $$\langle c^{\dagger}_{x,y}c_{x',y}\rangle\stackrel{\text{near the b.}} {\longrightarrow}e^{-i z_d(y) (x-x')} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{d\tilde{k}d\tilde{q}}{(2\pi)^2} ~\frac{e^{-i\left[\tilde{k}\left(x-x_{{\rm a.c.}}(y)\right)+ \frac{\tilde{k}^3C(y)}{3}\right]+i\left[\tilde{q}\left(x'-x_{{\rm a.c.}}(y)\right)+ \frac{\tilde{q}^3 C(y)}{3}\right]}}{i(\tilde{k}-\tilde{q}-i0)} \, ,$$ where $\tilde{k}=k-z_d(y)$ and $\tilde{q}=q-z_d(y)$. Finally, with the same rescaled variables $X$ and $K$ (and $X'$ and $Q$) as above, we obtain $$\label{airy_stat} \langle c^{\dagger}_{x,y}c_{x',y}\rangle\stackrel{\text{near the b.}} {\longrightarrow } e^{-i z_d(y) (x-x')} C(y)^{-1/3} \, K(X, X') \, ,$$ where $K(X,X')$ is the Airy kernel [80] $$\label{Airy_k} K(X,X')= \int_{0}^{\infty} d\lambda \operatorname{Ai}(X+\lambda)\operatorname{Ai}(X'+\lambda).$$ The result \eqref{airy_stat} can be easily proved by observing that $[i(\tilde{k}-\tilde{q}-i0)]^{-1}=\int_{0}^{\infty}d\lambda e^{-i\lambda(\tilde{k}-\tilde{q}-i0)}$ and recalling the definition of the Airy function given above. In passing, notice that one equivalent and somewhat more explicit way of writing the Airy kernel is $K(X,X')=[\operatorname{Ai}(X)\operatorname{Ai}'(X')-\operatorname{Ai}(X')\operatorname{Ai}'(X)]/(X-X')$. For later purposes though, the integral form \eqref{Airy_k} is slightly more convenient.

In summary, the appearance of the Airy kernel close to the arctic curve is very generic: it comes from a degenerate critical point when one performs the stationary phase approximation [81]. Once the propagator is expressed in the form of the Airy kernel (\ref{airy_stat}), the Tracy-Widom distribution follows straightforwardly, as we review briefly now.

### Tracy-Widom distribution

To understand the emergence of the Tracy-Widom distribution and the Airy process we adapt the original arguments of [51] (see also the more recent Ref. [82]). Let us consider the interval of the $x$-axis $E_s$ defined by the points $x$ such that $s\leq X< \infty$, with $X$ the scaling edge variable and $y$ fixed ($y^2\not=1$). We then ask the question: what is the generating function $\chi(\lambda,s)$ for the particle number $N(E_s)$ in the region $E_s$ when $x$ is near the boundary (in the sense discussed in the previous paragraph) ? Formally, $$\chi(\lambda,s)\,=\,\mathinner{\langle{\Psi_0|e^{i\lambda N(E_s)}|\Psi_0}\rangle}\,=\,\sum_{n=0}^{\infty}e^{i\lambda n}~\mathbb P[N(E_s)=n],$$ where $\mathbb P[N(E_s)=n]$ is the probability that exactly $n$ particles are in the interval $E_s$. By standard free fermion manipulations, one can recast the function $\chi(\lambda,s)$ in the form of a determinant, $\det \left[ \left( \delta_{x,x'} + (e^{i \lambda} -1) \left< c^\dagger_x c_{x'} \right> \right) \right]$, where $x$ and $x'$ are integers both restricted to the semi-infinite interval where $X, X' \in E_s$. Then, replacing the lattice propagator $\left< c^\dagger_x c_{x'} \right>$ by its asymptotic form (\ref{airy_stat}), and considering $x$, $x'$ (or equivalently $X$, $X'$) as continuous variables, $\chi (\lambda,s)$ becomes the Fredholm determinant $$\chi(\lambda,s)=\det[1+(e^{i\lambda}-1)K(X,X')],$$ where it is understood that $X, X'$ both belong to the interval $E_s$. To see this, notice that the phase in \eqref{airy_stat} may be removed by a gauge transformation $c^\dagger_x \rightarrow e^{i z_d x} c^\dagger_x$, which leaves the determinant invariant. Also, notice that the Jacobian $\frac{\partial x}{\partial X}= C(y)^{1/3}$ exactly compensates the prefactor in \eqref{airy_stat}. It follows that all probabilities $\mathbb P[N(E_s)=n]$ can be expressed as $$\mathbb P[N(E_s)=n]=\left.\frac{(-1)^n}{n!}\frac{d^n}{du^n}\det(1-uK)\right|_{u=1}.$$ Now let us define $F_2(X)dX$ as the probability that the rightmost particle is in the interval $[X,X+dX]$. Then the cumulative distribution of $F_2(X)$ is $\mathbb P[N(E_s)=0]=\int_{-\infty}^s d X~F_2(X)$, which is equal to $\det(1-K)$, where $1-K$ is an operator on the interval $E_s$ with $s=X$, as above. This is the gaussian unitary Tracy-Widom distribution.

### Airy process

The last step is to introduce the stochastic process [50] associated with the motion of the right-most particle. Since we are dealing with free fermion problems, we know already that all multipoint correlation functions may be obtained from the fermion correlators $\langle c^{\dagger}_{x_1,y_1}\dots c^{\dagger}_{x_n,y_n}c_{x_1',y_1'}\dots c_{x_n',y_n'} \rangle$, which themselves are obtained from Wick's theorem. The only relevant question is: what is the propagator? In mathematical language, this means that the process one is trying to define must be a determinantal point process. The question is: what is its kernel ?

As usual in this paper, the propagator (or kernel) is given by (\ref{cf_exact}), and we need to analyze its scaling behavior. Here, the scaling behavior we want to focus on is, following [50], $$\label{eq:airy_process_scaling} \frac{y_0}{R} \; {\rm fixed} , \qquad \frac{x-x_{{\rm a.c.}}(y_0)}{R^{1/3}} \quad {\rm fixed} , \qquad \frac{x'-x_{{\rm a.c.}}(y_0)}{R^{1/3}} \quad {\rm fixed}, \qquad \frac{y-y_{0}}{R^{2/3}} \quad {\rm fixed} , \qquad \frac{y'-y_{0}}{R^{2/3}} \quad {\rm fixed} , \qquad {\rm and} \; R \rightarrow \infty .$$ That is, we zoom in on a region around the point $(x_{{\rm a.c.}}(y_0), y_0)$ on the arctic curve. The reason why the powers $R^{1/3}$ and $R^{2/3}$ show up is again related to the way we perform the stationary phase approximation. This time, in (\ref{cf_exact}), we are going to integrate on $k$ around the stationary point $z_d (y)$, which is itself close to $z_d(y_0)$, so we write \begin{eqnarray*} k -z_d(y_0) & = & (k-z_d(y)) + (z_d(y)-z_d(y_0)) \, \simeq \, \frac{1}{C(y_0)^{1/3}} \left( K - i \tau \right) \, , \end{eqnarray*} with $K = (k-z_d(y)) C(y_0)^{1/3}$ and $\tau \, = \, i (y-y_0) \, C(y_0)^{1/3} \, \left(\left.\frac{dz_d(y)}{dy}\right|_{y=y_0}\right)$. [Notice that $\tau$ is of order $O(1)$ in the regime (\ref{eq:airy_process_scaling}); one can also check that $\tau$ is always real.] The same relation holds for the rescaled variables $Q$ and $\tau'$ in terms of $q$ and $y'$. Then, expanding the argument of the exponential in (\ref{cf_exact}) as previously, one finds that, \begin{eqnarray} \label{extended_airy} \nonumber \langle c^{\dagger}_{x,y}c_{x',y'}\rangle & \stackrel{\text{near the b.}}{\longrightarrow} & e^{-i z_d(y) (x-x_{{\rm a.c.}}(y))- i\varphi(z_d(y) )} e^{i z_d(y') (x'-x_{{\rm a.c.}}(y'))+ i\varphi(z_d(y') )} \, C(y_0)^{-1/3} \\ && \quad \times \quad \int_{-\infty}^{\infty}\frac{d K}{2\pi}\int_{-\infty}^{\infty}\frac{dK'}{2\pi} \frac{e^{-i K X-i\frac{K^3}{3}+iX'K'+i\frac{K'^3}{3}}}{i(K-i\tau-K'+i\tau'-i0)}, \end{eqnarray} where $X= C(y_0)^{-1/3} (x-x_{\rm a.c.} (y))$, and the same relation holds for $X'$, $\tau'$ in terms of $x'$, $y'$. The phases may be removed by a gauge transformation as above, and one is left with $\langle c^{\dagger}_{x,y}c_{x',y'}\rangle\rightarrow C(y_0)^{-1/3} K_{Ai}(X,\tau; X', \tau')$, where $K_{Ai}$ is the extended Airy kernel, $$\label{ex_Airy} K_{Ai}(X,\tau; X' \tau') \; = \; \left\{ \begin{array}{lcl} \displaystyle \int_{0}^{\infty}d\lambda \ e^{-\lambda(\tau-\tau')}\operatorname{Ai}(X+\lambda)\operatorname{Ai}(X'+\lambda) && {\rm if} \quad \tau\geq\tau' \\ \displaystyle -\int_{\infty}^{0}d\lambda~e^{-\lambda(\tau-\tau')} \operatorname{Ai}(X+\lambda)\operatorname{Ai}(X'+\lambda) && {\rm if} \quad \tau < \tau' \, . \end{array} \right.$$ The determinantal point process defined by this extended Airy kernel is called the Airy process [50,51,17].

## Useful Hilbert transforms

In this appendix we compute the Hilbert transforms needed in the main text. We always work on the space of periodic, real, analytic functions on $\mathbb{R}$, with period $2\pi$, and with zero average: $\int_{-\pi}^\pi dk f(k) \, = \, 0$. On this space, the Hilbert transform $\tilde{f}$ of the function $f$ may be defined in (at least) three equivalent ways:

1. as the principal value of the integral $$\tilde{f}(k) = {\rm p.v.} \int_{-\pi}^\pi \frac{dk'}{2\pi} f(k') {\rm cot} \left( \frac{k-k'}{2} \right)$$
2. using the Fourier series of $f(k)$ to define the one of $\tilde{f}(k)$: $$f(k) = \sum_{p \geq 1} \left[ a_p \cos (p k) + b_p \sin (p k) \right] \Rightarrow \tilde{f}(k) = \sum_{p \geq 1} \left[ a_p \sin (p k) - b_p \cos (p k) \right]$$
3. as the unique periodic, real-valued function on $\mathbb{R}$, which is such that there exists an analytic continuation of $k \mapsto f(k) + i \tilde{f}(k)$ to a complex neighborhood of $\mathbb{R}$ that includes the upper half-plane ${\rm Im}\, k >0$.

We are free to pick any of those three equivalent definitions when we perform concrete calculations. While the Hilbert transform appeared in part 2. in the first form (the principal value), the two other definitions are also very useful. For example the Hilbert transform in the toy model immediately follows from the second definition.

All other models may be treated by making use of the third definition. The basic ingredient is the complex function (for some real parameter $u$ with $0<u<1$): $$k \mapsto \log (1+ u e^{i k}) .$$ This function is analytic in the upper half-plane ${\rm Im}\, k >0$. Therefore, according to the third definition above, if we look at the pair of real-valued functions $f$ and $\tilde{f}$ defined for $k \in \mathbb{R}$ as $$f(k) = {\rm Re} \left[ \log (1+ u e^{i k}) \right] \tilde{f}(k) = {\rm Im} \left[ \log (1+ u e^{i k}) \right] ,$$ then $\tilde{f}$ is the Hilbert transform of $f$. This is all we need for the dimer model on the honeycomb lattice, since the dispersion relation in this model turns out to be exactly $\varepsilon_{{\rm honey.}}(k) \, =\, -f(k)$. In the six-vertex and square dimer models, the situation is similar (we use the notations of parts 4., 5., and write $\kappa$ instead of $k$ for the argument of the functions): \begin{eqnarray*} 2\varepsilon_{{\rm 6v.}} (\kappa) & = & - \log \left( \frac{c + b \cos \kappa} {c - b \cos \kappa} \right) \, = \, -2\, {\rm Re} \left[ \log ( 1 + u e^{i \kappa} ) \right]\, + \,2 \,{\rm Re} \left[ \log ( 1 - u e^{i \kappa} ) \right] \, , \end{eqnarray*} with $\frac{b}{c} = \frac{2u}{1+u^2}$. Then the Hilbert transform is \begin{eqnarray*} 2\tilde{\varepsilon}_{{\rm 6v.}} (\kappa) & = & -2\, {\rm Im} \left[ \log ( 1 +u e^{i \kappa} ) \right]\, + \,2 \,{\rm Im} \left[ \log ( 1 - u e^{i \kappa} ) \right] \\ &=& i \, \log \left(\frac{ 1+ i \frac{ 2u}{1-u^2} \sin \kappa} { 1 - i \frac{2u}{1-u^2} \sin \kappa} \right)\\ & =& i \, \log \left(\frac{ a+ i b \sin \kappa} { a - i b \sin \kappa} \right) \, , \end{eqnarray*} as claimed in the main text. In the last line, we have used the assumptions of part 5., namely $a^2+b^2=c^2$, and $a>0$, $b>0$, $c>0$. The Hilbert transform for the dimer model is exactly the same, provided one makes the substitutions $a=v$, $b=u$, and $c=w$.

## Partition functions from bosonization

In this appendix, we recover the partition functions of all the models studied in this paper, using the bosonization procedure presented in section 2.. We start by introducing the general formalism of bosonization and applying it to our toy model, i.e. the XX chain. In the rest of the appendix, we shall see how this procedure can be modified to compute the partition function of the honeycomb and square lattice dimer model as well as the 6-vertex model.

### Bosonization and the XX chain

The idea is to think for a moment of the real time analog of the partition function, the return probability $\tilde{Z}(t)=\mathinner{\langle{\Psi_0|e^{-i H t}|\Psi_0}\rangle}$. In bosonized form, the following substitution holds $$e^{-i t H} \mathinner{|{\Psi_0}\rangle} \; \rightarrow \; \mathcal{N}: e^{-i t \int \frac{dk}{2\pi} \varepsilon (k) \partial \varphi (k)} : \mathinner{|{\Psi_0}\rangle} \,,$$ up to some unknown prefactor $\mathcal{N}$. It is obtained by looking at the norm of that state, which is one due to the unitarity of the real time evolution: \begin{eqnarray*} \mathinner{\langle{\Psi_0}|} e^{i t H} e^{-i t H} \mathinner{|{\Psi_0}\rangle} &=&1\\ &=&\mathcal{N}^2 \mathinner{\langle{\Psi_0}|} : e^{i t \int \frac{dq'}{2\pi} \varepsilon (q') \partial \varphi (q'- i \epsilon)} : : e^{-i t \int \frac{dq}{2\pi} \varepsilon (q) \partial \varphi (q)} : \mathinner{|{\Psi_0}\rangle} \\ &=&\mathcal{N}^2\; \exp \left( t^2 \int \frac{dq}{2\pi} \int \frac{dq'}{2\pi} \varepsilon (q) \varepsilon (q') \partial_q \partial_{q'} \mathinner{\langle{\Psi_0}|} \varphi (q' -i \epsilon)\varphi (q) \mathinner{|{\Psi_0}\rangle} \right)\,, \end{eqnarray*} By definition the expectation value of a single time-ordered exponential is one, so the return probability is given by $Z(t)=\mathcal{N}$. Therefore, we have \begin{equation*} \tilde{Z}(t)=\exp \left( -\frac{t^2}{2} \int \frac{dq}{2\pi} \int \frac{dq'}{2\pi} \varepsilon (q) \varepsilon (q') \partial_q \partial_{q'} \mathinner{\langle{\Psi_0}|} \varphi (q' -i \epsilon)\varphi (q) \mathinner{|{\Psi_0}\rangle} \right). \end{equation*} The partition function we are looking for, $Z(R)=\mathinner{\langle{\Psi_0 |e^{-2RH}|\Psi_0}\rangle}$, may be obtained by Wick-rotating back $t\to 2iR$: $$\label{eq:prescription} Z(R)=\tilde{Z}(2iR)=\exp \left( 2R^2 \int \frac{dq}{2\pi} \int \frac{dq'}{2\pi} \varepsilon (q) \varepsilon (q') \partial_q \partial_{q'} \mathinner{\langle{\Psi_0}|} \varphi (q' -i \epsilon)\varphi (q) \mathinner{|{\Psi_0}\rangle} \right) .$$ As we shall see, all the partition functions may be obtained by applying the formula (\ref{eq:prescription}). Before doing so, let us however mention that this formula can be put on a more rigorous footing using the so-called infinite wedge formalism (see [21,83]). We explain the method on the toy model Hamiltonian \begin{equation*} H=\int_{-\pi}^{\pi} \frac{dk}{2\pi}\varepsilon(k) c^\dagger(k)c(k) \end{equation*} for some dispersion relation $\varepsilon(k)$ that has zero mean value. What we really have in mind is the particular case $\varepsilon(k)=-\cos k$, but the following argument is general. $\varepsilon(k)$ may be written in Fourier series: $$\varepsilon(k)=\sum_{p\geq 1}\, [\varepsilon]_{p} e^{ipk}+\sum_{p\geq 1}\, [\varepsilon]_{-p} e^{-ipk}={\varepsilon}_r(k)+{\varepsilon}_l(k),$$ where we have separated the positive (right) Fourier modes from the negative (left) ones. Similarly, we also write $H={H}_r +{H}_l$, and make the observation that ${H}_l\mathinner{|{\Psi_0}\rangle}=0$ and $\mathinner{\langle{\Psi_0}|} {H}_r=0$. As is also explained in Ref. [21], the result is a priori sentitive to the choice regularization. Here we choose to consider a system of large even size $L$ with open boundary conditions, with sites running from $-L/2+1/2$ to $L/2-1/2$, and only then take the limit $L\to \infty$. This is clearly the most natural choice for the XX chain; for the dimer and vertex models, it will actually be imposed.

Let us specify a cosine dispersion relation for simplicity. The idea now is to apply the Baker-Campbell-Hausdorff (BCH) formula. The commutator between ${H}_l$ and ${H}_r$ is given by $$\label{eq:somecommutator} [{H}_r,{H}_l]=\frac{1}{4}\sum_{x=-L/2+1/2}^{L/2-3/2}\left(c_{x+1}^\dagger c_{x+1}-c_{x}^\dagger c_x\right).$$ The crucial point is that $\mathinner{|{\Psi_0}\rangle}$ is an eigenstate of the commutator with eigenvalue $-1/4$, independent on $L$. One can also check that all higher order nested commutators up to order $L/2$ have expectation value zero in the domain wall initial state. Therefore, if $L$ is sent to infinity first, the commutator (\ref{eq:somecommutator}) may be considered as a scalar when taking expectation values in the DWIS. Using the BCH formula we obtain \begin{eqnarray} Z(R)&=&\mathinner{\langle{\Psi_0 |e^{-2R {H}_r} e^{-2R^2[{H}_r,{H}_l]} e^{-2R {H}_l}|\Psi_0}\rangle}\\ &=&e^{-2 R^2 [{H}_r,{H}_l]}\\ &=&e^{R^2/2}. \end{eqnarray} The above can be reformulated in bosonic language. To do that, we introduce a set of operators $a_n$, defined as $$a_n= \sum_{x\in \mathbb{Z}+1/2} c^\dagger_{x-n} c_x$$ These operators satisfy the Heisenberg commutation relations $$\label{eq:comm} [a_n,a_m]=n \delta_{n,-m},$$ in the sense described above, that is, when taking expectation values in the DWIS. The $a_n$ are nothing but the modes of the bosonic field $\varphi$. Using (\ref{eq:comm}), computing the commutator between ${H}_r$ and ${H}_l$ trivialises, and we find $$\label{eq:genpartitionfunction} Z(R)=\exp\left(2R^2 \sum_{p\geq 1} p [\varepsilon]_p [\varepsilon]_{-p}\right).$$ The equality between (\ref{eq:genpartitionfunction}) and (\ref{eq:prescription}) can easily be shown by expanding $\varepsilon(k)$ in Fourier series. With this at hand, we are able to compute all the partition functions of interest in the present paper. The only complication will be that the mean energy $\mathinner{\langle{\Psi_0|H|\Psi_0}\rangle}$ will not necessary be zero.

### Dimers on the honeycomb lattice

We discuss now dimers on the honeycomb lattice, and regularize by looking at a finite system of size $L$ finite and even, and then take the limit $L\to \infty$. The diagonalization of the transfer matrix is similar to the text, the only difference being that the momenta are now discrete. For simplicity we label the sites from $1$ to $L$ here, so that the DWIS is $\mathinner{|{\Psi_0}\rangle}=\prod_{j=1}^{L/2} c_j^\dagger \mathinner{|{0}\rangle}$. We still have $T=e^{-H}$, with $$H=-\frac{1}{2}\sum_{m=1}^L \log \left[1+u^2+2u \cos \left(\frac{m}{L+1}\right) \right] c_{m}^\dagger c_m \qquad,\qquad c_m^\dagger=\sqrt{\frac{2}{L+1}}\sum_{j=1}^{L} \sin\left( \frac{m\pi j}{L+1} \right)c_{j}^\dagger$$ Let us now show that the expectation value of $H$ in the DWIS is not zero. A simple calculation shows $\mathinner{\langle{c_m^\dagger c_m}\rangle}=1/2$ for all $m$. Using this and expanding the logarithm in Fourier series, we obtain \begin{eqnarray}\nonumber \mathinner{\langle{H}\rangle}&=& \sum_{p=1}^{\infty} \frac{(-u)^p}{2p} \sum_{m=1}^{L}\cos \left(\frac{p m\pi}{L+1}\right)\\ \nonumber &=& \sum_{p=1}^{\infty} \frac{(-u)^p}{2p} \left[-\frac{1+(-1)^p}{2}+(L+1)\delta_{p,2L+2}\right]\\ \label{eq:energy0} &=&\frac{1}{4}\log (1-u^2)-\frac{1}{4}\log \left(1-u^{2L+2}\right) \end{eqnarray} Since $u<1$, the second term can be safely ignored, but not the first one, which is finite. We now use the following trick. We write the partition function as \begin{equation*} Z=\mathinner{\langle{T^{N}}\rangle}=\mathinner{\langle{e^{-2N H}}\rangle}=e^{-2N \mathinner{\langle{H}\rangle}}\mathinner{\langle{e^{-2N \left(H-\mathinner{\langle{H}\rangle}\right)}}\rangle} \end{equation*} and bosonize $H-\mathinner{\langle{H}\rangle}$ – which has average exactly zero – according to (\ref{eq:prescription}). Combining (\ref{eq:energy0}) with (\ref{eq:genpartitionfunction}), we obtain \begin{eqnarray}\nonumber Z&=& e^{-\frac{N}{2}\log (1-u^2)}\times e^{-\frac{2N^2}{4}\log (1-u^2)}\\ &=&\left(\frac{1}{1-u^2}\right)^{\frac{N}{2}(N+1)} \end{eqnarray} This result is exact in the limit $N \to \infty$ and in agreement with the results of [51] for the partition function of the multilayer PNG droplet model.

### Dimers on the square lattice and 6-vertex model

We apply exactly the same method to the dimer model. For convenience, we set the vertical weights to $v=1$; this keeps the partition function finite. The Hamiltonian for a finite size system is $$H=-\sum_{m=1}^{L} {\rm arcsinh}\left(u \cos \frac{m\pi}{L+1}\right) c_{+,m}^\dagger c_{+,m}.$$ Note that we have used the fact that the two bands are the analytic continuations of each other. Expanding again in Fourier modes, we find $$\mathinner{\langle{H}\rangle}=-\frac{1}{4}\log \left(1+u^2\right).$$ Using the same two band bosonization method explained in Section. 2. and 4., we find that the result (\ref{eq:genpartitionfunction}) can still be used, albeit with the Fourier coefficients in terms of the $\kappa$ variable. The Fourier coefficients of the dispersion relation $\varepsilon(\kappa)=-\frac{1}{2}\log\left(\frac{w+u\cos \kappa}{w-u \cos \kappa}\right)$ are given by $[\varepsilon]_{2p+1}=(\frac{w-v}{u})^{2p+1}/(2p+1)$ and $[\varepsilon]_{2p}=0$ otherwise. We finally obtain $\sum_{p\geq 1}p[\varepsilon]_p [\varepsilon]_{-p}={\rm arctanh}\left(\left[\frac{w-v}{u}\right]^2\right)=\frac{1}{2}\log(w/v)= \frac{1}{4}\log(1+u^2)$. Hence the partition function reads \begin{eqnarray}\nonumber Z&=&e^{\frac{2N}{4}\log (1+u^2)}e^{\frac{2N^2}{4}\log(1+u^2)}\\ &=&(1+u^2)^{\frac{N}{2}(N+1)} \end{eqnarray} which is the celebrated Aztec diamond partition function [2]. We switch now to the 6-vertex model with domain-wall boundary conditions and start by setting $a=1$, which makes the partition function finite. The calculation is very similar to that of the dimer model, the only real difference being that $\mathinner{\langle{H}\rangle}=0$. The contribution coming from the bosonized part can be deduced from that of the dimers, we simply need to replace $u$ by $b$ and $w$ by $c$. We recover the well known result [10,9,20] $$Z_{\rm dwbc}=c^{N^2}.$$

The purpose of this appendix is to provide some additional technical details regarding the derivation of the coefficients $r_{x,y}(k)$ and $s_{x,y}(k)$ (\ref{eq:ident_square}) and (\ref{eq:ident_6v}) needed to obtain exact finite-size correlations for models with two bands. We treat separately the dimer model on the square lattice and the six vertex model with domain wall boundary conditions.

### Dimers on the square lattice.

Here we derive the result \begin{eqnarray}\label{eq:odd_action_2bis} &&T_{\rm o}\ c_{2j+1/2}\ T_{\rm o}^{-1}=\frac{1}{v}\, c_{2j+1/2},\nonumber \\ &&T_{\rm o}\ c_{2j-1/2}\ T_{\rm o}^{-1}=\frac{1}{v^2}\left(-u c_{2j-3/2}+ v c_{2j-1/2}-u c_{2j+1/2}\right), \end{eqnarray} for the action of the odd transfer matrix $T_{\rm o}$ on the annihilation operators. The easiest is to go to momentum space. Recall the action on the creation operators is given by $$T_{\rm o} \left(\begin{array}{c}c^\dagger(k)\\ c^\dagger(k+\pi)\end{array}\right) T_{\rm o}^{-1}=M_{\rm o}(k)\qquad,\qquad M_{\rm o}(k)= \left(\begin{array}{cc} v+u\cos k & iu \cos k \\ iu \cos k & v-u\cos k \end{array}\right)\left(\begin{array}{c}c^\dagger(k)\\ c^\dagger(k+\pi)\end{array}\right).$$ Now diagonalizing the matrix $M_{\rm o}(k)$, one can define new creation operators $\tilde{c}^\dagger_+(k), \tilde{c}_+(k)$, so that the transfer matrix reads $T_{\rm o}=\exp\left(-\int_{-\pi}^{\pi}\frac{dk}{2\pi}\tilde{\varepsilon}(k)\tilde{c}^\dagger_+(k) \tilde{c}_+(k)\right)$. Hence $T_{\rm o} \tilde{c}_+^\dagger (k)T_{\rm o}^{-1}=e^{-\tilde{\varepsilon}(k)}\tilde{c}_+^\dagger(k)$ and also $T_{\rm o} \tilde{c}_+ (k)T_{\rm o}^{-1}=e^{\tilde{\varepsilon}(k)}\tilde{c}_+(k)$. Said differently, $T_{\rm o}$ acts on the creation operators as a two by two matrix $M_{\rm o}(k)$, and acts on the annihilation operators as the inverse of $M_{\rm o}(k)$: $$\label{eq:someinversestuff} T_{\rm o} \left(\begin{array}{cc} c(k)&c(k+\pi) \end{array} \right) T_{\rm o}^{-1}=\left(\begin{array}{cc} c(k)&c(k+\pi) \end{array} \right) M_{\rm o}^{-1}(k)\qquad,\qquad M_{\rm o}^{-1}(k)= \frac{1}{v^2}\left(\begin{array}{cc} v-u\cos k & -iu \cos k \\ -iu \cos k & v+u\cos k \end{array}\right)$$ Going back to real space, the equation (\ref{eq:someinversestuff}) becomes (\ref{eq:odd_action_2bis}), which was used in the main text.

### Six-vertex model.

We now focus on the six vertex model, for which the determination of the coefficient $r_{x,y}(k)$ and $s_{x,y}(k)$ is slightly more involved. The propagator of the modes $c^\dagger_+(k)$, $c_+(k')$ is $$\label{eq:prop_pm_6v_a} \left< c^\dagger_{+,y}(k) c_{+,y'}(k') \right> \, = \, e^{y \varepsilon (\kappa) - i N \tilde{\varepsilon}(\kappa)} e^{-y' \varepsilon (\kappa') + i N \tilde{\varepsilon}(\kappa')} \mathinner{\langle{\Psi_0}|} c^\dagger_+(k) c_+(k') \mathinner{|{\Psi_0}\rangle} \, ,$$ again with $\kappa = \kappa(k)$ and $\kappa' = \kappa(k')$. Below, we use once again the change of basis between the modes $c^\dagger (k)$, $c^\dagger (k+\pi)$, and $c^\dagger_+ (k)$, $c^\dagger_- (k)$, given by $U(k)$. It will then be convenient to use (\ref{eq:prop_pm_6v_a}) in the following matrix form \begin{eqnarray} \label{eq:prop_pm_6v_b} \nonumber && \left( \begin{array}{cc} \left< c^\dagger_{+,y} (k) c_{+,y'} (k') \right> & \left< c^\dagger_{+,y} (k) c_{-,y'} (k') \right> \\ \\ \left< c^\dagger_{-,y} (k) c_{+,y'} (k') \right> & \left< c^\dagger_{-,y} (k) c_{-,y'} (k') \right> \end{array} \right) \, = \, \left( \begin{array}{cc} \left< c^\dagger_{+,y} (k) c_{+,y'} (k') \right> & \left< c^\dagger_{+,y} (k) c_{+,y'} (k'+\pi) \right> \\ \left< c^\dagger_{+,y} (k+\pi) c_{+,y'} (k') \right> & \left< c^\dagger_{+,y} (k+\pi) c_{+,y'+\pi} (k') \right> \end{array} \right) \\ && = \, \mathinner{\langle{\Psi_0}|} \left( \begin{array}{c} e^{E(\kappa)} c^\dagger_+(k) \\ e^{E(\kappa+\pi)} c^\dagger_+(k+\pi) \end{array} \right) \left( \begin{array}{cc} e^{-E(\kappa')} c_+(k') & e^{-E(\kappa'+\pi)} c_+(k'+\pi) \end{array} \right) \mathinner{|{\Psi_0}\rangle} \end{eqnarray} where $E(\kappa) \, \equiv \, y \varepsilon (\kappa) - i N \tilde{\varepsilon}(\kappa)$; notice that we have used $c^\dagger_-(k) = c^\dagger_+(k+\pi)$.

Thus, focusing on the fermion creation/annihilation operators, it is clear that one needs to be able to take into account the conjugation by $T_{\rm e}^{\pm 1/2}$. A little calculation shows that $$\begin{array}{ccc} T_{\rm e}^{1/2} c^\dagger(k) T_{\rm e}^{-1/2} \, = \, v^\dagger(k) \left( \begin{array}{c} c^\dagger(k) \\ c^\dagger(k+\pi) \end{array} \right) & & T_{\rm e}^{-1/2} c^\dagger(k) T_{\rm e}^{1/2} \, = \, v^\dagger(k+\pi) \left( \begin{array}{c} c^\dagger(k) \\ c^\dagger(k+\pi) \end{array} \right) \\ T_{\rm e}^{1/2} c(k) T_{\rm e}^{-1/2} \, = \, \left( \begin{array}{cc} c(k) & c(k+\pi) \end{array} \right) v(k+\pi) & & T_{\rm e}^{-1/2} c(k) T_{\rm e}^{1/2} \, = \, \left( \begin{array}{cc} c(k) & c(k+\pi) \end{array} \right) v(k) \\ \\ {\rm with} \quad v(k) \, = \, \left( \begin{array}{c} \frac{a+b+c}{2 \sqrt{a(b+c)}} +\frac{c+b-a}{2 \sqrt{a(b+c)}} \cos k \\ \frac{c+b-a}{2 \sqrt{a(b+c)}} \sin k \end{array} \right) \, . \end{array}$$ Let us focus on an example: the case $y \in 2\mathbb{Z} + \frac{1}{2}$, $y' \in 2\mathbb{Z} + \frac{1}{2}$ with $y \geq y'$. The propagator in $k$-space is \begin{eqnarray} \nonumber \left< c^\dagger_{y} (k) c^\dagger_{y'} (k') \right> & = & \frac{ \mathinner{\langle{\Psi_0}|} T_{\rm e}^{1/2} T^{2(\frac{N-y-1/2}{2})} \, T_{\rm e}^{1/2} c^\dagger(k) T_{\rm e}^{-1/2} \, T^{2(\frac{y-y'}{2})} \, T_{\rm e}^{1/2} c(k') T_{\rm e}^{-1/2} \, T^{2(\frac{N+y'+1/2}{2})} T_{\rm e}^{1/2} \mathinner{|{\Psi_0}\rangle} } {\mathinner{\langle{\Psi_0}|} T_{\rm e}^{1/2} T^{2N} T_{\rm e}^{1/2} \mathinner{|{\Psi_0}\rangle}} \\ \nonumber &=& v^\dagger (k) U^\dagger(k) \left( \begin{array}{cc} \left< c_{+,y}^\dagger(k) c_{+,y'}(k') \right> & \left< c_{+,y}^\dagger(k) c_{-,y'}(k') \right> \\ \\ \left< c_{-,y}^\dagger(k) c_{+,y'}(k') \right> & \left< c_{-,y}^\dagger(k) c_{-,y'}(k') \right> \end{array} \right) U(k') v(k'+\pi) \\ \nonumber &=& \mathinner{\langle{\Psi_0}|} v^\dagger (k) U^\dagger(k) \left( \begin{array}{c} e^{E(\kappa)} c^\dagger_+(k) \\ e^{E(\kappa+\pi)} c^\dagger_+(k+\pi) \end{array} \right) \left( \begin{array}{cc} e^{-E(\kappa')} c_+(k') & e^{-E(\kappa'+\pi)} c_+(k'+\pi) \end{array} \right)U(k') v(k'+\pi) \mathinner{|{\Psi_0}\rangle} \, . \end{eqnarray} In the second line, we used the change of basis given by $U(k)$ in Eq. (\ref{eq:diag_6v}), and in the third line we used (\ref{eq:prop_pm_6v_b}), with $E(\kappa) \equiv (y+\frac{1}{2}) \varepsilon(\kappa) - i N \tilde{\varepsilon}(\kappa)$ (the shift $+1/2$ comes from the power of the transfer matrix, which involves $y+\frac{1}{2}$ instead of $y$). Note also that $U(k+\pi)=U(k)$ for the six vertex model. We then have the following 'identity'. Whenever the fields $c^\dagger_y(k)$ and $c_y(k)$ appear in a bracket $\left< . \right>$, we may replace it by the combination $$\left\{ \begin{array}{rcl} c^\dagger_y (k) & \rightarrow & v^\dagger (k) U^\dagger(k) \left( \begin{array}{c} e^{E(\kappa)} c^\dagger_+(k) \\ e^{E(\kappa+\pi)} c^\dagger_+(k+\pi) \end{array} \right) \\ \\ c_y (k') & \rightarrow & \left( \begin{array}{cc} e^{-E(\kappa')} c_+(k') & e^{-E(\kappa'+\pi)} c_+(k'+\pi) \end{array} \right)U(k') v(k'+\pi) \end{array} \right.$$ and evaluate the r.h.s in the DWIS $\mathinner{|{\Psi_0}\rangle}$. When we come back to real-space, the 'identity' can be simplified, \begin{eqnarray} \nonumber (y \in 2 \mathbb{Z}+ \frac{1}{2}) \qquad c^\dagger_{x,y} & \rightarrow & \int_{-\pi}^\pi \frac{dk}{2\pi} e^{-i k x} v^\dagger (k) U^\dagger(k) \left( \begin{array}{c} e^{E(\kappa)} c^\dagger_+(k) \\ e^{E(\kappa+\pi)} c^\dagger_+(k+\pi) \end{array} \right) \\ \nonumber &&= \, \int_{-\pi}^\pi \frac{dk}{2\pi} e^{- i k x} \, v^\dagger(k) U^\dagger (k) \left( \begin{array}{c} e^{E(k)} c^\dagger_+(k) \\ 0 \end{array} \right) + \int_{-\pi}^\pi \frac{dk}{2\pi} e^{- i k x} \, v^\dagger(k) U^\dagger (k) \left( \begin{array}{c} 0 \\ e^{E(k+\pi)} c^\dagger_+(k+\pi) \end{array} \right) \\ \nonumber &&= \, \int_{-\pi}^\pi \frac{dk}{2\pi} e^{- i k x} \, v^\dagger(k) U^\dagger (k) \left( \begin{array}{c} e^{E(k)} c^\dagger_+(k) \\ 0 \end{array} \right) + \int_{-\pi}^\pi \frac{dk}{2\pi} e^{- i (k-\pi) x} \, v^\dagger(k-\pi) U^\dagger (k) \left( \begin{array}{c} 0 \\ e^{E(k)} c^\dagger_+(k) \end{array} \right) \\ &&= \, \int_{-\pi}^\pi \frac{dk}{2\pi} e^{- i k x} \, \left[ v^\dagger(k) U^\dagger (k) \left( \begin{array}{c} 1 \\ 0 \end{array} \right) + v^\dagger(k-\pi) U^\dagger (k) \left( \begin{array}{c} 0 \\ e^{i \pi x} \end{array} \right) \right] e^{E(\kappa)} c^\dagger_+(k) \, . \end{eqnarray} It is convenient to distinguish the cases $x \in 2 \mathbb{Z}+ \frac{1}{2}$ and $x \in 2 \mathbb{Z}- \frac{1}{2}$, for which $e^{i \pi x} = \pm i$. All the possible cases are summarized in table III.

 $r_{x,y}(k)$ $y \in 2\mathbb{Z} + \frac{1}{2}$ $y \in 2\mathbb{Z} - \frac{1}{2}$ $x \in 2\mathbb{Z} + \frac{1}{2}$ $e^{\frac{\varepsilon (\kappa)}{2}} \left[ v^\dagger(k) U^\dagger (k) \left( \begin{array}{c} 1 \\ 0 \end{array} \right) + v^\dagger(k-\pi) U^\dagger (k) \left( \begin{array}{c} 0 \\ e^{i \pi x} \end{array} \right) \right]$ $e^{-\frac{\varepsilon (\kappa)}{2}} \left[ v^\dagger(k+\pi) U^\dagger (k) \left( \begin{array}{c} 1 \\ 0 \end{array} \right) + v^\dagger(k) U^\dagger (k) \left( \begin{array}{c} 0 \\ e^{i \pi x} \end{array} \right) \right]$ $\displaystyle = \, e^{\frac{\varepsilon (\kappa)}{2}} \left[ \frac{a+b+c}{2 \sqrt{a(b+c)}} e^{ -i \frac{k-\kappa}{2}} + \frac{c+b-a}{2 \sqrt{a(b+c)}} e^{-i \frac{k+\kappa}{2}} \right]$ $\displaystyle = \, e^{-\frac{\varepsilon (\kappa)}{2}} \left[ \frac{a+b+c}{2 \sqrt{a(b+c)}} e^{ -i \frac{k-\kappa}{2}} - \frac{c+b-a}{2 \sqrt{a(b+c)}} e^{-i \frac{k+\kappa}{2}} \right]$ $x \in 2\mathbb{Z} - \frac{1}{2}$ $e^{\frac{\varepsilon (\kappa)}{2}} \left[ v^\dagger(k) U^\dagger (k) \left( \begin{array}{c} 1 \\ 0 \end{array} \right) + v^\dagger(k-\pi) U^\dagger (k) \left( \begin{array}{c} 0 \\ e^{i \pi x} \end{array} \right) \right]$ $e^{-\frac{\varepsilon (\kappa)}{2}} \left[ v^\dagger(k+\pi) U^\dagger (k) \left( \begin{array}{c} 1 \\ 0 \end{array} \right) + v^\dagger(k) U^\dagger (k) \left( \begin{array}{c} 0 \\ e^{i \pi x} \end{array} \right) \right]$ $\displaystyle = \,e^{\frac{\varepsilon (\kappa)}{2}} \left[ \frac{a+b+c}{2 \sqrt{a(b+c)}} e^{ i \frac{k-\kappa}{2}} + \frac{c+b-a}{2 \sqrt{a(b+c)}} e^{i \frac{k+\kappa}{2}} \right]$ $\displaystyle = \, e^{-\frac{\varepsilon (\kappa)}{2}} \left[ \frac{a+b+c}{2 \sqrt{a(b+c)}} e^{ i \frac{k-\kappa}{2}} - \frac{c+b-a}{2 \sqrt{a(b+c)}} e^{i \frac{k+\kappa}{2}} \right]$
 $s_{x,y}(k)$ $y \in 2\mathbb{Z} + \frac{1}{2}$ $y \in 2\mathbb{Z} - \frac{1}{2}$ $x \in 2\mathbb{Z} + \frac{1}{2}$ $e^{-\frac{\varepsilon (\kappa)}{2}} \left[ \left( \begin{array}{cc} 1 & 0 \end{array} \right) U(k) v(k+\pi) + \left( \begin{array}{cc} 0 & e^{-i \pi x} \end{array} \right) U(k) v(k) \right]$ $e^{\frac{\varepsilon (\kappa)}{2}} \left[ \left( \begin{array}{cc} 1 & 0 \end{array} \right) U(k) v(k) + \left( \begin{array}{cc} 0 & e^{-i \pi x} \end{array} \right) U(k) v(k-\pi) \right]$ $\displaystyle = \, e^{-\frac{\varepsilon (\kappa)}{2}} \left[ \frac{a+b+c}{2 \sqrt{a(b+c)}} e^{ i \frac{k-\kappa}{2}} - \frac{c+b-a}{2 \sqrt{a(b+c)}} e^{i \frac{k+\kappa}{2}} \right]$ $\displaystyle = \, e^{\frac{\varepsilon (\kappa)}{2}} \left[ \frac{a+b+c}{2 \sqrt{a(b+c)}} e^{ i \frac{k-\kappa}{2}} + \frac{c+b-a}{2 \sqrt{a(b+c)}} e^{i \frac{k+\kappa}{2}} \right]$ $x \in 2\mathbb{Z} - \frac{1}{2}$ $e^{-\frac{\varepsilon (\kappa)}{2}} \left[ \left( \begin{array}{cc} 1 & 0 \end{array} \right) U(k) v(k+\pi) + \left( \begin{array}{cc} 0 & e^{-i \pi x} \end{array} \right) U(k) v(k) \right]$ $e^{\frac{\varepsilon (\kappa)}{2}} \left[ \left( \begin{array}{cc} 1 & 0 \end{array} \right) U(k) v(k) + \left( \begin{array}{cc} 0 & e^{-i \pi x} \end{array} \right) U(k) v(k-\pi) \right]$ $\displaystyle = \,e^{-\frac{\varepsilon (\kappa)}{2}} \left[ \frac{a+b+c}{2 \sqrt{a(b+c)}} e^{ -i \frac{k-\kappa}{2}} - \frac{c+b-a}{2 \sqrt{a(b+c)}} e^{-i \frac{k+\kappa}{2}} \right]$ $\displaystyle = \, e^{\frac{\varepsilon (\kappa)}{2}} \left[ \frac{a+b+c}{2 \sqrt{a(b+c)}} e^{- i \frac{k-\kappa}{2}} + \frac{c+b-a}{2 \sqrt{a(b+c)}} e^{-i \frac{k+\kappa}{2}} \right]$
Table III. Functions $r_{x,y}(k)$ and $s_{x,y}(k)$ appearing in Eq. (\ref{eq:ident_6v}).
It turns out the formulas can be simplified even further. Let us demonstrate how so for $r_{x,y}$ in the case $x,y\in 2\mathbb{Z}+1/2$ (all other cases will be slight variants of this). First, we rewrite $r_{x,y}(k)$ as \begin{equation*} r_{x,y}(k)=\frac{1}{\sqrt{a}} e^{\frac{\varepsilon(\kappa)-ik}{2}}\left[\cos \frac{\kappa}{2}+i \frac{a}{\sqrt{b+c}}\sin \frac{\kappa}{2}\right], \end{equation*} and then compute its square. Using $e^{\varepsilon(\kappa)-ik}=\frac{a\cos \kappa-i c \sin \kappa}{c+b\cos \kappa}$, we obtain \begin{eqnarray*} \left[r_{x,y}(k)\right]^2&=&\frac{a\cos \kappa-i c \sin \kappa}{c+b\cos \kappa} \frac{\left[(b+c)\cos \frac{\kappa}{2}+i a \sin \frac{\kappa}{2}\right]^2}{a(b+c)}\\ &=&\frac{a\cos \kappa-i c \sin \kappa}{c+b\cos \kappa} \frac{b+c\cos \kappa+i a \sin \kappa}{a}\\ &=&1-i\frac{b}{a}\sin \kappa \end{eqnarray*} It is then easy to see that the correct square root is $$r_{x,y}(k)=\sqrt{1-i\frac{b}{a}\sin \kappa}.$$ Repeating the same procedure on all possible cases, we find in the end the simple result $$r_{x,y}(k)=s_{x,y}(k)=\sqrt{1-i \sigma \frac{b}{a}\sin \kappa}\qquad,\qquad \sigma=(-1)^{x-y}.$$

## References

1. William Jockusch, James Propp, and Peter Shor. Random domino tilings and the arctic circle theorem. arXiv preprint math/9801068, 1998.
2. Noam Elkies, Greg Kuperberg, Michael Larsen, and James Propp. Alternating-sign matrices and domino tilings (part i). Journal of Algebraic Combinatorics, 1(2):111–132, 1992.
3. Noam Elkies, Greg Kuperberg, Michael Larsen, and James Propp. Alternating-sign matrices and domino tilings (part ii). Journal of Algebraic Combinatorics, 1(3):219–234, 1992.
4. Henry Cohn, Michael Larsen, and James Propp. The shape of a typical boxed plane partition. New York J. Math, 4(137):165, 1998.
5. Vladimir Korepin and Paul Zinn-Justin. Thermodynamic limit of the six-vertex model with domain wall boundary conditions. Journal of Physics A: Mathematical and General, 33(40):7053, 2000.
6. KONSTANTIN Palamarchuk and NICOLAI Reshetikhin. The 6-vertex model with fixed boundary conditions. arXiv preprint arXiv:1010.5011, 2010.
7. F Colomo and AG Pronko. The arctic circle revisited. Contemporary Mathematics, 458:361, 2008.
8. Filippo Colomo and AG Pronko. The arctic curve of the domain-wall six-vertex model. Journal of Statistical Physics, 138(4-5):662–700, 2010.
9. Vladimir E Korepin. Calculation of norms of bethe wave functions. Communications in Mathematical Physics, 86(3):391–418, 1982.
10. Anatoli G Izergin. Partition function of the six-vertex model in a finite volume. In Soviet Physics Doklady, volume 32, page 878, 1987.
11. Anatoli G Izergin, David A Coker, and Vladimir E Korepin. Determinant formula for the six-vertex model. Journal of Physics A: Mathematical and General, 25(16):4315, 1992.
12. Anatoly M Vershik and Sergei V Kerov. Asymptotics of plancherel measure of symmetrical group and limit form of young tables. Doklady Akademii Nauk SSSR, 233(6):1024–1027, 1977.
13. Benjamin F Logan and Larry A Shepp. A variational problem for random young tableaux. Advances in mathematics, 26(2):206–222, 1977.
14. Craig Rottman and Michael Wortis. Statistical mechanics of equilibrium crystal shapes: Interfacial phase diagrams and phase transitions. Physics Reports, 103(1):59–79, 1984.
15. Bernard Nienhuis, Hendrik Jan Hilhorst, and HWJ Blöte. Triangular sos models and cubic-crystal shapes. Journal of Physics A: Mathematical and General, 17(18):3559, 1984.
16. Henry Cohn, Noam Elkies, James Propp, et al. {Local statistics for random domino tilings of the Aztec diamond}. Duke Mathematical Journal, 85(1):117–166, 1996.
17. Kurt Johansson. The arctic circle boundary and the airy process. Annals of probability, pages 1–30, 2005.
18. Richard Kenyon. Lectures on dimers. arXiv preprint arXiv:0910.3129, 2009.
19. Richard Kenyon, Andrei Okounkov, and Scott Sheffield. Dimers and amoebae. Annals of mathematics, pages 1019–1056, 2006.
20. NM Bogoliubov, AG Pronko, and MB Zvonarev. Boundary correlation functions of the six-vertex model. Journal of Physics A: Mathematical and General, 35(27):5525, 2002.
21. Andrei Okounkov and Nikolai Reshetikhin. Correlation function of schur process with application to local geometry of a random 3-dimensional young diagram. Journal of the American Mathematical Society, 16(3):pp. 581–603, 2003.
22. Andrei Okounkov and Nicolai Reshetikhin. Random skew plane partitions and the pearcey process. Communications in mathematical physics, 269(3):571–609, 2007.
23. F Colomo and AG Pronko. On two-point boundary correlations in the six-vertex model with domain wall boundary conditions. Journal of Statistical Mechanics: Theory and Experiment, 2005(05):P05010, 2005.
24. Fabien Alet, Jesper Lykke Jacobsen, Grégoire Misguich, Vincent Pasquier, Frédéric Mila, and Matthias Troyer. Interacting classical dimers on the square lattice. Physical Review Letters, 94(23):235702, 2005.
25. Olav F. Syljuå sen and M. B. Zvonarev. Directed-loop monte carlo simulations of vertex models. Phys. Rev. E, 70:016118, Jul 2004.
26. David Allison and Nicolai Reshetikhin. Numerical study of the 6-vertex model with domain wall boundary conditions. In Annales de l'institut Fourier, volume 55, pages 1847–1869, 2005.
27. Leticia F Cugliandolo, Giuseppe Gonnella, and Alessandro Pelizzola. Six vertex model with domain-wall boundary conditions in the bethe-peierls approximation. arXiv preprint arXiv:1501.00883, 2014.
28. Anthony Perret, Zoran Ristivojevic, Pierre Le Doussal, Grégory Schehr, and Kay J Wiese. Super-rough glassy phase of the random field x y model in two dimensions. Physical review letters, 109(15):157205, 2012.
29. Pasquale Calabrese and John Cardy. Evolution of entanglement entropy in one-dimensional systems. Journal of Statistical Mechanics: Theory and Experiment, 2005(04):P04010, 2005.
30. Pasquale Calabrese and John Cardy. Time dependence of correlation functions following a quantum quench. Physical review letters, 96(13):136801, 2006.
31. Pasquale Calabrese and John Cardy. Quantum quenches in extended systems. Journal of Statistical Mechanics: Theory and Experiment, 2007(06):P06008, 2007.
32. Pasquale Calabrese and John Cardy. Entanglement and correlation functions following a local quench: a conformal field theory approach. Journal of Statistical Mechanics: Theory and Experiment, 2007(10):P10004, 2007.
33. Spyros Sotiriadis and John Cardy. Inhomogeneous quantum quenches. Journal of Statistical Mechanics: Theory and Experiment, 2008(11):P11003, 2008.
34. Spyros Sotiriadis, Pasquale Calabrese, and John Cardy. Quantum quench from a thermal initial state. EPL (Europhysics Letters), 87(2):20002, 2009.
35. Spyros Sotiriadis and John Cardy. Quantum quench in interacting field theory: a self-consistent approximation. Physical Review B, 81(13):134305, 2010.
36. Andrea Gambassi and Pasquale Calabrese. Quantum quenches as classical critical films. EPL (Europhysics Letters), 95(6):66007, 2011.
37. Jérôme Dubail and Jean-Marie Stéphan. Universal behavior of a bipartite fidelity at quantum criticality. Journal of Statistical Mechanics: Theory and Experiment, 2011(03):L03002, 2011.
38. Jean-Marie Stéphan and Jérôme Dubail. Local quantum quenches in critical one-dimensional systems: entanglement, the loschmidt echo, and light-cone effects. Journal of Statistical Mechanics: Theory and Experiment, 2011(08):P08019, 2011.
39. Curtis T Asplund and Steven G Avery. Evolution of entanglement entropy in the d1-d5 brane system. Physical Review D, 84(12):124053, 2011.
40. Jean-Marie Stéphan and Jérôme Dubail. Logarithmic corrections to the free energy from sharp corners with angle 2Ï. Journal of Statistical Mechanics: Theory and Experiment, 2013(09):P09002, 2013.
41. Curtis T Asplund and Alice Bernamonti. Mutual information after a local quench in conformal field theory. Physical Review D, 89(6):066015, 2014.
42. DM Kennes, V Meden, and R Vasseur. Universal quench dynamics of interacting quantum impurity systems. Physical Review B, 90(11):115101, 2014.
43. Kevin Kuns and Donald Marolf. Non-thermal behavior in conformal boundary states. Journal of High Energy Physics, 2014(9):1–19, 2014.
44. John Cardy. Quantum quenches to a critical point in one dimension: some further results. arXiv preprint arXiv:1507.07266, 2015.
45. Dennis Bernard and Benjamin Doyon. A hydrodynamic approach to non-equilibrium conformal field theories. arXiv preprint arXiv:1507.07474, 2015.
46. Jean-Marie Stéphan. Emptiness formation probability, toeplitz determinants, and conformal field theory. Journal of Statistical Mechanics: Theory and Experiment, 2014(5):P05010, 2014.
47. Pasquale Calabrese, Christian Hagendorf, and Pierre Le Doussal. Time evolution of one-dimensional gapless models from a domain wall initial state: stochastic loewner evolution continued? Journal of Statistical Mechanics: Theory and Experiment, 2008(07):P07013, 2008.
48. Gesualdo Delfino and Jacopo Viti. Phase separation and interface structure in two dimensions from field theory. Journal of Statistical Mechanics: Theory and Experiment, 2012(10):P10009, 2012.
49. T Antal, Z Rácz, A Rákos, and GM Schütz. Transport in the xx chain at zero temperature: Emergence of flat magnetization profiles. Physical Review E, 59(5):4912, 1999.
50. Michael Prähofer and Herbert Spohn. Universal distributions for growth processes in 1+ 1 dimensions and random matrices. Physical review letters, 84(21):4882, 2000.
51. Michael Prähofer and Herbert Spohn. Scale invariance of the png droplet and the airy process. Journal of statistical physics, 108(5-6):1071–1106, 2002.
52. Herbert Spohn. Exact solutions for kpz-type growth processes, random matrices, and equilibrium shapes of crystals. Physica A: Statistical Mechanics and its Applications, 369(1):71–99, 2006.
53. Mikio Nakahara. Geometry, topology and physics. CRC Press, 2003.
54. Elliott H Lieb. Solution of the dimer problem by the transfer matrix method. Journal of Mathematical Physics, 8(12):2339–2341, 1967.
55. Fabien Alet, Yacine Ikhlef, Jesper Lykke Jacobsen, Grégoire Misguich, and Vincent Pasquier. Classical dimers with aligning interactions on the square lattice. Physical Review E, 74(4):041124, 2006.
56. Jean-Marie Stéphan, Shunsuke Furukawa, Grégoire Misguich, and Vincent Pasquier. Shannon and entanglement entropies of one-and two-dimensional critical wave functions. Physical Review B, 80(18):184421, 2009.
57. Nicolas Allegra. Exact solution of the 2d dimer model: Corner free energy, correlation functions and combinatorics. Nuclear Physics B, 894:685–732, 2015.
58. Cédric Boutillier, Jérémie Bouttier, Guillaume Chapuy, Sylvie Corteel, and Sanjay Ramassamy. Dimers on rail yard graphs. arXiv preprint arXiv:1504.05176, 2015.
59. P Zinn-Justin. Six-vertex model with domain wall boundary conditions and one-matrix model. Physical Review E, 62(3):3411, 2000.
60. Scott Sheffield. Gaussian free fields for mathematicians. Probability theory and related fields, 139(3-4):521–541, 2007.
61. Berbard Nienhuis. Coulomb gas formulations of two-dimensional phase transitions. In Phase Transitions and Critical Phenomena, volume 11. Academic press, 1987.
62. HWJ Blöte and HJ Hilhorst. Roughening transitions and the zero-temperature triangular ising antiferromagnet. Journal of Physics A: Mathematical and General, 15(11):L631, 1982.
63. Alexei Borodin and Patrik L Ferrari. Anisotropic growth of random surfaces in 2+ 1 dimensions. Citeseer, 2008.
64. Leonid Petrov et al. Asymptotics of uniformly random lozenge tilings of polygons. gaussian free field. The Annals of Probability, 43(1):1–43, 2015.
65. Jeffrey Kuan. The gaussian free field in interlacing particle systems. Electron. J. Probab, 19(72):1–31, 2014.
66. Maurice Duits. On global fluctuations for non-colliding processes. arXiv preprint arXiv:1510.08248, 2015.
67. Richard Kenyon, Andrei Okounkov, and Scott Sheffield. Dimers and amoebae. Annals of mathematics, pages 1019–1056, 2006.
68. Jarrett Lancaster and Aditi Mitra. Quantum quenches in an x x z spin chain from a spatially inhomogeneous initial state. Physical Review E, 81(6):061134, 2010.
69. Jorn Mossel and Jean-Sébastien Caux. Relaxation dynamics in the gapped xxz spin-1/2 chain. New Journal of Physics, 12(5):055028, 2010.
70. Thiago Sabetta and Grégoire Misguich. Nonequilibrium steady states in the quantum xxz spin chain. Physical Review B, 88(24):245114, 2013.
71. David Senechal. An introduction to bosonization. In Theoretical Methods for Strongly Correlated Electrons, pages 139–186. Springer, 2004.
72. W_P Su, JR Schrieffer, and Ao J Heeger. Solitons in polyacetylene. Physical Review Letters, 42(25):1698, 1979.
73. Sunil Chhita and Kurt Johansson. Domino statistics of the two-periodic aztec diamond. arXiv preprint arXiv:1410.2385, 2014.
74. Philippe Di Francesco and Rodrigo Soto-Garrido. Arctic curves of the octahedron equation. Journal of Physics A: Mathematical and Theoretical, 47(28):285204, 2014.
75. Eldad Bettelheim and Leonid Glazman. Quantum ripples over a semiclassical shock. Physical review letters, 109(26):260602, 2012.
76. T. Antal, P. L. Krapivsky, and A. Rákos. Logarithmic current fluctuations in nonequilibrium quantum spin chains. Phys. Rev. E, 78:061115, Dec 2008.
77. Jacopo Viti, Jean-Marie Stéphan, Jerome Dubail, and Masud Haque. Inhomogeneous quenches in a fermionic chain: exact results. arXiv:1507.08132, 2015.
78. A. G. Abanov. Hydrodynamics of correlated systems. Emptiness Formation Probability and Random Matrices. Nato Science Series II: Mathematics, Physics and Chemistry, vol. 221 Springer, 2005.
79. Eldad Bettelheim and Paul B. Wiegmann. Universal fermi distribution of semiclassical nonequilibrium fermi states. Phys. Rev. B, 84:085102, Aug 2011.
80. Craig A. Tracy and Harold Widom. Level-spacing distributions and the airy kernel. Comm. Math. Phys., 159(1):151–174, 1994.
81. Arthur Erdélyi. Asymptotic Expansions. Dover Publications, 1956.
82. Viktor Eisler and Zoltán Rácz. Full counting statistics in a propagating quantum front and random matrix spectra. Phys. Rev. Lett., 110:060602, Feb 2013.
83. Victor G. Kac. Infinite dimensional Lie algebras. Cambridge University Press, 1994.
Filters:

## {{childPeer.user.name}}

{{getCommentDate(childPeer.date) | amDateFormat:'MMM Do YYYY, HH:mm'}}

## {{childPeer.user.name}}

{{getCommentDate(childPeer.date) | amDateFormat:'MMM Do YYYY, HH:mm'}}