Finite-size scaling of free energy in the dimer model on a hexagonal domain

TL;DR

This paper analyzes the finite-size scaling behavior of free energy in the dimer model on a hexagonal lattice, deriving asymptotic expansions and discussing their physical significance and universality.

Contribution

It provides the first detailed asymptotic expansion of free energy for the dimer model on a hexagonal domain, including various boundary conditions and weight configurations.

Findings

Derived asymptotic expansion of free energy as lattice mesh tends to zero.

Identified universal features of the free energy expansion coefficients.

Discussed physical interpretation and universality of the results.

Abstract

We consider dimer model on a hexagonal lattice. This model can be seen as a "pile of cubes in the box". The energy of configuration is given by the volume of the pile and the partition function is computed by the classical MacMahon formula or, more formally, by the determinant of Kasteleyn matrix. We use the expression for the partition function to derive the scaling behavior of free energy in the limit of lattice mesh tending to zero and temperature tending to infinity. We consider the cases of finite hexagonal domain, of infinite height box and coordinate-dependent Boltzmann weights. We obtain asymptotic expansion of free energy and discuss the universality and physical meaning of the expansion coefficients.

Tables

Table 1. Table 1: Values of the fitted parameters and analytically computed expansion coefficients for various weight functions φ 𝜑 \varphi and values of a , b 𝑎 𝑏 a,b .

	$\phi \left(z\right)$	$a,b$	$f_0$	$f_1$	$12ab{f}_2$	$f_3$
analytic	$\left(2+\cos \left(z\right)\right)/3$	1,3	$0.472206688$	$0$	$1$	$-0.043883000$
numeric	$\left(2+\cos \left(z\right)\right)/3$	1,3	$0.472206693$	$0.000000730$	$1.000730689$	$-0.043827958$
analytic	$\left(2+\cos \left(z\right)\right)/3$	2,3	$0.235922467$	$0$	$1$	$-0.030398000$
numeric	$\left(2+\cos \left(z\right)\right)/3$	2,3	$0.235922467$	$0.000000044$	$1.000088528$	$-0.030394694$
analytic	$a+z/2$	1,3	$0.097288596$	$0$	$1$	$-0.033688300$
numeric	$a+z/2$	1,3	$0.097288593$	$0.000000848$	$1.000848267$	$-0.033624441$
analytic	$a+z/2$	2,3	$0.032804447$	$0$	$1$	$-0.015094000$
numeric	$a+z/2$	2,3	$0.032804447$	$-0.000000002$	$0.999996808$	$-0.015094162$

Equations

$$Z=\sum _{\mathrm{conf}}\prod _{e\in \mathrm{conf}}w(e)=|\det \mathcal{K}|.$$

$$(Df)(v)=\sum _u\mathcal{K}(v,u)f(u).$$

$$h_{ij}\le h_{i-1,j},h_{ij}\le h_{i,j-1}.$$

$$E[conf]=\sum _{i,j}{h}_{ij}=\mathrm{Volume}$$

$$Z=\sum _{conf}{e}^{-\frac{E[conf]}{T}}=\sum _{conf}{q}^{\mathrm{Volume}[conf]},$$

$$Z[M,N,K,q]=\prod _{i=1}^M\prod _{j=1}^N\prod _{k=1}^K\frac{1-q^{i+j+k-1}}{1-q^{i+j+k-2}}.$$

$$\mathcal{K}(w,b)=q^{\Re w+\Im w}\mathrm{if}\Im w=\Im b$$

$$\mathcal{K}(w,b)=1\mathrm{if}\Im w\ne \Im b.$$

$$f=-\frac{1}{V}\ln Z(M,N,K,q).$$

$$V=2(MN+NK+MK)=2(ab+bc+ca){\epsilon }^{-2}.$$

$$\begin{gathered} f=-\frac{1}{V}\ln Z= \\ =-\sum _{i=1}^M\sum _{j=1}^N\sum _{k=1}^K\frac{1}{V}\left[\ln \left(1-e^{-\epsilon (i+j+k-3)}{e}^{-2\epsilon }\right)-\ln \left(1-e^{-\epsilon (i+j+k-3)}{e}^{-\epsilon }\right)\right]. \end{gathered}$$

$$f=-\sum _{i=1}^M\sum _{j=1}^N\sum _{k=1}^K\sum _{n=1}^{\infty }\frac{1}{V}\frac{1}{n}\left[-e^{-n\epsilon (i+j+k-3)}{e}^{-2n\epsilon }+e^{-n\epsilon (i+j+k-3)}{e}^{-n\epsilon }\right].$$

