Arctic octahedron in three-dimensional rhombus tilings and related integer solid partitions

M. Widom^{1}, R. Mosseri^{2}, N. Destainville^{3,†}, F. Bailly^{4}

1Carnegie Mellon U. Department of Physics, Pittsburgh, PA 15213, USA. 2Groupe de Physique des Solides, Universités Paris 6 et 7, 2 place Jussieu, 75251 Paris Cedex 05, France.

3Laboratoire de Physique Quantique, UMR CNRS-UPS 5626, Université Paul Sabatier, 31062 Toulouse Cedex 04, France.

4Laboratoire de Physique du Solide-CNRS, 92195 Meudon Cedex, France.

                     †Corresponding author; e-mail: Nicolas.Destainville@irsamc.ups-tlse.fr

(January 18, 2002)

Three-dimensional integer partitions provide a convenient representation of codimension-one three-dimensional random rhombus tilings. Calculating the entropy for such a model is a notoriously difficult problem. We apply transition matrix Monte Carlo simulations to evaluate their entropy with high precision. We consider both free- and fixed-boundary tilings. Our results suggest that the ratio of free- and fixed-boundary entropies is σ_{free}/σ_{fixed} = 3/2, and can be interpreted as the ratio of the volumes of two simple, nested, polyhedra. This finding supports a conjecture by Propp concerning the “arctic octahedron phenomenon” in three-dimensional random tilings.

Key-words: Random tilings; Integer partitions; Configurational entropy; Boundary effects; Transition matrix Monte Carlo algorithms.

I. INTRODUCTION

Since the discovery of quasicrystals in 1984^{1}, quasiperiodic tilings (such as Penrose tilings) and random rhombus tilings^{2,3} have been extensively studied as paradigms of quasicrystal structure. Quasicrystals are metallic alloys exhibiting exotic symmetries that are forbidden by usual crystallographic rules: octagonal, decagonal, dodecagonal or icosahedral symmetries. Nonetheless, the existence of sharp Bragg peaks in their diffraction patterns demonstrates a long-range translational and orientational order. As compared to perfect quasiperiodic tilings, specific degrees of freedom called “phason flips” are active in random tilings. Despite their random character, the latter still display the required long-range quasiperiodic structure.

When tiles are appropriately decorated with atoms, random tilings become excellent candidates for modeling real quasicrystalline materials^{4}. Therefore the statistical mechanics of random tilings is of fundamental interest for quasicrystal science. In particular, the configurational entropy of their phason fluctuations contributes to the free energy of quasicrystal, and might be a key ingredient in order to understand quasicrystal stability^{5}.

In parallel, random tilings have become an active field of research in discrete mathematics and computer science (see ^{6,7} for instance), and many challenging questions remain open for investigation in this field.

The relation between random tilings and integer partitions provides an important tool for the calculation of random tilings entropy^{8–12}. Integer partitions are arrays of integers, together with suitable inequalities between these integers. One-to-one correspondences can be established between integer partitions and tilings of rhombi filling specified polyhedra. However, it seems that such strictly controlled “fixed” boundary conditions inflict a non-trivial macroscopic effect on random tilings^{8}, even in the thermodynamic limit, lowering the entropy per tile below the entropy with free or periodic boundary conditions. This effect a priori makes difficult a calculation of free-boundary entropies via the partition method.

This boundary sensitivity is well described, for the simple case of hexagonal tilings^{6,13}, in terms of a spectacular effect known as the “arctic circle phenomenon”: the constraint imposed by the boundary effectively freezes macroscopic regions near the boundary, where the tiling is periodic and has a vanishing entropy density. Outside these “frozen” regions the entropy density is finite and we call the tiling “unfrozen”. The boundary of the unfrozen region appears to be a perfect circle inscribed in the hexagonal boundary. The entropy density varies smoothly within the unfrozen region, reaching a maximum equal to the free boundary entropy density at the center.

This quantitative result has never been generalized to higher dimension or codimension tilings, because its generalization requires the knowledge of the free boundary entropy if one wants to use the same proof as in the hexagonal case. However, it has recently been conjectured by James Propp (unpublished), thanks to numerical simulations, that in dimensions higher than 2, the corresponding arctic region should be a polytope itself, with flat boundaries. This conjecture renews the interest for the partition method since the relation between both entropies becomes amazingly simple in this case.

