Domain wall fluctuations of the six-vertex model at the ice point

We report on Monte-Carlo simulations of the six-vertex model with domain wall boundary conditions. In thermal equilibrium such boundary conditions force a fluctuating line separating the disordered region from the perfectly ordered ones. Specifically we study the ice point at which all vertex weights are equal. With high precision the one-point fluctuations of the line are confirmed to be of order N13 and governed by the Tracy–Widom distribution. Furthermore, the non-universal scaling coefficients are computed for a wide range of interaction strengths.


Introduction
The six-vertex model is one of the most famous Bethe ansatz solvable models of Statistical Mechanics in two dimensions [1,2].While various physical interpretations are used, in our contribution the model will be viewed as a statistical ensemble of possibly touching, but non-crossing up-right paths, see Figure 1.In Figure 1 we imposed domain wall boundary conditions (DWBC), which means that the lines are kept fixed at the boundary of the square.The lattice paths will be viewed as level lines of a height function, which takes the value 0 in the South-West corner and increases by 1 whenever crossing a level line.Thus the height function takes the value 0 at the South (S) and East (E) borders, while increasing as a slope 1 staircase at the West (W) border and decreasing as a slope −1 N E S W Figure 1: One realization of the level lines at ∆ = 1  2 for N = 64 with DWBC.Even for such a fairly small system one clearly distinguishes a disordered zone located in the center from the perfectly ordered facets situated in the four corners.Our focus is the facet edge, more specifically its SE section.staircase at the North (N) border.In the interior of the square the height is allowed to fluctuate.Studied will be a random surface, which is defined through the graph of an integer-valued random height function over the basic domain [0, N + 2] 2 increasing in the East-to-West and South-to-North directions.
Physically one could think of a soap bubble, i.e. an elastic membrane with non-zero surface tension which is pinned at a given frame.Another interpretation is a threedimensional crystal, whose shape is modelled by a height function satisfying the monotonicity conditions mentioned above.
We are interested in the large N limit.As suggested by Figure 1, one expects a deterministic macroscopic shape.It consists of an interior zone surrounded by four perfectly ordered regions.In the interior domain the macroscopic shape is rounded with nonzero curvature.Microscopically the height is disordered with fluctuations of order log(N ), more precisely the height statistics are those of a massless Gaussian field with a covariance depending on the specific macroscopic location in the interior zone [3].The height gradients are critical as signalled by a slow decay of their correlations.On the other hand, in the four corners the height function is flat with no fluctuations.
In addition to the non-crossing constraint, the lines have a nearest neighbor interaction with a strength usually denoted by ∆.Its precise definition will be given below.At ∆ = 0, the model can be handled using free fermion techniques.In this case an explicit formula for the macroscopic shape is available [4,5,6,7,8], which indeed has four perfectly ordered facets located in the four corners.The facet edge is an ellipse reflection symmetric with respect to the diagonals and touching each one of the four bounding edges.For general ∆ the information is less complete.For −1 < ∆ < 1 the macroscopic shape is qualitatively similar to the one with free fermion parameters.For ∆ > 1 the disordered region degenerates into a line with either diagonal or anti-diagonal orientation and perfect facets in the lower respectively upper triangle.For ∆ < −1, minimizing the surface free energy an inner facet appears.Inside this facet the level lines are purely zig-zag with distance 1.A defect costs only an energy of order one.Hence statistical errors can be accommodated with nonzero density, decreasing to zero as ∆ → −∞.The outer facet edge gets further deformed, but in essence does not notice the birth of the inner facet.The disordered zone has the shape of an annulus.We refer to [9] for typical line configurations at various choices of ∆.
The macroscopic shape has a universal feature known as Pokrovsky-Talapov (PT) law, which refers to macroscopic crystal shapes [10].Close to the facet edge the deviation from the facet plane varies with the position x relative to the facet edge as s(x) = 0 for x ≤ 0 and s(x) = ±x δ for x > 0, the left describing the facet and the right the rounded part.The PT law predicts in our context the value of δ to be δ = 3/2.This implies that close to the facet the line density has a square root singularity, which in turn translates into a distance of order N 1 3 between neighboring level lines [28].In case of DWBC, except for the free fermion parameters, the analytical formula for the macroscopic shape is not sufficiently explicit to deduce PT.For the outer facet edge, however, a parametric representation has been accomplished [11,12,13,14].The representation is based on exact hole probabilities, i.e. the probability that the rectangle [n, N + 2] × [0, m] is free of lines.
The focus of our contribution are fluctuations of the facet edge.Based on a random walk analogy one might expect that the fluctuations are Gaussian with size N 1 2 , but this would miss the non-crossing constraint between neighboring level lines.In fact, assuming the validity of PT, the facet edge fluctuations are expected to be of order N 1 3 .Finer details can be understood through the connection with stochastic growth processes.For this purpose one regards N as a time parameter.Increasing the system size from N to N + 1 induces a random growth of the facet edge.At the free fermion point the update turns out to be governed by a Markov chain with a structure very similar to the discrete time asymmetric simple exclusion process starting from step initial conditions.With such input, it is ensured that the Kardar-Parisi-Zhang (KPZ) theory of growing interfaces [15] applies.In particular, as established in [16], in the limit of large N the facet edge statistics at a single reference point is determined by the Tracy-Widom distribution of the largest eigenvalue of a N × N GUE random matrix.
With this background one might ask whether the GUE Tracy-Widom distribution is special for the free fermion parameters or possibly a generic feature of the six-vertex model.Amongst experts, the latter version seems to be favored.But when pressed the actual evidence is scarce.One natural approach is to test through Monte Carlo simulations [17,18,19].However size and number of samples are not large enough for our purposes.More extensive Monte Carlo simulations are reported in [20].We try to improve the situation through Monte-Carlo simulations of the six-vertex model at the diagonal point with ∆ = 1  2 .This specific choice of parameters is known as ice point for which all admissible height configurations have the same weight.The corresponding growth process is no longer a Markov chain.But KPZ scaling theory is still applicable.We then follow the standard route and first determine the non-universal coefficients of the KPZ scaling theory [21].They are expressed in terms of the macroscopic edge function and the derived formulas are valid for the outer edge at any ∆.Thereby we arrive at a parameter-free numerical fit of our numerical data.
From the probabilistic side indirect but very strong evidence has been obtained.According to KPZ theory, relative to a coordinate system rotated by π/4, the minimum of the SE facet edge should be distributed according to the GOE Tracy-Widom distribution [22].For parameters at the ice point this property has been proved recently in [23].
The ice point is also of interest in combinatorics, more specifically for random alternating sign matrices (ASM).A matrix A ij , i, j = 1, ..., N + 1, is ASM if A ij takes values 0, ±1 and, along each row and column, the +1 and −1 entries alternate and sum up to 1. Investigated are statistical properties of ASM under the uniform distribution [24,25].For DWBC the connection is based on a mapping between ASM and vertex configurations under which the uniform distribution of ASM matrices corresponds to the ice point.Thus our results also bear on fluctuations of random ASMs.In the following section we define the six-vertex model and discuss its phase diagram for the case of DWBC.As input to the KPZ scaling theory the non-universal coefficients are computed for the outer facet edge.We discuss Monte Carlo (MC) schemes.The numerical simulations have been carried out for the ice point labeled as "ASM" and for the symmetric free fermion point tagged as "Aztec" for reasons explained below.Our results are displayed in Section 5, including a detailed comparison with the predictions of the theory.

