Finite-size scaling of free energy in the dimer model on a hexagonal domain
Anton Nazarov Sergey Paston





Finite-size scaling of free energy in the dimer model on a hexagonal domain
    A.A. Nazarov^1,2, S.A. Paston^1

^1Department of Physics, St. Petersburg State University, 

Ulyanovskaya 1, 198504 St. Petersburg, Russia

^2email:antonnaz@gmail.com

    
===========================================================================================================================================================================




  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.





§ 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
<cit.>. 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 <cit.>, Temperley and Fisher <cit.>
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 <cit.> and by Fan
and Wu <cit.> 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
<cit.>. Later, the well-known limit shape
phenomenon <cit.> was discovered for dimer models. First, the “Arctic circle”
theorem was proven for domino tilings of the domain in the form of “Aztec diamond”
 <cit.>. Then similar result was obtained for a hexagonal domain on the
hexagonal lattice <cit.>. Soon the connection of these results with the theory of
random matrices was established <cit.>. The papers
<cit.> 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 <cit.> and numerical investigation <cit.>. Computation of
correlation functions is a common problem for all vertex models <cit.> as well as
for dimer models. Another problem of great interest is the study of limit shapes in various cases
<cit.>.

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 <cit.>.

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
<cit.>. Scaling limit of dimer model is proven to be
described by a Gaussian free-field theory <cit.>, 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
<cit.>.



§ MODEL DEFINITION