On the other hand, except an early Ansatz^{9} and some exact numerical results for small tilings^{11} (up to about 300 tiles), almost nothing is known about the entropy of codimension-one tilings of dimension larger than 2. Note also that some numerical Monte Carlo simulations provided an estimate of the entropy of (codimension 3) 3-dimensional random tilings with icosahedral symmetry^{14}.

The present paper is devoted to a numerical investigation of codimension-one three-dimensional tilings. Thanks to a powerful transition matrix Monte Carlo algorithm, we achieve precise estimates of both fixed- and free-boundary entropies. The latter is calculated via a modified partition method, which produces tilings with fixed boundaries that do not impose any strain to the tilings, thus generalizing a former two-dimensional approach^{13}. Comparing both entropies, we support Propp’s conjecture with good confidence.

The paper is organized as follows: section II reviews the relation between random tilings and integer partitions, introduces the different boundary conditions considered in this paper, and describes the arctic region phenomenon conjecture. Section III describes our Monte Carlo method and our numerical results. Discussions and conclusions are displayed in the last section.

II. PARTITIONS, TILINGS AND BOUNDARY CONDITIONS

A. Generalities

In this paper we consider three-dimensional tilings of rhombohedra which tile a region of Euclidean space without gaps or overlaps. A standard method^{2,15} to generate tilings of rhombohedra (or of rhombi in two dimensions) consists of selecting sites and tiles in a D-dimensional cubic lattice, then projecting them into a d-dimensional subspace with D > d. The difference D − d is called the codimension of the tilings. The class of symmetry of a tiling is determined by both its dimension and its codimension. We consider in this paper codimension-one, three-dimensional random tilings (i.e. D = 4 and d = 3). These tilings are composed of four different rhombohedra. All rhombohedra have identical shapes but they occur in four possible different orientations. Each rhombohedron is the projections of one of the four different three-dimensional faces of a four-dimensional hypercube. The interested reader can refer to [^{11}] for a review on codimension-one tilings.

The correspondence between codimension-one tilings and solid partitions, a reformulation of the cut-and-project method, is analyzed in detail in reference [^{11}], and generalized to higher codimensions in reference [^{12}].

We now define three-dimensional solid partitions. Consider a three-dimensional array of sides k_{1} × k_{2} × k_{3}. Fix an integer p > 0, called the height of the partition problem. Put non-negative integers in the array, smaller than p, with the constraint that these integers decrease in each of the three directions of space. More precisely, if i_{1}, i_{2} and i_{3} are indices attached to the boxes of the array (1 ≤ i_{α} ≤ k_{α}), we denote by n_{i_{1},i_{2},i_{3}} the integral variables attached to these boxes (the parts). Our partition array contains N_{p} = k_{1}k_{2}k_{3} parts. The partition constraint is

0 ≤ n_{i_{1},i_{2},i_{3}} ≤ p  (1)

and

n_{i_{1},i_{2},i_{3}} ≥ n_{j_{1},j_{2},j_{3}}  (2)

whenever i_{1} ≤ j_{1}, i_{2} ≤ j_{2} and i_{3} ≤ j_{3}.

FIG. 1. One-to-one correspondence between plane partitions (left), hexagonal tilings (center), and de Bruijn lines (right) in two dimensions^{9,11}; De Bruijn lines of two families (among three) have been represented.

Figure 1 displays correspondence in two dimensions for easy visualization. Plane partitions are k × l arrays of non-negative integers smaller than a height p and decreasing in each row and column (left). This correspondence is constructed as follows: stack 3-dimensional cubes above the boxes of the partition so that the height of each stack (center) equals the value of the corresponding partition box (left). Then project this stacking along the (1,1,1) direction of the cubic lattice. Faces of the cubes project to rhombi. The so-obtained rhombus tiling fills a hexagon of sides k, l and p (center).

