Power series expansions for the planar monomer-dimer problem

TL;DR

This paper develops a computational method to generate extensive power series expansions for the free energy of the planar monomer-dimer model, surpassing previous results despite high complexity.

Contribution

It introduces a new approach to compute longer power series expansions for the monomer-dimer problem, enabling more accurate bounds and numerical estimates.

Findings

Computed nearly three times more terms than previous work.

Provided tighter bounds and more precise numerical values for free energy.

Demonstrated potential applicability to other complex models.

Abstract

We compute the free energy of the planar monomer-dimer model. Unlike the classical planar dimer model, an exact solution is not known in this case. Even the computation of the low-density power series expansion requires heavy and nontrivial computations. Despite of the exponential computational complexity, we compute almost three times more terms than were previously known. Such an expansion provides both lower and upper bound for the free energy, and allows to obtain more accurate numerical values than previously possible. We expect that our methods can be applied to other similar problems.

Tables

Table 1. Table 1 . Values of a k subscript 𝑎 𝑘 a_{k}

$k$	$a_k$	$k$	$a_k$	$k$	$a_k$
$2$	$\frac{1}{16}$	$23$	$\frac{4312434281365}{17803292276948992}$	$44$	$\frac{18487601206244410582171859}{292772819290992435013642878976}$
$3$	$\frac{1}{192}$	$24$	$\frac{5789230773063}{25895697857380352}$	$45$	$\frac{74150661042096992710148129}{1225560638892526472150132981760}$
$4$	$\frac{7}{1536}$	$25$	$\frac{69044819053441}{337769972052787200}$	$46$	$\frac{297604910587450946018199331}{5125071762641474338082374287360}$
$5$	$\frac{41}{10240}$	$26$	$\frac{272097812497681}{1463669878895411200}$	$47$	$\frac{1194303993371769853836734501}{21411410919479937234655252578304}$
$6$	$\frac{181}{61440}$	$27$	$\frac{1068966474984721}{6323053876828176384}$	$48$	$\frac{4789513328571295127284133845}{89369367316090172805517575979008}$
$7$	$\frac{757}{344064}$	$28$	$\frac{601281977474899}{3891110078048108544}$	$49$	$\frac{19188774086998950351884051009}{372689276467099444040030742380544}$
$8$	$\frac{3291}{1835008}$	$29$	$\frac{16672616519735441}{117021532717594968064}$	$50$	$\frac{76803645872757902332578961121}{1552871985279581016833461426585600}$
$9$	$\frac{14689}{9437184}$	$30$	$\frac{66545602395606901}{501520854503978434560}$	$51$	$\frac{307176141884436645170078617001}{6465018061163969947633186347417600}$
$10$	$\frac{64771}{47185920}$	$31$	$\frac{267471214350929957}{2144433998568735375360}$	$52$	$\frac{1228026136368811312663436458705}{26894475134442114982154055205257216}$
$11$	$\frac{276101}{230686720}$	$32$	$\frac{1080431496491179115}{9149585060559937601536}$	$53$	$\frac{4909003176336757275553467075425}{111796641735328007376797249088520192}$
$12$	$\frac{1132693}{1107296256}$	$33$	$\frac{4374403039126240385}{38959523483674573012992}$	$54$	$\frac{19627584575160129028816787257753}{464386050285208646026696265444622336}$
$13$	$\frac{4490513}{5234491392}$	$34$	$\frac{17705045340400677607}{165577974805616935305216}$	$55$	$\frac{78505240133588264624896189049521}{1927640208731054757091946762222960640}$
$14$	$\frac{17337685}{24427626496}$	$35$	$\frac{71484177460946258777}{702452014326859725537280}$	$56$	$\frac{314123632091141305526902518303973}{7996137162143634547936964346998947840}$
$15$	$\frac{65867621}{112742891520}$	$36$	$\frac{287529593953850293471}{2975090884207876484628480}$	$57$	$\frac{1257288843192384664389299749835521}{33147623144886339580538688565741092864}$
$16$	$\frac{249437227}{515396075520}$	$37$	$\frac{1151710503160001680385}{12580384310364734849286144}$	$58$	$\frac{1677695623930304081656255827713551}{45775289104843040373124855638404366336}$
$17$	$\frac{955110593}{2336462209024}$	$38$	$\frac{4596336312298962012663}{53117178199317769363652608}$	$59$	$\frac{20147683002193594117896886735926057}{568577275196997764634603470034917392384}$
$18$	$\frac{3740591431}{10514079940608}$	$39$	$\frac{18298456303802689186745}{223953508083610054614319104}$	$60$	$\frac{26879884904186172110556704720248631}{784244517513100365013246165565403299840}$
$19$	$\frac{15039656569}{47004122087424}$	$40$	$\frac{72784234597284215364691}{942962139299410756270817280}$	$61$	$\frac{322682332818808295011085893297500673}{9729948929145584189655867681252122296320}$
$20$	$\frac{61727254227}{208907209277440}$	$41$	$\frac{289698911730110389042529}{3965276688335983693036257280}$	$62$	$\frac{1290942327848947576849492154270349133}{40217122240468414650577586415842105491456}$
$21$	$\frac{255640084561}{923589767331840}$	$42$	$\frac{1155125274097244765650075}{16654162091011131510752280576}$	$63$	$\frac{5163832046366445947035366917883833877}{166142865649148204785992652078560829243392}$
$22$	$\frac{50273131919}{193514046488576}$	$43$	$\frac{4616317010648384103125561}{69866240967168649264619323392}$	$64$	$\frac{983546099095446058993477411998292607}{32667107224410092492483962313449748299776}$