Domain wall boundary conditions
The most direct route to the six-vertex model is through a tiling problem.Given are six distinct tiles of side-length 1, where we follow the convention of [7,13], see Figure 2.An admissible tiling is defined by having no broken lines.For DWBC we consider the square, Λ N = [0, N + 2] 2 ⊂ R 2 , of side length N + 2, with the origin at the SW corner point, compare with Figure 3.The boundary region is The West boundary has horizontal line segments placed at integers and the North boundary has vertical line segments also placed at integers.The square [ 1 2 , N + 3 2 ] 2 is tiled under the constraint that the thereby generated lines end only at the West and North edge.The six tiles have Boltzmann weights ω α > 0, α = 1, ..., 6.The weight of an admissible tiling is given by the product where α(i, j) is the index of the tile centered at (i, j) and i, j = 1, ..., N + 1.An admissible tiling corresponds to a collection of up-right lattice paths, whose possible turns are only at the grid points of Z 2 .The paths may touch, but they never cross.These paths are the level lines of a height function ϕ N (x 1 , x 2 ), 0 ≤ x 1 , x 2 ≤ N + 2, which is constant on each unit square, increases by one at each NW oriented crossing of a level line, and is normalized as The conventional parametrization of the Boltzmann weights is given by The free energy depends on a, b, c only through a/c and b/c, which defines the interaction strength H, V control the slopes of ϕ N .For DWBC it suffices to set H = 0 = V , which leaves a/c, b/c as free parameters.Symmetric case means a = b.The free fermion parameters form the quarter circle defined by ∆ = 0, see Figure 4. Ice point refers to a = 1, b = 1, c = 1, which implies a uniform distribution over all admissible line configurations.In the following we use the acronym ASM for the ice point.The bijection from admissible DWBC configurations to (N + 1) × (N + 1) matrices is given by Using DWBC it follows that the so defined matrix is ASM.The model at the symmetric free fermion point, a = b = 1, c = √ 2, will be labelled as "Aztec" for the following reason: Due to DWBC one can equivalently assign the weights corresponding to the so called 2-enumeration of ASMs, where each −1 in a given ASM has weight 2. As explained in [7] in a more general context, there is a natural one-to-one correspondence between domino tilings of an Aztec diamond with size N and 2-enumerated ASMs of size N + 1 (to account for the 2-enumeration one introduces a 7th tile of type, say, 5 ′ , which is a twin of tile 5. Domino tilings can then be translated into vertex configurations via a local rule).In this case the macroscopic facet edges form a circle which has been termed arctic circle, since it separates the four frozen zones at the corners from the disordered zone in the interior.For general weights, one expects that there is a deterministic macroscopic limiting shape function ϕ ma such that lim ⌊•⌋ denoting integer part.In Figure 4 we display the phase diagram for ϕ ma .For the domain D (disordered), ϕ ma has perfect facets in the four corners of [0, 1] 2 , clockwise with tiles of type 1,4,2,3, compare with Figure 1 and 2, and interpolates smoothly with a non-zero curvature in-between.In the domain AF (anti-ferroelectric) in addition there appears a central facet, isolated from the outer facets, which has a pattern consisting of alternating tiles of type 5,6 with Poisson like random errors.In the domains F and F (both ferroelectric) the interior zone degenerates to a line, either diagonal or anti-diagonal.The facets consist either of type 1,2 (diagonal interface) or type 3,4 (anti-diagonal interface).
For finite N there is the well-known Izergin-Korepin formula [26] for the free energy.Its asymptotics, including subleading terms, is scattered in various contributions and summarized in [27].Some qualitative features of the phase diagram can be recovered.The features for interest to us concern the local free energy, which is beyond the Izergin-Korepin formula.