Following the same construction in three dimensions, four-dimensional hypercubes are stacked above the threedimensional partition array, with the heights of the stacks equal to the corresponding parts. Then project into three dimensions along the (1,1,1,1) direction of the hypercubic lattice. Like in the hexagonal case, the so-obtained tilings fill a polyhedron, a “rhombic dodecahedron” (RD) of integral sides k_{1}, k_{2}, k_{3} and p (see the outer frame in figure 2). We call these tilings “RD-tilings” and denote their configurational entropy per tile (^{1}) by σ_{fixed}.

An alternative description of the same tilings and partitions is in terms of de Bruijn lines^{16} and membranes. Figure 1 (right) shows the de Bruijn lines of the tiling (center). These are formed by connecting centers of parallel edges on every rhombus. There are three families of de Bruijn lines, each family with an average orientation perpendicular to the tile edges. Although de Bruijn lines within a family never cross, de Bruijn line crossings among different families occur at the center of each rhombus. For three-dimensional tilings, we have instead de Bruijn membranes.

Finally, we note that every tiling can be represented as a d-dimensional directed hypersurface embedded in a Ddimensional space. Beginning with d = 2 tilings, this hypersurface consists of the faces of the stacking of cubes visible from the (1,1,1) direction (figure 1 (center)). When projected along the same direction on the “real” space E, it becomes a plane tiling. It is “directed” because no tile overlaps occur during the projection. The same relationship holds for d = 3 tilings: the directed hypersurface consists of the three-dimensional faces of the stacking of hypercubes viewed along the (1,1,1,1) direction. Since the hypersurface is directed, it can be seen as a single-valued, continuous, piecewise linear function φ from the real space E = ℝ^{3} to ℝ. The value of ϕ is just the height along the (1,1,1,1) axis. In the thermodynamic limit, these piecewise linear functions φ can be coarse-grained to obtain smooth functions ϕ : ℝ^{3} → ℝ which only contain large-scale fluctuations of the original corrugated membranes^{3,13}.

B. Boundary conditions

Polyhedral boundary conditions such as the rhombic dodecahedron bounding RD tilings have macroscopic effects on random tilings. In the “thermodynamic limit” of large system size, the statistic ensemble is dominated by tilings which are fully random only inside a finite fraction of the full volume and are frozen in macroscopic domains. By frozen, we mean they exhibit simple periodic tilings in these domains with a vanishing contribution to the entropy. In two dimensions, this is known as the “arctic circle phenomenon”, as described in introduction^{6,13}.

Such boundary conditions are not very physical. Indeed, even if one can imagine situations where the quasicrystal is constrained by a flat interface (e.g. growth experiments on a crystalline substrate), the previous considerations rely on the assumption that, in the physical quasicrystalline material, the tiles are elementary unbreakable structures. However, the system could possibly lower its free energy by breaking some tiles (a line of tiles in 2D or a surface in 3D), in order to adapt to the constraint, and thereby free the remainder the tiling from the constraint. For sufficiently large tilings a net lowering of free energy results. For example, the energy cost of a broken line of tiles grows linearly with the length of this line, whereas the free energy difference between a constrained tiling and a free one grows like the number of tiles and therefore like length squared (in two dimensions). In the thermodynamic limit, even if it costs a great amount of energy to break a tile, the system will eventually prefer to pay this cost. Therefore, the “true” thermodynamic entropy density is the free-boundary entropy^{17}.

Consequently, it is desirable to relate fixed boundary condition entropies to the more physical free boundary ones. Fortunately, there exists an exact formal relation between these entropies^{13}. A quantitative relation has even been calculated in the hexagonal case^{6,13}, but it has not been possible (so far) to extend quantitative relationships to more general classes of tilings. However, a conjecture by James Propp has recently brought some new hope (see section II C).

To exploit the calculational advantages of a partition representation, while achieving the physical free-boundary entropy in the thermodynamic limit, we adapt the partition method so that the corresponding tilings exhibit no frozen

1Recall that the configurational entropy per tile is the logarithm of the number of configurations, divided by the number of tiles.

regions. The new boundary, even though fixed, has no macroscopic effect on tiling entropy in the thermodynamic limit. The tilings become homogeneous, displaying the free-boundary entropy density throughout.

