# Fluctuation of the phase boundary in the six-vertex model with Domain Wall
Boundary Conditions: a Monte Carlo study

Ivar Lyberg<EMAIL_ADDRESS>PhD Stony Brook University  Vladimir
Korepin  C.N. Yang Institute for Theoretical Physics, Stony Brook University,
Stony Brook, New York, 11794-3840, USA  Jacopo Viti<EMAIL_ADDRESS>(corresponding author) INFN, sezione di Firenze, via G. Sansone 1, 50019 Sesto
Fiorentino, Italy

###### Abstract

We consider the six-vertex model with Domain Wall Boundary Conditions in a
square lattice of dimension $N\times N$. Our main interest is the study of the
fluctuations of the extremal lattice path about the arctic curves. We address
the problem through Monte Carlo simulations. At $\Delta=0$, the fluctuations
of the extremal path along any line parallel to the square diagonal were
rigorously proven to follow the Tracy-Widom distribution. We provide strong
numerical evidence that this is true also for other values of the anisotropy
parameter $\Delta$ ($0\leq\Delta<1$). We argue that the typical width of the
fluctuations of the extremal path about the arctic curves scales as $N^{1/3}$
and provide a numerical estimate for the parameters of the scaling random
variable.

##  1 Introduction

In probability theory [1], the sum of a large number of random variables once
divided by the number of summands might approach a deterministic value. For
instance, the law of large numbers implies that the sample mean of $N$
independent identically distributed random variables converges to their
distribution average for large $N$. Moreover, if those random variables have
finite variance, fluctuations of their sum are normally distributed with a
variance of order $\sqrt{N}$. The convergence toward smooth functions of
combinations of a large number of random variables is at the root of the
effectiveness of the methods of statistical mechanics and also provides a
solid ground for numerical techniques such as Monte Carlo. This phenomenon is
not restricted to local observables but extends to geometrical quantities such
as random curves. An instructive example discussed in [2] (but see also [3])
is a one-dimensional random walk starting from the origin and conditioned in
$T$ steps to reach the point $X$. For $X,~{}T\rightarrow\infty$ with $X/T=v$
finite, the random paths performed by the walker are with probability one
within any neighborhood of the classical trajectory $X=vT$. Limit shapes [4]
are smooth nonrandom curves that describe the thermodynamic limit of random
lattice paths. Historically, one of the most studied examples is the arctic
circle [5] that appears in dimer coverings of the Aztec diamond. This domino
tiling problem is equivalent [6] to the six vertex model [7, 8] with Domain
Wall Boundary Conditions (DWBC) at its free fermionic point [9, 10, 11].

Let $N$ be a positive integer, a configuration of the six vertex model on a
$N\times N$ square lattice with DWBC [12, 13, 14, 15, 16, 17] is specified by
$N$ non-intersecting paths traveling on the edges of the lattice. The lines
enter from above and exit on the left. On each vertex, the path can have the
six possible configurations given in Fig. 1; the Boltzmann weights associated
with those are $a$ for the vertices of type $a_{1}$ and $a_{2}$, $b$ for the
vertices of type $b_{1}$ and $b_{2}$ and $c$ for the vertices of type $c_{1}$
and $c_{2}$. The phase diagram of the model is characterized by the ratio
$b/a$ and the parameter $\Delta$ [9], defined as

$\Delta:=\frac{a^{2}+b^{2}-c^{2}}{2ab}.$

In this work, we will focus on the domain $0\leq\Delta<1$, which is included
in the so-called disordered phase $|\Delta|<1$ see [18] for a comprehensive
study of the phase diagram. The case $\Delta=0$ is the free fermionic point
mentioned a few lines above.

Figure 1: The six possible configurations of a lattice path in the six-vertex
model are denoted by thick lines. The continuous lines are lattice edges while
dashed lines are edges of the dual lattice.

The $N$ random paths cannot intersect and for large enough $N$ they can
fluctuate only within a region whose boundary, after rescaling all the
distances by $N$, is arbitrarily close to a curve, dubbed the arctic curve
[5]. For $\Delta=0$ and $a=b$, such a curve is the circle
$(x/R)^{2}+(y/R)^{2}=1$, depicted in Fig. 2 and inscribed in the square
lattice. Inside the arctic curve the model is in a disordered phase, outside
is in a trivial ferromagnetic phase. The Cartesian coordinates are also
defined in Fig. 2 and lengths, such as the circle radius $R$, are measured in
suitable lattice units. In this paper, we follow the convention to divide the
square diagonal into $N$ intervals, therefore $x,y$ are semi-integers in the
interval $[-N/2,N/2]$; for simplicity, we assume $N$ a power of two from now
on.

The Cartesian equation of the arctic curve depends on the parameter $\Delta$
and was calculated in [19, 20], see also [3]. A rigorous derivation at
$\Delta=1/2$ has been given in [21].

The limit shape, however, is a property of the model in the thermodynamic
limit and it is relevant to understand how the interface between the disorder
and the ferromagnetic phase fluctuates due to finite-size effects. In this
case, the phases are not separated by a smooth curve; rather, the distance
from a corner to the disordered phase varies according to a probability
distribution. Consider, for instance, a quarter of the original lattice: the
southeast (SE) quadrant on the right of Fig. 2. When $N$ is large, there are
no paths in a neighborhood of the SE corner and all the vertices are of type
$a_{2}$. For later purposes, introduce the parameter $\eta$ as

$\eta:=2y/N.$ (1)

Take $|\eta|<1/2$; as we move along a constant-$y$ line toward the boundary
(for $x>0$), we cross $N/4$ lattice paths, see Fig. 2. Denote then by
$X_{ext}(\eta)$ the $x$-coordinate of the intersection of a constant-$y$ line
with the furthest path from the origin. Refs. [22, 23] tackled rigorously the
study of the fluctuations of the extremal path in the context of dimer
coverings. The same conclusions apply to the six-vertex with DWBC at
$\Delta=0,~{}$[6]. In particular, it can be proven that for $N\gg 1$,

$X_{ext}(\eta)\stackrel{{\scriptstyle N\gg
1}}{{\longrightarrow}}\Lambda(\eta)N\pm\Gamma(\eta)N^{1/3}\mathcal{Z},$ (2)

where $\mathcal{Z}$ follows a Tracy-Widom (TW) distribution [24] for any
$|\eta|<1/2$. The sign in front of Eq. (2), is the one of $X_{ext}(\eta)$ in
the coordinate system defined in Fig. 2. Moreover, the parameters $\Lambda$
and $\Gamma$ are $O(1)$: $\Lambda N$ is the $x$-coordinate of the intersection
of the arctic curve (an ellipse if $a\not=b$) with the constant-$y$ line while
$\Gamma$ can be derived from Ref. [25], see Sec. 3.

