Finite size scaling in the dimer and six-vertex model

We present results of the Monte-Carlo simulations for scaling of the free energy in dimers on the hexagonal lattice. The traditional Markov-chain Metropolis algorithm and more novel non-Markov Wang-Landau algorithm are applied. We compare the calculated results with the theoretical prediction for the equilateral hexagon and show that the latter algorithm gives more precise results for the dimer model. For a non-hexagonal domain the theoretical results are not available, so we present the numerical results for a certain geometry of the domain. We also study the two-point correlation function in simulations of dimers and the six-vertex model. The logarithmic dependence of the correlation function on the distance, which is in accordance with the Gaussian free field description of fluctuations, is obtained.


Introduction
The six vertex model was introduced as the model of two-dimensional ice on a square lattice [1]. The oxygen atoms occupy the vertices of the lattice, the hydrogen atoms can be located on the edges close the the oxygen atom, but there can be only one hydrogen atom on each edge. The position of the hydrogen atom on the edge is then indicated by an arrow leading to the oxygen atom of the molecule. Thus there are exactly two incoming and two outgoing arrows at each vertex. The model was proposed as a toy model to study ice melting and crystallization, then became an important example of strong dependence of the scaling behavior on the boundary conditions [2,3] and a test-bed for the transfer matrix method for the exact solution of lattice models [4].
The dimer model appeared as a very simplified model of solutions [5]. The molecules are represented by the rigid tiles on a lattice. The partition function was computed on various lattices and with different boundary conditions [6,7,8].
The six vertex and dimer models are the integrable lattice models of statistical physics. Since they allow exact solutions only for a few particular cases [4,9], both the models and their generalizations are under an active theoretical [10,11,12] and numerical investigation [13,14].
Specifically, the well-known limit shape phenomenon [15] was discovered for these models [16] and a connection with the theory of random matrices was established [17].
For the studied models, the height function can be defined under some conditions. In the scaling limit, this function is fixed ("frozen") outside of an analytical curve which is usually called "the Arctic circle" [18], though its shape depends on the particular model. Inside this curve, the height function weakly converges to a certain surface. The fluctuations around this surface are described by the Gaussian free field. Although the limit shape phenomenon [16] is studied by the advanced methods of algebraic geometry, the analytical description of the limit curve is obtained only for a few particular cases, such as domino tiling of the so-called "Aztec diamond" [17,8]. From the point of view of the probability theory, limit shapes in the dimers and the six-vertex model are related to the Tracy-Widom distribution [19] of the random matrix eigenvalues.
Scaling behavior of the free energy depends on the geometry of the model and the area of the unfrozen domain. In this report, we present results of the Monte-Carlo simulations for scaling of the free energy in dimers on the hexagonal lattice. The traditional Markov-chain Metropolis algorithm [20] and more novel non-Markov Wang-Landau algorithm [21,22] are applied for the study. We compare the calculated results with the theoretical prediction for the equilateral hexagon and show that the latter algorithm gives more precise results for the dimer model. For a non-hexagonal domain the theoretical results are not available, so we present the numerical results for a certain geometry of the domain.
We also study the two-point correlation function in simulations of dimers and the sixvertex model. We observe the appearance of the limit shape. The logarithmic dependence of the correlation function on the distance, which is in accordance with the Gaussian free field description of fluctuations, is obtained.

Lattice models 2.1. The dimer model
The configurations of the dimer model are perfect matchings of a bipartite graph with some choice of weights. We consider coverings of the hexagonal lattice consisting of the subsets of lattice edges such that every vertex is the endpoint of exactly one edge.
We can draw a rhombus on a dual lattice around each edge in the configuration. The picture of "cubes in the corner" presented in the Fig. 1 is obtained. Let us write on the top of each uppermost cube the height of its column of cubes. Looking at this picture from the top, we obtain a height function defined on the rectangular domain of the square lattice. Let us define the sizes M , N , and K of the sides of the hexagon. The above description can be formalized by setting the non-negative numbers up to K in the boxes of the rectangular M × N table so that a value in each box is not greater than values in the adjacent upper and left boxes The weight of a particular configuration is given by the exponent of the volume of all cubes or by a sum of the height function values: We set Boltzman constant equal to 1, then the partition function is where q = exp (−1/T ) .
In general the partition function of the dimer model is given by the determinant of the Kasteleyn matrix [7,8,6], that can be seen as a discrete Dirac operator on a lattice [23]. For this particular case, the partition function is given by the classical Macmahon combinatorial formula [24] The scaling limit is achieved on the infinite lattice as T → ∞. To study the scaling behavior using Monte-Carlo simulations on the finite lattices we will consider temperatures proportional to the lattice size and, then, extrapolate them to infinity.