The idea is to consider tilings (we focus on “diagonal” tilings with k_{1} = k_{2} = k_{3} = p) that fill a regular octahedron O instead of the rhombic dodecahedron RD. We call these tilings O-tilings. Such an octahedron is displayed in figure 2. It is inscribed in an RD and has puckered boundaries instead of flat ones. Despite this puckering, the boundaries are effectively flat in the thermodynamic limit. Eight vertices of the RD must be truncated to produce the O.

FIG. 2. The puckered octahedral boundary conditions (k_{1} = k_{2} = k_{3} = p). The octahedron O is inscribed in the rhombic dodecahedron RD of side p coming from the solid partition method. The volume ratio |RD|/|O| = 3/2.

We now explain why O-tilings have a free boundary entropy even though their boundary is fixed. The tiling entropy is contained in fluctuations of the associated hypersurface φ. Small-scale fluctuations are integrated in an entropy functional s[ϕ] which takes into account the total number of possible piecewise linear hypersurfaces φ that are close to the smooth one ϕ (see reference [^{13}] for discussion). Functions maximizing s[ϕ] represent the dominant macroscopic states of the system. Note that s[ϕ] can be written as a functional of the gradients of ϕ, known as the phason gradient or phason strain in the quasicrystal community. And s[ϕ] is maximum and equal to the free-boundary entropy σ_{free} when this gradient vanishes everywhere.

Fixed boundaries on tilings translate into fixed boundary conditions for the functions ϕ. Therefore s[ϕ] must be maximized on a restricted set of functions, F. For rhombic dodecahedral boundaries on RD tilings, the boundaries of functions ϕ ∈ F are non-flat polyhedra, and the phason gradient cannot vanish everywhere. Therefore their entropy density σ_{fixed} is bounded below σ_{free}. For octahedral boundaries on O tilings, the functions ϕ are also constrained by a fixed boundary condition. But in this case, the boundary is flat and strain-free. It does not impose any phason gradient on the functions ϕ. The phason gradient can vanish everywhere and the maximum of s[ϕ] equals σ_{free}.

FIG. 3. Boundary conditions for the partition problem associated with the strain-free tiling problem with boundary conditions of Fig. 2. Two views of the truncated cube (of side p) are provided. The marked vertex is the same in both views.

Figure 3 illustrates changes in the partition array boundary conditions needed to achieve O tilings instead of RD tilings. The partition array is no longer a cube, but two opposite pyramidal corners have been truncated, leaving a slab of D_{3d} symmetry. This slab contains N_{p} = p^{3} − (p − 1)p(p + 1)/3 parts.

We indicate minimal and maximal values of the parts: the three faces in the topmost figure bear minimal values of their adjacent boxes, whereas the three hidden faces bear maximal values. These values are visible on the bottom in the view in the bottommost figure. Individual parts are bounded by integers depending on the position of the part in the array. Since interior parts are bounded by their neighbors, it is sufficient to impose these bounds on surface parts, as indicated in figure 3. The lower bounds range from 0 to p − 1, with constant values on lines parallel to diagonal of cube faces. The upper bounds are read on the opposite faces, and range from 1 to p.

Using the methods discussed in section II A, such partitions generate tilings filling the octahedron O. For an original cube of side k_{1} = k_{2} = k_{3} = p, the octahedron O contains N_{t} = 4p^{3} − 4(p − 1)p(p + 1)/3 tiles. These partitions are used in section III to estimate numerically σ_{free}.

C. Arctic octahedron conjecture (Propp)

The relation between free and fixed boundary entropies in not known in general, because the generalization of the proof of the hexagonal case requires the knowledge of the free boundary entropy. In two dimensions, this relation is non-trivial because the unique function ϕ_{max} maximizing s[ϕ] has a complex expression^{6,13}, leading to the arctic circle phenomenon: the statistically dominating tilings are frozen (and periodic) outside a circle inscribed in the boundary hexagon.

In three dimensions, James Propp has recently explored numerically the typical shape of a RD-tiling and has conjectures that the two-dimensional circle then becomes not a sphere but a regular octahedron, inscribed in RD like O. More precisely, the unfrozen region is not exactly O but tends towards O at the large size limit.

This conjecture has a crucial consequence: the statistical ensemble of RD-tilings is dominated by tilings periodic outside O and random inside O, equivalently O-tilings completed by eight periodically tiled pyramids to fill RD. Taking into account the tiles in the frozen regions to calculate an entropy per tile, one finally gets