Scaling theory
We turn to the finer structure of the facet edge in case of DWBC and parameters in the domain D. By symmetry, it suffices to consider only a quarter section of the edge, for which we choose the SE facet, i.e. {(x 1 , x 2 ) ∈ [0, N +2] 2 | ϕ N (x 1 , x 2 ) = 0}.It will be convenient to view its NW boundary as the graph of a function h N in a coordinate system centered at (N + 1, 1) and rotated by π/4, compare with Figure 3.We refer to h N restricted to [−N, ..., N ] as the edge function.One always has h N (−N ) = N = h N (N ).Furthermore h N is piecewise linear between integers and takes integer values on {−N, . . ., N } with increments being ±1.For example, the edge function h N corresponding to Figure 3, where N = 2, is specified by the values Attaining a macroscopic shape as in ( 5) strongly indicates the existence of a limiting edge function h ma such that lim A parametric representation of h ma is computed in [11,12], see also the recent contributions [13,14].
Assuming that the edge fluctuations are governed by the KPZ universality class, one expects that to leading order in N [28].Here ξ TW is distributed according to the Tracy-Widom distribution from GUE random matrix theory with distribution function F 2 (s) = Prob(ξ TW ≤ s) and Γ(x) is a model dependent parameter.The Tracy-Widom distribution has a negative mean.Therefore Eq. (7) indicates that the actual first level line is slightly above the macroscopic edge profile.One knows that for |x| = 1 2 the limit (7) fails [29,30] and convergence non-uniform in x is the rule.The most rapid convergence is expected to be close to x = 0.The KPZ scaling theory relates Γ(x) to more directly accessible properties of the model [31].For DWBC this amounts to linking Γ(x) with h ma (x).In more complicated models such input might not be available and has then to be extracted from the numerical data (or from an actual experiment [32]).
The KPZ scaling theory starts from a growth model, which on the macroscopic scale has a growth velocity v(∂ x h) depending on the local slope ∂ x h.The shape of the growing cluster is hence governed by For the six-vertex model with DWBC the size N plays the role of the time parameter.We do not attempt to write down a stochastic dynamics linking N to N + 1.But one can still expect the validity of (8) for a suitable choice of v.If (8) holds, then in our context h ma must be determined through a self-similar solution of the form h(x, t) = th ma (t −1 x) and one concludes that h ma (x) is the Legendre transform of v(u), u = h ′ ma .In particular In the next few lines we suppress the x-dependence for conciseness.
The KPZ scaling theory [21] asserts that The coefficient A is determined by the local roughness of the interface.More precisely, for large N , j → h N ⌊N x⌋ + j looks like a random walk with drift, implying the variance Note that here one requires 1 ≪ |j| ≪ N 2 3 , the latter being the scale of mesoscopic fluctuations.
It remains to compute A. For large N the edge h N is separated from the neighboring level line by N 1 3 .Thus, for the computation of A, contacts between the two lines can be ignored and it suffices to only consider a single up-right lattice path.The number of tiles 1 and 2 is not modified.The up-and right-unit segments lines have weights be λ and be −λ , where c = 1 without loss of generality and we introduced the parameter λ to control the bias.We use the diagonal transfer matrix, consistent with the frame introduced for h N .The two-step transfer matrix T is obtained as for i even, with all other matrix elements vanishing.T is a two-periodic Toeplitz matrix.Summing over all 2n-step walks, (X i ) i=0,...,2n with σ i := X i − X i−1 ∈ {±1}, i = 1, . . ., 2n, starting with X 0 = 0 or X 0 = 1 and allowing for arbitrary end points, yields the partition function with ψ = (1, 1).Here the 2 × 2 matrix L(λ) is given by For large n the partition function is dominated by the largest eigenvalue of L(λ), which is determined to Asymptotically the mean of the random walk is given by (Z ′ 2n /Z 2n ) ≃ n(E ′ /E) = 2nu, where u is equated with h ′ ma .By definition the variance (Z ′ 2n /Z 2n ) ′ ≃ n(E ′ /E) ′ = 2nA.Somewhat unexpectedly, the asymmetry parameter λ can be eliminated explicitly from the two equations Reintroducing c the result reads Thus, according to ( 9), (10), and ( 19), the scale coefficient Γ(x) is determined through Since there is no restriction on b/c, our argument remains valid for the outer facet for all ∆ < 1.The classic result of Colomo and Pronko [11] yields h ma explicitly for ∆ = 1 2 and ∆ = 0. Our result (19) allows quantitative prediction of the scaling form (7) along each direction x with |x| < 1 2 to be compared with numerical results, see ( 24), ( 26) in Section 5.
In principle the KPZ scaling theory should also apply to the inner facet in the case ∆ < −1.However the required input is not yet available.Firstly the position of the macroscopic facet edge is not known sufficiently explicit.Furthermore, since the facet has Poisson type defects the microscopic facet edge is only fuzzily defined.To obtain the coefficient A one had to rely on a more sophisticated reasoning.