Table 2. Table 2 . Comparison with  [ 28 ] . Digits in square brackets mean the corresponding digit in lower and upper bounds.

$p$	[28]	Our estimate
$10/20$	$0.633195588930\left[4-5\right]$	$0.6331955889305251415416\left[5-6\right]$
$11/20$	$0.650499726669\left[5-8\right]$	$0.6504997266695759205\left[7-8\right]$
$12/20$	$0.66044120984\left[2-5\right]$	$0.66044120984322136\left[2-4\right]$
$13/20$	$0.6625636470\left[2-4\right]$	$0.66256364703101\left[3-4\right]$
$14/20$	$0.65620036\left[0-1\right]$	$0.656200361027\left[4-5\right]$
$15/20$	$0.64039026\left[3-5\right]$	$0.6403902642\left[8-9\right]$
$16/20$	$0.6137181\left[3-4\right]$	$0.613718137\left[2-7\right]$
$17/20$	$0.573983\left[2-3\right]$	$0.573983\left[2-3\right]$
$18/20$	$0.51739\left[1-2\right]$	$0.51739\left[1-3\right]$
$19/20$	$0.435\left[8-9\right]$	$0.435\left[8-9\right]$

Equations

$$f_2(p)=\underset{n,m\to \infty }{\lim }\frac{\ln a_p(mn)}{mn}.$$

$$f_2(p)=\frac{1}{2}\left((2\ln 2-1)p-p\ln p\right)-(1-p)\ln (1-p)+\sum _{j=2}^{\infty }{a}_j{p}^j.$$

$$f_2(p)=\frac{1}{2}\left((2\ln 2+1)p-p\ln p\right)+\sum _{j=2}^{\infty }{b}_j{p}^j.$$

$${\Theta }_{m,n}(z)=\sum _{s=0}^{\lfloor mn/2\rfloor }{a}_{2s/mn}(m,n)z^s.$$

$$\Theta (z)=\underset{m,n\to \infty }{\lim }{\left({\Theta }_{m,n}(z)\right)}^{\frac{1}{mn}},\text{ and }\ln \Theta (z)=\underset{m,n\to \infty }{\lim }\frac{\ln {\Theta }_{m,n}(z)}{mn}.$$

$$f_2(p)=\underset{z\in {\mathbb{R}}_{+}}{\inf }\left\{-\frac{p}{2}\ln z+\ln \Theta (z)\right\}.$$

$$\frac{d}{\mathrm{dz}}\left(-\frac{p}{2}\ln z+\ln \Theta (z)\right)=-\frac{p}{2z}+{\left(\ln \Theta (z)\right)}^{\prime }=0,$$

$$p(z)=2z{\left(\ln \Theta (z)\right)}^{\prime }.$$

$$\ln \Theta (z)=\underset{p\in \mathbb{R}}{\sup }\left\{\frac{p}{2}\ln z+f_2(p)\right\},$$

$$\frac{d}{\mathrm{dp}}\left(\frac{p}{2}\ln z+f_2(p)\right)=\frac{\ln z}{2}+f_2^{\prime }(p)=0.$$

$$f_2(p)=\frac{1}{2}\underset{0}{\overset{p}{\int }}z(p)\mathrm{dp},\text{ where }p(z)=2z{\left(\ln \Theta (z)\right)}^{\prime }.$$

$$f_2(p)=\frac{1}{2}\underset{0}{\overset{p}{\int }}\ln \left(\frac{p}{4}+O(p^2)\right)\mathrm{dp}=\frac{1}{2}\underset{0}{\overset{p}{\int }}\left(\ln p-2\ln 2+O(p)\right)\mathrm{dp}=\frac{1}{2}\ln p-\frac{1}{2}p-p\ln 2+\sum _{k=2}^{\infty }{a}_k{p}^k.$$

$$\ln \Theta (z)-\left(S_N-3S_{N-1}+3S_{N-2}-S_{N-3}\right)=O(z^{N-1}),$$