σ_{fixed} = \frac{N_{O}}{N_{RD}} σ_{free} = \frac{2}{3}σ_{free},  (3)

since the ratio of the numbers of tiles in RD and O is 3/2. The next section is dedicated to the calculation of this ratio via Monte Carlo simulations.

III. MONTE CARLO

Conventional Monte Carlo simulations using Metropolis sampling^{19} are useful for generating typical tiling configurations. They can reveal temperature driven phase transitions^{20} 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^{21} and transfer matrices^{22} yield reasonably accurate free energies and absolute entropies.

We employ a variant of the transition matrix method^{24} 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 RD-tilings, as well as strong support for Propp’s arctic circle conjecture.

A. Transition Matrix

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

E(P) = ∑_{ijk} n_{ijk},  (4)

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 (see Fig. 4) a single partition element increases or decreases by one unit resulting in an energy change ΔE = ±1. The ground state of this model is the lowest height legal partition. For the boundary conditions employed here, the ground state is unique and we denote its energy as E_{min}.

FIG. 4. Example of single vertex flip: 4 rhombohedra tiling a rhombic dodecahedral region inside of a larger tiling can be “rotated” thus leading to a new tiling. The counterpart of such a flip in the solid partition representation increases or decreases a single part n_{ijk} by one unit.

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 becomes large, 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) = 2N_{p}  (5)

holds for every partition P. This value is twice the number of parts because each part can be raised or lowered. The calculation of n_{±,0} is easy (though not fast) and exact.

In a single Monte Carlo step, the probability for upwards or downwards transitions at infinite temperature is n_{±}(P)/2N_{p}, and the probability a move will be rejected is n_{0}(P)/2N_{p}. 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)/2N_{p}  (6)

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)/2N_{p} over those partitions P occurring during a Monte Carlo simulation. The acceptance ratios obey the sum rule

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

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)/Z, where

Z = ∑_{E} W (E)  (8)

is the total number of legal partitions (the “partition function”). Multiplying the probability W (E)/Z 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)/Z = ω_{−}(E + 1)W (E + 1)/Z.  (9)

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

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

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 (9) to extract the density of states even though our simulation is performed entirely at finite temperatures.

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. We store the matrix as the integer valued sum of the integers n_{±,0}(P), and impose the normalization (7) 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, ordinary 4 byte integers can 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} = k_{1}k_{2}k_{3}p. 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 W (E) 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. For a given configuration (partition P), the time required to calculate n_{±,0}(P) is proportional to the partition size N_{p}, and hence to the time required to attempt N_{p} flips. 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_{F L} = N_{p}/5 single vertex flips, 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_{F L}/100 steps. The total number of attempted flips in an N_{MC} = 10^{6} run on a 10x10x10 partition is thus in excess of 4x201x10^{6}x10^{3}/5=1.6 × 10^{11} 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 states at the energy extremes, 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. Fig. 5 plots the number of partitions sampled as a function of energy during the longest runs (length 10^{6}) for the p = 4 partition. Note that the number of hits exceeds 10^{6} uniformly for each energy in the range 0 ≤ E ≤ 256.

FIG. 5. Simulation data for p = 4 fixed boundary tiling. Top panel: Number of configurations sampled (units of 10^{6}). Center panel: Density of states W (E) (units of 10^{14} per energy). Bottom panel: Microcanonical entropy S(E) = ln W (E).

The density of states, W (E) shown in figure 5 (center), is nearly a Gaussian. W (E) reaches a peak value of 1.7 × 10^{15} states at energy E = 128, and has a half width at half maximum of ΔE = 20. As the system size grows this width grows more slowly than the number of tiles, so the density of states asymptotically approaches a delta function. Although we sampled a total of 8 × 10^{6} configurations at E = 128, this represents a fraction of only about 4.7 × 10^{−9} of the total number that exist at this energy.

The microcanonical entropy, defined as log W (E), is plotted in the lower panel of Fig. 5. This plot reveals the expected symmetry around E_{mid}. The degree to which symmetry is broken can be used as an indicator of errors accumulated during our iterative calculation (10), because sampling errors at low and high energies need not cancel. By inspection it can be seen that this error is quite low, since the entropy returns essentially to zero at high energy. The total entropy S = ln Z, and the entropy density σ = S/N_{t}.