Monte Carlo simulations
We construct a Markov chain in such a way that its unique stationary measure is the normalized Gibbs distribution (1) of the six-vertex model with DWBC.First we explain the standard detailed balance approach, discussing further options at the end of the section.For the Monte Carlo scheme used here, it is convenient to view the volume under the height function ϕ N as made up of unit cubes stacked on top of each other.An allowed MC move corresponds to adding or removing a single cube in such a way, that the sixvertex constraint is maintained.The MC moves are conveniently encoded by introducing plaquettes, each consisting of four tiles arranged as a square, see Figure 5.There are 6 4 distinct plaquette states.The total number of admissible plaquette states turns out to be 82.Out of these only 32 are allowed to undergo a transition.To enumerate all possible transitions one groups them into 16 pairs of plaquette states having identical outward lines.Let us denote by (S, S ′ ) such a pair of plaquette states, such that S ′ arises from S by adding a unit cube in the height function picture as in Figure 5.The probability to jump from S to S ′ is denoted by p + (S → S ′ ).With probability p + (S → S) = 1 − p + (S → S ′ ) there is no jump.Correspondingly we denote the jump from S ′ to S by p − (S ′ → S) and the probability of no jump is p − (S ′ → S ′ ) = 1 − p − (S ′ → S).The condition of detailed balance is imposed as where w(S) is the weight of the plaquette state S as in (1).To write down more explicitly the ratio w(S ′ )/w(S) we introduce four binary variables η NE , η SW , η NW , η SE taking values 0, 1 and referring to the respective tile of plaquette S, η NE = 0 for tile 4 at NE, 1 for tile 6 at NE, η SW = 0 for tile 3 at SW, 1 for tile 6 at SW, η NW = 0 for tile 2 at NW, 1 for tile 5 at NW, η SE = 0 for tile 1 at SE, 1 for tile 5 at SE.
Then, with these conventions (recall the weights a, a, b, b, c, c for tiles 1 through 6, respectively), .
The detailed balance condition determines only the ratio between two corresponding transition probabilities.For parallel update schemes a conventional choice is At the ASM point w(S ′ )/w(S) = 1 and all transition probabilities equal 1 2 .For the actual simulation we note that disjoint plaquettes can be updated in parallel and independently.Their corresponding transition probabilities have to be multiplied.There are four distinct ways to group all tiles into disjoint plaquettes.Hence a complete Monte Carlo step consists of four consecutive sub-routines with parallel update of the distinct plaquette coverings.In the case of ASM this procedure can be contracted to two half-steps.All configurations reachable during a half-step have the same transition probability.Denoting by k the number of flippable plaquettes in a half-step, one of these 2 k configurations can be chosen by consuming only k random bits from a bitwise reliable random number generator.For such parallel MC update scheme any other choice of parameters complicates the algorithm and can lead to bottlenecks, i.e. very small transition rates, which would slow down convergence.
A generic difficulty of Markov Chain Monte Carlo simulations is to ensure that the chain equilibrates.In the community close to probability theory a popular choice is the coupling from the past algorithm [33], which ensures perfect equilibration.Early coupling from the past simulations of the DWBC six-vertex model can be found in [34].But the so obtained number of samples is too small for a reliable check on KPZ universality.
In our simulations at the ASM point we also use the coupling from the past algorithm.The height functions, ϕ N , carry a natural partial ordering with a maximal (largest) and minimal (smallest) element.First one has to check that the stochastic dynamics respects this partial order.With the choice (23), such monotonicity is preserved provided a ≤ c, b ≤ c as can be checked by walking through all 16 possible local transitions S ↔ S ′ [35].Particularly monotonicity holds for ASM.Given such input one evolves the maximal and minimal height configuration under the same space-time noise from time −T , T > 0, to time 0. If at time 0 both height functions agree, one has obtained a valid sample of the equilibrium distribution.If the height functions disagree, one has to rerun, starting for example from time −2T , thereby retaining the noise during [−T, 0] from the prior run.If there is still disagreement the procedure has to be repeated.The prescribed time T has to be chosen with care.If too large, the computing time after coalescence is wasted.If too small, one has to restart the simulation leading to a considerable overhead in computing time.To make a reliable guess, we use the statistics of the time of coalescense for moderate system sizes by following the volume difference between the two extremal configurations.Empirically, when doubling the system size the number of steps for coalescence has to be increased by a factor of 4.5 approximately.For example, it takes roughly 10 6 Monte-Carlo steps (consisting of two half steps for the even and odd sublattice each) for system size N = 510 to have coalescence with probability > 99.999% rendering a negligible overhead for extending and repeating an almost identical simulation.
There are other schemes, which we have not tried.One option is called Gibbsian resampling.One observes that conditioned on its two neighboring level lines, a given level line has an explicit Gibbsian distribution.Thus one sequentially equilibrates randomly chosen level lines with respect to a frozen background.Superficially similar approaches  are directed-loop Monte Carlo [36], loop-cluster update [17] and a full lattice multi-cluster algorithm [37].

