Entropy of three-dimensional integer partitions and related three-dimensional rhombus tilings

M. Widom^{1}, R. Mosseri^{2}, N. Destainville^{3}

1Carnegie Mellon U. Department of Physics, Pittsburgh, PA 15213, USA 2GPS, Université Paris 7 et 6, 2 place Jussieu, 75251 Paris Cedex 05, France 3Laboratoire de Physique Quantique, Université Paul Sabatier, 31062 Toulouse Cedex 04, France.

Random tiling models show that quasiperiodicity can arise spontaneously as the maximum entropy state of a random covering of space by rigid objects. Actually calculating the entropy for a given model is a notoriously difficult problem. We develop “transition matrix” Monte Carlo methods to calculate random tiling entropies, and apply them to the study of codimension 1 random tilings which are equivalent to integer partitions in three dimensions. We compare the entropies of “fixed” and “free” boundary tilings and show this ratio equals the ratio of volumes of two simple polyhedra.

I. INTRODUCTION

Since the discovery of quasicrystals in 1984 [1], many models for quasicrystal structure invoke tilings (coverings of space without gaps or overlaps) by rhombahedra, presuming these tiles to be decorated in some manner by atoms. Some proposals for the formation and stability of quasicrystals presume energetic interactions among tiles that force quasiperiodicity, while others note the configurational entropy available to noninteracting tilings and note that high entropy also favors quasiperiodicity. In both models, the diffraction pattern contains delta function Bragg peaks.

The latter class of models are known as random tiling models. Proof of the concept depends on the ability to calculate the entropy, a well known difficult problem in statistical mechanics. This paper presents a new method for entropy calculation based on transition matrix Monte Carlo. We apply the method to three dimensional, codimension 1 tilings that can be described as projections from four dimensional space into three dimensional space. They can also be considered as integer partitions. Integer partitions are sets of integers arranged in more or less complex arrays, together with suitable inequalities between these integers. The basic idea is that a one-to-one correspondance can be established between such partitions and tilings of rhombus filling some well-defined polytopes.

Unfortunately, it was rapidly understood that such exotic boundaries have a non-trivial macroscopic effect on random tilings [2], leading to a lower entropy per tile at the large-size limit than in the free or periodic boundary case. This phenomenon was later well described, at least in the simple case of hexagonal tilings [3,4]. It is known as the “arctic circle phenomenon”: the constraint imposed by the boundary is so strong that there are macroscopic frozen regions near the boundary. It is only near its very center that the tiling has “forgotten” the presence of the boundary. The unfrozen region appears to be a circle inscribed in the hexagonal boundary in this case.

However, this quantitative result has never been generalized to higher dimension or codimension tilings, because its generalization requires the knowledge of the free boundary entropy. Recently, it has been conjecture by Propp [3], on the basis of small numerical simulations that in dimensions higher than 2 the corresponding arctic region should be a polytope itself, with flat boundaries.

The present paper is devoted to a numerical inspection of codimension one tilings. On the one hand, thanks to a powerful transition matrix Monte Carlo algorithm, we are able to have precise estimates of both fixed boundary and free boundary entropies. The latter is calculated via a modified partition method, which produces tilings with fixed boundaries but which do not impose any strain to the tilings, thus generalizing a former two-dimensional approach [4]. On the other hand, comparing both entropies, we are able to validate Propp’s conjecture with good confidence.

II. PARTITIONS, TILINGS AND BOUNDARY CONDITIONS

The correspondance between codimension-one tilings and solid partitions is a reformulation of the cut-and-project method. It was analysed in detail in reference [ [5]]. There is a generalization of this correspondance to higher codimensions (see [ [6]] and references therein), but we will focus here on codimension one.

Solid partitions are defined as follows [5]: let us consider a three-dimensional array of sides k × l × m. Fix an integer p > 0, called the height of the partition problem. Put non-negative integers in that array, smaller than p, with the following constraints: these integers are decreasing in each of the 3 directions of space. Each integer partition is equivalent to one tiling of the rhombic dodecahedron “RD” (outer frame in Fig. II) by rhombahedra. Because of the fixed boundaries we call the entropy of this tiling σ_{fixed}.

FIG. 1. The puckered octahedral boundary conditions. The octahedron is inscribed in the rhombic dodecahedron coming from the solid partition method.

We also consider “free” boundary tilings. The boundaries of these tilings actually are fixed, but are inside a frame with zero phason strain. It is presumed that the entropy per vertex of a free boundary tiling equals the entropy per vertex of any tiling with vanishing global phason strain in the thermodynamic limit. These boundary conditions are illustrated in figure II, and they fill up an octahedral “O” region (filled region in Fig. II). We call the entropy of this tiling σ_{free}. Note that the volume ratio V_{RD}/V_{O} = 3/2.