Table I shows the convergence of entropy data as run length grows. We show the convergence only for the largest system size, which represents our worst case. Table II summarizes our data for all system sizes studied. We report the entropy resulting from the run of length N_{MC} = 10^{6}, and use the difference between that run and the length N_{MC} = 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 of order 1/p and log(p)/p because of boundary effects, which have the same order of magnitude for the values of p considered here. Indeed, these two terms are the first correction terms in the exactly solvable one- and two-dimensional random tilings. This logarithmic form fits the data better than a simple correction of order 1/p. For free boundary tilings the data fit well to

σ_{free}(p) ≃ 0.214 − \frac{0.046}{p} − 0.052\frac{log(p)}{p},  (11)

from which we conclude that σ_{free} = 0.214(2). The value is obtained from a fit excluding the data point with p = 1, while the uncertainty estimate comes from instead excluding the data point with p = 10. For the fixed boundary data we find

σ_{fixed}(p) ≃ 0.145 + \frac{0.034}{p} − 0.0049\frac{log(p)}{p},  (12)

from which we conclude that σ_{fixed} = 0.145(3). This value is very close to a conjectured limit^{9} of 0.139, but given our small uncertainty, we believe the conjectured value is not exact. For the ratio, we fit to

\frac{σ_{free}(p)}{σ_{fixed}(p)} ≃ 1.48 − \frac{0.54}{p} − 0.37\frac{log(p)}{p},  (13)

and we conclude the ratio σ_{free}/σ_{fixed} = 1.48(3). This value equals 3/2 within our uncertainty.

It would be desirable to achieve higher precision in this ratio in the future. At present, we are limited to system sizes of p = 10 or less because for larger systems the density of states W (E) exceeds 10^{308}, the limiting floating point number on our Pentium 4 processor. Either specialized floating point arithmetic or a 64-bit processor will be required to treat systems of size p = 11 and above.

IV. DISCUSSION

Returning to Propp’s conjecture, we recall (see equation (3)) we expected a ratio σ_{free}/σ_{fixed} = 3/2. In fact, with few extra hypotheses, we now demonstrate that our numerical results provide a strong support to Propp’s conjecture. Suppose that, as in the two-dimensional case, the function ϕ_{max} maximizing the entropy functional s[ϕ] is unique. This is a usual hypothesis of the random tiling model (RTM)^{3}. If ϕ_{0} is the piecewise linear function which vanishes identically in the central octahedral region O and has maximum strain (and 0 entropy) in the eight pyramidal regions comprising RD − O, then s[ϕ_{0}] = 2/3 σ_{free}, and thus s[ϕ_{0}] = σ_{fixed}. Therefore ϕ_{0} = ϕ_{max}, by uniqueness of ϕ_{max}.

As a conclusion, if σ_{free}/σ_{fixed} = 3/2 and ϕ_{max} is unique, then Propp’s conjecture is true. Note that this conjecture can be generalized to higher dimensions, making in principle possible a calculation of the above ratio for any dimension, at least as far as codimension one problems are concerned.

In terms of de Bruijn membranes, the above result means that all de Bruijn membranes are straight (or flat), at least at large scales. Indeed if de Bruijn membranes are flat in a RD-tiling, the 4 de Bruijn families intersect in the central octahedron O, but only intersect 3 by 3 in the 8 pyramidal corners, leading to a frozen tiling in these 8 regions and a strain-free tiling in O.

By comparison, one-dimensional de Bruijn lines in two dimensional hexagonal tilings are not straight since that would lead to a hexagonal central arctic region with uniform entropy density. If σ_{free}^{3→2} denotes the free boundary entropy of hexagonal tilings, one would get a fixed boundary diagonal entropy per tile of 3 σ_{free}^{3→2}/4 = 0.242 instead of 0.261 (using similar arguments as above). Fluctuations of the de Bruijn lines raise the entropy from 0.242 to 0.261. Entropic repulsion of the lines causes them to bend.