$$f=-\sum _{n=1}^{\infty }\frac{1}{V}\frac{1}{n}{e}^{-n\epsilon }\left(1-e^{-n\epsilon }\right)\sum _{i=1}^M\sum _{j=1}^N\sum _{k=1}^K{e}^{-n\epsilon (i+j+k-3)}.$$

$$f=-\sum _{n=1}^{\infty }\frac{1}{V}\frac{1}{n}\frac{e^{-n\epsilon }}{{\left(1-e^{-n\epsilon }\right)}^2}\left(1-e^{-na}\right)\left(1-e^{-nb}\right)\left(1-e^{-nb}\right).$$

$$\chi (z)=e^{-z}{\left(\frac{z}{1-e^{-z}}\right)}^2,$$

$$H_n=\left(1-e^{-na}\right)\left(1-e^{-nb}\right)\left(1-e^{-nc}\right).$$

$$\chi (z)=1-\frac{z^2}{12}+\frac{z^4}{240}+\mathcal{O}(z^6).$$

$$f=-\frac{1}{2(ab+bc+ca)}\sum _{n=1}^{\infty }\frac{\chi (n\epsilon )}{n^3}{H}_n.$$

$$\begin{gathered} {f|}_{\epsilon \to 0}=\frac{1}{2(ab+bc+ca)}[{\mathrm{Li}}_3(e^{-a})+{\mathrm{Li}}_3(e^{-b})+{\mathrm{Li}}_3(e^{-c})-{\mathrm{Li}}_3(e^{-a-b}) \\ -{\mathrm{Li}}_3(e^{-b-c})-{\mathrm{Li}}_3(e^{-a-c})+{\mathrm{Li}}_3(e^{-a-b-c})-\zeta (3)]. \end{gathered}$$

$${\frac{\partial f}{\partial \epsilon }|}_{\epsilon \to 0}=0.$$

$$\frac{{\partial }^{2}f}{\partial {\epsilon }^2}=-\frac{1}{2(ab+bc+ca)}\sum _{n=1}^{\infty }\frac{{\chi }^{\prime \prime }(n\epsilon )}{n}{H}_n.$$

$$\frac{{\partial }^{2}f}{\partial {\epsilon }^2}=-\frac{1}{2(ab+bc+ca)}\left(\sum _{n=1}^{\infty }\frac{{\chi }^{\prime \prime }(n\epsilon )}{n}+\sum _{n=1}^{\infty }\frac{{\chi }^{\prime \prime }(n\epsilon )}{n}(H_n-1)\right).$$

$$\begin{gathered} {\chi }^{\prime \prime }(0)\sum _{n=1}^{\infty }\frac{H_n-1}{n}= \\ ={\chi }^{\prime \prime }(0)\sum _{n=1}^{\infty }\frac{1}{n}\left(-e^{-an-bn-cn}+e^{-an-bn}+e^{-an-cn}+e^{-bn-cn}-e^{-an}-e^{-bn}-e^{-cn}\right)= \\ =\frac{1}{6}\ln \left(\frac{(e^a-1)(e^b-1)(e^c-1)(e^{a+b+c}-1)}{(e^{a+b}-1)(e^{b+c}-1)(e^{a+c}-1)}\right). \end{gathered}$$

$${\chi }^{\prime \prime }(z)=e^{-z}\frac{e^{2z}\left(z^2+4e^z\left(z^2-1\right)+e^{2z}\left(z^2-4z+2\right)+4z+2\right)}{{\left(e^z-1\right)}^4}=e^{-z}\xi (z).$$

$$\sum _{n=1}^{\infty }\frac{{\chi }^{\prime \prime }(n\epsilon )}{n}=\sum _{n=1}^{\infty }\frac{e^{-n\epsilon }}{n}\xi (n\epsilon )=\sum _{n=1}^{\infty }\frac{e^{-n\epsilon }}{n}\xi (0)+\epsilon \sum _{n=1}^{\infty }{e}^{-n\epsilon }\frac{\xi (n\epsilon )-\xi (0)}{n\epsilon }.$$

$$Q(z)=\frac{\xi (z)-\xi (0)}{z}.$$

$$\epsilon \sum _{n=1}^{\infty }{e}^{-n\epsilon }Q(n\epsilon )=\epsilon \sum _{n=1}^{\infty }{e}^{-n\epsilon }\frac{\xi (n\epsilon )-\xi (0)}{n\epsilon }\approx {\int }_0^{\infty }{e}^{-z}\frac{\xi (z)-\xi (0)}{z}dz={\int }_0^{\infty }{e}^{-z}Q(z)\mathrm{dz}.$$