FIG. 2. Boundary conditions for the partition problem associated with the strain free tiling problem with boundary conditions of Fig. II. Two views of the truncated cube are provided: two “corners” have been cut. The marked vertex is the same in both views. On the other hand, we have indicated the minimal and maximal values of parts.

III. MONTE CARLO

Conventional Monte Carlo simulations using Metropolis sampling [9] are useful for generating typical configurations. They can reveal temperature driven phase transitions [13] and yield quantitative evaluations of phason elastic constants from fluctuations. Evaluation of entropy by Monte Carlo simulation demands specialized techniques. Advanced methods using specialized dynamics based on entropic sampling [12] and transfer matrices [11] yield reasonably accurate free energies and absolute entropies.

We employ a variant of the transition matrix method [7] that couples a conventional Metropolis Monte Carlo simulation with a data collection and analysis scheme that constructs an approximate transition matrix. The density of states is an eigenvector of the exact transition matrix, and its sum equals the total number of states. This method yields highly accurate absolute entropies with impressive efficiency. The method itself is described in section III A. Results described in section III B include the disproof of a conjectured value of the entropy of fixed boundary 4 → 3 tilings, as well as strong support for a conjecture on the ratio of free boundary to fixed boundar entropies for the same tiling model.

A. Transition Matrix

For any legal partition P = {n_{ijk}} we define its “energy” as its total height

E(P) = ∑_{ijk} n_{ijk}  (1)

and incorporate the constraint of legality by defining E = ∞ for any partition P that violates the ordering conditions inside the partition or at the boundaries. Under single vertex flip dynamics a single partition element increases or decreases by one unit resulting in an energy change ΔE = ±1. The gound state of this model is the lowest height legal partition. For the boundary conditions employed here, the ground state is unique and has energy E_{min}.

In a Metropolis Monte Carlo simulation of the partition, a randomly chosen partition element is randomly increased or decreased by one. The change is accepted if it lowers the energy (i.e. we attempted to decrease it, and the result is a legal partition). If the change raises the energy, it is accepted with probability exp (−E/T). Note that illegal partitions are never accepted. Low temperatures bias the partitions towards low heights, while high temperatures remove any bias according to height. As T → ∞ the Metropolis simulation faithfully reproduces the random tiling ensemble in which all configurations have equal weight. If we choose, we may take negative values of temperature T. Large negative temperatures again reproduce the random tiling ensemble, while small negative temperatures favor partitions of maximal height.

For a given partition P, a number of partitions n_{±}(P) can be reached by single upwards or downwards steps. Additionally, a certain number n_{0}(P) of steps are forbidden due to the partition constraints. The sum rule

n_{−}(P) + n_{0}(P) + n_{+}(P) = N  (2)

holds for every partition P. In a single Monte Carlo step, the probability for upwards or downwards transitions at infinite temperature is n_{±}(P)/N, and the probabilty a move will be rejected is n_{0}(P)/N. At finite temperatures the transition probabilities must be multiplied by appropriate Boltzmann factors.

During a Monte Carlo simulation, many different partitions with a given energy may be sampled. We can calculate exactly the infinite temperature acceptance ratios of transitions by averaging over states with energy E,

ω_{±,0}(E) = \frac{1}{W (E)} ∑_{P (E)} n_{±,0}(P)/N  (3)

where the sum is over all partitions of energy E and the normalization W (E) is the total number of partitions with energy E. In general the set P (E) and the value W (E) are not known. However, we may calculate the rates approximately by averaging n_{±,0}(P)/N over those partitions P occuring during a Monte Carlo simulation. The acceptance ratios obey the sum rule

ω_{−}(E) + ω_{0}(E) + ω_{+}(E) = 1  (4)

for any energy E.

Consider the infinite temperature detailed balance condition. At infinite temperature the probability that the energy of partition has value E is W (E)/W, where W is the total number of legal partitions. Multiplying this probability by the acceptance rate of upwards transitions ω_{+}(E) yields the total forwards transition rate. The backwards rate is obtained in similar fashion. Detailed balance requires that the total rate of forward transitions from energy E to energy E + 1 must equal the total rate of backwards transitions, hence

ω_{+}(E)W (E)/W = ω_{−}(E + 1)W (E + 1)/W.  (5)

It is useful to rearrange the detailed balance eq. (5) to find

W (E + 1) = \frac{ω_{−}(E + 1)}{ω_{+}(E)}W (E)  (6)

which allows us to extract the full density of states function W (E) using uniqueness of the ground state, W (E_{min}) = 1.