Figure 2: Left. Lattice path configuration in the six vertex model with DWBC
for $N=64$, and $a=b$ at $\Delta=0$. The possible vertex configurations have
been colored differently: vertices of type $a_{1}$ are colored in red, of type
$a_{2}$ in white, of type $b_{1}$ in yellow, of type $b_{2}$ in blue while of
type $c_{1}$ and $c_{2}$ in green. The black circle is the arctic curve for
the given choice of the parameters. Right. The SE corner of the lattice path
configuration on the left for $N=64$. The furthest path from the center of the
lattice intersects a constant $y$ line in a point with abscissa
$X_{ext}(\eta)$. The picture shows $X_{ext}(\eta=0)$.

The TW distribution characterizes fluctuations of extrema of random variables
in several different contexts, among them: the spectrum of Hermitian random
matrices [24], ordered sequences in random permutations [26], random growth
models [22, 27] and their generalizations [28, 31], the KPZ [32] equation
[33], quantum dynamics [34, 25, 35]. Many of these problems can be mapped into
each other [36] and one might wonder [37] how universal is the TW
distribution. Ref. [29] proved for instance that the maximum reached by the
extremal path at $\Delta=1/2$ follows a TW distribution (although with
$\beta=1$).

In this paper, we will study numerically the probability distribution of
$X_{ext}(\eta)$ for values of $\Delta\not=0$ in the interval $[0,1)$. We will
argue that after a rescaling analogous to Eq. (2), this random variable
follows, for large enough $N$, again a TW distribution. This result is not
unexpected and can be justified heuristically. Non-intersecting lattice paths
can be mapped to fermionic trajectories in a Euclidean space time [25]. The
fermionic particles are interacting if $\Delta\not=0$, however the particle
furthest from the origin, which is also the fastest, behaves as if were free
[36, 38, 39]. This observation is also at the root of the so-called tangent
method [3].

The rest of the paper is organized as follows. In Sec 2, we review briefly the
numerical method. In Section 3 we discuss the comparison with the case
$\Delta=0$ and further analyze the lattice path fluctuations for
$\Delta\not=0$ in Sec. 4. The analysis in Sec. 4 focuses on the lattice
diagonal ($\eta=0$), while some comments on results at $\eta\not=0$ are
contained in the two conclusive appendices. Finally, in Sec. 5 we summarize
our conclusions.

##  2 Monte Carlo study of the six-vertex model with DWBC

Method.—We modify a Monte Carlo (MC) algorithm for simulating the six-vertex
model with DWBC that was originally proposed in [40]. The MC algorithm
involves two moves called first flip and second flip; see Fig. 3. The moves
modify the configuration of the paths along the edges of a plaquette of the
square lattice. The state $S$ of a plaquette is a list that contains the
configurations of its four vertices: $S=\\{v_{1},v_{2},v_{3},v_{4}\\}$, with
$v_{i}$ given in Fig. 1. The Boltzmann weight of $v_{i}$ will be denoted by
$w(v_{i})=a,b,~{}c$. The weight of the state $S$ is then
$W(S)=\prod_{i=1}^{4}w(v_{i})$. An MC move takes then a state $S$ to a new
state $S^{\prime}$. Since in a flippable plaquette, each vertex can have two
possible configurations, there are sixteen distinct plaquette states which are
respectively first or second flippable. They have up to nine distinct weights.

Figure 3: Left. The two MC moves. Right. The six-vertex model with DWBC with
$N=4$. The picture also shows an allowed state. This state has four flippable
plaquettes. The plaquettes indicated in blue, are second flippable while the
plaquettes indicated in red are first flippable.

The computer code has two lists, each containing the plaquettes flippable in
one of the two different ways. We can formulate the algorithm as follows.
First, randomly search among all plaquettes on the lattice, until a plaquette
contained in one of the two lists is found. A flip will then be attempted, and
if the attempt is successful, then both lists will be updated: at the
plaquette itself and at the four plaquettes sharing a side with it. The
probability $P$ of an attempted flip that modifies the state of the plaquette
$S$ into $S^{\prime}$ is the Boltzmann weight of the new state, namely

$P=W(S^{\prime}).$ (3)

In the computer program, the largest of the three weights $a$, $b$, and $c$ is
1 and therefore Eq. (3) is well defined. The above conditions satisfy detailed
balance and ensure ergodicity. An equivalent MC algorithm with Glauber
dynamics has been employed in [41] and [42] to obtain numerical estimates for
the arctic curves and one-point functions of the model with several values of
$\Delta$ in the disorder phase. To ensure the system has thermalized we rely
then on the results of those papers. We define a MC sweep as $N^{2}$ accepted
flips. The largest value of $N$ used in this paper is 256 whereas the longest
side of the lattice used in the mentioned papers was as long as 500 units. We,
therefore, consider it known at this time how long this system needs to
equilibrate. From a rigorous perspective, Ref. [43] shows that the Glauber
dynamics is rapidly mixing at $\Delta=1/2$, while Refs. [44, 45] proved that
the Markov chain exponentially slows down for large and negative $\Delta$
111The arguments of [45] break down in the whole disorder phase and although
not rigorously proven, it is expected that Glauber dynamics is rapidly mixing
for all $|\Delta|<1$.

A reliable test of the TW distribution is challenging in a statistical model,
see for instance [46] and requires a much larger number of sweeps than that
employed previously in Refs. [41, 42]. In this paper, we will construct
samples with $10^{6}$ MC sweeps.

In the disorder phase, analogous MC algorithms are nowadays commonly adopted
in numerical experiments with the six vertex model [47, 48, 49, 50, 51, 52].
For earlier numerical approaches to ice models, see [53, 54]. Finally, an
efficient algorithm to simulate interacting dimers has been proposed in [55]
but so far it has not been implemented in the six-vertex model.

$X_{ext}(0)$$X_{ext}(\eta)$ Figure 4: Left. We show the numerical data for
distributions of the random variables $X_{ext}(\eta=0)$ and
$X_{ext}(\eta=1/4)$ along the lines $y=0$ and $y=N/8$, see also Eq. (1), at
$\Delta=1/2$ and $N=256$. Right. A realization of the random variable
$X_{ext}(\eta)$ for two different values of $\eta$ at $\Delta=1/2$ and
$N=256$.

The extremal lattice path.—In the coordinate system 222Notice that in Ref.
[25], the square diagonal is divided into $2N$ intervals. Therefore distances
measured along the diagonal in this paper are twice as big as the ones
considered in [25]. This choice also implies that the constant $\Gamma$ in Eq.
(2) is half of the one obtained analytically in [25] at $\Delta=0$, see Sec.
3. defined in Fig. 1 with origin in the center of the lattice, the random
variable $X_{ext}(\eta)$ is the position, along a line at constant $y$, of the
last occurrence of a vertex of different type from those in the frozen corner.
The right panel of Fig. 4 illustrates a realization of $X_{ext}(\eta)$ for
$\eta=0$ and $\eta=1/4$ on the NW corner of lattice with $N=256$ and
$\Delta=1/2$. Nearby this corner the vertices are in the state $a_{1}$, see
Fig. 1 and Fig. 2; the black dashed line is the arctic curve [19].