In contrast, the repulsion between de Bruijn membranes is sufficiently weak that they are not forced away from their flat configuration. This is possibly related to the fact that the fluctuations of a free 2-dimensional directed membrane in 3 dimensions are of order \sqrt{log(L)} where L is its linear size^{3}, whereas these fluctuations are larger and of order \sqrt{L} in the corresponding 2-dimensional case^{3}. Therefore it is natural to suppose that the flatness of the arctic region will persist in dimensions higher than 3 where the fluctuations are even smaller, since they are bounded^{3}.

Our result emphasizes important dimensional dependence of the transition between the frozen and unfrozen regions. Indeed, in 2 dimensions, the transition is continuous, since the entropy density is 0 by the arctic circle and then continuously varies to reach its maximum value near the center of the hexagon, with a non-zero gradient everywhere except near the center. By contrast, the situation seems to be radically new and different in 3 dimensions, since the entropy density appears to be constant in the arctic octahedron O, with a vanishing gradient everywhere and a discontinuous transition at the boundary of O.

This result (as well as its possible generalization to higher dimensions) is a strong support in favor of the partition method, of which it was formerly believed that it could not easily provide relevant results about free boundary entropies. Indeed, provided Propp’s conjecture is correct and the boundary of the arctic region is strain-free, the ratio of both entropies is nothing but a ratio of volumes of suitable polytopes. The latter two conditions should be fulfilled as soon as the entropic repulsion between de Bruijn membranes is sufficiently weak, that is as soon as the spatial dimension is 3 or greater.

ACKNOWLEDGMENTS

We acknowledge useful conversations with James Propp and Cristopher Moore. We thank Robert Swendsen for suggesting the computational approach. This work was supported in part by NSF grants DMR-0111198 and INT-9603372.

1 D. Shechtman, et al., Metallic Phase with Long-Range Orientational Order and No Translational Symmetry, Phys. Rev. Lett. 53, 1951 (1984). 2 V. Elser, Comment on “Quasicrystals: a New Class of Ordered Structures”, Phys. Rev. Lett. 54, 1730 (1985). 3 C.L. Henley, Random Tiling Models, in Quasicrystals, the State of the Art, Ed. D.P. Di Vincenzo, P.J. Steinhart (World Scientific, 1991), 429. 4 M. Mihalkovic, et al., Total-energy-based prediction of a quasicrystal structure, Phys. Rev. B (to appear, 2002). 5 M. Widom, K.J. Strandburg and R.H. Swendsen, Quasicrystal equilibrium state, Phys. Rev. Lett. , 58, 706 (1987). 6 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). 7 M. Luby, D. Randall and A. Sinclair, Markov Chain Algorithms for Planar Lattice Structures, Proc. 36th Ann. Symp. on the Foundations of Comp. Sci. (IEEE, 1995) p. 150 8 V. Elser, Solution of the Dimer Problem on an Hexagonal Lattice with Boundary, J. Phys. A: Math. Gen. 17, 1509 (1984). 9 R. Mosseri, F. Bailly, C. Sire, Configurational Entropy in Random Tiling Models, J. Non-Cryst. Solids, 153&154, 201 (1993).

10 R. Mosseri, F. Bailly, Configurational Entropy in Octagonal Tiling Models, Int. J. Mod. Phys. B, Vol 7, 6&7, 1427 (1993). 11 N. Destainville, R. Mosseri, F. Bailly, Configurational Entropy of Codimension-One Tilings and Directed Membranes, J.

Stat. Phys. 87, Nos 3/4, 697 (1997). 12 N. Destainville, R. Mosseri, F. Bailly, Fixed-boundary octagonal random tilings: a combinatorial approach, J. Stat. Phys.

102, 147 (2001). 13 N. Destainville, Entropy and boundary conditions in random rhombus tilings, J. Phys. A: Math. Gen 31, 6123 (1998). 14 K.J. Strandburg, Entropy of a three-dimensional random-tiling quasicrystal, Phys. Rev. B 44, 4644 (1991). 15 M. Duneau, A. Katz, Quasiperiodic patterns, Phys. Rev. Lett. 54, 2688 (1985); A.P. Kalugin, A.Y. Kitaev and L.S. Levitov,