$$W(d_1,d_2)=\{\begin{array}{l}1,\text{if }{d}_1\text{ and }{d}_2\text{ do not overlap,}\\ 0,\text{if }{d}_1\text{ and }{d}_2\text{ overlap.}\end{array}$$

$${\Theta }_{m,n}(z)=\sum _{s=0}^{\infty }\frac{z^s}{s!}\sum _{(d_1,\dots ,d_s)\in D_{m,n}^s}\left(\prod _{1⩽i<j⩽s}W(d_i,d_j)\right),$$

$${\Theta }_{m,n}(z)=\sum _{s=0}^{\infty }\frac{z^s}{s!}\sum _{\mathbf{d}=(d_1,\dots ,d_s)\in D_{m,n}^s}\sum _{G\in {\mathcal{G}}_s}F(\mathbf{d},G),\text{ where }F(\mathbf{d},G)=\prod _{(ij)\in E(G)}F(d_i,d_j),$$

$${\mathcal{W}}_{m,n}(G)=\sum _{\mathbf{d}\in D_{m,n}}F(\mathbf{d},G).$$

$$\ln {\Theta }_{m,n}(z)=\sum _{s=0}^{\infty }\frac{z^s}{s!}\sum _{G\in {\mathcal{C}}_s}{\mathcal{W}}_{m,n}(G),$$

$${\left(m-h(\mathbf{d})+1\right)}_{+}\cdot {\left(n-w(\mathbf{d})+1\right)}_{+},$$

$${\mathcal{W}}_{m,n}(G)$$

$$\stackrel{(ii)}{=}\sum _{\mathbf{d}\in {\mathcal{T}}_s}\left(\sum _{{\mathbf{d}}^{\prime }\in D_{m,n}^s,{\mathbf{d}}^{\prime }\sim \mathbf{d}}F({\mathbf{d}}^{\prime },G)\right)$$

$$\stackrel{(iii)}{=}\sum _{\mathbf{d}\in {\mathcal{T}}_s}{\left(m-h(T)+1\right)}_{+}\cdot {\left(n-w(T)+1\right)}_{+}\cdot F(\mathbf{d},G).$$

$$\ln {\Theta }_{m,n}(z)=\sum _{s=0}^{\infty }\frac{z^s}{s!}\sum _{G\in {\mathcal{C}}_s}{\mathcal{W}}_{m,n}(G)$$

$$=\sum _{s=0}^{\infty }\frac{z^s}{s!}\sum _{G\in {\mathcal{C}}_s}\left(\sum _{\mathbf{d}\in {\mathcal{T}}_s}{\left(m-h(\mathbf{d})+1\right)}_{+}\cdot {\left(n-w(\mathbf{d})+1\right)}_{+}\cdot F(\mathbf{d},G)\right)$$

$$=\sum _{s=0}^{\infty }\frac{z^s}{s!}\sum _{\mathbf{d}\in {\mathcal{T}}_s}{\left(m-h(\mathbf{d})+1\right)}_{+}\cdot {\left(n-w(\mathbf{d})+1\right)}_{+}\cdot \mathcal{W}(\mathbf{d}).$$

$$\ln {\Theta }_{m,n}(z)=\sum _{s=0}^{\infty }\frac{z^s}{s!}\sum _{\mathbf{d}\in {\mathcal{T}}_s}{\left(m-h(\mathbf{d})+1\right)}_{+}\cdot {\left(n-w(\mathbf{d})+1\right)}_{+}\cdot \mathcal{W}(\mathbf{d}).$$

$$\ln \Theta (z)=\underset{m,n\to \infty }{\lim }\frac{\ln {\Theta }_{m,n}(z)}{mn}=\sum _{s=0}^{\infty }\frac{z^s}{s!}\sum _{\mathbf{d}\in {\mathcal{T}}_s}\underset{m,n\to \infty }{\lim }\left(\frac{{\left(m-h(\mathbf{d})+1\right)}_{+}{\left(n-w(\mathbf{d})+1\right)}_{+}}{mn}\right)\mathcal{W}(\mathbf{d}).$$

$$\ln \Theta (z)=\sum _{s=0}^{\infty }\frac{z^s}{s!}\sum _{\mathbf{d}\in {\mathcal{T}}_s}\mathcal{W}(\mathbf{d}).$$

$$S_N=\sum _{s=0}^{\infty }\frac{z^s}{s!}\sum _{\mathbf{d}\in {\mathcal{T}}_s}\left(\genfrac{}{}{0ex}{}{(N-w(\mathbf{d})-h(\mathbf{d})+3{)}_{+}}{3}\right)\mathcal{W}(\mathbf{d}).$$

$$\sum _{m+n=N}(m-h(\mathbf{d})+1{)}_{+}\cdot (n-w(\mathbf{d})+1{)}_{+}.$$

$$\sum _{k=1}^{N-p+1}k\left(N-p+2-k\right)=\left(N-p+2\right)\left(\sum _{k=1}^{N-p+1}k\right)-\left(\sum _{k=1}^{N-p+1}{k}^2\right).$$

$$1-\left(\left(\genfrac{}{}{0ex}{}{N-p+3}{3}\right)-3\left(\genfrac{}{}{0ex}{}{N-p+2}{3}\right)+3\left(\genfrac{}{}{0ex}{}{N-p+1}{3}\right)-\left(\genfrac{}{}{0ex}{}{N-p}{3}\right)\right).$$

Full text

This work was supported by the Austrian Science Fund FWF grant Y464-N18.

1. Introduction

The exact solution of the closed-packed dimer plane model obtained in [1, 2, 3] is a fundamental result in statistical mechanics and combinatorics. In particular, it implies that the number of tilings of $m\times n$ rectangle using dimers grows as $e^{\frac{G}{\pi}mn}$, where $G\approx 0.916$ is the Catalan constant. Similar results were later obtained for other shapes (see [4] and references therein). Applications to physics suggest two natural further questions: what if the dimension of the lattice is higher (i.e. we compute the number of tilings of a hyperrectangle), and what if we consider tilings using both dimers and monomers with a fixed proportion. For the first question, we refer the reader to [5, 6, 7] and references therein. The second questions originates from the study of liquid mixtures on crystal surfaces in [8], see also [9] for comparison with experimental data. Monomer-dimer systems also arise in connection with the Ising model and the Heisenberg model, see [10, Sect. 5]. For both these questions, the exact solution is out of reach so far. However, even finding the answer numerically leads to very challenging computational problems, because the underlying combinatorial counting problems are very hard, and even a small change of the parameters of the problem makes computations much harder or even unfeasible. For example, in [11] it is proved that the monomer-dimer tilings counting problem is #P-complete in the sense of theoretical computer science.

In this paper we will focus on the second question, namely on the case of planar monomer-dimer tilings with fixed dimer density. Let us state the problem precisely. We denote by $a_{p}(m,n)$ the number of tilings of an $m\times n$ rectangle by monomers and dimers with exactly $\lfloor pmn/2\rfloor$ dimers, where $m,n$ are positive integers, and $p\in[0,1]$. Then $p$ is roughly the fraction of the area covered by dimers. We are interested in the limit

[TABLE]

In other words, we want to determine the constant $\lambda$ such that $a_{p}(m,n)\sim e^{\lambda mn}$ as a function of $p$. From the point of view of statistical mechanics, $f_{2}(p)$ is equal to the negative of the Helmholtz free energy per lattice site expressed in units of the thermal energy $k_{B}T$. Some lower and upper bounds for $f_{2}(p)$ were rigorously proved in [12, 13]. However, these bound are not very tight.

Another approach, taken in a series of papers [14, 15, 16] and independently in [17, IV.A] is to expand this function as a power series in $p$ (in some of these papers also expansion with respect to the dimension was discussed). In the former papers, the authors look for a representation of $f_{2}(p)$ of the form

[TABLE]

In [17], the representation is of the form (see details in Section 2)

[TABLE]

Expanding $(1-p)\ln(1-p)$ into a Taylor series in $p$, it is easy to move back and forth between (1) and (2) (see (18)). An important observation is that $a_{j}>0$ and $b_{j}<0$ for all known $a_{j}$ and $b_{j}$. Under the assumption that this pattern holds for all $j$, truncations of (1) and (2) provide lower and upper bound for $f_{2}(p)$, respectively. Thus, computing more terms of these series would result in tighter bounds for $f_{2}(p)$.

Previously, the record in the number of computed terms in (1) and (2) was $23$ (i.e. from $a_{2}$ till $a_{24}$) obtained in [16, Table II]. This result is highly nontrivial, because the underlying combinatorial problem has exponential complexity, so the cost of every next term is usually higher by some factor. In this paper, we compute $63$ terms (from $a_{2}$ to $a_{64}$, see Table LABEL:table:ak) i.e., is almost three times more than what was previously possible.

The contribution of our paper is two-fold. First, our approach allows us to compute significantly more terms for both series (1) and (2) than were previously known and, combining them, we obtain very accurate values of $f_{2}(p)$. Moreover, we provide additional support for the important conjecture that $a_{j}>0$ and $b_{j}<0$ for all $j\geqslant 2$. Second, we show how methods of computer algebra (guessing, modular computation, etc.) can be applied to study models in statistical mechanics. We expect that our methods can be used in other problems of this type (monomer-polymer mixtures, other types of lattices, etc.) in order to push computational limits further.

The rest of the paper is organized as follows. In Section 2 we collect some known results and approaches that connect power series expansion of $f_{2}(p)$ to the combinatorial data. Section 3 contains Theorem 1, a main combinatorial ingredient of our computation. Section 4 contains the description of our algorithm together with all computer algebra machinery used to speed it up. Finally, in Section 5 we describe our implementation, provide numerical results, and compare them to the previous work.

The author is grateful to Manuel Kauers and Doron Zeilberger for introducing him to the subject and for helpful discussions, and to the referees for the comments on the manuscript.

2. Reduction to grand-canonical partition function

For every fixed $m$ and $n$, consider the grand-canonical partition function

[TABLE]

We consider a thermodynamic limit of $\Theta_{m,n}(z)$ (its existence is proved in [10, VIII])

[TABLE]

Since $\Theta_{m,n}^{\prime}(0)=2mn-m-n$, $\ln\Theta(z)=2z+O(z^{2})$. Theorems [10, 8.8A, 8.8B] rewritten in our notation (i.e., replacing $\mu$ with $\ln z$, $g(\mu)$ with $\ln\Theta(z)$, $\rho$ with $\frac{p}{2}$, and $h(\rho)$ with $f_{2}(p)$) state that

[TABLE]

We compute $z$, where the expression $-\frac{p}{2}\ln z+\ln\Theta(z)$ is minimal

[TABLE]

hence

[TABLE]

Since $\Theta(z)=1+2z+O(z^{2})$, then $2z\left(\ln\Theta(z)\right)^{\prime}=4z+O(z^{2})$. Hence, there exists a unique compositional inverse $z(p)$ of $p(z)$ and $z(p)$ is a formal power series (see [18, Th. 1.8]). Moreover, $z(p)=\frac{p}{4}+O(p^{2})$.

Due to [10, Eq. 8.24], $\ln\Theta(z)$ is a convex function in $\ln z$ for all $z\in\mathbb{R}_{+}$. Then equation (5) can be seen as a fact that $-f_{2}(p)$ as a function of $\frac{p}{2}$ is a Legendre transform of $\ln\Theta(z)$ as a function of $\ln z$ (for introduction to Legendre transform, see [19] and [20, §14]). Involutivity of Legendre transform (see [20, §14.C]) implies that

[TABLE]

and the supremum on the right-hand side is reached at $p=p(z)$. On the other hand, we can find this $p$ also by differentiation

[TABLE]

Hence, $f_{2}^{\prime}(p(z))=-\frac{\ln z}{2}$. Substituting $z(p)$ into $z$, we obtain $f_{2}^{\prime}(p)=-\frac{1}{2}\ln z(p)$. Integrating with respect to $p$, and using the initial condition $f_{2}(0)=0$, we conclude that

[TABLE]

The same formula is obtained in [17, IV.A] using another argument.

Finally, we deduce expansion (2) from (8). Using $z(p)=\frac{p}{4}+O(p^{2})$, we conclude that

[TABLE]

The latter expression is exactly of the same form as the right-hand side in (2).

3. Computation of $\Theta(z)$ using $\Theta_{m,n}(z)$

The goal of the present section is to prove the following theorem, which provides a way to compute the thermodynamical limit $\ln\Theta(z)$.

Theorem 1**.**

For every integer $N\geqslant 4$,

[TABLE]

where $S_{M}=\sum\limits_{m+n=M}\ln\Theta_{m,n}(z)$.

In what follows, we will use some properties of the Mayer expansion following [21, Section 2.2].

Let $R_{\infty,\infty}$ be the first quadrant of the plane. We denote the $m\times n$ rectangle whose lower-left corner is the origin by $R_{m,n}$. By definition, $R_{m,n}\subset R_{\infty,\infty}$ for all $m$ and $n$. We denote the set of all dimers in $R_{m,n}$ by $D_{m,n}$. By definition, the cardinality of $D_{m,n}$ is $2mn-m-n$, and $D_{m,n}\subset D_{\infty,\infty}$ for every $m,n\in\mathbb{Z}_{>0}$. For $d_{1},d_{2}\in D_{\infty,\infty}$, we introduce $W(d_{1},d_{2})$ by

[TABLE]

Using this notation, the grand-canonical partition function introduced in (3) can be written as (see formula (1.1a) in [21])

[TABLE]

where $D_{m,n}^{s}$ stands for the set of all ordered $s$-tuples of elements of $D_{m,n}$. Unlike (3), formula (10) includes an infinite sum. However, since among every $\lfloor mn/2\rfloor+1$ dimers there exists at least one pair of overlapping dimers, all terms with $s>\lfloor mn/2\rfloor$ vanish. We introduce $F(d_{1},d_{2})=W(d_{1},d_{2})-1$ for every $d_{1},d_{2}\in D_{\infty,\infty}$. Then (10) can be rewritten as (see formula (2.7) in [21])

[TABLE]

and $\mathcal{G}_{s}$ denotes the set of all graphs on $\{1,\ldots,s\}$, and $E(G)$ is the set of edges of a graph $G$. Changing the order of summation, we obtain $\Theta_{m,n}(z)=\sum\limits_{s=0}^{\infty}\frac{z^{s}}{s!}\sum\limits_{G\in\mathcal{G}_{s}}\mathcal{W}_{m,n}(G)$, where

[TABLE]

In [21, p. 1161] it is shown that (see formula (2.11a))

[TABLE]

where $\mathcal{C}_{s}$ is the set of all connected graphs on $\{1,\ldots,s\}$. The above calculations work in quite general context and do not exploit the structure of $D_{m,n}$. Now we will perform a more careful analysis of (12) and (13) in our setting.

For a tuple $\mathbf{d}=(d_{1},\ldots,d_{s})\in(D_{m,n})^{s}$, we construct a graph with vertices labeled $1,\ldots,s$ such that there is an edge between $i$ and $j$ if and only if $d_{i}$ and $d_{j}$ overlap. We call a tuple $\mathbf{d}$ connected if the corresponding graph is connected. The set of connected tuples in $D_{m,n}$ of length $s$ is denoted by $\left(D_{m,n}\right)^{s}_{c}$. For $\mathbf{d}\in\left(D_{m,n}\right)^{s}_{c}$, we define the height (resp., the width) of $\mathbf{d}$ to be the number of rows (resp., columns) having nontrivial intersection with at least one of the dimers in $\mathbf{d}$. We denote it by $h(\mathbf{d})$ (resp., $w(\mathbf{d})$). Two tuples $\mathbf{d}_{1}=(d_{1}^{1},\ldots,d_{s}^{1})$ and $\mathbf{d}_{2}=(d_{1}^{2},\ldots,d_{s}^{2})$ are said to be translation-equivalent if there exists a translation $\pi$ of the plane by some vector such that $\pi(d_{i}^{1})=d_{i}^{2}$ for every $1\leqslant i\leqslant s$. This is an equivalence relation, and we write it as $\mathbf{d}_{1}\sim\mathbf{d}_{2}$.

The following facts follow straightforwardly from the definitions

Lemma 1**.**

1. (i)

For every tuple $(d_{1},\ldots,d_{s})\in\left(D_{m,n}\right)^{s}\setminus\left(D_{m,n}\right)^{s}_{c}$, the corresponding summand in (12) vanishes. 2. (ii)

If $\mathbf{d}_{1}=(d_{1}^{1},\ldots,d_{s}^{1})$ and $\mathbf{d}_{2}=(d_{1}^{2},\ldots,d_{s}^{2})$ are translation-equivalent, then $F(\mathbf{d}_{1},G)=F(\mathbf{d}_{2},G)$ for every graph $G\in\mathcal{G}_{s}$. 3. (iii)

For every connected tuple $\mathbf{d}\in\left(D_{\infty,\infty}\right)^{s}_{c}$, the number of tuples $\mathbf{d}^{\prime}\in\left(D_{m,n}\right)^{s}_{c}$ such that $\mathbf{d}\sim\mathbf{d}^{\prime}$ is exactly

[TABLE]

where $(x)_{+}:=\max(x,0)$.

We denote by $\mathcal{T}_{s}$ a set of tuples in $\left(D_{\infty,\infty}\right)^{s}_{c}$ that contains exactly one representative of every equivalence class of translation-equivalent connected tuples. Due to Lemma 1 we can rewrite (12) as

[TABLE]

For $\mathbf{d}\in\left(D_{m,n}\right)_{c}^{s}$, we define $\mathcal{W}(\mathbf{d})=\sum\limits_{G\in\mathcal{C}_{s}}F(\mathbf{d},G)$. Using this notation and (14), we can rewrite (13) as

[TABLE]

Hence

[TABLE]

Now we want to obtain a similar expression for $\ln\Theta(z)$ defined in (4)

[TABLE]

Since $\lim\limits_{m,n\to\infty}\frac{\left(m-h(\mathbf{d})+1\right)_{+}\cdot\left(n-w(\mathbf{d})+1\right)_{+}}{mn}=1$, we obtain

[TABLE]

We are now ready to deduce Theorem 1 from (15) and (16).

Lemma 2**.**

For every $N\in\mathbb{Z}_{>0}$

[TABLE]

Proof.

By Lemma 1, the coefficient of $\frac{z^{s}}{s!}\mathcal{W}(\mathbf{d})$ is equal to

[TABLE]

If $N<p:=w(\mathbf{d})+h(\mathbf{d})$, the above expression is equal to $0=\binom{(N-w(\mathbf{d})-h(\mathbf{d})+3)_{+}}{3}$. Otherwise, it is equal to

[TABLE]

It can be verified by direct computation using the formula for the sum of squares that the latter expression is equal to $\binom{N-w(\mathbf{d})-h(\mathbf{d})+3}{3}$. This proves the lemma. ∎

Fix some $s\leqslant N-2$ and $\mathbf{d}\in\mathcal{T}_{s}$. We will prove that all summands of the form $\frac{z^{s}}{s!}\mathcal{W}(\mathbf{d})$ on the left-hand side of (9) cancel. Since $\mathbf{d}$ is connected, it contains at least $w(\mathbf{d})-1$ horizontal dimers and at least $h(\mathbf{d})-1$ vertical dimers. Hence, $w(\mathbf{d})+h(\mathbf{d})-2\leqslant s\leqslant N-2$, so $p:=w(\mathbf{d})+h(\mathbf{d})\leqslant N$. This inequality together with Lemma 2 and (16) implies that the coefficient of $\frac{z^{s}}{s!}\mathcal{W}(\mathbf{d})$ on the left-hand side of (9) is equal to

[TABLE]

Expanding the brackets, we verify that this expression is zero for every $N-p\geqslant 0$. This concludes the proof of Theorem 1.

4. Description of the algorithm

4.1. General algorithm.

Combining (8) and Theorem 1, we obtain Algorithm 1, the first version of an algorithm for computing the first $n$ terms of $f_{2}(p)$. Note that

• •

line 1 is correct due to Theorem 1;

• •

line 1 is correct due to (8);

• •

procedure $\operatorname{ComputeTheta}$ is described in subsection 4.2;

• •

procedure $\operatorname{InversePowerSeries}(a(z))$ computes a power series $z(p)$ given a power series $p(z)$ (see [18, Th. 1.8]).

Several improvements can be made:

• •

Computation of $\operatorname{ComputeTheta}(i,j)$ deals with a very long vector of possibly very large numbers (see Section 4.2). In order to fit into the memory, we perform computation modulo several primes and use chinese remaindering and rational reconstruction to obtain the final result (see Section 4.4).

• •

The output of Algorithm 1 with input $n$, let us call it $\tilde{g}_{n}(p)$, coincides with $f_{2}(p)$ only modulo $O(z^{n})$. Nevertheless, the first few nonzero coefficients of $f_{2}(p)-\tilde{g}_{n}(p)$ turn out to satisfy linear recurrence relations with respect to $n$, so they can be computed easily. This allows us to “correct” these terms and obtain a more precise result. See Section 4.3 for further details.

• •

Since we need only the first $n$ terms of $\Theta(z)$, it is sufficient to compute only the first $n$ terms for every computed $\Theta_{i,j}(z)$. Therefore, all intermediate polynomials can also be truncated.

With these improvements, we obtain the final version of our algorithm. For more details, see the source code (see Section 5.1).

4.2. Computation of $\Theta_{m,n}(z)$.

We will compute $\Theta_{m,n}(z)$ using an optimization of the transfer–matrix method (see [22, §4.7]). Fix a positive integer $m$. Let $n$ be a nonnegative integer, and $0\leqslant N<2^{m}$. Viewing $N$ as a vector of $m$ bits, we denote the $i$-th bit of $N$ by $N[i]$. We denote by $F_{N}^{(m,n)}$ the polygon obtained form the $m\times n$ rectangle by adding one additional cell (we will call it an external cell) to the end of every row such that $N[i]=1$, where $i$ is the index of the row. For example, $F_{5}^{(4,6)}$ is shown on Figure 1. In particular, $F_{0}^{(m,n+1)}$ is the same as $F_{2^{m}-1}^{(m,n)}$.

We introduce polynomial $P_{N}^{(m,n)}(z)$ to be a generating function for the number of tilings of $F_{N}^{(m,n)}$ such that every external cell is covered by a horizontal dimer, i.e. $P_{N}^{(m,n)}(z)=\sum\limits_{j=0}^{m(n+1)}a_{N,j}^{(m,n)}z^{j}$, where $a_{N,j}^{(m,n)}$ is the number of monomer-dimer tilings of $F_{N}^{(m,n)}$ with exactly $j$ dimers such that every external cell is covered by a horizontal dimer. We will call such tilings rigid. In particular, $\Theta_{m,n}(z)=P_{0}^{(m,n)}(z)$. We denote by $P^{(m,n)}$ the vector $(P_{0}^{(m,n)}(z),\ldots,P_{2^{m}-1}^{(m,n)}(z))$.

Remark 1**.**

It can be shown (using techniques from [23, §V.6.]), that there exists a matrix $M$ with entries in $\mathbb{Z}[z]$ such that $P^{(m,n+1)}=MP^{(m,n)}$. Hence, $\Theta_{m,n}(z)$ can be computed as the first coordinate of $M^{n}P^{(m,0)}$. However, in our computations $m$ can be any natural number up to $30$, so $M$ can have $2^{30}\times 2^{30}=2^{60}\approx 10^{18}$ entries. Luckily, the matrix $M$ is highly structured (see [24]), so there exists a faster algorithm for computing $P^{(m,n+1)}$ from $P^{(m,n)}$.

We present an algorithm (Algorithm 2) that computes $P^{(m,n+1)}$ from $P^{(m,n)}$ in-place (i.e. with $O(1)$ additional space) using $O(m2^{m})$ arithmetic operations. We denote the number of ones in the binary representation of $N$ by $\operatorname{BinDig}(N)$.

Proposition 1**.**

Algorithm 2 is correct.

Proof.

We will prove by induction on $j$ that after the $j$-th iteration of the loop in line 2 (for $j=0$ it means the moment just before the first iteration) $\widetilde{P}_{N}:=z^{\operatorname{BinDig}(N)}\cdot P[N]$ is the generating polynomial for the number of monomer-dimer tilings of $F_{N}^{(m,n)}$ satisfying the following $A_{j}$ property:

$A_{j}$ property:* the tiling is rigid, and the right-most cell in rows with the number greater than $j$ is covered by a horizontal dimer. *

First we prove the base case, where $j=0$. Due to the loop in line 2, $P[N]=P_{2^{m}-1-N}^{(m,n)}(z)$. Since the binary representation of $2^{m}-1-N$ can be obtained from the binary representation of $N$ by inverting all $m$ bits, adding a horizontal dimer to the end of every row without an external cell provides us a bijection between the set of rigid tilings of $F_{2^{m}-1-N}^{(m,n)}$ and the set of tilings of $F_{N}^{(m,n+1)}$ with $A_{0}$ property (see Fig. 2). This map adds $\operatorname{BinDig}(N)$ new dimers, so the corresponding generating polynomials differ by the factor $z^{\operatorname{BinDig}(N)}$.

Assume now that $j>0$. For $N$ such that $N[j]=1$, properties $A_{j-1}$ and $A_{j}$ are the same, so the corresponding component of vector $P$ should not be changed. Assume that $N[j]=0$. We denote the last cell of the $j$-th row in $F_{N}^{(m,n+1)}$ by $c$. Consider an arbitrary monomer-dimer tilings of $F_{N}^{(m,n+1)}$ with property $A_{j}$. There are three options for $c$:

1. (1)

Cell $c$ is covered by a horizontal dimer. Then this tiling has also property $A_{j-1}$ and is already counted in $z^{\operatorname{BinDig}(N)}P[N]$. 2. (2)

Cell $c$ is covered by a monomer. Replacing this monomer by a horizontal dimer, we establish a bijection between such tilings of $F_{N}^{(m,n)}$ and tilings of $F_{N+2^{m-j}}^{(m,n+1)}$ with property $A_{j-1}$ (see Figure 3). Due to the induction hypothesis, the generating polynomial for the latter is $\widetilde{P}_{N+2^{m-j}}$. Hence, in order to take into account tilings where $c$ is covered by a monomer, we should add $\frac{1}{z}\widetilde{P}_{N+2^{m-j}}$ to $\widetilde{P}_{N}$. This is equivalent to $P[N]\mathrel{+}=P[N+2^{m-j}]$ in line 2.

3. (3)

Cell $c$ is covered by a vertical dimer. This dimer cannot cover also the cell below $c$ due to $A_{j}$ property. Hence, it covers $c$ and the cell above, say $d$, so $N[j-1]=0$. Replacing this dimer with two horizontal dimers, we establish a bijection between such tilings of $F_{N}^{(m,n)}$ and tilings of $F_{N+2^{m-j}+2^{m-j+1}}^{(m,n+1)}$ with property $A_{j-1}$. These cases are counted in line 2.

Since the $A_{n}$ property is just rigidness, after multiplication by an appropriate degree of $z$ in line 2 we obtain the vector $P^{(m,n+1)}$. ∎

Remark 2**.**

Algorithm 2 can be parallelized. Consider an iteration of the loop in line 2 with $j>0$. Then, during the iteration, coordinates of $P[N]$ with different $N[0]$ do not interact, so the whole vector can be divided into two halves (depending on $N[0]$), and these halves can be processed by separate threads. Taking into accout $N[1]$, we can divide the work between four threads, and so on. In our computation, we used $32$ threads (so, we divided the work based on $N[0],\ldots,N[4]$).

Finally, using Algorithm 2, we can write a pseudocode for procedure $\operatorname{ComputeTheta}(m,n)$, see Algorithm 3.

4.3. Correction terms

We can compute more terms of $\Theta(z)$ and, consequently, of $f_{2}(p)$ if we examine carefully the right-hand side of (9). Below we write down the first nonzero term of the right-hand side of (9) for $N=4,5,\ldots$

[TABLE]

Denote the sequence of coefficients by $\{a_{n}\}_{n=1}^{\infty}$. Using Guess package ([25], for introduction to guessing, see [26, §4]) we find that this sequence (we computed first $50$ terms) satisfies the following recurrence relation

[TABLE]

Using (17), we can compute $a_{n}$ easily, so we get one more correct term of $\Theta(z)$. Instead of giving a rigorous proof of (17) which is long and involved, we would like to explain informally why it is natural to expect such a relation.

Formula (16) shows that the coefficient of $z^{s}$ in $\Theta(z)$ is a sum of weights of all connected polyominos constructed from $s$ overlapping dimers. On the other hand, the argument after Lemma 2 shows that the coefficient of $z^{s}$ in

[TABLE]

is a sum of weights over all connected polyominos constructed from $s$ overlapping dimers with the sum of height and width at most $N-2$. Hence, the coefficient of $z^{N-1}$ in their difference is a sum of weights of all connected polyominos constructed from $N$ overlapping dimers with the sum of the height and the width exactly $N-1$ (the sum can not be larger for a connected polyomino). These requirements on a polyomino are quite restrictive, by a combinatorial argument one can see that all such polyominos are “of a similar shape” as those in Figure 5. More precisely, there exist two cells ($a$ and $b$ in the figure), maybe coinciding, such that each of them is connected to two sides of an $m\times n$ ($m+n=N-1$) rectangle by straight lines, and $a$ and $b$ are connected by a path such that at each step the path becomes closer to $b$ (all such paths have the same length). Counting such polyominos is a standard combinatorial problem (similar counting problems for polyominos are discussed in [22, §4.7.5]), that is very likely to result in a formula satisfying a linear recurrence.

Moreover, the same argument shows that there also should be a combinatorial description and a similar recurrence for the second nonzero term in the left-hand side of (9), the third, the fourth and so on. Our data was enough to discover and verify five formulas of this type (from the first until the fifth nonzero term in (9)). This is the recurrence for the second nonzero coefficient

[TABLE]

We omit the others, because they are too large. However, in our program we do not use recurrences themselves, but the closed form expression for their solutions. This allows us to compute five more terms of $\Theta(z)$ and, consequently, of $f_{2}(p)$.

4.4. Modular computation

The largest $n$ we used as an input of the algorithm in our computation was $65$. Taking into account correction terms, this means that $\operatorname{ComputeTheta}$ is invoked with parameters $30$ and $31$. Hence, the vector $P$ in Algorithm 2 will have $2^{30}\approx 10^{9}$ entries. Every entry is a polynomial (in our computations it is a truncated polynomial with only $70$ terms), hence in total we have $7.5\cdot 10^{10}$ integers at every moment. Since these integers represent the number of tilings of a rectangle, they grow fast, so storing them all exactly would require at least several terabytes of memory. However, the final result is a list of coefficients of a power series for $f_{2}(p)$, that is just $65$ rational numbers. A standard way to deal with such situation (see [26, §4.2]) is to use computations modulo prime $p$ for intermediate steps. If $p\leqslant 2^{31}-1$, then all numbers will fit into $32$ bits, and the whole vector $P$ will occupy just $270$ gigabytes. Repeating this computation for different primes, we can reconstruct the coefficients of $f_{2}(p)$ using the chinese remaindering (see [27, §5.4]) and the rational reconstruction procedure (see [27, §5.10]).

The question is how many primes we should take. We start with $2^{31}-1$, and add new prime numbers until the result of the reconstruction stabilizes. It turned out that fifteen prime numbers (from $2^{31}-1=2147483647$ down to $2147483269$) are enough, however we computed several more in order to make sure that the result is correct. The correctness of the result is further justified by the comparison in Section 5.

5. Numerical results and implementation

5.1. Implementation.

We implemented most of our algorithm in Sage except the function $\operatorname{ComputeTheta}$, which was implemented in C. See the source code in github.com/pogudingleb/monomer˙dimer˙tilings. Computation modulo one prime with $n=5$ took about two days using $32$ cores and $270$ gigabytes of memory. Since we need fifteen primes, the whole computation took about one month.

5.2. Numerical results.

Table LABEL:table:ak contains $a_{k}$’s (defied in (1)) obtain by our computation. Expanding $(1-p)\ln(1-p)$ into Taylor series at $p=0$, we obtain the following formula expressing $b_{k}$ defined in (2) via $a_{k}$:

[TABLE]

We introduce following truncated versions of (1) and (2)

[TABLE]

All computed $63$ values $a_{k}$ are positive, all computed $63$ values $b_{k}$ are negative. Assuming that this pattern persists, we can write

[TABLE]

This provides us with lower and upper bound for $f_{2}(p)$. We plot both $L_{64}(p)$ and $U_{64}(p)$ together for $p\in[0,1]$ on Figure 6(a). Dashed curve on this plot is $-p\ln p-(1-p)\ln(1-p)$, that is the negative value of the free energy for monomer-monomer problem with two different types of monomers. We also plot both $L_{64}(p)$ and $U_{64}(p)$ for $p\in[0.9,1]$ on Figure 6(b).

Plots of $L_{64}(p)$ and $U_{64}(p)$ in Figure 6(a) are indistinguishable, the difference between them in Figure 6(b) is visible only very close to $p=1$. On Figure 6(b) we also see that lower bound is much more accurate at $p=1$. The difference $U_{64}(p)-L_{64}(p)$ does not exceed $2.3\cdot 10^{-16}$ for $p\in[0,0.5]$ and $2.1\cdot 10^{-6}$ for $p\in[0,0.9]$. Note that for $U_{24}(p)-L_{24}(p)$ (these two bounds could be computed using results of [16]) these numbers are $9.3\cdot 10^{-11}$ and $7.5\cdot 10^{-4}$, respectively, so our bound reduces the error by several orders of magnitude.

5.3. Comparison with [28].

We already compared our result to the previously known best bound used power series expansion from [16]. However, another method of computing lower and upper bounds for $f_{2}(p)$ based on the empirically observed inequality [28, Eq. 16] for strips was proposed in [28]. In this paper bounds for $p=\frac{1}{20},\ldots,\frac{20}{20}$ were computed (see [28, Table II]). We compare our results with this computation in Table 2. The table shows that for $p$ close to one Kong’s results may be more accurate. On the other hand, our bound is much more precise for $p\leqslant\frac{17}{20}$.