The probability distribution, $\text{Prob}(X_{ext}(\eta)=x)$, will be
approximated numerically at finite $N$ by a normalized histogram which we
denote by $h_{\eta}$; see Appendix. A for more detail about its construction.
The left panel of Fig. 4 shows two examples at $\eta=0$ and $\eta=1/4$ for
$\Delta=1/2$ and $N=256$. The two histograms are not symmetric (the right tail
is longer than the left) and obviously have different variance and mean. Their
properties will be studied in the next two sections.

##  3 Extremal lattice path fluctuations at $\Delta=0$ for
$|\eta|<\frac{1}{2}$

Special attention must be paid to the case $\Delta=0$, since it is the only
point at which there are known exact results for the random process that
describes the fluctuations of the coordinate $X_{ext}(\eta)$. Let us take for
simplicity $a=b$; in [23, 22], the random fluctuations were proven for large
$N$ to converge to the so-called Airy process for any $|\eta|<1/2$. By
following [25], we can derive analytically the parameters $\Gamma$ and
$\Lambda$ in Eq. (2). It turns out that

$\displaystyle\Lambda(\eta)=\pm\frac{1}{2}\sqrt{\frac{1}{2}-\eta^{2}}$ (4)
$\displaystyle[\Gamma(\eta)]^{3}=\frac{(1-4\eta^{2})^{2}}{16(2-4\eta^{2})^{3/2}},$
(5)

where the two branches of Eq. (4) depend on the sign of $X_{ext}(\eta)$.
Notice that $\Gamma(\eta)$ is positive and monotonically decreasing with
$\Gamma(1/2)=0$. When the right-hand-side of Eq. (5) vanishes, the
intersection point $X_{ext}(1/2)=N/4$ lies on the boundary of the square-
lattice and its fluctuations are no longer described by the TW distribution.
In this case, Ref.[56], see also [57] for a related discussion, proved that
the latter are Gaussian with an horizontal width $O(\sqrt{N})$.

We then exploit Eqs. (4) and (5) to benchmark the validity of the MC
simulation. In particular, according to Eq. (2), the probability distribution
($\text{Prob}(\mathcal{Z}(\eta)=\chi)$) of the random variable $\mathcal{Z}$
can be evaluated numerically as 333The notation $\lfloor x\rfloor$ indicates
the integer part of $x$.

Figure 5: Left. Distribution of the random variable $\mathcal{Z}(\eta)$
defined in Eq. (2) calculated from the MC simulation at $\eta=1/4$ and $N=256$
for $\Delta=0$. The red curve is the TW distribution
$\mathcal{F}_{2}^{\prime}(\chi)$. There are no free parameters in the
comparison. Right. The empirical cumulative distribution of the random
variable $\mathcal{Z}_{0}(\eta)$ is plotted together with the cumulative
distribution $\mathcal{F}_{2}$ at $\eta=1/4$ and $N=256$ for $\Delta=0$.

$f_{\eta}(\chi)=\Gamma(\eta)N^{1/3}h_{\eta}\bigl{(}\lfloor\Gamma(\eta)N^{1/3}\chi+N\Lambda(\eta)\rfloor\bigr{)}.$
(6)

Upon substituting Eqs. (4) and (5) into Eq. (6), one compares the MC data at
finite $N$ with the TW distribution $\mathcal{F}_{2}^{\prime}(\chi)$. The
result of this comparison is displayed on the left of Fig. 5 for $\eta=1/4$
and $N=256$. It shows a convincing agreement already for such a value of the
lattice size; it should be remarked that there are no free parameters.

From a statistical viewpoint, it is possible to verify if the random variable
$\mathcal{Z}$ is distributed according to $\mathcal{F}_{2}^{\prime}$, by
performing a Kolmogorov-Smirnov test [58]. The latter requires to construct
the empirical cumulative distribution of $\mathcal{Z}$, which will be denoted
by $F_{\eta}(\chi)$. On the right panel of Fig. 5 is plotted $F_{\eta}(\chi)$
at $\eta=1/4$ and $N=256$ together with the cumulative distribution of the TW
distribution, denoted by $\mathcal{F}_{2}(\chi)$. The test statistics [59, 60]
is
$D_{n}:=\sqrt{n}~{}\text{sup}_{\chi}|F_{\eta}(\chi)-\mathcal{F}_{2}(\chi)|$,
where $n$ is the number of points in Fig. 5. For $n=14$, we obtain
$D_{14}=1.021\dots$, which corresponds to a $p$-value of about $0.2$ [61],
well above the conventionally accepted threshold of $0.1$.

##  4 Numerical study of the cumulants for $0\leq\Delta<1$ and $\eta=0$

Figure 6: Left. Numerical data for the skewness of the empirical distribution
$X_{ext}(\eta=0)$ as a function of $N^{-2/3}$ for $N=64,128,256$. The
intercept of the regression line is the skewness in the thermodynamic limit.
The dashed red line is the TW value at $\beta=2$ which equals $0.224\dots$.
Right. Numerical data for the excess kurtosis of the empirical distribution
$X_{ext}(\eta=0)$ as a function of $N^{-2/3}$ for $N=64,128,256$. The
intercept of the regression line is the excess kurtosis in the thermodynamic
limit. The dashed red line is the TW value at $\beta=2$ which equals
$0.093\dots$.

In this section we report numerical results for the fluctuations of the random
variable $X_{ext}(\eta=0)$, i.e. along the diagonal of the square lattice for
values of the anisotropy in the domain $0\leq\Delta<1$. Details about the
construction of the empirical probability distribution of fluctuations of the
extremal lattice path are given in the Appendix A. To lighten the notation we
will drop the $\eta$-dependence in the remaining part of the discussion and
always assume $\eta=0$. At large $N$, and for any $|\Delta|<1$ the
fluctuations of $X_{\text{ext}}$ shall obey (see Eq. (2))

$X_{ext}\stackrel{{\scriptstyle N\gg 1}}{{\longrightarrow}}N\Lambda+\Gamma
N^{1/3}\mathcal{Z},$ (7)

where $\mathcal{Z}$ is a TW random variable and the parameters $\Lambda$ and
$\Gamma$ are non-universal and $O(1)$. In particular $\Lambda$ can be
estimated from the mean of $X_{ext}$ and $\Gamma$ from its variance.