Note that these infinite temperature transition rates may be calculated from finite temperature Metropolis simulations. Indeed, certain low energies E may occur very rarely in high temperature simulations so we must perform low temperature simulations to generate partitions with energy E, from which the infinite temperature transition rates may be calculated. By sweeping over a range of temperatures we obtain accurate values of ω_{±}(E) for all energies. Then, we use the infinite temperature detailed balance condition to extract the density of states even though all of our calculations are performed at finite temperature.

The accuracy of our result is controlled by the accuracy with which we determine ω_{±}(E). For each partition, the transition numbers n_{±,0}(P) are calculated exactly and added into the row of the transition matrix with energy E(P). Since there may be many partitions with the same energy E, each having different transition numbers, the accuracy with which a row of the transition matrix is determined is limited by our ability to generate a representative sample of partitions. The time required to calculate n_{±,0}(P) is proportional to the partition size N., and hence to the time required to attempt N flips.

We store the matrix as the integer valued sum of the integers n_{±,0}(P), and impose the normalization (4) after data collection is complete. Because the values of n_{±,0}(P) can be quite large, and we visit each energy a very large number of times, we ordinary 4 byte integers could not hold the data. We implemented special procedures to handle storage and algebraic manipulations of large integers.

B. Numerical Data

Fixed boundary partitions are initialized at zero height (and thus E = E_{min} = 0). Their maximum energy E_{max} = n_{x}n_{y}n_{z}h. We accumulate data in the transition matrix through four sweeps over temperature: from T_{min} to T_{max} during which time the mean energy grows from E_{min} to about E_{max}/2; from −T_{max} to −T_{min} during which time the mean energy grows from E_{max}/2 to E_{max}; from −T_{min} back to −T_{max}; from T_{max} back down to T_{min}. Each sweep visits 201 temperatures in a geometric sequence from the initial to the final temperature. Free boundary tilings are initialized with a flat partition at energy close to (E_{min} + E_{max})/2. We again perform four sweeps: from T_{max} down to T_{min}; from T_{min} back up to T_{max}; from −T_{max} to −T_{min}; from −T_{min} back to −T_{max}.

For both fixed and free boundary tilings, symmetries of the model dictate that the density of states is symmetric about the midpoint energy E_{mid} = (E_{min} + E_{max})/2. The density of states has a strong maximum at the midpoint.

By including both positive and negative temperature sweeps we sample both low and high energy states. By reversing our sweeps we mitigate possible systematic sampling errors associated with the direction of the sweep. At each temperature during a sweep, we accumulate data on N_{MC} sample configurations. Because the data accumulation is time consuming, we can perform many single vertex flips between sample without significantly slowing down the simulation. We take N_{fl} = N/5 which yields rough equality between time spent flipping and collecting data. Prior to data collection at any temperature, we anneal at fixed temperature for N_{MC} ∗ N_{fl}/100 steps. The total number of attempted flips in an N_{MC} = 10^{6} run on a 10x10x10 partition is thus in excess of 4 × 201 × 10^{6} × 10^{3}/5 = 1.6 × 10^{12} attempted flips. A run of this length takes 5 days on a 1.7GHz Pentium 4 processor.

This protocol was chosen to ensure approximately uniform coverage of energies from E_{min} to E_{max}. The energy distribution of partitions visited at each temperature overlaps strongly the energy distribution of partitions visited at the previous and following temperatures. The value of T_{min} is chosen sufficiently low to guarantee some coverage of the extreme energy states. Because there are few sttes at the energy exremes there is an entropic barrier to reaching these extremes. This causes a spike in the coverage close to the extreme energies that cannot be avoided using Metropolis sampling. The value of T_{max} is chosen sufficiently high that the midpoint energy is a local maximum in the coverage.

Densities of states for fixed and free partitions reveal the expected symmetry and strong maximum around E_{mid}. The degree to which symmetry is broken can be used as an indicator of errors errors accumulated during our iterative calculation (6), because sampling errors at low and high energies need not cancel. By inspection it can be seen that this error is quite low.

Table I shows the convergence of entropy data as run length grows. We show the convergence only for the largest system sizes, which represent our worst case. Table II summarizes our data for all system sizes studied. We report the entropy resulting from the run of length 10^{6}, and use the difference between that run and the length 10^{5} run as an estimate of the uncertainty.

To extrapolate our data to infinity, we fit to functional forms. For the free boundary entropy, we expect finite-size corrections to the entropy of order 1/L, because of boundary effects. In the corresponding two-dimensional case, this form appeared to perfectly describe finite-size corrections. In our case, the data fit to σ_{free} ≃ 0.20. For the fixed boundary data we expect two leading correction terms, of order log(L)/L and 1/L respectively, which have the same order of magnitude for the values of L considered here. Indeed, these two terms are the first correction terms in the excatly solvable one- and two-dimensional random tilings and we extrapolate them to the three-dimensional ones. In particular, this form happens to fit the data better than a simple correction of order 1/L. We find σ_{fixed} ≃ 0.145. For the ratio, we fit to \frac{σ_{free}}{σ_{fixed}} ≃ 1.49.