Asymptotics for the ASM facet
We generated samples for system sizes N = 2, 6, 14, 30, 62, 126, 254, 510, 1022, 2046 at the ice point.Our multi-spin implementation with 64 bit words produces 64 independent height functions in a single run.We recorded the shape of the disordered region for all four corners.For large N these samples are essentially independent.In our implementation one run consumed roughly (N/64) 4 seconds on a single core (as of 2018) leading to 256 almost independent corner samples (Mersenne twister was used as pseudorandom number generator).In our numerical data analysis, for comparison we also display results for the free fermion point a 2 = b 2 = 1 2 , which allows for well-known simplifications.One first maps the DWBC six-vertex model to the Aztec diamond, see for example the discussion in [7], which is a dimer tiling of a square rotated by π/4 and width 2 at the corners.The shuffling algorithm ensures that the steady state at size N is transformed to the steady state at size N + 1. Secondly under the shuffling algorithm the microscopic facet edge evolves under its own MC dynamics known also as a particular corner growth model, which can be viewed as a discrete time version of TASEP [16].Thus, instead of simulating a two-dimensional system, one can directly update h N from h N −1 , with an effort only proportional to N .Hence compared to the ASM model much larger system sizes can be reached for Aztec.Sample sizes for the corner growth model are N = 2 k , k = 2, . . ., 16, and complexity is optimally O(N 2 ).Our straightforward multi-spin implementation took roughly 3 × 10 −5 N 2.5 seconds for 10 6 realizations with 64 samples each.Figure 6 displays the number of samples used for the two scenarios.Total single-core computation time is more than 5 years for ASM and less than a year for Aztec.In actual fact, up to 100 CPU cores worked in parallel.
As a first test, the dynamical exponent 1/3, theoretically predicted in (7), is determined numerically from the variance of the edge function.For the Aztec diamond the corresponding logarithmic slope 2/3 has been rigorously proven in [16].For ASM the same exponent is confirmed in Figure 7 with high precision.Using the known macroscopic circle shape h ma , first determined in [38], together with the results from ( 10) and ( 19), the model dependent parameters for the macroscopic shape and local fluctuations for the Aztec point are obtained as Figure 8 displays the MC result for the probability distributions of the random variable h N (0) in the Aztec case for different values of N , transformed to the KPZ scale and compared to the limiting GUE Tracy-Widom distribution F ′ 2 (s) [39].At the ASM point the macroscopic shape determined in [11] and the fluctuations derived with (10) and ( 19) are Figure 9 shows the collapse of probability distributions for h N (0) in the ASM case, this time on the scale as predicted by KPZ scaling theory, and again compared to the Tracy-Widom density F ′ 2 (s).Note that in the plots of Figure 8 and 9 there is no free fitting parameter.So far we took into account only the fluctuations of the facet edge along the diagonal.More generally, in the limit of large N , the scaling theory predicts for each tile the probability to be in the frozen region of the lower right corner, compare with Figure 1.One expects the deviation from the properly rescaled distribution function F 2 to be small, In Figures 10 and 11 we present for Aztec and ASM, respectively, the numerical value of the difference in (28) as a function of k and h.For various values of N plotted are the results in the range |k| ≤ 0.9h.Note that for the extreme values of k, corresponding to | k N | close to 1 2 , the convergence to zero is particularly slow.For Aztec the convergence to 0 as N → ∞ has been established.The error is expected to be of order N −1/3 .Although the qualitative plot in Figure 11 is less conclusive, we still conjecture a similar behaviour for the case of ASM.
Fixing the origin of k, h, N with precision higher than order N 1/3 is somewhat arbitrary for h and N , while the origin of k can be determined by symmetry.We hope that by appropriately choosing these offsets or a more complicated correction term of order 1 (size of a tile) the convergence rate can be improved.As known from other models [40,41] convergence might be even of order N −2/3 .

2 Figure 3 :
Figure 3: One realization of the level lines for N = 2 for prescribed boundary tiles.For better visualization the tiles themselves are not shown. AF

Figure 4 :
Figure 4: Phase diagram of the six-vertex model with DWBC.

Figure 5 :
Figure 5: An example of a plaquette S with update S ′ .

Figure 6 :
Figure 6: ASM-Aztec number of samples vs size.

Figure 7 :
Figure 7: Dynamical scaling of the variance.

Figure 8 :
Figure 8: Data collapse of rescaled distributions for Aztec diamond.

Figure 9 :
Figure 9: Data collapse of rescaled distributions for ASM.