It should be understood that Eq. (7) is an asymptotic expansion valid for
large enough $N$. In particular, the coefficients $\Lambda$ and $\Gamma$ are
the $N\rightarrow\infty$ limit of non-trivial functions of the lattice size,
which are not known analytically for $\Delta\not=0$. The numerical data will
be then extrapolated to $N\rightarrow\infty$ by assuming the same functional
form for the finite size-corrections of the distribution cumulants discussed
in [38]. A numerical study of $\Lambda$ has been perfomed in detail in Ref.
[41] and will be not repeated here.

| $\Delta=0$ | $\Delta=3/8$ | $\Delta=1/2$ | $\Delta=5/8$
---|---|---|---|---
$s$ | $0.221(5)$ | $0.217(4)$ | $0.220(4)$ | $0.214(4)$
$k$ | $0.089(2)$ | $0.092(2)$ | $0.093(3)$ | $0.092(2)$
Table 1: The regression estimates for the parameter $s$ and $k$. The number in
the bracket is the estimated error in the last digit. The TW values for $s$
and $k$ are $0.224\dots$ and $0.093\dots$ respectively.

Skewness and Excess Kurtosis.—Since the random variable $X_{ext}$ is related
to the TW random variable $\mathcal{Z}$ in Eq. (7) by a linear (affine)
transformation, cumulants of order larger than two such as the skewness and
excess kurtosis, are the same for the two distributions. Their study [30]
provides a convenient parameter free test of the universality conjecture and
will be considered first. We denote then by $s_{N}$ and $k_{N}$ the empirical
skewness and excess kurtosis obtained from the Monte Carlo data for $X_{ext}$
at finite $N$. Following Ref. [38] we assume that the leading finite-size
correction for large $N$ to the cumulants is $O(N^{-2/3})$, namely

$s_{N}\stackrel{{\scriptstyle N\gg
1}}{{\longrightarrow}}s+aN^{-2/3}+o(N^{-2/3});~{}~{}k_{N}\stackrel{{\scriptstyle
N\gg 1}}{{\longrightarrow}}k+bN^{-2/3}+o(N^{-2/3}),$ (8)

where $o(N^{-2/3})$ represents quantities that vanish faster than $N^{-2/3}$
for $N\rightarrow\infty$. The exponent $2/3$ of the first subleading
correction in Eq. (8) was proposed in Ref. [38] in the context of the quantum
XXZ spin chain when analyzing fluctuations of the density profile of the
interacting fermions close to the edge. Albeit the numerical methods (Ref.
[38] employed DMRG techniques) and the systems are different, the convergence
patterns of cumulants toward the thermodynamic limit are resemblant.

A linear fit of the numerical data for $s_{N}$ and $k_{N}$ against $N^{-2/3}$
is shown in Fig. 6 by taking into account points at $N=64,128,256$ for
$\Delta=0,\frac{3}{8},\frac{1}{2}$ and $\frac{5}{8}$. The numerical estimates
for the coefficients $s$ and $k$ in Eq. (2) correspond then to $y$-intercept
of the regression lines; they are collected in Tab. 1 together with their
statistical uncertainties. The red dashed line in Fig. 6 is the TW value. The
agreement is convincing although not prefect for the skewness at
$\Delta\not=0$ and $\Delta\not=1/2$. The data points in Fig. 6 and their error
bars are obtained by averaging over several Monte Carlo histograms (5-10
histograms).
Variance at the ice point.— Let us denote by $\sigma_{N}^{2}$ the variance of
the emprical distribution $X_{ext}$ at finite $N$ and let $\sigma^{2}$ the
variance of the TW distribution. From Eq. (7) and the same scaling conjecture
proposed by Ref.[38] for the higher cumulants, it follows

$\left(\frac{\sigma_{N}^{2}}{\sigma^{2}}\right)^{3/2}\stackrel{{\scriptstyle
N\gg 1}}{{\longrightarrow}}N\Gamma^{3}+cN^{1/3}+o(1).$ (9)

In Sec. 3, we recalled the expression for the parameter $\Gamma^{3}$ along the
lattice diagonal at $\Delta=0$, namely
$\left.\Gamma^{3}\right|_{\Delta=0}=\frac{1}{32\sqrt{2}}$ see Eq. (5) for
$\eta=0$. More recently, Ref. [63] proposed a phenomenological (and not
rigorous) formula for $\Gamma^{3}$ at arbitrary values of $\Delta$, which
involves the second derivative of the arctic curve evaluated at the
intersection point with the lattice diagonal. At $\Delta=1/2$, the so-called
ice point, the arctic curve is an ellipse [19] and the parameter $\Gamma^{3}$
was conjectured to be [63]:
$\left.\Gamma^{3}\right|_{\Delta=1/2}=\frac{1}{16\sqrt{3}}$. The Ansatz
elaborated in [63] relies on the KPZ [32] scaling theory applied to the
fluctuations of the extremal lattice path. Fig. 7 shows a numerical estimation
of the parameter $\Gamma^{3}$ obtained from a linear regression analogous to
that performed for the skweness and excess kurtosis in the previous paragraph.
The agreement with the exact [25] and conjectured [63] expressions at
$\Delta=0$ and $\Delta=1/2$ is prefect. An equivalent result can be obtained
by averaging the distribution of the extremal lattice path along lines
parallel to the lattice diagonal and requiring best fit of the Monte Carlo
data with the TW distribution, see Appendix B.

Figure 7: Numerical estimation of the coefficient $\Gamma^{3}$ in Eq. (7) from
the variance of the empirical distribution of $X_{ext}$ along the lattice
diagonal ($\eta=0$). The predictions at $\Delta=0$ and $\Delta=1/2$ are the
$y$-intercepts of the regression lines, see Eq. (9) and perfectly agree with
the theoretical values represented by the dashed lines, see main text.

##  5 Conclusions

In this paper, we analyzed numerically, with Monte Carlo techniques, the
fluctuations of the extremal lattice path in the six vertex model with DWBC
for values of the anisotropy parameter in the domain $0\leq\Delta<1$. Exact
results at the free fermion point ($\Delta=0$) show that the intersection of
the furthest lattice path with the lattice diagonal follows the Tracy-Widom
distribution. The same is true along any line parallel to the diagonal, as
long as the contact point [56] of the arctic curve with the lattice is
avoided. Our numerical simulations were tested thoroughly against these
analytic results, showing convincing agreement already for lattices with side-
length $N=256$ and without fitting parameters.

Further, we provided strong numerical evidence that the fluctuations of the
extremal lattice path follow the Tracy-Widom distribution also for values of
the anisotropy parameter in the domain $0\leq\Delta<1$. In particular, we
verified the consistency of the scaling hypothesis in Eq. (2) by determining
the skewness and excess kurtosis of the distribution of $X_{ext}$ along the
lattice diagonal. We also verfied, see Eq. (9), that the width of the
fluctuations about the arctic curve is of order $O(N^{1/3})$ with a prefactor
that could be extrapolated from the variance of the Monte Carlo data and
agrees perfectly with the exact value at $\Delta=0$ and the recent conjecture
in [63] at $\Delta=1/2$.