$$\frac{{\partial }^{2}f}{\partial {\epsilon }^2}\approx -\frac{1}{2(ab+bc+ca)}\left[\frac{1}{6}\ln \left(\frac{(e^a-1)(e^b-1)(e^c-1)(e^{a+b+c}-1)}{(e^{a+b}-1)(e^{b+c}-1)(e^{a+c}-1)}\right)-\frac{1}{6}\ln \epsilon +{\int }_0^{\infty }{e}^{-z}\frac{\xi (z)-\xi (0)}{z}dz\right].$$

$$f\approx f_0+f_1\epsilon +f_2{\epsilon }^2\ln \epsilon +f_3{\epsilon }^2+\mathcal{O}({\epsilon }^3),$$

$$\begin{gathered} f_0=\frac{1}{2(ab+bc+ca)}[{\mathrm{Li}}_3(e^{-a})+{\mathrm{Li}}_3(e^{-b})+{\mathrm{Li}}_3(e^{-c})-{\mathrm{Li}}_3(e^{-a-b}) \\ -{\mathrm{Li}}_3(e^{-b-c})-{\mathrm{Li}}_3(e^{-a-c})+{\mathrm{Li}}_3(e^{-a-b-c})-\zeta (3)], \end{gathered}$$

$$f_1=0,$$

Full text

Introduction

The dimer model appeared in 1937 in an attempt to extend a statistical theory of perfect solutions in chemistry to the case of liquid mixtures with molecules of two very distinct sizes [1]. The molecules were represented by the rigid tiles on a lattice and the number of tilings was approximately estimated. But an exact computation was not accessible at the time.

In 1961, Kasteleyn [2], Temperley and Fisher [3] represented the partition function of the dimer model as a Pfaffian of the signed adjacency matrix (‘‘Kasteleyn matrix’’) thus allowing the computation of the number of tilings and of free energy scaling limit. This result was very elegantly used by Fisher to solve Ising model [4] and by Fan and Wu [5] to compute free energy for a certain case of eight-vertex model.

Further studies of dimer models revealed the connection to the theory of alternating matrices [6, 7]. Later, the well-known limit shape phenomenon [8] was discovered for dimer models. First, the ‘‘Arctic circle’’ theorem was proven for domino tilings of the domain in the form of ‘‘Aztec diamond’’  [9]. Then similar result was obtained for a hexagonal domain on the hexagonal lattice [10]. Soon the connection of these results with the theory of random matrices was established [11]. The papers [12, 13] present a detailed exposition of the limit shape phenomenon in the dimer models.

Dimer models are the integrable lattice models of statistical physics that are now under an active theoretical [14, 15] and numerical investigation [16]. Computation of correlation functions is a common problem for all vertex models [17] as well as for dimer models. Another problem of great interest is the study of limit shapes in various cases [18, 19].

Configurations of dimer model on a hexagonal lattice are in one-to-one correspondence with the configurations of five-vertex model that appears for certain choice of parameters in the well-known six-vertex model [20, 21].

Study of dimer models on various lattices and domains led to interesting connections with the geometry of curved manifolds and with spectra of discrete and continuous Dirac and Laplace operators [22, 23]. Scaling limit of dimer model is proven to be described by a Gaussian free-field theory [24], but finite-size corrections were not considered previously. These corrections are important to close the gap between numerical simulations and theoretical results.

We consider a particular case of the dimer model on a hexagonal domain of hexagonal lattice, that can be seen as a ‘‘pile of cubes in the corner of the box’’. The energy of configuration is the total number of cubes. For this particular case we use an exact combinatorial formula for the partition function to derive the expressions for scaling limit of free energy and first three terms of the finite-size corrections. We show that the first term is identically zero. The second term which encapsulates logarithmic dependence on the mesh size is connected with the central charge of the effective field theory and with the geometry of limit shape. The third term has a geometry-dependent contribution that is written explicitly in terms of elementary functions and a universal contribution, which is constant in all the cases that we consider.

We then consider the box of infinite height and the case of coordinate-dependent Boltzmann weights. In these cases we are also able to write explicit expressions for finite-size corrections. We demonstrate the universality of the logarithmic term and the constant contribution to the third term.

Our results are supported by numeric checks and Monte-Carlo simulations presented in our previous publication [25].

1 Model definition