The six vertex model
The six vertex model is introduced as the model of two-dimensional ice on a square lattice [1,4]. Each edge of the lattice obtains an orientation, an arrow, in such a way that for each vertex there are exactly two incoming and two outgoing edges, so-called "ice rule". The six types of possible vertices are presented in Fig. 2.
conf for each vertex takes a predefine value according to Fig. 2. These Boltzmann weights are the matrix elements of the R-matrix, which is the key ingredient for the integrability and satisfy the Yang-Baxter equation [4,9].
The behavior of the six-vertex model strongly depends on the boundary conditions [2,3]. We consider the model with the so-called "domain wall boundary conditions". Considering the domain as a vertical rectangle, it means that the oriented edges are entering the domain on the top and bottom borders and leaving the domain on the left and right borders.
We can represent the model with such conditions as a set of non-intersecting paths as follows. The paths start at the top border, can go along the arrows only downward or rightward. Each edge can belong to at most one path. The paths end on the right border. These paths are regarded as level lines of the height function. So we can represent the model in simulations by the rectangular table of height function values h ij .
The universality class of the scaling behavior is determined by the value of the combination We study the case ∆ = 0 when the limit shape phenomenon takes place and fluctuations are described by the Gaussian free field [25,26].

Metropolis algorithm
The classical Metropolis algorithm based on the Markov chains can be applied to simulate both the dimer and the six vertex models [13,14]. The Markov chain is used to generate a random sequence of the model states. The average of some observable over these states tends to the expectation value as a number of states grows. The crucial requirement for the correctness of the algorithm is the detailed balance condition for the transition probabilities.
The following naive algorithm was implemented for the dimer model: (i) choose a random box (i, j) in a table and change the height function value ∆h = ±1 randomly; (ii) if such a change does not break the condition (1) for the height function, we change the height by ∆h with the probability exp (−∆h/T ), otherwise we stay in the same configuration; (iii) turn to the first step. Averaging over the configurations obtained on the second step of the algorithm, we obtain estimates for the thermodynamic observables.
Note that the detailed balance condition holds since the probability to choose the transition from a configuration to another one with one more cube is the same as the probability to choose a transition to a configuration with one less cube. The described algorithm can be optimized by saving a list of boxes where an addition or decrease of height function is allowed. However, one must preserve the detailed balance condition and should choose a box on the first step with the probability equal to the ratio of possible changes in current and new configurations. Note also that the Metropolis algorithm suffers from the critical slowdown. This problem is usually avoided by the use of Wolff [27] or Swendsen-Wang [28] cluster algorithms. In our case, we cannot apply cluster algorithms, since we need to preserve the condition (1) and similar one for the six vertex model. The free energy, F = −T log Z, is unobtainable for any particular temperature because the Metropolis algorithm does not calculate the partition function, Z. We have to use the relation The entropy S can be obtained in simulations by a numerical integration of the heat capacity: Therefore, we applied the Metropolis algorithm for various temperatures T i and obtained the values of observables C(T i ), E(T i ) . We calculated the entropy as a sum S(T i ) = i j=1 C(T j )(T j − T j−1 )/T j and, as a result, computed a list of free energies F (T i ). The Metropolis algorithm can be applied to the six vertex model [13] in a similar way.

The Wang-Landau algorithm
We use the Wang-Landau algorithm [21,22] to simulate the energy distribution ρ(E) = e g(E) of the dimer model. The realization of this algorith is as follows [29,30]. The energy range is split into some number of intervals, which can coincide with the number of discrete energies. The algorithm starts with a random configuration (random height function which satisfies the condition (1)), an empty array of logarithms of energy densities g(E 1 ), . . . , g(E s ), an empty visitation histogram n(E 1 ), . . . , n(E s ) and some initial value (usually 1) of constant a. On each step of the algorithm, a box (i, j) and a change of the value of the height function ∆h at this box are chosen randomly. A new configuration with the changed energy value h ij + ∆h is accepted with the probability e g(Enew)−g(E old ) if the condition (1) is satisfied. The visitation number n(E) is increased by 1 and g(E) is increased by a. Otherwise, the current configuration is kept, n(E) and g(E) remain unchanged. This procedure is repeated until the visitation histogram becomes relatively flat, i.e. min(n(E))/ max(n(E)) > 0.8. Then, the value of a is divided by 2, the histogram is emptied and the next step of the algorithm begins.
In the limit, the distribution ρ(E) becomes stationary. In practice, about 25 such steps are enough to do in simulations. Having the energy distribution ρ(E), the partition function is given by Calculation of the partition function is an important advantage of the Wang-Landau algorithm over the Metropolis [20] and Wolff algorithms [27]. Having the partition function, we obtain the free energy F = −T log Z. An example of simulation results for the dimer model on the equilateral hexagonal lattice with the side length equal to 10 is presented in Fig. 3.2 (left panel).