Our finite size scaling analysis requires relatively small lattices and
samples of $O(10^{6})$. Therefore it is computationally less demanding than
other Monte Carlo attempts to uncover the Tracy-Widom distribution [46]. The
results are fully consistent with those that Ref. [38] obtained in quantum
spin chains through DMRG simulations.

Our conclusions should extend also to cases $-1<\Delta\leq 0$, which we did
not analyze in this paper but belong to the same disordered phase [18] of the
model. Finally, it would be interesting to study fluctuations of the boundary
of the inner phase separation curve in the antiferromagnetic regime, see for
instance [64] or understand the statistics of interior lattice paths [65]. To
circumvent the problem of exponentially slow mixing of the Markov chain Monte
Carlo dynamics for large and negative $\Delta$, one might implement the
algorithm proposed recently in [66].

## Acknowledgments

We thank Alexander Abanov, Filippo Colomo, Giacomo Gori, Andrei Pronko,
Herbert Spohn and Jean-Marie Stéphan for enlightening discussions and interest
in this work. JV is especially grateful to Filippo Colomo and Jean-Marie
Stéphan for a careful reading of the manuscript and Herbert Spohn for sharing
his unpublished draft on the subject.

## Appendices

##  Appendix A Numerical evaluation of the distribution of $X_{ext}(\eta)$ for
$|\eta|<\frac{1}{2}$

In this Appendix, we discuss in detail the numerical algorithm that we used to
construct the empirical distribution $X_{ext}(\eta)$ for $|\eta|<\frac{1}{2}$.

For simplicity, we will consider only the isotropic case $a=b$. Consider the
$1+k_{\max}$ parallel lines $y=y_{k}$ with $y_{k}=k/2$ and
$k=0,1,\dots,k_{\max}$, see Fig. 2, and define, see Eq. (1),

$\eta_{k}:=2y_{k}/N=k/N.$ (10)

As discussed before, the parameter $k_{\max}$ must be below $N/2$, in this
paper we chose $k_{\max}=N/4$. We begin by describing how the distribution of
the random variable $X_{ext}$, along the line $y=0$, that is the lattice
diagonal, will be constructed. We start with a random lattice state $S^{(0)}$
and an $N/2$-dimensional vector $v_{0}^{(0)}$:

$v_{0}^{(0)}=(0,~{}0,~{}0,~{}\dots,~{}0).$ (11)

The state $S^{(0)}$ will be used to obtain a new vector $v_{0}^{(1)}$ in the
following way. We move from each of the four corners towards the center of the
lattice until the first vertex with a state different from the state at the
corner is met. The integer variable $x$ will be used to count the steps. For
example, beginning at the NW corner at $x=0$, we walk along the diagonal
until, at step $x=x_{{\rm NW}}$, we reach a vertex which has a state different
from $a_{1}$. Beginning instead at the SW corner, we continue until we reach a
vertex which does not have state $b_{2}$; and so on. The vector $v_{0}^{(1)}$
will then have the entries

$v_{0}^{(1)}(x)=v_{0}^{(0)}(x)+\sum_{\iota\in\\{{\rm NE},~{}{\rm NW},~{}{\rm
SW},~{}{\rm SE}\\}}\delta_{x,~{}x_{\iota}};~{}x=0,~{}1,~{}\dots,~{}N/2-1.$
(12)

Next, we perform one MC sweep, see Sec. 2. This will generate a new,
independent random lattice state $S^{(1)}$. The process is then repeated and
we will obtain a new vector $v_{0}^{(2)}$ from the vector $v_{0}^{(1)}$ and
the state $S^{(1)}$. Altogether $T-1$ sweeps will be made. The final vector
$v_{0}^{(T-1)}$ will thus be obtained from $T$ states
$S^{(0)},~{}S^{(1)},~{}\dots,~{}S^{(T-1)}$. It is convenient to use the
opposite direction, from the center of the lattice towards the corner. To this
end, we define the normalized histogram $h_{0}$, see Sec. 2, as follows

$h_{0}(x):=v_{0}^{(T-1)}\bigl{(}N/2-1-x\bigr{)}/4T;~{}x=0,~{}1,~{}\dots,~{}N/2-1.$
(13)

The distribution obtained from Eq. (13) is illustrated on the upper left side
of Fig. 4 for $\Delta=1/2$ and $N=256$. Here $x$ is the distance from the
center along the diagonal, which in the coordinate systems of Fig. 2 coincides
indeed with the value of the $x$-coordinate in the NW corner. As mentioned in
Sec. 1, the diagonal distance between two neighboring vertices is $1$. In this
paper, $T$ will have the value $T=10^{6}$ throughout.

In order to provide evidence of the validity of the scaling hypothesis in Eq.
(7) also at $\eta\not=0$ and away from the free fermionic point, we construct
$k_{\max}$ histrograms of $X_{ext}(\eta\not=0)$, i.e. along lines parallel to
the lattice diagonal. The vectors
$\\{v_{\eta_{k}}^{(T-1)}\\}_{k=1}^{k_{\max}}$ are built taking into account
that the process is begun on side of the lattice and not in the corner.
Therefore, there are eight simultaneous processes and not four. We thus obtain
another set of vectors $\\{v_{\eta_{k}}^{(T-1)}\\}_{k=1}^{k_{\max}}$. From
these, we will obtain normalized, reversed histograms $\\{h_{\eta_{k}}\\}$ if,
for each $k=1,\dots,k_{\max}$, we define $h_{\eta_{k}}$ as

$h_{\eta_{k}}(x):=v_{\eta_{k}}^{(T-1)}(N/2-1-x+y_{k})/8T;\quad
x=-y_{k},\dots,~{}N/2-1-y_{k}.$ (14)

The histogram $h_{1/4}(x)$ obtained when $k=k_{\max}=N/4$ is displayed for
$\Delta=1/2$ and $N=256$ on the bottom left of Fig. 4.

##  Appendix B Distribution average over lines parallel to the lattice
diagonal

$N$ | $\Delta=0$ | $\Delta=3/8$ | $\Delta=1/2$ | $\Delta=5/8$
---|---|---|---|---
$64$ | $1.822(13)$ | $2.465(12)$ | $2.849(14)$ | $3.47(3)$
$128$ | $3.348(12)$ | $4.557(22)$ | $5.32(4)$ | $6.47(3)$
$256$ | $6.29(7)$ | $8.60(14)$ | $10.05(16)$ | $12.46(24)$
Table 2: The parameter $[\Gamma_{N}]^{3}$ as a function of $N$ and $\Delta$.
The values of $N$ are contained in the leftmost column and the values of
$\Delta$ are shown in the uppermost row. The number in the bracket is the
estimated error in the last digit.