The configurations of the dimer model are perfect matchings (sets of non-touching edges, covering all the vertices) on some graph ${\cal G}$ with some choice of weights $\omega(e)$ on the edges. The model is solvable on the bipartite graphs, i.e. the partition function can be computed if the weights are introduced in such a way that for each face bounded by 0 mod 4 edges there is an odd number of negative edge weights and each face with 2 mod 4 edges has an even number of of negative edge weights. Then the signs of the edge weights form a so-called ‘‘Kasteleyn orientation’’ on graph, the weighting is called ‘‘Kasteleyn weighting’’ [24, 13].

For a bipartite graph ${\cal G}$, color the vertices black and white in such a way that all the vertices adjacent to the black one are white. Denote by $B,W$ the sets of black and white vertices and by $b,w$ the elements of these sets.

The weights can be encoded as the ‘‘Kasteleyn matrix’’ – weighted, signed adjacency matrix ${\cal K}$ with the matrix elements ${\cal K}(w,b)$ equal to the weight of the edge $w\to b$: ${\cal K}(w,b)=\omega(w\to b)$.

Then the partition function is equal to the absolute value of the determinant of the Kasteleyn matrix[2, 3]:

[TABLE]

Kasteleyn matrix defines a discrete Dirac operator $D$, the action of $D$ on a function $f$ defined on vertices is given by:

[TABLE]

Kenyon [22, 23] and others [26] considered asymptotics of the determinants of the discrete Dirac and Laplace operators, the problem that, as can be seen from the above, is very close to the scaling of the free energy. But the finite size corrections to the free energy scaling were not computed.

We consider coverings of the hexagonal domain on 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\times N$ table so that a value in each box is not greater than values in the adjacent upper and left boxes

[TABLE]

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:

[TABLE]

We set Boltzmann constant equal to 1 and choose system of units in such a way that there is no coupling constant in the expression for energy. Then the partition function is

[TABLE]

where $q=\exp\left(-1/T\right)$ .

For this particular case, the partition function is given by the classical Macmahon combinatorial formula [27]

[TABLE]

MacMahon formula is obtained for the following definition of the Kasteleyn matrix. Embed the hexagonal lattice in the complex plane $\mathbb{C}$ in such a way that some edges are parallel to the real line with and corresponding vertices have coordinates with integer real and imaginary parts. Then take

[TABLE]

The free energy per site is defined as*∗* 111For convenience we omit the factor $\frac{1}{T}$ in the usual definition of the free energy.

[TABLE]

Here $V$ is the number of vertices, it is twice the number of dimers and twice the number of cube faces:

[TABLE]

We are interested in the scaling limit, combined with the thermodynamic limit, when $T\to\infty,$ and $M,N,K\to\infty$, such that ratios $\frac{M}{T}=a,\quad\frac{N}{T}=b,\quad\frac{K}{T}=c$ remain fixed. In what follows we use $\varepsilon=\frac{1}{T}$, which can be seen as the scale of the model, e.g. mesh size due to our choice of the system of units.

In the next section we compute the asymptotic expansion of the free energy $f$ in $\varepsilon$ and derive $\varepsilon$-independent closed expressions for the first several coefficients in this expansion.

2 The computation of the free energy asymptotic expansion

First we substitute MacMahon formula (4) into the free energy definition (1) and obtain

[TABLE]

We need to notice that for a finite $\varepsilon$ and each value of $i,j,k$ the logarithms have the form $\ln(1-x)$, where $0<x<1$. Thus we can represent the logarithms as a power series of $x$ and obtain

[TABLE]

Here we have an absolutely convergent series, so we can change the order of summation and rewrite the free energy as

[TABLE]

The triple sum factorizes into the product of three sums of the type $\sum_{i=0}^{M-1}e^{-n\varepsilon i}$, that are just the sums of geometric progression. So we have $\sum_{i=0}^{M-1}e^{-n\varepsilon i}=\frac{1-e^{-Mn\varepsilon}}{1-e^{-n\varepsilon}}$, and obtain

[TABLE]

Denote by $\chi(z)$ the function

[TABLE]

and by $H_{n}$ the product of $\varepsilon$-independent exponents

[TABLE]

The function $\chi(z)$ is smooth in 0, even and has the following Taylor series:

[TABLE]

Using the definition of the volume (7) the free energy is written as

[TABLE]

We are interested in the scaling behavior in $\varepsilon$ when $\varepsilon\to 0$. The function $f$ is regular in $\varepsilon$ for $\varepsilon>0$. So we need to compute $\left.f\right|_{\varepsilon\to 0}$ and the derivatives $\left.\frac{\partial f}{\partial\varepsilon}\right|_{\varepsilon\to 0}$, $\left.\frac{\partial^{2}f}{\partial\varepsilon^{2}}\right|_{\varepsilon\to 0}$.