Now let us return to Propp’s conjection [3]. If the artic circle phenomenon indeed causes a total freezing in the thermodynamic limit, then the total entropies of fixed and free tilings should be equal, and the ratio of entropies per vertex should equal the ratio of polytope volumes V_{RD}/V_{O} = 3/2. From the data, we see this conjecture is strongly supported.

[1] D. Shechtman, et al. Phys. Rev. Lett. 53,1951 (1984).

[2] V. Elser. – Solution of the Dimer Problem on an Hexagonal Lattice with Boundary, J. Phys. A: Math. Gen. 17, 1509 (1984).

[3] H. Cohn, M. Larsen, J. Propp, The Shape of a Typical Boxed Plane Partition, New York J. of Math. 4, 137 (1998); H. Cohn, R. Kenyon, J. Propp, A Variational Principle for Domino Tilings, Journal of the AMS 14, 297 (2001).

[4] N. Destainville, Entropy and boundary conditions in random rhombus tilings, J. Phys. A: Math. Gen 31, 6123 (1998).

[5] N. Destainville, R. Mosseri, F. Bailly, Configurational Entropy of Codimension-One Tilings and Directed Membranes, J. Stat. Phys. 87, Nos 3/4, 697 (1997).

[6] N. Destainville, R. Mosseri, F. Bailly, Fixed-boundary octagonal random tilings: a combinatorial approach, J. Stat. Phys. 102, 147 (2001).

[7] J.S. Wang, T.K. Tay and R.H. Swendsen, Phys. Rev. Lett. 82, 476 (1999); R.H. Swendsen, B. Diggs, J.-S. Wang, S.-T. Li and J.B. Kadane, Int. J. Mod. Phys. C10, 1563 (1999);

[8] M. Widom Phys. Rev. Lett. 70, 2094 (1993); P. Kalugin J. Phys. A 27, 3599 (1994); J. de Gier and B. Nienhuis, Phys. Rev. Lett. 76, 2918 (1996); J. de Gier and B. Nienhuis, J. Phys. A 31, 2141 (1998);

[9] K.J. Strandburg, L.-H. Tang and M.V. Jaric Phys. Rev. Lett. 63, 314 (1989); L.-H. Tang, Phys. Rev. Lett. 64, 2390 (1990); D. Joseph and M. Baake, J. Phys. A 29, 6709 (1996) M. Oxborrow and C.L. Henley, Phys. Rev. B 48, 6966 (1993); M. Luby, D. Randall and A. Sinclair, Proc. 36th Ann. Symp. on the Foundations of Comp. Sci. (IEEE, 1995) p. 150

[10] M. Widom, D.P. Deng and C.L. Henley Phys. Rev. Lett. 63, 310 (1989); W. Li, H. Park and M. Widom, J. Stat. Phys. 66, 1 (1992); M.E.J. Newman and C.L. Henley, Phys. Rev. B 52, 6386 (1995)

[11] L.J. Shaw and C.L. Henley, J. Phys. A 24, 4129 (1991);

[12] F. Gahler, Proc. ICQ5, 236 (1995);

[13] H.-C. Jeong and P.J. Steinhardt, Phys. Rev. B48, 9394 (1993);

TABLE I. Convergence of entropy data for increasing run length. The residual R = log W (E_{max}) measures the cumulative error in W (E).

N  10^{4}  10^{5}  10^{6}
σ_{fixed}  0.147385  0.147400  0.147349
R_{rixed}  0.180  0.429  0.002
σ_{free}  0.197300  0.197234  0.197273
R_{free}  0.143  -0.202  -0.036

TABLE II. Size-dependent entropies. Values in parentheses are uncertainties in final digit. Values without uncertainties are exact.

L  σ_{free}  σ_{fixed}  σ_{free}/σ_{fixed}
1  0.1732868  0.1732868  1.000
2  0.1732868  0.1601239  1.080
3  0.17947(2)  0.1545769  1.161
4  0.18455(6)  0.1516677  1.217
5  0.18829(3)  0.1502(1)  1.254
6  0.19108(4)  0.1492(1)  1.281
7  0.193(1)  0.1486(1)  1.299
8  0.19486(2)  0.147997(3)  1.317
9  0.196183(1)  0.14762(2)  1.329
10  0.19727(4)  0.14735(5)  1.339
∞  0.215(5)  0.145(1)  1.48(3)