We assume that all the histograms $h_{\eta_{k}}$ obtained in Appendix A once
shifted and rescaled to have the same mean and variance converge for large $N$
to the same distribution. To test the scaling hypothesis in Eq. (2), we then
average over them as described in the following.

First one calculates the mean $\mu_{k}$ and the variance $\sigma_{k}$ of the
distributions $h_{\eta_{k}}$ for all $k=0,\dots,k_{\max}$. We have

$\mu_{k}=\sum_{x=-y_{k}}^{N/2-1-y_{k}}x~{}h_{\eta_{k}}(x),\quad\sigma_{k}^{2}=\sum_{x=-y_{k}}^{N/2-1-y_{k}}(x-\mu_{k})^{2}~{}h_{\eta_{k}}(x).$
(15)

At this point, we will leave the actual lattice and define our distributions
in a new, one-dimensional space. To perform with greater numerical accuracy
the averaging procedure, it is convenient to magnify the bin width by a
constant factor $\theta$, such that $1<\theta<k_{\max}$. In terms of the new
variable $t=\theta x$, the probability distributions in Eq. (14) are then
444The entries of histograms are now indexed by $\theta k_{\max}+1$ integers
and, provided they have the same mean and variance, can be summed without the
need of defining a piecewise function on the reals. The optimal value of
$\theta$ is chosen in order to to render the averaged histogram as smooth as
possible. Its effect on the fluctuations of the random variable $X_{ext}$ can
be trivially traced back, see Eq. (21).

$h^{\sharp}_{\eta_{k}}(t)=\frac{1}{\theta}h_{\eta_{k}}\left(\left\lfloor\frac{t}{\theta}\right\rfloor\right).$
(16)

The value $\theta=8$ will be used in this paper. By following the assumption
stated at the beginning of this section, we further perform for all
$k=0,\dots,k_{\max}$ the linear change of variable

$\xi=\frac{\sigma_{0}}{\sigma_{k}}(t-\theta\mu_{k}),$ (17)

which transforms the probability distributions of Eq. (16) into

$h_{\eta_{k}}^{\flat}(\xi)=\frac{\sigma_{k}}{\theta\sigma_{0}}h_{\eta_{k}}\left(\left\lfloor\frac{\sigma_{k}}{\sigma_{0}\theta}\xi+\mu_{k}\right\rfloor\right).$
(18)

It is easy to show, by recalling Eq. (15), that all the probability
distributions in Eq. (18) have zero mean and variance
$(\theta\sigma_{0})^{2}$. According to our working hypothesis, they are then
different histogram representations of the same probability distribution of
the random variable

$\Xi:=\theta X_{ext}(\eta=0)+C.$ (19)

In Eq. (19), $C$ is a constant that ensures that $\Xi$ has zero mean. We
estimate the probability distribution of the random variable $\Xi$ with the
average

$\bar{h}_{0}(\xi):=\frac{1}{1+k_{\max}}\sum_{k=0}^{k_{\max}}h^{\flat}_{\eta_{k}}(\xi);$
(20)

in Eq. (20) the variable $\xi$ belongs to the union of all the support of the
histograms $h_{\eta_{k}}$ in Eq. (18).

Instead of calculating the cumulants as in Sec. 4, one could test the scaling
hypothesis in Eq. (2) (and Eq. (7)) by directly fitting Eq. (20) with the TW
distribution. In particular, recalling Eq. (19), we analyze weather

$\bar{h}_{0}(\xi)\stackrel{{\scriptstyle N\gg
1}}{{\longrightarrow}}\frac{1}{\Gamma_{N}\theta}\mathcal{F}_{2}^{\prime}\left(\frac{\xi-C-\theta\Lambda_{N}}{\Gamma_{N}\theta}\right),$
(21)

with $\mathcal{F}_{2}^{\prime}$, the TW distribution, by fitting at finite $N$
the two parameters $\Lambda_{N}$ and $\Gamma_{N}$ above. As also discussed
already in Sec. 4, $\Lambda_{N}$ and $\Gamma_{N}$ are non-trivial and unknown
functions of the lattice size $N$ such that (see Eq. (2))

$\lim_{N\rightarrow\infty}\frac{\Lambda_{N}}{N}=\Lambda(\eta=0),\quad\lim_{N\rightarrow\infty}\frac{\Gamma_{N}}{N^{1/3}}=\Gamma(\eta=0).$
(22)

Working with $N$ up to $256$, it becomes crucial to evaluate the leading
finite size corrections to those parameters, in particular following [38] one
has, see Sec. 4,

$\frac{[\Gamma_{N}(0)]^{3}}{N}\stackrel{{\scriptstyle N\gg
1}}{{\longrightarrow}}[\Gamma(0)]^{3}+cN^{-2/3}+o(N^{-2/3}).$ (23)

In Tab. 2 are collected the numerical estimates for $[\Gamma_{N}(0)]^{3}$
obtained by requiring best fit of the Monte Carlo data for the averaged
distribution of the extremal lattice path with the TW distribution (see Eq.
(21)). A linear fit of the numerical data for $[\Gamma_{N}]^{3}/N$ given in
Tab. 2 against $N^{-2/3}$ is shown on the right panel of Fig. 8 by taking into
account points at $N=64,128,256$. The numerical prediction for the coefficient
$[\Gamma(0)]^{3}$ in Eq. (2) corresponds then to $y$-intercept of the
regression lines. As in Sec. 4, see Fig. 7, the red dashed lines are the exact
[25] and conjectured values [63] at $\Delta=0$ and $\Delta=1/2$ respectively.
The agreement is perfect also with this approach.

Figure 8: Left. Numerical data for the distribution $\bar{f}_{0}(\chi)$ in Eq.
(24) at different values of $\Delta$ and $N=256$. The parameters $\Gamma_{N}$
and $\Lambda_{\beta}$ are obtained by requiring a best fit of the data with
the TW distribution, which is the red curve underneath. Right.  Numerical data
for $[\Gamma_{N}]^{3}/N$ from Tab. 2 together with the regression lines, see
Eq. (23). The red dashed lines are the analytical results for the
$y$-intercept at $\Delta=0$ and $\Delta=1/2$, Sec. 4.

Finally, Eq. (21) is equivalent to the statement that for large enough $N$

$\bar{f}_{0}(\chi):=\Gamma_{N}\theta~{}\bar{h}_{0}\bigl{(}\Gamma_{N}\theta~{}\chi+\theta\Lambda_{N}+C\bigr{)}$
(24)

converges to the TW distribution for any $0\leq\Delta<1$. On the left of Fig.
8, we show on the same plot the numerical results for $\bar{f}_{0}(\chi)$ at
$N=256$ and $\Delta=0,~{}\frac{3}{8},~{}\frac{1}{2},~{}\frac{5}{8}$. The
continuous red curve underneath, which is barely visible, is the TW
distribution $\mathcal{F}_{2}^{\prime}(\chi)$.