For $\left.f\right|_{\varepsilon\to 0}$ we have $\chi(0)=1$ and the sum $\sum_{n=1}^{\infty}\frac{H_{n}}{n^{3}}$ is just a sum of polylogarithms $\mathrm{Li}_{s}(z)=\sum_{n=1}^{\infty}\frac{z^{n}}{n^{s}}$ of third order:

[TABLE]

Here Riemann zeta function appears as a particular value of polylogarithm $\mathrm{Li}_{3}(1)=\zeta(3)$.

The derivative of the function $\chi(z)$ is zero, $\chi^{\prime}(0)=0$, and the series $\sum_{n=1}^{\infty}\frac{\chi^{\prime}(n\varepsilon)}{n^{2}}H_{n}$ are convergent for a finite $\varepsilon$, thus we have for a derivative

[TABLE]

Now let us compute the second derivative

[TABLE]

Second derivative of $\chi$ is finite $\chi^{\prime\prime}(0)=-\frac{1}{6}$. But there is a difficulty with the limit $\varepsilon\to 0$ in this expression due to the poor convergence of the series. Let us rewrite the second derivative as a sum of two series:

[TABLE]

The second sum is now convergent, since $H_{n}-1$ decays exponentially in $n$. It can be presented as a combination of series for logarithms:

[TABLE]

First sum in the expression (19) can be expressed as follows. First, consider second derivative of $\chi(z)$, separate an exponent and denote what is left by $\xi(z)$:

[TABLE]

The sum is then rewritten

[TABLE]

The first term is series for the logarithm $\ln\left(1-e^{-\varepsilon}\right)\approx\ln\varepsilon$. Denote by $Q(z)$ the difference of the function $\xi$ value in $z$ and zero, divided by z:

[TABLE]

The function $Q(z)$ is analytic in $z$.

Then the second term can approximated by the integral:

[TABLE]

[TABLE]

Since we have logarithmic term in the second derivative, we have the behavior of free energy of the form

[TABLE]

as expected in general in two-dimensional critical systems [28]. Taking the derivatives of such expression by $\varepsilon$, we get $f_{0}=f(0)$, $f_{1}=\left.\frac{\partial f}{\partial\varepsilon}\right|_{\varepsilon\to 0}$, $\frac{\partial^{2}f}{\partial\varepsilon^{2}}=f_{2}(2\ln\varepsilon+3)+f_{3}$.

Thus we have the coefficients

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

3 Limit shape and physical meaning of the expansion coefficients

In the paper [25] we have presented results of Monte Carlo simulations using Metropolis and Wang-Landau algorithms that support the form of expansion (26). In particular we got $f_{1}=0.04\pm 0.04$ which is consistent with $f_{1}=0$.

The dimers in the scaling limit are in general described by a free-fermion theory [29], but on a hexagonal lattice (and on the graphs that admit Kasteleyn orientation) the model is reformulated in terms of height function, that can be seen as a bosonization of the theory [30].

In the scaling limit $\varepsilon\to 0$ so called ‘‘limit shape phenomenon’’ [9, 10] appears in the dimer model. The areas around the corners of the domain are ‘‘frozen’’ with height function values being fixed. An analytical ‘‘Arctic curve’’ delimits frozen regions from the region where the behavior is described by the effective free field theory [13, 31, 12]. Moreover, the height function in the scaling limit tends to certain analytical surface and fluctuations around this ‘‘limit shape’’ are described by a gaussian free field theory in a curved background of ‘‘limit shape’’. A confomal map can be used to send the internal region of the ‘‘Arctic curve’’ onto the upper half-plane.

The scaling behavior (26) of the logarithm of the partition function is generic in two-dimensional models [28].

Thus we can interpret first two terms $f_{0}$ and $f_{1}$ as a bulk and boundary free energies in the corresponding field theory. The term $f_{1}$ corresponds to the contribution of the ‘‘limit shape’’. Since we have a deterministic frozen boundary, we naturally have $f_{1}=0$.

The term proportional to the logarithm of the scale $\varepsilon$ is also universal [28] and should appear in all two-dimensional theories with boundary. In the paper [28] Cardy and Peschel argued that on a smooth manifold of a characteristic length $L$ with a smooth boundary such a term must have the following form:

[TABLE]

where ${\bf c}$ is central charge of the effective field theory and $\chi$ is the Euler characteristic of the manifold $\chi=2-2h-b$, where $h$ is the number of handles and $b$ is the number of boundaries. Naive use of the formula (31) with our result $\delta F=\frac{1}{12}$ would give ${\bf c}=\frac{1}{2}$, but it is not correct, since we need to consider the field theory in the curved background.