Al_{0.86}Mn_{0.14}: a Six-Dimensional Crystal, JETP Lett. 41, 145 (1985); A.P. Kalugin, A.Y. Kitaev and L.S. Levitov, 6-Dimensional Properties of Al_{0.86}Mn_{0.14}, J. Phys. Lett. France 46, L601 (1985).

16 N.G. de Bruijn, Algebraic theory of Penrose’s non-periodic tilings of the plane, Ned. Akad. Wetensch. Proc. A84, 39 (1981);

Dualization of multigrids, J. Phys. France 47, C3-9 (1986) 17 To make the “free-boundary” entropy well defined, we need to add some weak constraint to ensure the tiling remains connected

and reasonably compact. We presume that the free-boundary entropy equals that for strain-free periodic boundary conditions. 18 M. Widom, Bethe ansatz solution of the square-triangle random tiling model, Phys. Rev. Lett. 70, 2094 (1993); P. Kalugin,

Square-triangle random-tiling model in the thermodynamic limit, J. Phys. A 27, 3599 (1994); J. de Gier and B. Nienhuis, Exact solution of an octagonal random tiling model, Phys. Rev. Lett. 76, 2918 (1996); J. de Gier and B. Nienhuis, Bethe ansatz solution of a decagonal rectangle-triangle random tiling, J. Phys. A 31, 2141 (1998);

19 K.J. Strandburg, L.-H. Tang and M.V. Jaric, Phason elasticity in entropic quasicrystals, Phys. Rev. Lett. 63, 314 (1989); L.-H. Tang, Random-tiling quasicrystal in three dimensions, Phys. Rev. Lett. 64, 2390 (1990); D. Joseph and M. Baake, boundary conditions, entropy and the signature of random tilings, J. Phys. A 29, 6709 (1996) M. Oxborrow and C.L. Henley, Random square-triangle tilings: a model for twelvefold-symmetric quasicrystals, Phys. Rev. B 48, 6966 (1993).

20 H.-C. Jeong and P.J. Steinhardt, Finite-temperature elasticity phase transition in decagonal quasicrystals, Phys. Rev. B 48,

9394 (1993) 21 F. Gähler, Thermodynamics of Random Tiling Quasicrystals, Proceedings of the 5th International Conference on Quasicrys-

tals, Ed. C. Janot, R. Mosseri, p. 236 (World Scientific, 1995); 22 L.J. Shaw and C.L. Henley, A transfer-matrix Monte Carlo study of random Penrose tilings, J. Phys. A 24, 4129 (1991); 23 M. Widom, D.P. Deng and C.L. Henley, Transfer-matrix analysis of a two-dimensional quasicrystal, Phys. Rev. Lett. 63, 310

(1989); W. Li, H. Park and M. Widom, Phase diagram of a random tiling quasicrystal, J. Stat. Phys. 66, 1 (1992); M.E.J. Newman and C.L. Henley, Phason elasticity of a three-dimensional quasicrystal: a transfer matrix study, Phys. Rev. B 52, 6386 (1995)

24 J.S. Wang, T.K. Tay and R.H. Swendsen, Transition matrix Monte Carlo reweighting and dynamics, Phys. Rev. Lett. 82, 476 (1999); R.H. Swendsen, B. Diggs, J.-S. Wang, S.-T. Li and J.B. Kadane, Transition matrix Monte Carlo, Int. J. Mod. Phys. C10, 1563 (1999);

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

N_{MC}  10^{3}  10^{4}  10^{5}  10^{6}
σ_{fixed}  0.147827  0.147385  0.147400  0.147349
R_{fxixed}  3.474  0.180  0.429  0.002
σ_{free}  0.197516  0.197300  0.197234  0.197273
R_{free}  1.214  0.143  -0.202  -0.036

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

p  σ_{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.158
4  0.18455(6)  0.1516677  1.213
5  0.18829(3)  0.15017(2)  1.252
6  0.19108(4)  0.14918(6)  1.281
7  0.19320(4)  0.14848(1)  1.299
8  0.19486(2)  0.14780(1)  1.317
9  0.19618(1)  0.14762(2)  1.329
10  0.19727(4)  0.14735(5)  1.339
∞  0.214(2)  0.145(3)  1.48(3)