## References

  * [1] W. Feller, An Introduction to Probability Theory and Its Applications vol. I, 3rd Edition, Wiley (1991), ISBN: 978-0-471-25708-0.
  * [2] A. Okunkov, Limit shapes, real and imaginary, AMS colloquium 2007, available at http://math.columbia.edu/~okounkov/AMScolloq.pdf.
  * [3] F. Colomo and A.Sportiello, Arctic curves in the six-vertex model on generic domains: The Tangent Method, J. Phys. A 164 (2016) 1488, https://arxiv.org/abs/1605.01388.
  * [4] R. Kenyon and A. Okounkov, Limit shapes and the complex Burgers equation, Acta Math. 199 (2007), 263–302, https://arxiv.org/abs/math-ph/0507007.
  * [5] W. Jockusch, J. Propp, and P. Shor, Random domino tilings and the arctic circle theorem, arXiv: math.CO/9801068, https://arxiv.org/abs/math/9801068.
  * [6] P. Ferrari and H. Spohn, Domino Tilings and the six-vertex model at its free fermion point, J. Phys. A: Math. Gen. 39 (2006), 10297-10306, https://arxiv.org/abs/cond-mat/0605406.
  * [7] L. Pauling, The Structure and Entropy of Ice and of Other Crystals with some Randomness of Atomic Arrangement, J. Am. Chem. Soc. 57 (12): 2680-2684, (1935).
  * [8] E. Lieb, Phys. Rev. 162, 162, (1967); E.H. Lieb, Phys. Rev. Lett. 18 (1967) 692; Phys. Rev. Lett. 18 (1967) 1046; Phys. Rev. Lett. 19 (1967) 108.
  * [9] R.J. Baxter, Exactly solved models in statistical mechanics, Dover Book on Physics, (2008), ISBN: 978-0486462714.
  * [10] V. Korepin, N. Bogoliubov and A. Izergin, Quantum inverse scattering method and correlation functions, Cambridge University Press (1993).
  * [11] N. Reshetikhin, Lectures on integrable models in statistical mechanics, In: “Exact methods in low-dimensional statistical physics and quantum computing”, Proceedings of Les Houches School in Theoretical Physics, Oxford University Press, (2010), https://arxiv.org/abs/1010.5031.
  * [12] V. Korepin, Calculation of Norms of Bethe Wavefunctions, Comm. Math. Phys. 86, 391-418 (1982).
  * [13] A. Izergin, Partition Function of the 6-Vertex Model in a Finite Volume, (Russian) Dokl. Akad. Nauk USRR, v. 297, 331-333, (1987).
  * [14] A. Izergin, D. Coker, and V. Korepin, Determinant formula for the six-vertex model, J. Phys. A 25 (1992), 4315–4334.
  * [15] P. Zinn-Justin, Six-vertex model with domain wall boundary conditions and one-matrix model, Phys. Rev. E 62 (2000), 3411–3418, https://arxiv.org/abs/math-ph/0005008.
  * [16] K. Eloranta. Diamond Ice. Journal of Statistical Physics 96 (1999), 1091-1109.
  * [17] V. Korepin and P. Zinn-Justin. Thermodynamic limit of the six-vertex model with domain wall boundary conditions. J. Phys. A 33 (2000), p. 7053-7066, https://arxiv.org/abs/cond-mat/0004250
  * [18] P. Bleher and K. Liechty. Random Matrices and the Six-Vertex Model. CRM Monograph Series Volume: 32, (2014), ISBN: 978-1-4704-0961-6.
  * [19] F. Colomo and A. Pronko, The arctic curve in the six vertex model, J. Stat. Phys. 138, 662-700 (2010), https://arxiv.org/abs/0907.1264.
  * [20] F. Colomo, A. Pronko and P. Zinn Justin, The arctic curve of the domain wall six vertex model in its anitferroelectric regime, J. Stat. Mech: Theor. Exp. (2010) L03002, https://arxiv.org/abs/1001.2189.
  * [21] A. Aggarwal, Arctic Boundaries of the Ice Model on Three-Bundle Domains, https://arxiv.org/abs/1812.03847.
  * [22] K. Johansson, Shape Fluctuations and Random Matrices, Commun.Math.Phys. 209, 437 - 476 (2000), https://arxiv.org/abs/math/9903134.
  * [23] K. Johansson, The arctic circle boundary and the Airy process. The Annals of Probability (2005), Vol. 33, No. 1, p. 1-30, https://arxiv.org/abs/math/0306216.
  * [24] H. Widom and C. Tracy, Level Spacing distributions and the Airy kernel, Commun. Math. Phys. 159 (1994) 151-174, https://arxiv.org/abs/hep-th/9211141.
  * [25] N. Allegra, J. Dubail, J-M. Stéphan, and J. Viti. Inhomogeneous field theory inside the arctic circle. J. Stat. Mech. (2016) 053108, https://arxiv.org/abs/1512.02872
  * [26] J. Baik, P. Deift and K. Johansson, On the Distribution of the Length of the Longest Increasing Subsequence of Random Permutations, J. Amer. Math. Soc. 12 (1999), no. 4, 1119–1178, https://arxiv.org/abs/math/9810105.
  * [27] M. Praehofer and H. Spohn, Scale Invariance of the PNG Droplet and the Airy Process, J. Stat. Phys. 108 (5-6): 1071-1106 (2002), https://arxiv.org/abs/math/0105240.
  * [28] C.Tracy and H. Widom, Asymptotics in asep with step initial condition, Communications in Mathematical Physics 290(1), 129 (2009), https://arxiv.org/abs/0807.1713.
  * [29] A. Ayyer, S. Chhita and K. Johansson, GOE fluctuations for the maximum of the top path in alternating sign matrices, https://arxiv.org/abs/2109.02422.
  * [30] K. Takeuchi and M. Sano, Universal fluctuations of growing interfaces: evidence in turbulent liquid crystals, Phys. Rev. Lett. 104, 230601 (2010). https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.104.230601.
  * [31] Z. Chen, J. de Gier, I. Hiki and T. Sasamoto, Exact confirmation of 1D nonlinear fluctuating hydrodynamics for a two-species exclusion process, Phys. Rev. Lett. 120, 240601 (2018), https://arxiv.org/abs/1803.06829.
  * [32] M. Kardar, G. Parisi, and Y.-C. Zhang, Dynamic scaling of growing interfaces, Phys. Rev. Lett. 56, 889–892 (1986).
  * [33] T. Sasamoto and H. Spohn, The one-dimensional KPZ equation: an exact solution and its universality, Phys. Rev. Lett. 104, 230602 (2010), https://arxiv.org/abs/1002.1883.
  * [34] V. Eisler and Z. Racz, Full counting statistics in a propagating quantum front and random matrix spectra, Phys. Rev. Lett. 110, 060602 (2013), https://arxiv.org/abs/1211.2321
  * [35] J. Viti, J-M. Stéphan, J. Dubail and M. Haque, Inhomogeneous quenches in a fermionic chain: exact results, EPL 115 (2016) 40011, https://doi.org/10.1209/0295-5075/115/40011.
  * [36] H. 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 (2006), https://arxiv.org/abs/cond-mat/0512011, Fundamental Problems in Statistical Physics.
  * [37] P. Deift, Universality for mathematical and physical systems, arXiv:math-ph/0603038 (2006), https://arxiv.org/abs/math-ph/0603038.
  * [38] J-M. Stéphan, Free fermions at the edge of interacting systems, SciPost Phys. 6, 057 (2019), https://arxiv.org/abs/1901.02770.
  * [39] J-M. Stéphan, Extreme boundary conditions and random tilings, SciPost Phys. Lect. Notes 26 (2021), https://arxiv.org/abs/2003.06339.
  * [40] D. Allison and N. Reshetikhin. Numerical study of the 6-vertex model with domain wall boundary conditions. Annales de l’Institut Fourier, Tome 55, no 6 (2005), p. 1847-1869, https://arxiv.org/abs/cond-mat/0502314.
  * [41] I. Lyberg, V. Korepin, and J. Viti, The density profile of the six vertex model with domain wall boundary conditions, J. Stat. Mech. (2017) 053103, https://arxiv.org/abs/1612.06758.
  * [42] I. Lyberg, V. Korepin, G. A. P. Ribeiro and J. Viti, Phase separation in the six-vertex model with a variety of boundary conditions, J. Math. Phys. 59, 053301 (2018), https://arxiv.org/abs/1711.07905.
  * [43] D. Randall and P. Tetali. Analyzing Glauber dynamics by comparison of Markov chains. Journal of Mathematical Physics, 41(3):1598–1615, 2000, https://doi.org/10.1063/1.533199.
  * [44] T. Liu, Torpid mixing of Markov chains for the six-vertex model on $\mathbb{Z}_{2}$. In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques( APPROX/RANDOM 2018), https://arxiv.org/abs/1809.02703.
  * [45] M. Fahrbach and D. Randall, Slow Mixing of Glauber Dynamics for the Six-Vertex Model in the Ordered Phases, Proceedings of the 23rd International Conference on Randomization and Computation (RANDOM 2019), https://arxiv.org/abs/1904.01495.
  * [46] C. Mendl and H. Spohn, Searching for the Tracy-Widom distribution in nonequilibrium processes, : Phys. Rev. E 93, 060101(R) (2016) https://arxiv.org/abs/1512.06292.
  * [47] K. Eloranta, The bounded 19-vertex model (2017) https://arxiv.org/abs/1710.03609; K. Eloranta, The bounded 15-vertex model (2018), https://arxiv.org/abs/1807.07567.
  * [48] R. Keesman and J. Lamers, A numerical study of the F-model with domain-wall boundaries, Phys. Rev. E 95, 052117 (2017), https://arxiv.org/abs/1702.05474.
  * [49] D. Keating and A. Sridhar, Random Domino Tilings with the GPU, J. Math. Phys. 59, 091420 (2018), https://arxiv.org/abs/1804.07250.
  * [50] B. Debin, P. Di Francesco and E. Guitter, Arctic curves of the twenty-vertex model with domain wall boundaries, J Stat Phys 179, 33-89 (2020), https://arxiv.org/abs/1910.06833.
  * [51] J-F. de Kemmeter, B, Debin and P. Ruelle, Arctic curves of the 6V model with partial DWBC and double Aztec rectangles, J. Phys. A: Math. Theor. 55 (2022) 305004, https://arxiv.org/abs/2203.08506.
  * [52] P. Belov and N. Reshetikhin, The two-point correlation function in the six-vertex model, J. Phys. A: Math. Theor. 55 155001 (2022), https://arxiv.org/abs/2012.05182.
  * [53] G. T. Barkema and M. E. J. Newman. Monte Carlo simulation of ice models. Phys. Rev. E 57 (1998), 1155, https://arxiv.org/abs/cond-mat/9706190.
  * [54] O. F. Syljuåsen and M. B. Zvonarev. Directed-loop Monte Carlo simulations of vertex models. Phys. Rev. E 70 (2004), 016118, https://arxiv.org/abs/cond-mat/0401491.
  * [55] F. Alet, Y. Ikhlef, J. Jacobsen, G. Misguich and V. Pasquier, Classical dimers with aligning interactions on the square lattice, Phys. Rev. E 74, 041124 (2006), https://arxiv.org/abs/cond-mat/0607747.
  * [56] V. Gorin, From Alternating Sign Matrices to the Gaussian Unitary Ensemble, Comm. Math. Phys., 332, no. 1 (2014), https://arxiv.org/abs/1306.6347.
  * [57] V. Kapitonov and A. Pronko, Six-Vertex Model as a Grassmann Integral, One-Point Function, and the Arctic Ellipse, Journal of Mathematical Sciences volume 264, pages 313–346 (2022), https://doi.org/10.1007/s10958-022-06000-w
  * [58] M. Hollander, D. Wolfe and E. Chicken, Nonparametric Statistical Methods, John Wiley and Sons, 3rd Edition, 2015.
  * [59] A. N. Kolmogorov, Sulla determinazione empirica di una legge di distribuzione. Giorn. Ist. Ital. Attuari 4, 83-91 (1933).
  * [60] W. Feller, On the Kolmogorov-Smirnov Limit Theorems for Empirical Dis- tributions. Annals of Mathematical Statistics, 19, 177-189 (1948) https://doi.org/10.1214/aoms/1177730243
  * [61] J. Vrbik, Deriving CDF of Kolmogorov-Smirnov Test Statistic. Applied Mathematics, 11, 227-246 (2020) https://doi.org/10.4236/am.2020.113018.
  * [62] https://reference.wolfram.com/language/ref/FindFit.html
  * [63] M. Praehofer and H. Spohn, Domain wall fluctuations of the six-vertex model at the ice point, https://arxiv.org/abs/2305.09502.
  * [64] S. Chhita and K. Johansson, Domino statistics of the two-periodic Aztec diamond, Advances in Mathematics 294 (2016) 37–149, https://arxiv.org/abs/1410.2385
  * [65] B. Debin, J-F. de Kemmeter and P. Ruelle, Fluctuations and arctic curve in the Aztec diamond, https://arxiv.org/abs/2301.00600
  * [66] J-Y. Cai and T. Liu, An FPTAS for the square lattice six-vertex and eight-vertex models at low temperatures, Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms SODA (2021), https://epubs.siam.org/doi/abs/10.1137/1.9781611976465.92