In the case of a corner singularity on the boundary with the angle $\Theta$ the logarithmic contribution to the free energy is changed to

[TABLE]

This formula was recently used by N. Allegra [32] to show that the central charge in the dimer model on the square lattice with a corner monomer is ${\bf c}=1$, contrary to some previous claims of ${\bf c}=-2$.

In the case of manifold with conic singularity with the semi-angle $\Theta$, the logarithmic contribution presented in [28] is:

[TABLE]

In our case we interpret $2(ab+bc+ca)$ as a volume of the domain, thus

[TABLE]

Using Cardy-Peschel formulas, presented above, $\delta F$ can be computed from value of the central charge for Gaussian free field ${\bf c}=1$, the curvature of metric, induced by the limit shape and conformal map from the internal region of the ‘‘Arctic curve’’ to the upper half-plane. Such a computation can be used to check our result on $f_{2}$ but is quite lengthy and will be presented in a separate publication.

The last term $f_{3}$ depends only on the shape of the domain through $a,b,c$ with a universal contribution that is equal to the integral of a function $e^{-z}\frac{\xi(z)-\xi(0)}{z}$. In the next sections we show that this constant appears due to the treatment of logarithmic divergency even for the coordinate-dependent Boltzmann weights.

4 Numeric checks

To check our result (26) numerically, we first need to evaluate the integral in $f_{3}$. The function $e^{-z}\frac{\xi(z)-\xi(0)}{z}$ is smooth and but some care is required when taking the limit $z\to 0$ numerically. Evaluating the integral numerically with Mathematica we get

[TABLE]

The asymptotic expansion (26) gives very good approximation for free energy. In Figure 2 we present exact and approximate values for $f$ and, on the right panel, demonstrate that inaccuracy is actually $\mathcal{O}(\varepsilon^{4})$.

5 Box of infinite height

Consider the case of cubes in the box of infinite height. In this case $K=\infty$ and $M,N$ are finite. Now the hexagonal domain of dimer configurations has side of infinite length and total number of dimers is infinite, so it is more natural to set the volume of the system $V=MN$.

Putting $K=\infty$ in the partition function, we see that all terms with $k$ are cancelled and we have

[TABLE]

The computation of free energy is almost exactly the same as in finite box case. We expand the logarithm in the free energy and then factorize the sum:

[TABLE]

Computing the sums of geometric progressions we obtain almost the same expression as (11), except that $H_{n}$ now has two multipliers $H_{n}=\left(1-e^{-na}\right)\left(1-e^{-nb}\right)$:

[TABLE]

The function $\chi(z)=e^{-z}\left(\frac{z}{1-e^{-z}}\right)^{2}$ is the same as in (12).

Doing the same computations as in section 2, we obtain following results for the coefficients $f_{0},f_{1},f_{2},f_{3}$:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Again in the Fig. 3 we see good numerical agreement between the asymptotic expansion (26) and the exact formula (37).

6 Box of infinite height with coordinate-dependent weight

Now let us consider the generalization to the case when $q$ depends upon the coordinate. It’s usually done by varying $q$ on diagonal slices numbered by $t$. We need proper setup to consider the scaling limit. Let $\varphi(t)$ be a continuous bounded smooth function for $t\in\left[-a,b\right]$. Denote by $q_{i}=e^{-\varepsilon\varphi\left(i\varepsilon\right)}$ for $i=(-M+1)\dots N-1$. We index diagonal slices in such a way, that top left corner of the $M\times N$ box has index $N-M$, top right corner – $N$, bottom left – $-M$ and bottom right – [math]. (See Fig. 4 ). The slices with the indexes $N$ and $-M$ are empty so we set $q_{N}=q_{-M}=0$. There can be at most $N+M-1$ non-empty slices.

The partition function for this and more general case of arbitrary inner shape of the box was presented in the paper [33]. In this case it contains products of weights $q_{i}$ for slices with negative and positive indexes:

[TABLE]

(see Theorem 2 in [33]). If we set $q_{i}=q$, we recover formula (36).

Now let us again consider free energy dependence on $\varepsilon$:

[TABLE]

The terms $f_{0}$ in the asymptotic expansion (26) can be easily obtained from this expression. To do so we need to approximate the sums by the integrals, using \sum_{k=0}^{i}\varepsilon\varphi(b-a-k\varepsilon)\approx\int_{0}^{i\varepsilon}\varphi(b-a-x)\mathrm{dx}\leavevmode\nobreak\and similar approximation for the outer sums:

[TABLE]

Substituting in this expression $\varphi(x)\equiv 1$ we immediately get formula (39), since polylogarithms are integrals of logarithms. Although the logarithm has a divergency in the point $y=z=0$, it disappears after double integration.

To compute the term $f_{1}$ we can use the following simple Lemma:

Lemma 1**.**

Let $a=M\varepsilon$ and $b=N\varepsilon$. For an integrable almost everywhere analytic function $f$ and analytic function $\varphi$ the sum is approximated by an integral with the corrections of order $o(\varepsilon)$:

[TABLE]

Proof.

Use Euler-Maclaurin formula to approximate the sums by the integrals with corrections. It is enough to consider linear in $\varepsilon$ terms. First use it for the sum in the argument of function $f$:

[TABLE]

For a double sum the following simple analogue of Euler-Maclaurin formula can be used. Let $G$ be an analytic function in the rectangle $[0,a]\times[0,b]$, then

[TABLE]

This formula can be easily derived by dividing the volume of integration into squares with side $\varepsilon$ and substituting Taylor series for $G$ into integrals.

Denote by $G_{\varepsilon}(y,z)$ the function of two arguments that appear in the double sum in (47):

[TABLE]

Use formula (48) to express the double sum in (47) as the integral plus correction:

[TABLE]

Now expand the function $G_{\varepsilon}(y,z)$ under integral into Taylor series in $\varepsilon$:

[TABLE]

On the other hand, the derivatives of $G_{\varepsilon}$ in $y,z$ are computed as

[TABLE]

Substituting equations (51),(52) into (50) and using formula (48) to change the double sum of the derivatives of $G_{\varepsilon}$ into double integral, we see that the linear correction in $\varepsilon$ is cancelled.

∎

Using this lemma for our free energy density $f$, which is analytic and $\varphi$ which is smooth, we instantly see that

[TABLE]

Next terms in the asymptotic expansion require more care, since we need to single out the contribution of the order $\varepsilon^{2}\ln\varepsilon$.

Thus, to derive the expressions for higher order terms of the expansion (26) we need to expand the logarithm, change the order of summation and rewrite the free energy in the same way as we did before:

[TABLE]

Now we again can factorize double sum into product of two sums:

[TABLE]

where we assume that $\sum_{k=1}^{0}\varphi(b-a-k\varepsilon)=0$ and $\sum_{l=1}^{0}\varphi(b-a+l\varepsilon)=0$. Note that for $\varphi(t)\equiv 1$ we recover formula (37).

Introduce following notations:

[TABLE]

The functions $r_{\mp}$ can be approximated using Euler-Maclaurin formula as:

[TABLE]

Denote by $\Psi_{n}^{\mp}$ the following functions:

[TABLE]

These functions are analytic in $\varepsilon$, moreover, the higher derivatives of $\Psi^{\mp}_{n}$ grow in $n$ not too rapidly. Thus, the sums of the series written in terms of these functions and their derivatives converge, as we will see below.

With these notations we can write the free energy as

[TABLE]

The first term $f^{(0)}$ can be computed explicitly as

[TABLE]

Denote by $H_{n}(\varepsilon)$ the numerator of this expression:

[TABLE]

Note that

[TABLE]

can be rewritten as

[TABLE]

where

[TABLE]

Now for $f^{(0)}$ we have

[TABLE]

where $\tilde{\chi}_{n}(z,\varepsilon)$ is defined as

[TABLE]

It is important to include $H_{n}(\varepsilon)$, since $H_{n}(\varepsilon)$ contains decaying exponents in $n$.

Taking the derivatives of $f^{(0)}$ in $\varepsilon$ we get

[TABLE]

where the function $\chi(z)$ is the same as in the equation (12)

[TABLE]

and

[TABLE]

The last sum in the last line of the equation (68) is the same as in the section 5, thus we can compute it:

[TABLE]

where $Q(z)$ was defined in the equation (23).

Thus the contribution of $f^{(0)}$ to the asymptotic expansion (26) consists of the term

[TABLE]

that is cumbersome but does not depend on $\varepsilon$ and can be computed explicitly and the expression (71), which gives us logarithmic behavior in $\varepsilon$.

Now we need to consider other terms in the expression (60). For $f^{+}$ we have

[TABLE]

The sum over $j$ can be approximated by an integral using Euler-Maclaurin formula as

[TABLE]

Since $r(x)$ is an analytic function and $r(0)=0$ we can write $r(x)=xr^{\prime}(0)+x^{2}\zeta(x)$, so