The configurations of the dimer model are perfect matchings (sets of non-touching edges, covering
all the vertices) on some graph G with some choice of weights ω(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” <cit.>.

For a bipartite graph 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 K with
the matrix elements K(w,b) equal to the weight of the edge w→ b: K(w,b)=ω(w→ b).

Then the partition function is equal to the absolute value of the determinant of the Kasteleyn
matrix<cit.>: 

    Z=∑_conf∏_e∈confw(e)=| K|.


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

    (Df)(v)=∑_u K(v,u) f(u).


Kenyon <cit.> and others <cit.>
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. <ref> 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

    h_ij≤ h_i-1,j,   h_ij≤ h_i,j-1.


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:

    E[conf]=∑_i,j h_ij=Volume

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

    Z=∑_conf e^-E[conf]/T=∑_confq^Volume[conf],

where q=exp(-1/T) .

For this particular case, the partition function is given by the classical Macmahon combinatorial
formula <cit.>

    Z[M,N,K,q]=∏_i=1^M∏_j=1^N∏_k=1^K1-q^i+j+k-1/1-q^i+j+k-2.



MacMahon formula is obtained for the following definition of the Kasteleyn matrix. Embed the
hexagonal lattice in the complex plane ℂ 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

    K(w,b)=q^ w+ w  if   w= b
     K(w,b)=1   if   w≠ b.



The   free energy per site is defined as^*
[For convenience we omit the factor 1/T in the usual definition of the free energy.]

    f=-1/Vln Z(M,N,K,q).

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

    V=2(MN+NK+MK)=2(ab+bc+ca) ε^-2.


We are interested in the scaling limit, combined with the thermodynamic limit, when T→∞,
and M,N,K→∞, such that ratios M/T=a,  N/T=b,   K/T=c
remain fixed. In what follows we use ε=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 ε and
derive ε-independent closed expressions for the first several coefficients in this
expansion.
  


§ THE COMPUTATION OF THE FREE ENERGY ASYMPTOTIC EXPANSION


First we substitute MacMahon formula (<ref>) into the free energy definition (<ref>) and
obtain

    f=-1/Vln Z =
    =- ∑_i=1^M∑_j=1^N∑_k=1^K1/V[
      ln(1-e^-ε (i+j+k-3) e^-2ε) -ln(1-e^-ε
          (i+j+k-3) e^-ε)].

We need to notice that for a finite ε 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

    f=-∑_i=1^M∑_j=1^N∑_k=1^K∑_n=1^∞1/V1/n[-e^-nε (i+j+k-3) e^-2nε+e^-nε
          (i+j+k-3) e^-nε].

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

    f=-∑_n=1^∞1/V1/ne^-nε(1-e^-nε)∑_i=1^M∑_j=1^N∑_k=1^Ke^-nε (i+j+k-3) .

The triple sum factorizes into the product of three sums of the type ∑_i=0^M-1e^-nε
  i, that are just the
 sums of geometric progression. So we have ∑_i=0^M-1e^-nε
  i=1-e^-Mnε/1-e^-nε, and obtain

    f=-∑_n=1^∞1/V1/ne^-nε/(1-e^-nε)^2(1-e^-na)(1-e^-nb)(1-e^-nb).

Denote by χ(z) the function

    χ(z)=e^-z(z/1-e^-z)^2,

and by H_n the product of ε-independent exponents

    H_n=(1-e^-na)(1-e^-nb)(1-e^-nc).


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

    χ(z)=1 - z^2/12 + z^4/240+𝒪(z^6).


Using the definition of the volume (<ref>) the free energy is written as

    f=-1/2(ab+bc+ca)∑ _n=1^∞χ(nε)/n^3 H_n.

We are interested in the scaling behavior in ε when ε→ 0. The function f
is regular in ε for ε>0. So we need to
compute  .f|_ε→ 0 and the derivatives .∂ f/∂ε|_ε→ 0, .∂^2 f/∂ε^2|_ε→ 0.

For .f|_ε→ 0 we have χ(0)=1 and the sum ∑_n=1^∞H_n/n^3 is just a sum of  polylogarithms
Li_s(z)=∑_n=1^∞z^n/n^s of third order:

    .f|_ε→ 0 =1/2(ab+bc+ca)[_3(e^-a)+_3(e^-b)+_3(e^-c)-
        _3(e^-a-b).
    .-_3(e^-b-c)-    _3(e^-a-c)+    _3(e^-a-b-c)-ζ(3)].

Here Riemann zeta function appears as a particular value of polylogarithm _3(1)=ζ(3).

The derivative of the function χ(z) is zero, χ'(0)=0, and the series
∑_n=1^∞χ'(nε)/n^2H_n are convergent for a finite
ε, thus we have for a derivative

    .∂ f/∂ε|_ε→ 0=0.


Now let us compute the second derivative

    ∂^2 f/∂ε^2=-1/2(ab+bc+ca)∑ _n=1^∞χ”(nε)/n H_n.

Second derivative of χ is finite χ”(0)=-1/6. But there is a difficulty with the
limit ε→ 0 in this expression due to the poor convergence of the series. Let us
rewrite the second derivative as a sum of two series:

    ∂^2 f/∂ε^2=-1/2(ab+bc+ca)(∑ _n=1^∞χ”(nε)/n+
      ∑_n=1^∞χ”(nε)/n(H_n-1)).

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:

    χ”(0)∑_n=1^∞H_n-1/n=
    =χ”(0)∑_n=1^∞1/n(-e^-a
          n-b n-c n+e^-a n-b n+e^-a n-c n+e^-b n-c n - e^-a n-e^-b n-e^-c n)=
    
      =1/6ln((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)).


First sum in the expression (<ref>) can be expressed as follows. First, consider second
derivative of χ(z), separate an exponent and denote what is left by ξ(z):

    χ”(z)=e^-ze^2 z(z^2+4 e^z
       (z^2-1)+e^2 z(z^2-4
       z+2)+4
       z+2)/(e^z-1)^4=e^-zξ(z).

The sum is then rewritten

    ∑ _n=1^∞χ”(nε)/n=
      ∑_n=1^∞e^-nε/nξ(nε) =
      ∑_n=1^∞e^-nε/nξ(0)+  ε∑_n=1^∞e^-nεξ(nε)-ξ(0)/nε.

The first term is series for the logarithm ln(1-e^-ε)≈lnε. 
Denote by Q(z) the difference of the function ξ value in z and zero, divided by z:

    Q(z)=ξ(z)-ξ(0)/z.

The function Q(z) is analytic in z.

Then the second term can approximated by the integral:

    ε∑_n=1^∞e^-nεQ(nε)=ε∑_n=1^∞e^-nεξ(nε)-ξ(0)/nε≈∫_0^∞ e^-zξ(z)-ξ(0)/z dz=∫_0^∞e^-zQ(z).


    ∂^2 f/∂ε^2≈-1/2(ab+bc+ca)[1/6ln((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))-1/6lnε+∫_0^∞ e^-zξ(z)-ξ(0)/z dz].

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

    f≈ f_0+f_1ε +f_2ε^2lnε+ f_3ε^2+𝒪(ε^3),

as expected in general in two-dimensional critical systems <cit.>. Taking the
derivatives of such expression by ε, we get f_0=f(0),
f_1=.∂ f/∂ε|_ε→ 0,
∂^2 f/∂ε^2=f_2(2lnε+3)+f_3. 

Thus we have the coefficients

    f_0=1/2(ab+bc+ca)[_3(e^-a)+_3(e^-b)+_3(e^-c)-
        _3(e^-a-b).
    .-_3(e^-b-c)-    _3(e^-a-c)+    _3(e^-a-b-c)-ζ(3)],


    f_1=0,


    f_2=-1/2(ab+bc+ca)1/12,

and

    f_3=-1/2(ab+bc+ca)[1/12ln((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))-1/8+ 1/2∫_0^∞ e^-zQ(z) dz
        ].






§ LIMIT SHAPE AND PHYSICAL MEANING OF THE EXPANSION COEFFICIENTS



In the paper <cit.> we have presented results of Monte Carlo simulations using
Metropolis and Wang-Landau algorithms that support the form of expansion (<ref>). In
particular we got f_1=0.04± 0.04 which is consistent with f_1=0. 

The dimers in the scaling limit are in general described by a
free-fermion theory <cit.>, 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 <cit.>.

In the scaling limit ε→ 0 so called “limit shape
phenomenon” <cit.> 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
<cit.>. 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 (<ref>) of the logarithm of the partition function is generic in two-dimensional
models <cit.>.

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 ε is also universal
<cit.> and should appear in all two-dimensional theories with boundary. In the paper
<cit.> 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:

    δ F = -1/6 cχln L,

where c is central charge of the effective field theory and
χ is the Euler characteristic of the manifold χ=2-2 h-b,
where h is the number of handles and b is the number of
boundaries. Naive use of the formula (<ref>) with our result
δ F=1/12 would give c=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
Θ the logarithmic contribution to the free energy is changed to

    δ F =  c/24(Θ/π-π/Θ) ln L.

This formula was recently used by N. Allegra <cit.> to show that the central charge
in the dimer model on the square lattice with a corner monomer is c=1, contrary to some
previous claims of c=-2.

In the case of manifold with conic singularity with the semi-angle
Θ, the logarithmic contribution presented in <cit.> is:

    δ F =  c/24(Θ/π-4π/Θ) ln L.



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

    δ F = (2(ab+bc+ca)) f_2lnε.


Using Cardy-Peschel formulas, presented above, δ F can be
computed from value of the central charge for Gaussian free field
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ξ(z)-ξ(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. 



§ NUMERIC CHECKS


To check our result (<ref>) numerically, we first need to evaluate the integral in f_3.
The function e^-zξ(z)-ξ(0)/z is smooth and but some care is required when taking the
limit z→ 0 numerically. Evaluating the integral numerically with Mathematica we get

    ∫_0^∞e^-zQ(z)=∫_0^∞e^-zξ(z)-ξ(0)/z dz = -0.080842.


The asymptotic expansion (<ref>) gives very good approximation for free energy. In Figure
<ref> we present exact and approximate values for f and, on the right panel,
demonstrate that inaccuracy is actually 𝒪(ε^4). 





§ BOX OF INFINITE HEIGHT



Consider the case of cubes in the box of infinite height. In this case K=∞ 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=∞ in the partition function, we see that all terms with k are cancelled and we have

    Z[M,N,q]=∏_i=1^M∏_j=1^N1/1-q^i+j-1.


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:

    f=-1/MN∑_i=1^M∑_j=1^Nln(1-e^-ε(i+j-2)e^-ε)=
    
      =∑_n=1^∞ε^2/ab1/ne^-nε(∑_i=1^Me^-nε(i-1))
      (∑_j=1^Ne^-nε(j-1)).

Computing the sums of geometric progressions we obtain almost the same expression as (<ref>),
except that H_n now has two multipliers H_n =(1-e^-na)(1-e^-nb):

    f=-1/ab∑_n=1^∞χ(nε)/n^3 H_n.

The function χ(z)=e^-z(z/1-e^-z)^2 is the same as in (<ref>).

Doing the same computations as in section <ref>, we obtain following results
for the coefficients f_0, f_1, f_2, f_3:

    f_0=1/ab[_3(e^-a)+_3(e^-b)-
        _3(e^-a-b)-ζ(3)],


    f_1=0,


    f_2=1/ab1/12,


    f_3=-1/ab[1/12ln((e^a-1)(e^b-1)/(e^a+b-1))-1/8+ 1/2∫_0^∞ e^-zξ(z)-ξ(0)/z dz
        ].


Again in the Fig. <ref> we see good numerical agreement between the asymptotic
expansion (<ref>) and the exact formula (<ref>).





§ 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 φ(t) be a continuous bounded smooth function for t∈[-a,b]. Denote by
q_i=e^-εφ(iε) for i=(-M+1)… N-1. We index diagonal slices in
such a way, that top left corner of the M× N box has index N-M, top right corner – N,
bottom left – -M and bottom right – 0. (See Fig. <ref> ). 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 <cit.>. In this case it contains products of weights
q_i for slices with negative and positive indexes:

    Z=∏_i=0^N-1∏_j=0^M-11/1-q_N-M^-1(∏_k=0^iq_N-M-k)(∏_l=0^jq_N-M+l).

(see Theorem 2 in <cit.>).
If we set q_i=q, we recover formula (<ref>).

Now let us again consider free energy dependence on ε:

    f=-ε^2/ab∑_i=0^N-1∑_j=0^M-1ln(1-e^εφ(b-a)e^-ε∑_k=0^iφ(b-a-kε)
        e^-ε∑_l=0^jφ(b-a+lε)).


The terms f_0  in the asymptotic expansion (<ref>) can be easily obtained from
this expression. To do so we need to approximate the sums by the integrals, using
∑_k=0^iεφ(b-a-kε)≈∫_0^iεφ(b-a-x)
and similar approximation for the outer sums:

    f_0=-1/ab∫_0^b∫_0^aln(1-e^-∫_0^yφ(b-a-x)
        e^-∫_0^zφ(b-a+x̅))=
    
      =-1/ab∫_0^b∫_0^aln(1-e^-∫_-y^zφ(b-a+x)).

Substituting in this expression φ(x)≡ 1 we immediately get formula (<ref>), 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:

  Let a=Mε and b=Nε.
  For an integrable almost everywhere analytic function f and analytic function φ the sum is approximated by
  an integral with the corrections  of order o(ε):
  
    ∑_i=0^M-1∑_j=0^N-1ε^2f(ε∑_k=-i^jφ(kε))=
        ∫_0^a∫_0^bf(∫_-y^zφ(x)) +o(ε).


  Use Euler-Maclaurin formula to approximate the sums by the integrals with corrections. It is
  enough to consider linear in ε terms. First use
  it for the sum in the argument of function f:
  
    ∑_i=0^M-1∑_j=0^N-1ε^2f(ε∑_k=-i^jφ(kε))=
        ∑_i=0^M-1∑_j=0^N-1ε^2f(∫_-iε^jεφ(x)+ε/2(φ(jε)-φ(-iε))+𝒪(ε^2)).

  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]×[0,b], then

    ∫_0^a∫_0^bG(y,z)   ≈∑_i=0^M-1∑_j=0^N-1{ε^2G(iε,jε)+
        ε^3/2((∂_y+∂_z)G)(iε,jε)}+o(ε^2).

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

Denote by G_ε(y,z) the function of two arguments that appear in the double sum in (<ref>):

    G_ε(y,z)=f(∫_-y^zφ(x)+ε/2(φ(z)-φ(-y))+o(ε)).

Use formula (<ref>) to express the double sum in (<ref>) as the integral plus correction:

    ∑_i=0^M-1∑_j=0^N-1ε^2G_ε(iε,jε)=∫_0^a∫_0^b
      G_ε(y,z) -
      ε^3/2∑_i=0^M-1∑_j=0^N-1((∂_y+∂_z)G_ε)(iε,jε)+o(ε^3).

Now expand the function G_ε(y,z) under integral into Taylor series in ε:

    ∫_0^a∫_0^bG_ε(y,z) =
      ∫_0^a∫_0^bG_0(y,z) +∫_0^a∫_0^bε/2(φ(z)-φ(-y))·
      f'(∫_-y^zφ(x)) +o(ε).

On the other hand, the derivatives of G_ε in y,z are computed as

    ∂_zG_ε(y,z)=f'(∫_-y^zφ(x)+ε/2(φ(z)-φ(-y))+o(ε))φ(z)+o(ε)
    ∂_yG_ε(y,z)=-f'(∫_-y^zφ(x)+ε/2(φ(z)-φ(-y))+o(ε))φ(-y)+o(ε).

Substituting equations (<ref>),(<ref>) into (<ref>) and using formula
(<ref>) to change the double sum of the derivatives of G_ε into double integral,
we see that the linear correction in ε is cancelled. 





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

    f_1=0.


Next terms in the asymptotic expansion require more care, since we need to single out the
contribution of the order ε^2lnε.

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

    f=∑_n=1^∞ε^2/ab1/ne^nεφ(b-a)∑_i=0^N-1∑_j=0^M-1 e^-nε∑_k=0^iφ(b-a-kε)
        e^-nε∑_l=0^jφ(b-a+lε).

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

    f=∑_n=1^∞ε^2/ab1/ne^-nεφ(b-a)(∑_i=0^N-1
         e^-nε∑_k=1^iφ(b-a-kε))
       (∑_j=0^M-1e^-nε∑_l=1^jφ(b-a+lε)),

where we assume that ∑_k=1^0φ(b-a-kε)=0 and
∑_l=1^0φ(b-a+lε)=0. Note that for φ(t)≡ 1 we recover formula
(<ref>).



Introduce following notations:

    r_∓(iε)≡ε∑_k=1^iφ(b-a∓ kε).

The functions r_∓ can be approximated using Euler-Maclaurin formula as:

    r_∓(x)≈∫_0^xφ(b-a∓ y)  +ε/2(φ(b-a∓
        x)-φ(b-a))∓ε^2/12(φ'(b-a∓ x)-φ'(b-a)).

Denote by Ψ_n^∓ the following functions:

    Ψ^∓_n(x)≡ e^-nr_∓(x)-e^-nxr_∓'(0).

These functions are analytic in ε, moreover, the higher derivatives of Ψ^∓_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

    f=1/ab∑_n=1^∞1/ne^-nεφ(b-a)[ε^2∑_i=0^N-1e^-nε
        i r_-'(0)∑_j=0^M-1e^-nε j
        r_+'(0)+ε∑_i=0^N-1e^-nε i
        r_-'(0)ε∑_j=0^M-1Ψ^+_n(jε)+.
    . +ε∑_i=0^N-1Ψ^-_n(iε) ·ε∑_j=0^M-1e^-nε
      j r_+'(0)+ε∑_i=0^N-1Ψ^-_n(iε)
    ·ε∑_j=0^M-1Ψ^+_n(jε)]=
    
    =f^(0)+f^++f^-+f^×.

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

    f^(0)=1/ab∑_n=1^∞1/n e^-nεφ(b-a)ε^2(1-e^-n b r_-'(0))(1-e^-n a
            r_+'(0))/(1-e^-nε r_-'(0))(1-e^-nε r_+'(0)).


Denote by H_n(ε) the numerator of this expression:

    H_n(ε)=(1-e^-n b r_-'(0))(1-e^-n a r_+'(0)).

Note that

    r_∓'(0)=φ(b-a)∓ε/2φ'(b-a)+ε^2/12φ”(b-a)+𝒪(ε^3).

can be rewritten as

    r_∓'(0)=φ(b-a)(1+α_∓(ε)),

where

    α_∓(ε)=∓ε/2φ'(b-a)/φ(b-a)+ε^2/12φ”(b-a)/φ(b-a)+𝒪(ε^3).


Now for f^(0) we have

    f^(0)=1/ab∑_n=1^∞χ̃_n(nεφ(b-a),ε) /n^3φ(b-a)^2,

where χ̃_n(z,ε)  is defined as

    χ̃_n(z,ε)≡z^2 e^-z
        H_n(ε)/(1-e^-z(1+α_-(ε)))
        (1-e^-z(1+α_+(ε))).

It is important to include H_n(ε), since H_n(ε) contains decaying
exponents in n.

Taking the derivatives of f^(0) in ε we get

    (d/dε)^2 f^(0)=1/abφ(b-a)∑_n=1^∞1/n^3[
        .χ̃”_n,εε(nεφ(b-a),ε)|_ε=0+
        2nφ(b-a).χ̃”_n,zε(nεφ(b-a),ε)|_ε=0+.
    .
          n^2(φ(b-a))^2χ̃”_n,zz(nεφ(b-a),ε)  ]+
        𝒪(εlnε)=
    
        =1/abφ(b-a)[∑_n=1^∞1/n^3χ̃”_n,εε(0,0)+2φ(b-a)∑_n=1^∞1/n^2χ̃”_n,zε(0,0)]+
        
    
        +1/ab∑_n=1^∞H_n/nχ”(nεφ(b-a))+𝒪(εlnε),

where the function χ(z) is the same as in the equation (<ref>)

    χ(z)=z^2 e^-z/(1-e^-z)^2

and

    H_n=(1-e^-nbφ(b-a)) (1-e^-naφ(b-a)).


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

    1/ab∑_n=1^∞H_n/nχ”(nεφ(b-a))=
      1/ab(1/6lnε +1/6lnφ(b-a)+ ∫_0^∞
        e^-zQ(z) +χ”(0)∑_n=1^∞H_n-1/n)=
    
      =1/ab(1/6lnε +1/6lnφ(b-a)+ ∫_0^∞
        e^-zQ(z) +1/6ln((1-e^aφ(b-a))(1-e^aφ(b-a))/1-e^(a+b)φ(b-a))),

where Q(z) was defined in the equation (<ref>).

Thus the contribution of f^(0) to the asymptotic expansion (<ref>) consists of the term 

    f^(0)_3=1/abφ(b-a)[∑_n=1^∞1/n^3χ̃”_n,εε(0,0)+2φ(b-a)∑_n=1^∞1/n^2χ̃”_n,zε(0,0)],

that is cumbersome but does not depend on ε and can be computed explicitly and the
expression (<ref>), which gives us logarithmic behavior in ε.

Now we need to consider other terms in the expression (<ref>). For f^+ we have

    f^+=1/ab∑_n=1^∞1/n e^-nεφ(b-a)ε1-e^-nbr_-'(0)/1-e^-nε r_-'(0)ε∑_j=0^M-1Ψ^+_n(jε).

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

    ε∑_j=0^M-1Ψ^+_n(jε) = ∫_0^aΨ^+_n(x)  -
      ε/2(Ψ^+_n(a)-Ψ^+_n(0))+
      ε^2/12(Ψ^+'_n(a)-Ψ^+_n(0))+
      𝒪(ε^mΨ^+(m-1)_n).

Since r(x) is an analytic function and r(0)=0 we can write r(x)=xr'(0)+x^2ζ(x), so 

    Ψ_n(x)=e^-nr'(0)x(e^-n(r(x)-xr'(0))-1)=e^-nxr'(0)(e^-nx^2ζ(x)-1)

and we see that Ψ_n(0)=0 and Ψ_n'(0)=0.
Than in the limit n→∞ we have

    ∫_0^cΨ_n(x)
      ∫_0^ce^-nxr'(0)(-nx^2ζ(x)+…)∼1/n^2.

The derivatives of Ψ behave as Ψ_n^(m)(0)∼ n^m-1, since the derivatives should
fall on the bracket in the expression (<ref>) 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
ε^mΨ_n^(m-1)(0) to f^+ gives us

    ∑_n=1^∞1/n e^-nεφ(b-a)1-e^-nbr_-'(0)/1-e^-nε r_-'(0)ε^m+1Ψ_n^(m-1)(0)∼∑_n=1^∞1/n^4 e^-nεφ(b-a)1-e^-nbr_-'(0)/1-e^-nε r_-'(0) (ε n )^m+1.

For the second derivative in ε we thus see that main contribution behaves as
∑_n=11/n^2e^-nεφ(b-a)
  (nε)^m+10
if m≥ 3. For the terms with m=0,1,2 we need to
take the limit ε→ 0. Doing so, we obtain

    f^+_3≡.(d/dε)^2 f^+|_ε=0 =
    .(d/dε)^21/ab∑_n=1^∞1/n
      e^-nεφ(b-a)1-e^-nbr_-'(0)/1-e^-nε r_-'(0)ε(∫_0^aΨ_n^+(x)  -εΨ_n^+(a) +ε^2Ψ^+'_n(0))|_ε=0.

Note that ε/1-exp(-nε r'_-(0)) is analytic in ε. 
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:

    f^-_3≡.(d/dε)^2 f^-|_ε=0 =
    .(d/dε)^21/ab∑_n=1^∞1/n
      e^-nεφ(b-a)1-e^-nar_+'(0)/1-e^-nε r_+'(0)ε(∫_0^bΨ_n^-(x)  -εΨ_n^-(a) +ε^2Ψ^-'_n(0))|_ε=0.

At last, we need to consider
f^×=1/ab∑_n=1^∞1/ne^-nεφ(b-a)ε∑_i=0^N-1Ψ_n^-(iε)
∑_j=0^M-1Ψ_n^+(jε). Analyzing the contributions of the derivatives
Ψ_n^(m)(0), similarly to the discussion before the formula (<ref>), we see that
ε^mΨ_n^(m-1)(0)∼1/n^2. Thus

    ∑_n=1^∞1/ne^-nεφ(b-a)1/n^2ε^mΨ_n^(m-1)(0)∼∑_n=1^∞1/n^5e^-nεφ(b-a)(ε n)^m,
 and these terms
give zero contribution after the derivative (d/dε)^2 for
ε→ 0 and m≥ 3. The same holds for all the contributions with the derivatives of
Ψ in the both sums. In the terms with m<3 we can just take a double derivative over
ε and limit ε→ 0 and obtain

    f^×_3≡.(d/dε)^2
        f^×|_ε=0=.1/ab∑_n=1^∞1/n(d/dε^2)^2(∫_0^bΨ^-_n(x)
        -ε/2Ψ^-_n(b)+ε^2/12Ψ^-'_n(b))·.
    ·.(∫_0^aΨ^+_n(x)-ε/2Ψ^+_n(a)+
        ε^2/12Ψ^+_n(a))|_ε=0.

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:

    f_0= -1/ab∫_0^b∫_0^aln(1-e^-∫_-y^zφ(b-a+x)),


    f_1=0,


    f_2=1/ab1/12

and

    f_3=f^(0)_3+f^+_3+f^-_3+f^×_3+
      1/ab1/12ln(φ(b-a)(1-e^aφ(b-a))(1-e^aφ(b-a))/1-e^(a+b)φ(b-a))+
    
      +1/ab(1/8 +1/2∫_0^∞
        e^-zq(z) ),

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



 §.§ 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
(<ref>),(<ref>),(<ref>),(<ref>) using Mathematica computer algebra system
with the results of fitting exact numeric values of free energy for 1/ε=2…
200
by the formula f=f_0+
f_1ε+
f_2ε^2lnε+f_3ε^2+f_4ε^3+f_5ε^4. 

In the table we present results for various functions φ(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
<cit.> 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.

utphys
 
