Monomer-Dimer Tilings
F. Cazals, December 1997.
A fundamental problem in lattice statistics is the monomer-dimer problem, in which the sites of a regular lattice are covered by non-overlapping monomers and dimers, that is squares and pairs of neighbor squares. An example of such a tiling for a
chessboard with
and
is depicted below. The relative number of monomers and dimers can be arbitrary or may be constrained to some density
, and the problem can be generalized to any fixed dimension
. This model was introduced long ago to investigate the properties of adsorbed diatomic molecules on a crystal surface [Rob35], and its three-dimensional version occurs in the theory of mixtures of molecules of different sizes [Gug52] as well as the cell cluster theory of the liquid state [CoAl55]. Practically, most of the thermodynamic properties of these physical systems can be derived from the number of ways a given lattice can be covered, so that a considerable attention has been devoted to this counting question. For any fixed dimension
and any monomer density
, a
provably good polynomial time approximation algorithm
is exposed in [KenAl95]. But exact counting results are still unknown even in dimension two.
The goal of this worksheet is to show that these questions are amenable to an
automated computer algebra treatment
which goes from the specifications of the coverings constructions in terms of Combstruct grammars, to the asymptotics using rational generating functions and the numeric-symbolic method exposed in [GoSa96]. In particular we shall be interested in enumerating the tilings for a vertical strip of constant width
in terms of multivariate rational generating functions, from which the average number of pieces or the expected proportions of the three types of pieces in a random tiling are easily derived.
This will also enable us to establish a
provably good
sequence of upper and lower bounds for the connectivity constant
where
counts the number of ways to tile an
cheesboard.
But before getting started, we need to load the Combstruct library, as well as the piece of code doing the asymptotics of rational fractions:
> with(combstruct): with(gfun): read `ratasympt.mpl`;read `./gfsolve.mpl`;
References
[CoAl55] E.G.D. Cohen et al., A cell-cluster theory for the liquid state II, Physica XXI, 1955.
[Fin97]S. Finch, Favorite Mathematical Constants, http://www.mathsoft.com/cgi-shl/constant.bat.
[GoSa96] X. Gourdon and B. Salvy, Effective Asymptotics of linear recurrences with rational coefficients, Discrete Mathematics, Vol. 153, 1996.
[Gug52] E.A. Guggenheim, Mixtures, Clarendon Press, 1952.
[Ken95] C. Kenyon et al., Approximating the number of Monomer-Dimer Coverings of a Lattice, Proc. of the 25th ACM STOC, 1993.
[Rob35]J.K. Robert, Some properties of adsorbed films of oxygen on tungsten, Proc. of the Royal Society of London, Vol. A 152, 1935.
[Sloa95] N.J.A. Sloane and S. Plouffe, The Encyclopedia of Integer Sequences, Academic Press, 1995.
A step-by-step example
We first observe that the number
counting the different tilings of a vertical slice of width 1 has a well known expression: since height
can be reached from height
by adding a monomer and from height
with a vertical dimer, we have
with
, that is the Fibonacci recurrence. This can be checked directly with Combstruct:
> TGr:={T=Sequence(Union(monomer,dimer)),monomer=Z,dimer=Prod(Z,Z)};
And we can retrieve the corresponding rational Generating Function with gfsolve:
> gfsolve(TGr,unlabelled, z);
More interesting is the case m=2 which we examine examine now.
Tiling a slice of width
An example covering of a 2x6 lattice is depicted below. If we draw a horizontal line at height 0, it turns out that we do not `cut' any piece, which we encode by MM. At height 1, we cut the lefmost vertical dimer but just touch the monomer topmost side, which we encode by PM. At height 2 the leftmost P turned into an M since we now touch the dimer boundary, while on the right side we added a dimer and have a P. More generally, we shall
assign to each height of the construction containing a monomer or dimer boundary
a word of length
on the alphabet
as follows: the
th digit of the word is
if an horizontal line at this particular height splits a vertical domino located in the
th column, and
otherwise. To summarize our example we therefore have MM, PM, MP, MM, MM,MM at the heights 0,1,2,3,4,6. (BTW, M stands for Minus and P for Plus!)
This encoding is not one-to-one since whenever we find two consecutive Ms, we do not know wether they are on top of two monomers or of a horizontal dimer. But it is sufficient to incrementally build all the possible configurations by recording the status of the
fringe.
If
, the possible fringes are MM,MP,PM and each of them can be derived from a combination of the others and of monomers and dimers. For example, the configuration MM can be reached in 5 different ways by:
-stacking a horizontal dimer H, two monomers C,C, or two vertical dimers V,V on top of a MM configuration,
-adding a monomer C to the right column of a PM configuration or to the left one of a MP.
The remaining transitions follow similar rules. And in order to characterize the ordinate reached by the construction, we can mark the height reached by the bottommost piece whose elevation gain is 1 or 2 at each step of the construction. Putting everything together and associating the symbols
and
to the number of horizontal dimers, vertical dimers, monomers and the height yields the following Combstruct grammar:
>
Gr2:={MM=Union(Epsilon,Prod(S, MM, H), Prod(S, MM, C,C),
Prod(S,PM, C),Prod(S,MP, C),Prod(S,S,MM,V,V)),
PM=Union(Prod(S,MM, V, C), Prod(S,MP,V)),
MP=Union(Prod(S,MM, C, V), Prod(S,PM,V)),
H=Epsilon,V=Epsilon,C=Epsilon,S=Atom}:
The ordinary generating functions can be derived by Combstruct[gfsolve]:
> GF2Sys:=gfsolve(Gr2, unlabelled, z, [[h,H], [v,V], [c,C]]);
Furthermore we can isolate the GF corresponding to the MM fringes; the coefficient of
in this GF counts the number of ways to tile a chessboard
with respectively
and
horizontal and vertical dimers and monomers:
> GF2:=subs(GF2Sys,MM(z,h,v,c));
The number of configurations up to a given height independently of the number and kind of pieces used can be retrieved by erasing the dimers and monomers markers followed by a Taylor expansion:
> GF2h:=subs([h=1,v=1,c=1],GF2);series(GF2h,z=0,11);
This sequence does not appear in [Sloa95]. It can be checked that these values match those computed directly from the grammer by Combstruct[count]:
> seq(count([MM,Gr2], size=i), i=0..10);
Another way to compute the exact number of tilings for large values of
is through the recurrence equation satisfied by the Taylor coefficients and computed by gfun[diffeqtorec]:
> diffeqtorec(y(z)-GF2h,y(z),u(n));
> p2:=rectoproc(",u(n)):
> for i from 1 to 10 do i,p2(i) od;
For example:
> p2(1000);evalf(");
But as we shall see now, asymptotic estimates can be derived much faster.
Asymptotic estimates of the number of tilings
We have just seen that the number of configurations is encoded by the rational generating function GF2h(z). An elegant way to access its Taylor coefficients is therefore through a full partial fraction decomposition yielding linear denominators:
> fpf:=convert(GF2h,fullparfrac,z);
The term in
comes from the contributions of the roots of
in the expansion of
> el:=op(1,fpf);
and since there are 3 singularities, the main asymptotic contribution comes from the one with smallest modulus:
> fsolve(-3*_Z+1+_Z^3-_Z^2,_Z);
>
root1:=RootOf(-3*_Z+1+_Z^3-_Z^2,.3111078175);
root2:=RootOf(-3*_Z+1+_Z^3-_Z^2,-1.481194304);
root3:=RootOf(-3*_Z+1+_Z^3-_Z^2,2.170086487);
On this example the dominant pole is clearly
so that the main contribution is encoded by:
> el1:=subs(_alpha=root1,el);
> evalf(");
Extracting the term in
in the previous expression produces the estimate:
> es2:=n->.2067595751*(1/.3111078175)^(n+1);
> seq(es2(i),i=1..10);
To sum up, from the rational generating function we have:
-performed a full partial fraction decomposition,
-computed the singularities and sorted them by increasing moduli,
-extracted the contribution of the singularity with smallest modulus.
The key step consists in deciding which are the singularity (ies) with smallest modulus (i), and can be performed numerically using properties of polynomials with integer coefficients --see [GoSa96]. This is implemented by the
ratasympt
function --whose optional
th argument corresponds to the number of singularity layers the user wants to take into account. In particular to retrieve the main contribution, one writes:
> layer1:=ratasympt(GF2h,z,n,1);nbCfs1:=evalf(layer1);
And to take into account all the layers:
> layers:=ratasympt(GF2h,z,n);nbCfs:=evalf(layers);
We can check that the second approximation is more accurate:
> evalf(seq(subs(n=i, layer1), i=1..10));
> evalf(seq(subs(n=i, layers), i=1..10));
> seq(p2(i),i=1..10);
The proportion of monomers and dimers
We now address the computation of the average number of pieces in a random tiling. From the multivariate generating function
we can merge the three types of pieces as follows:
> GF2;stij:=subs([h=t,v=t,c=t], GF2);
The coefficient of
in
counts the number of tilings at height
with exactly
pieces of any type. To get the total number of pieces we just have to compute the derivative with respect to
and substitute
:
> sstij:=subs(t=1, diff(stij,t));
For example, the total number of dimers and monomers used in all the configurations tilling the square
is 20:
> series(",z=0,5);
As before, we can compute an estimate of the total number of pieces in all the configurations at a given height:
> ratasympt(sstij,z,n,1);
> nBDPiecesN:=evalf(");
So that the average number of pieces is asymptotically equivalent to:
> avNbD:=expand(nBDPiecesN/nbCfs1);
And the average number of pieces per layer in a tiling of height
is therefore:
> asympt("/n,n);
The number of occurrences and the proportions of dimers and monomers can be computed in the same way by erasing the irrelevant indeterminates:
>
pieceProportion:=proc(MGF, keptPiece)
local forSubs, stij, sstij, nbp;
forSubs:={h=1,v=1,c=1} minus {keptPiece=1};
stij:=subs([op(forSubs)], MGF);
sstij:=subs(keptPiece=1, diff(stij,keptPiece));
nbp:=evalf(ratasympt(sstij,z,n,1));
asympt(nbp/nBDPiecesN,n,2)
end:
And we end up with:
>
nbh:=pieceProportion(GF2,h);
nbv:=pieceProportion(GF2,v);
nbc:=pieceProportion(GF2,c);
Plotting routines archive
The figures above were plotted with the following functions:
>
dominoH:=proc(x,y) [[x,y], [x+2,y], [x+2,y+1], [x,y+1], [x,y]] end:
dominoV:=proc(x,y) [[x,y], [x+1,y], [x+1,y+2], [x,y+2], [x,y]] end:
dominoC:=proc(x,y) [[x,y], [x+1,y], [x+1,y+1], [x,y+1], [x,y]] end:
>
plot([dominoV(0,0), dominoC(0,2),dominoC(0,3),
dominoC(1,0),dominoV(1,1),dominoH(1,3),
dominoV(2,0),dominoC(2,2),
dominoC(3,0),dominoC(3,1),dominoV(3,2),
dominoH(0,4),dominoH(2,4)],scaling=constrained,color=blue);
>
> plot([dominoV(0,0), dominoC(1,0),dominoV(1,1),dominoC(0,2),dominoH(0,3),dominoV(0,4),dominoV(1,4)], scaling=constrained,color=blue);
Automatic counting in a slice of width
Computing the generating functions
We now show how to automate the previous computations for any integer
. The first task consists in generating the
words on the binary alphabet
, and this is easily done with a Combstruct grammar as follows:
>
allMPWords:=proc(m::integer)
local i, MPGr, mps1, mps2, Pm;
MPGr:={AllMP=Sequence(MP), MP=Union(M,P), M=Atom, P=Atom};
mps1:=allstructs([AllMP, MPGr], size=m);
mps2:=convert(map(proc(x) cat(op(x)) end, mps1), set);
Pm:=cat(seq(P,i=1..m));
[op(mps2 minus {Pm})]
end:
For example if
:
> allMPWords(3);
More interesting is the generation of the transitions between these words. Let
be one of them and suppose we want to figure out all the fringes
can be derived from. Suppose for example the
th letter of
is a
; this means that the
th letter of the fringe
was derived from was
and that a vertical dimer was put on top of this
. Similar rules applies if the
th digit is a
. And since the letter of a given fringe are independent --except for two consecutive
s that may come from an horizontal dimer, it suffices to recursively examine the digits from left to right as follows:
>
#--pattern is the fringe to be built, e.g. MMPMM
recComesFrom:=proc(pattern::string, idx::integer, prefix::string, mul::list, result::table)
local prodRes, m, Mm, Pm;
if (idx>length(pattern)) then #--stores the result into an indexed table
prodRes:=Prod(S,prefix, op(mul));
if not assigned(result[pattern]) then result[pattern]:={prodRes}
else result[pattern]:=result[pattern] union {prodRes}
fi
else
#--we examine the idx^{th} letter of the target
if substring(pattern,idx)=P then
recComesFrom(pattern, idx+1, cat(prefix,M), [op(mul), V], result)
else #target=M
recComesFrom(pattern, idx+1, cat(prefix,P), mul, result);
recComesFrom(pattern, idx+1, cat(prefix,M), [op(mul), C], result);
#--we may have MM=Prod(MM,H)
if (length(pattern)>idx) and (substring(pattern,idx+1)=M) then
recComesFrom(pattern, idx+2, cat(prefix,M,M), [op(mul), H], result)
fi
fi
fi;
#--some extra work for M^m
m:=length(pattern);
Mm:=cat(seq(M,i=1..m));
if pattern=Mm then
Pm:=cat(seq(P,i=1..m));
result[Mm]:=result[Mm] minus {Prod(S,Pm)}
union {Epsilon,Prod(S,S,Mm,seq(V,i=1..m))}
fi
end:
Here is the table for
:
> table3:=table():for i in allMPWords(3) do recComesFrom(i, 1, ``, [], table3) od:print(table3);
The tables entries are merged as follows:
>
setGrammarFromTable:=proc(aTable)
local aList, transitions, x;
aList:=op(op(aTable));#--[a={Prod(...), Prod(...)}, ...]
transitions:=seq(op(1,x)=Union(op(op(2,x))), x=aList);
{transitions} union {H=Epsilon,V=Epsilon,C=Epsilon,S=Atom}
end:
This yields the grammar:
> Gr3:=setGrammarFromTable(table3);
This is solved as usual:
>
MM3GFSys:=gfsolve(Gr3, unlabelled, z, [[h,H], [v,V], [c,C]]);
MM3GF:=subs(MM3GFSys,MMM(z,h,v,c));
Putting everything together, we end up with a procedure which takes
as entry and returns the grammar:
>
getGrammar:=proc(m::integer)
local i, MPTable;
MPTable:=table();
for i in allMPWords(m) do recComesFrom(i,1,``,[],MPTable) od;
setGrammarFromTable(MPTable)
end:
getMmGFun:=proc(m::integer)
local i, MPTable,Grm,MMmGFSys;
Grm:=getGrammar(m);
MMmGFSys:=gfsolve(Grm, unlabelled, z, [[h,H], [v,V], [c,C]]);
subs(MMmGFSys,cat(seq(M,i=1..m))(z,h,v,c));
end:
The computation to be carried out being quite heavy for 4-variate generating functions, we can alleviate it be keeping only the markers for the total number of pieces and the height:
>
getMmGFunZ:=proc(m::integer)
local i, MPTable,Grm,GrmM,MMmGFSys;
MPTable:=table();
for i in allMPWords(m) do recComesFrom(i,1,``,[],MPTable) od;
Grm:=setGrammarFromTable(MPTable);
MMmGFSys:=gfsolve(Grm, unlabelled, z);
subs(MMmGFSys,cat(seq(M,i=1..m))(z));
end:
Asymptotics
We can now compute the generating functions for small values of
:
> gf:='gf':
> for i from 1 to 5 do i,time(assign(gf[i],getMmGFunZ(i))),gf[i] od;
For bigger ones, the grammar size, that is
, inherently yields a linear system
with large coefficients whose resolution is very much time consuming. So that for
, a better alternative to running getMmGFunZ(m) is to retrieve the result in the archive below!
From these generating functions we can easily isolate the main contribution to the asymptotic equivalent with the ratasympt procedure:
> asGf:='asGf':
> for i from 1 to 6 do assign(asGf[i],ratasympt(gf[i],z,n,1)),evalf(asGf[i]) od;
It should be observed that these estimates correspond to huge expressions. For
for example:
> asGf[5];
As observed in [Fin97], if
denotes the number of tilings of a
chessboard, an interesting value for the physical applications is
.
No exact expression for this limit is known, although the approximation 1.940215531 is generally agreed on. The first terms of the sequence can be computed from the previous approximations and are consistent with 1.94:
> nn:='nn':
> for i from 1 to 6 do assign(nn[i],coeff(series(gf[i],z=0,i+1),z,i)),evalf((nn[i])^(1/(i*i))) od;
But more interesting is the following observation. Suppose for example
is a multiple of 6. To tile a
chessboard we can put side by side
slices of width 6. In this case
with
the singularity of smallest modulus of the denominator of
. If
is not a multiple of 6, it suffices to complete with at most 5 vertical stripes of width 1, but this does not change the limit. The interest in using as many slices of maximal width is to minimize the number of joints where the overlaps are not taken into account. The sequence
therefore provides lower bounds for the constant
. An upper bound can be obtained in the same way by having slices of width 6 overlap on a position, and the corresponding sequence is
.
> for i from 2 to 6 do i,(1/op(1,denom(evalf(asGf[i]))))^(1./i),(1/op(1,denom(evalf(asGf[i]))))^(1./(i-1)) od;
At last a trick we can use to try to guess the value of
is Romberg's convergence acceleration. Let
be a sequence known to converge to
. If the rate of convergence is of the form
, then
is
. On our example, although the upper bound does not make sense due to too erroneous initial values, after a single step the lower bound gets close to the commonly accepted value:
>
u[2]:=1.792852404:u[4]:=1.863497010:
v[2]:=3.214319743:v[4]:=2.293180643:
2*u[4]-u[2],2*v[4]-v[2];
u[3]:=1.838281935:u[6]:=1.888704987:
v[3]:=2.492402505:v[6]:=2.144850135:
2*u[6]-u[3],2*v[6]-v[3];
Generating functions archive
>
gf[1]:=-1/(-1+z^2+z);
> gf[2]:=-1/(-3*z+1-z^2+z^3)*(z-1);
> gf[3]:=-(z^4+z^3-4*z^2-z+1)/(14*z^2-1+4*z+z^6-10*z^4);
> gf[4]:=-(1-4*z-15*z^2+20*z^3+z^7-11*z^5-2*z^6+10*z^4)/(z^9-z^8-23*z^7+29*z^6+91*z^5-111*z^4-41*z^3+41*z^2+9*z-1);
> gf[5]:=-(z^18+2*z^17-45*z^16-68*z^15+654*z^14+870*z^13-3820*z^12-4700*z^11+9255*z^10+9448*z^9-11175*z^8-7532*z^7+6956*z^6+1994*z^5-1794*z^4-88*z^3+113*z^2+6*z-1)/(z^20+2*z^19-65*z^18-140*z^17+1281*z^16+2538*z^15-10366*z^14-17604*z^13+38553*z^12+50158*z^11-73623*z^10-60482*z^9+74665*z^8+26564*z^7-35106*z^6-898*z^5+4757*z^4+16*z^3-229*z^2-14*z+1);
> gf[6]:=-(-1+311*z^2-3891*z^3-12057*z^4-315889*z^6-2997721*z^7+218447*z^5+13467571*z^9+8754480*z^8+23*z-458919487*z^18-303976032*z^17+612805499*z^16+207743591*z^15-496137395*z^14-56233657*z^13+240612231*z^12-14684235*z^11-66016499*z^10+206819317*z^20+249194245*z^19-109*z^32-36273*z^29+861*z^31+7443809*z^24+37223601*z^23-123372421*z^21-54160427*z^22-6708699*z^25+z^34-29377*z^28+686517*z^27-338040*z^26+3521*z^30-7*z^33)/(1-576*z^2+6080*z^3+42422*z^4-443404*z^6+12931566*z^7-453004*z^5-83558644*z^9-25517604*z^8-36*z+4169343006*z^18+2978277152*z^17-4669345206*z^16-1630080704*z^15+3235975264*z^14+274712602*z^13-1335612340*z^12+154307596*z^11+295510396*z^10-2310327672*z^20-2919950172*z^19+5736*z^32+1503868*z^29-62874*z^31-149620588*z^24-626694028*z^23+1717916424*z^21+777289050*z^22+141424642*z^25-8*z^35-138*z^34+z^36-94620*z^28-19237868*z^27+13835164*z^26-81796*z^30+1224*z^33);
>
>
>
Conclusion
We showed that various parameters related to dimer-monomer tilings such as the average number of pieces or the relative numbers of horizontal dimers and monomers in a random tiling of height
in a strip of width
can be computed very easily using Combstruct and ratasympt. More precisely Combstruct is used to define the grammars the tilings are derived from, and ratasympt is used to perform asymptotic expansions on rational fractions with rational coeficients.
About the number
of different tilings of a
chessboard, altough the method presented here is limited due to the exponential growth of the grammar describing these tilings, the very first terms computed provide
provably good upper and lower bounds
for the connectivity constant
. More precisely:
Theorem . The connectvity constant for two dimensional monomer-dimer tilings satisfies
and