[TABLE]

and we see that $\Psi_{n}(0)=0$ and $\Psi_{n}^{\prime}(0)=0$. Than in the limit $n\to\infty$ we have

[TABLE]

The derivatives of $\Psi$ behave as $\Psi_{n}^{(m)}(0)\sim n^{m-1}$, since the derivatives should fall on the bracket in the expression (75) for the value to be non-zero, and then extra power of $x$ is produced, which in turn should be differentiated. So the contribution from $\varepsilon^{m}\Psi_{n}^{(m-1)}(0)$ to $f^{+}$ gives us

[TABLE]

For the second derivative in $\varepsilon$ we thus see that main contribution behaves as $\sum_{n=1}\frac{1}{n^{2}}e^{-n\varepsilon\varphi(b-a)}(n\varepsilon)^{m+1}\xrightarrow[\varepsilon\to 0]{}0$ if $m\geq 3$. For the terms with $m=0,1,2$ we need to take the limit $\varepsilon\to 0$. Doing so, we obtain

[TABLE]

Note that $\frac{\varepsilon}{1-\exp(-n\varepsilon r^{\prime}_{-}(0))}$ is analytic in $\varepsilon$. The derivative can be computed explicitly, but the expression is too long to present it here.

The contribution of $f^{-}$ is computed in exactly the same way:

[TABLE]

At last, we need to consider $f^{\times}=\frac{1}{ab}\sum_{n=1}^{\infty}\frac{1}{n}e^{-n\varepsilon\varphi(b-a)}\varepsilon\sum_{i=0}^{N-1}\Psi_{n}^{-}(i\varepsilon)\sum_{j=0}^{M-1}\Psi_{n}^{+}(j\varepsilon)$. Analyzing the contributions of the derivatives $\Psi_{n}^{(m)}(0)$, similarly to the discussion before the formula (77), we see that $\varepsilon^{m}\Psi_{n}^{(m-1)}(0)\sim\frac{1}{n^{2}}$. Thus

[TABLE]

and these terms give zero contribution after the derivative $\left(\frac{d}{d\varepsilon}\right)^{2}$ for $\varepsilon\to 0$ and $m\geq 3$. The same holds for all the contributions with the derivatives of $\Psi$ in the both sums. In the terms with $m<3$ we can just take a double derivative over $\varepsilon$ and limit $\varepsilon\to 0$ and obtain

[TABLE]

Again, this expression is easily computed explicitly, but is too cumbersome to present it here.

At last we can write the asymptotic expansion coefficients in the present case:

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where in the first line we have the contributions that depend on the choice of the weight function $\varphi$ and in the second line – the universal contribution that appears in the case $\varphi\equiv 1$ and in the box of finite height.

6.1 Numerical checks

Here we present results of comparison of numerical values of the expansion coefficients $f_{0},f_{1},f_{2},f_{3}$ that were obtained from analytical formulas (81),(82),(83),(84) using Mathematica computer algebra system with the results of fitting exact numeric values of free energy for $\frac{1}{\varepsilon}=2\dots 200$ by the formula $f=f_{0}+f_{1}\varepsilon+f_{2}\varepsilon^{2}\ln\varepsilon+f_{3}\varepsilon^{2}+f_{4}\varepsilon^{3}+f_{5}\varepsilon^{4}$.

In the table we present results for various functions $\varphi(z)$ and values of $a,b$.

We see very good coincidence between values, obtained numerically and analytically.

Conclusion and outlook

In the present paper we computed the asymptotic expansion of the free energy in the dimer model on a hexagonal domain of the hexagonal lattice, on a hexagon with a side of infinite length and in the model with coordinate-dependent Boltzmann weights. We presented numerical results supporting our computations. We’ve discussed the physical meaning of the expansion coefficients and argued that our results support the identification of the scaling behavior of the dimer model with the Gaussian free field theory.

In further work we will study the connection of the expansion coefficients with the spectral properties of Dirac operator on the non-frozen domain and study the universality of our results. We will also present the computations of the logarithmic correction to the free energy from the geometry of the ‘‘limit shape’’.

Acknowledgments

We are grateful to professor Nikolai Reshetikhin for his attention to this work. We thank Pavel Belov for useful discussions and general support. We thank Mikhail Vyazovsky for pointing out a mistake in the first version of the paper. We thank Istvan Prause for the reference [32] and Vadim Gorin for a discussion on geometry of the limit shape.

We thank the organizers and participants of the conference MQFT-2018 for the opportunity to present our results and useful discussions.

This research is supported by RFBR grant No. 18-01-00916.