Free energy scaling in the dimer model
Let us fix a domain D and cover it by a lattice with the lattice constant ε. Compute the logarithm of the partition function for the dimer model and consider its dependence on ε.
We made preliminary computations, expanding formula (2) in ε. There are also results in the literature for scaling of the determinant of a discrete Laplace operator [31,32] and partial results on Dirac operator (Kasteleyn matrix) [23]. From these results we conjecture the scaling behavior ln(Z D,ε ) = 1 where f i are the fitted parameters. We will present our computations in a separate publication.
Here we demonstrate the numerical evidence and simulation results that support this conjecture. Consider the equilateral hexagon with a side M as the domain, then ε = 1 M . The scaling limit takes place on an infinite lattice with q = exp (−1/T ) tending to one, i.e. T → ∞. We need to consider temperatures proportional to the lattice size, T = τ M , where τ is a dimensionless factor. So in the extrapolation to infinity, we keep τ to be constant. We divide the formula (7) by the volume, M 2 , and multiply by −T = −τ M . Scaling behavior of the free energy density becomes as following: First, we test the expansion (7) by the Macmahon formula (2). Using it for lattice sizes ranging from 1 to 20 and doing a fit by Eq. (8), we obtain the scaling parameters The uncertainty of the last decimal digit is given in brackets. The fit is shown in Fig. 4 and the results are statistically significant for f 0 , f 2 , f 3 . Therefore, our scaling behaviour, Eq. (7), is supported by the numerical data.
The uncertainty of the last decimal digit is given in brackets. Although the calculated results are less precise, the Monte-Carlo results contain the theoretical predictions within the uncertainties. The corresponding fit of the free energy density is shown in the right panel of Fig. 4.
Using the Metropolis algorithm on the lattices of size from 3 to 18, we obtained similar but even less precise values This loss of precision is due to the critical slowdown that cannot be avoided since we are interested in the critical scaling behavior. From the described comparisons, we see that the Wang-Landau algorithm is more efficient for study scaling behavior. We can apply it to study the scaling behavior (7) on domains with more complex geometry, since there is no known formula for the partition function on such domains.
So, the region i < M/2 for j < M/2 is forbidden. The limit shape and the two-point correlation function for a domain with a cutted square region is shown in Fig. 4. A dependence of the free energy density on the lattice constant is presented in the left panel of Fig. 4. On the right panel, we show a behavior of the heat capacity on the temperature.
We see that f 0 strongly depends upon the domain geometry. To draw some conclusion about f 1 , f 2 , f 3 , more extensive simulations are required.
5. Two-point correlation functions and the limit shape in the dimer model Scaling behavior of the dimer model inside the "the Arctic circle" is described by an effective free field theory -Gaussian free field (massless free boson) [33]. In the scaling limit height function converges to a free bosonic field. Thus, two-point height-height correlation function h(x)h(y) in the dimer model inside "the Arctic circle" demonstrates logarithmic behavior for x = y [8]: This logarithmic behavior corresponds to the liquid phase of the model. Phase diagram and limit shape phenomenon is discussed in the paper [16]. The simulated data is presented in Fig. 5. Fixing, for example, y and fitting the simulated data by a particular formula a log x−y M , as shown in Fig. 5 we obtain that a = −0.0567(3). If we now choose the normalization of the correlation function in such a way that h(y) 2 = 1, we obtain the coefficient −0.147 (8) which is sufficiently close to − 1 2π . In the paper [34] it was proven that height fluctuations around the limit shape in the scaling limit are described by the Gaussian free field. Our simulation result is in agreement with this theoretical observation.  6. Two-point correlation functions and the limit shape in the six vertex model A logarithmic behavior, analogous to the dimer model, is observed for the six vertex model in case of ∆ = 0. Fitting the simulated data by the formula a log |x − y| + b with fixed y and a,b as the fitted parameters, as shown in Fig. 8, we obtain the fitted values listed in Tab b has logarithmic dependence on M . Fitting b(M ) = −a log M , we obtain a = −0.165 ± 0.003, which is in agreement with the values of a in Tab. 1. Thus, we confirm a logarithmic behaviour of the correlation function inside the non-frozen region, that is predicted by the Gaussian free field description.

Conclusion
In this paper, we apply the Wang-Landau algorithm as well as the Metropolis algorithm to study the scaling behavior of the free energy in the dimer model. We compare a precision of the results obtained by these methods, as well as confront them to the theoretical prediction based on the combinatorial formula. We conclude that Wang-Landau algorithm gives more precise results and can be used for extracting the scaling parameters of the free energy for a domain with a complex geometry. We also present the two-point correlation function and illustrate the limit shape phenomenon for the dimer and six vertex model. A logarithmic dependence of the correlation function is in agreement with the Gaussian free field description of fluctuations around the limit shape.