Random tilings of high symmetry: II. Boundary conditions and numerical studies

M. Widom

Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA.

N. Destainville

Laboratoire de Physique Théorique, UMR CNRS-UPS 5152, Université Paul Sabatier, 31062 Toulouse Cedex 04, France.

R. Mosseri

Groupe de Physique des Solides, CNRS et Université 6, Campus Boucicaut, 140 Rue de Lourmel, 75015 Paris, France.

F. Bailly

LPSB-CNRS, 92195 Meudon Cedex, France. (October 21, 2003)

We perform numerical studies including Monte Carlo simulations of high rotational symmetry random tilings. For computational convenience, our tilings obey fixed boundary conditions in regular polygons. We obtain statistics on path counting and vertex coordination which compare well with predictions of mean-field theory and allow estimation of the configurational entropy, which tends to the value 0.568 per vertex in the limit of continuous symmetry. Tilings with phason strain appear to share the same entropy as unstrained tilings, as predicted by mean-field theory. We consider the thermodynamic limit and argue that the limiting fixed boundary entropy equals the limiting free boundary entropy, although these differ for finite rotational symmetry.

I. INTRODUCTION

A tiling is a filling, without gaps or overlaps, of a given region of a d-dimensional Euclidean space, with tiles which differ according to their shapes, sizes, and orientations. In the present paper, the tiles are d-dimensional rhombohedra, which we will generically call “rhombi” in the following. Each tile is the projection of a d-dimensional face of a D dimensional hypercube (D > d) into d-dimensional space. The difference D − d is known as the codimension of the tiling. Details of this construction can be found in a previous paper^{1}, hereafter referred to as “paper I”. We mainly focus on the case d = 2. A random rhombus tiling can be viewed as a fluctuating membrane in this higher dimensional space^{2}, the membrane being the union of hypercube faces (see paper I^{1}, section II.1).

Rhombus tilings are dual to de Bruijn grids^{1,3}. In two dimensions, these grid lines pass through the midpoints of parallel rhombus edges. Every rhombus edge orientation defines a family of effectively parallel de Bruijn grid lines. De Bruijn grid lines within a family never cross. In contrast, lines of different families do cross, and their crossing defines a rhombus of the tiling. There are D rhombus edge orientations and hence D families of de Bruijn grid lines.

Rhombus tilings provide simple models for quasicrystals^{4,5}, metal alloys that exhibit rotational symmetry forbidden by conventional crystallography. One of the key properties of random tilings^{6} is their configurational entropy that may play a role in stabilizing the quasicrystal state. The source of entropy is localized tile rearrangements known as phason flips. Groups of three adjacent rhombi may be permuted so that their perimeter remains fixed while their shared vertex moves to a nearby point. Our present focus, as in paper I, is the limit of high rotational symmetry obtained as D → ∞ at fixed d. We examine this limit because of its intrinsic interest and because it deepens our understanding of finite D tiling models relevant for real alloy systems.

FIG. 1. Fixed boundary D → 2 tilings with D = 3, 5 and 7. The side lengths of the polygons are all p = 8.

As usual in statistical mechanics, tiling systems can have different boundary conditions, such as free, fixed or periodic. For example, the tilings in figures 1, 2 or 12, have fixed polygonal boundaries. It is believed that free and periodic boundary tilings reach equal entropies at the thermodynamic limit. Fixed boundary tilings, in contrast, exhibit entropies that are strictly smaller than the free boundary entropy. For example, the fixed hexagonal-boundary entropy of 3 → 2 tilings (see Fig. 1, left) equals 0.261 (^{7}), while the free boundary entropy equals 0.323 (^{8}). This phenomenon can be understood^{7,9–11} by inspection of Fig. 1. The local entropy density displays a gradient between crystalline regions near the boundary, where the entropy density vanishes, and the central region, where the entropy density reaches a maximum value equal to the free boundary quasicrystalline value. Only at the very center of the tiling does the influence of the boundary disappear.

Paper I^{1} developed a mean-field theory for free boundary tilings, applicable to the limit of high rotational symmetry. Two earlier papers^{12,13} presented initial studies of this problem. The first paper^{12} proposed an upper bond on the entropy in the limit of large D and discussed problems associated with the thermodynamic limit of tiling models. The second one^{13} presented a preliminary mean-field approach of the entropy calculation. In the present paper we focus on fixed boundary tilings because they are easy to simulate numerically. We present results of our simulations, including an accurate estimate of the limiting entropy density. We investigate the role of fixed boundary conditions and argue for a boundary-condition-independent thermodynamic limit in the limit D → ∞. Similarly, we show that phason strain does not influence the entropy for large D.

The organization of this paper is as follows. We start in section II with a discussion of fixed boundary tilings, and describe their relationship to interesting problems of pure mathematics. Next, we review the problem of the thermodynamic limit in section III where we argue that free and fixed boundary tilings attain the same thermodynamic limit as D → ∞. Then, section IV describes our Monte Carlo simulations. In that section we explore the entropy, path counting statistics, vertex coordination statistics and the role of phason strain. The paper is written so that the reader does not need the notations of paper I, except in the appendices.

II. FIXED BOUNDARY TILINGS

Fixed boundary tilings such as those of figure 1 possess the most natural boundary conditions for tilings coded by combinatorial objects such as generalized partitions or sorting algorithms. Even though these boundary conditions are unnatural for bulk quasicrystals, this inconvenience is counter-balanced by the powerful tools provided by such codings in terms of understanding of tiling set structures and enumeration of tilings. Furthermore, we will demonstrate in this paper that fixed and free boundary tilings become equivalent in the large codimension limit, which will justify a posteriori the use of fixed boundary conditions.

As already described in the introduction, such fixed boundary tilings lack a proper thermodynamic limit for finite D because the boundary has a strong macroscopic effect on the whole tiling^{7,9–11,14}. This leads to a spectacular effect known as the “arctic circle phenomenon” in hexagonal (D = 3) tilings, where the tiling is periodic (and “frozen”) outside a perfect circle at the large size limit and random inside this circle. More generally and for larger D, fixed boundary tilings have an effective codimension that is smaller near the boundary than in the bulk. Very near the boundary the effective codimension vanishes. The tiling becomes a crystalline domain comprising a single tile type. The local entropy density vanishes because no phason flips are possible. In the membrane picture, this region corresponds to a flat area with large tilt. Further from the boundary the effective codimension grows, reaching a maximum at the center of the tiling. The variable codimension results in an entropy density gradient, growing from zero at the boundary to the maximal, free boundary, entropy density at the tiling center. Consequently, fixed boundary tilings have a smaller total entropy per tile than free boundary tilings.

In contrast, we will show in section III that when D becomes large (figure 2), this heterogeneity diminishes and a thermodynamic limit is restored. Since fixed boundary tilings are easier to specify and manipulate for arbitrary D, they are the appropriate tool to tackle large D entropies or other statistical properties. In the following, the fixed boundary entropies per tile will be denoted by σ̄_{D} and their limiting value by σ̄_{∞}. Corresponding values without the over-bar refer to free boundary tilings.

In terms of de Bruijn dualization, any two lines of two different families intersect inside a fixed boundary tiling^{15} (the grid is said to be “complete”). Duality defines associates tile with the intersection of two lines. Therefore if there are k_{i} lines in each family F_{i}, the number of tiles is

N_{T} = ∑_{1≤i<j≤D} k_{i}k_{j}.  (1)

If all k_{i} are set to 1, N_{T} = D(D − 1)/2.

FIG. 2. A 40 → 2 fixed boundary tiling with one de Bruijn line in each family.

The relationship between tilings and partitions has been widely explored in references^{7,15–17}. The idea is to code a random rhombus tiling by an array of integers satisfying certain ordering constraints. Figure 3 illustrates this point in the simple hexagonal case^{7}: there is a one-to-one correspondence between hexagonal tilings filling a centrally symmetric hexagon of sides k, l and p on the one hand and sets of integers arranged in a rectangular array k × l, decreasing in each row and each column, on the other hand. These integers are bounded between zero and the side length p. They decrease monotonically. Such a set of integers is called a plane partition.

FIG. 3. Three-dimensional representation of an hexagonal tiling filling a centrally symmetric hexagon (left). This object can be seen as a staking of cubes. The height of the different stacks can be arranged in an array (right). They are decreasing in rows and columns. This array together with these order relations constitutes a plane partition problem.

This point of view can be extended to any D → d problem^{15–17}: any D → d tiling can be coded with a generalized partition on an array related to a D − 1 → d tiling. This correspondence is also one-to-one. However, tilings coded by partitions have fixed polyhedral boundary conditions. In two dimensions, the boundaries are centrally symmetric 2D-gons (see figure 1).

We now present an analogy between fixed boundary tilings and algorithms for sorting lists^{15}. Although Computer Science is not our motivation, we are interested in the same enumeration problem^{18,19}. More generally, the rich topological and combinatorial properties of random tilings make them an active field of research in pure mathematics, combinatorics and computer science^{10,20–22}. We shall use results from these fields throughout this paper.

In the sorting language, a comparator [i; j] acts on a list (x_{1}, x_{2}, …, x_{D}) of numbers as follows: x_{i} and x_{j} are respectively replaced by min(x_{i}, x_{j}) and max(x_{i}, x_{j}). Following Knuth^{18}, we call a complete sorting algorithm a sequence of such comparators which sorts in the increasing order any list of real numbers (x_{1}, x_{2}, …, x_{D}). This sorting algorithm will be called primitive if each comparator can be written [i; i + 1]. We also suppose that this algorithm is not redundant, that is to say it does not contain any comparator [i; j] that could be suppressed because previous comparators already insure that x_{i} ≤ x_{j}. Knuth shows that a sequence of comparators is a sorting algorithm if it correctly sorts the completely reversed list (D, D − 1, …, 1). This means that a complete primitive sorting algorithm is a sequence of comparators [i; i + 1] that transforms the list (D, D − 1, …, 1) into the list (1, 2, …, D).

Such an algorithm has a diagrammatical representation in which the D variables x_{i} are represented by D horizontal lines. Each comparator [i; i + 1] is represented by a crossing between lines i and i + 1. Figure 4 illustrates this construction. A continuous line follows a number during the sorting process. Since every number must be compared to every other number, and since there is not any redundancy, every line crosses every other line, only once. There are \binom{D}{2} crossings.

FIG. 4. A diagram associated with a sorting algorithm acting on five element lists. A line follows a number during the sorting process. Each pair of lines cross, only once. Each crossing represents a comparator [i; i + 1].

We now establish the link between sorting algorithms and the de Bruijn representations^{1,3} of D → 2 tilings. Each line of the diagram 4 represents a de Bruijn line and crosses every other line exactly once. However, different sorting algorithms sometimes represent the same de Bruijn grid, since only the crossing topology is meaningful. For example, in figure 4, the fourth and the fifth comparator (i.e. [4;5] and [2;3]) are applied in this order (these comparator are circled in the figure). If they were applied in the reverse order, the algorithm would be different whereas the de Bruijn grid would be the same.

Therefore we define equivalence classes of sorting algorithms^{18,21}. We say that two successive comparators [i; i + 1] and [j; j + 1] commute if |i − j| > 1. Two algorithms are equivalent if they differ by a finite number of comparator commutations. Equivalence classes of D-element sorting algorithms are in one-to-one correspondence with D-family grids with one line per family, and therefore with tilings inscribed in polygons of side 1. Following Knuth, we denote by B_{D} this number of equivalence classes, whereas the number of algorithms is denoted by A_{D} (A_{D} ≥ B_{D}). Since each pattern of crossings (equivalence class of algorithms) defines a tiling, the tiling entropy density per tile is

σ̄_{D} = \frac{1}{N_{T}} log B_{D}.  (2)

FIG. 5. Entropy per tile σ̄_{D} of fixed boundary tilings (p = 1). We indicate the limiting values 0.568 (full line) and the close upper bound log 2 ≃ 0.693 (dashed line).

Sorting algorithms are easy to manipulate and provide efficient enumeration numerical tools. This analogy provides us with the number B_{D} of D → 2 tilings for small values of D ≤ 10 from the work of Knuth^{18} as extended by us^{12} (see Table I). Unfortunately it is not possible to reach the next value (D = 11) for the foreseeable future using our current algorithms. Figure 5 suggests that the entropy σ̄_{∞}(p = 1) has a well defined limit when D tends to infinity. Indeed, Knuth^{18} gives lower and upper bounds on B_{D} from which we get

\frac{1}{3} log 2 ≤ lim_{D→∞} σ̄_{D} ≤ 2 log 2.  (3)

Moreover, Björner (see reference^{19}, p. 270) derives a better upper bound:

lim_{D→∞} σ̄_{D} ≤ 1.44 log 2.  (4)

In the following section, we argue that large D, fixed boundary tilings, have the same entropy per tile as free boundary ones. Since we demonstrated in paper I^{1} that for free boundary entropies lim_{D→∞} σ_{D} ≤ log 2, we finally get:

lim_{D→∞} σ̄_{D} ≤ log 2.  (5)

This point of view generalizes to systems with more than one line per de Bruijn family^{15}, leading to the definition of partial sorting algorithms for merging pre-ordered lists of numbers. Suppose that we have D families of k_{i} numbers each (i = 1, …, D), and that in each family the numbers are presorted in increasing order. Algorithms which order the union of these sets of numbers are called partial sorting algorithms. The ideas are essentially the same as in the previous case, except that, since the numbers of a given family are already ordered, the corresponding lines do not cross within the family. The corresponding diagram is similar to a de Bruijn grid with D families of lines, k_{i} lines in each family. The tilings are equivalence classes of such algorithms. They are inscribed in polygons of sides k_{1}, …, k_{D}. In reference (^{15}), we derive an analytic expression for A_{D}(k_{1}, …, k_{D}). However, it proves impossible to calculate B_{D} analytically.

The diagonal case, with k_{1} = k_{2} = ⋯ = k_{D} = p will be of special interest below. Let B_{D}(p) denote the number of such tilings. Using the sorting algorithm analogy, we computed the numbers of tilings inscribed in polygons of side 2 (B_{D}(2)) for small values of D (see Table III).

III. THERMODYNAMIC LIMIT

In this section, we argue that the fixed boundary entropy in the large D limit equals the free boundary one and is independent of the phason strain. The (technical) proofs are provided in appendices A and B. Free boundary and fixed boundary entropies can be compared when both are known. For hexagonal tilings with D = 3, the free and fixed entropies are respectively 0.323^{8} and 0.261^{7}. For octagonal tilings with D = 4 the values are 0.434^{23} and 0.36(1)^{15}. Therefore the relative difference (σ − σ̄)/σ decreases between D = 3 and D = 4, consistent with a vanishing difference at large D. Data on p = 1 and p = 2 tilings (Tables I and III) suggest that when D becomes large, the entropy becomes independent of p.

Furthermore, a naïve observation of figures 1 leads to the following qualitative remarks: when D becomes large, the faceted regions of fixed boundary tilings, which are caused by the strong influence of the boundary^{17}, seem to be repelled towards the tiling boundaries. In these regions, few families of de Bruijn lines cross, so only a small number of tile types are present. The central region of the tiling, where a large number of de Bruijn line families cross, becomes dominant as D increases. The entropy density in the central region approaches the large D free boundary entropy σ_{∞}.

More precisely, let (x_{i}) be any set of numbers 0 < α ≤ x_{i} ≤ β, where α and β are fixed positive real numbers. Consider fixed boundary tilings of a 2D-gon with edge lengths k_{i} = x_{i}p as p → ∞. We impose such bounds on the side lengths in order to be sure that our tilings are effectively large codimension ones (see also paper I^{1} (section II.3)). The tile fractions are

n̄_{ij} = \frac{k_{i}k_{j}}{N_{T}} = \frac{k_{i}k_{j}}{k_{1}k_{2} + … + k_{D−1}k_{D}} ,  (6)

owing to eq. (1). In an unstrained tiling, where all k_{i} are equal,

n̄_{ij}^{∗} = \frac{2}{D(D − 1)}  (7)

whereas the condition 0 < α ≤ x_{i} ≤ β imposes

n̄_{ij} ≥ \frac{α^{2}}{β^{2}} n̄_{ij}^{∗} = \frac{α^{2}}{β^{2}} \frac{2}{D(D − 1)}.  (8)

This lower bound ensures that no tile fraction is vanishingly small and therefore that the tiling is really a large codimension one rather than a small codimension one with few defect lines added in.

Under this condition, it is proven in appendix A that the sequence σ̄_{D} reaches a thermodynamic limit σ̄_{∞} that does not depend on the set (x_{i}). Furthermore, these tilings are locally equivalent to free boundary high symmetry tilings. The local entropy density thus approaches σ_{∞}, and lim_{D→∞} σ̄_{D} = σ_{∞}. Moreover, this result remains valid if D becomes large while p is held fixed (even in the extreme case k_{i} = 1). In particular, the limits D → ∞ and p → ∞ commute with each other. The proofs of all these results are given in appendix A.

In conclusion, the local structure, and therefore the entropy per tile, of large D tilings is independent of their size, shape, tile fractions and boundary conditions. In all these cases, the local structure of these tilings is similar to free boundary large D tilings, where one encounters only one de Bruijn line per family in large tiling patches^{1}. In this sense, it is a true thermodynamic limit.

Henceforth, we concentrate on the study of tilings inscribed in 2D-gons of sides k_{i} = p = 1 as representative of the whole class of high symmetry tilings. In particular, we focus on p = 1 tilings to extract (numerically) the thermodynamic limit of the entropy density.

IV. SIMULATIONS

Monte Carlo numerical simulations^{24} are widespread in the fields of random tilings and quasicrystals where there exist many difficult unsolved theoretical questions^{25}. To check our claims of a thermodynamic limit, to obtain B_{D} for D > 10, and to get more precise numerical values than those obtained via mean-field arguments^{1}, we perform Monte Carlo calculations on large D tilings. The D → 2 configuration space is sampled for large D systems via single-vertex flip dynamics. This method is validated by the connectivity of the space of configurations in two dimensions^{20,15}.

These simulations utilize fixed boundary tilings for the sake of technical convenience. It is easier to code a fixed boundary tiling in the memory of a computer, and the entropy is easier to estimate in this case. We have established that the central region of a large D tiling is close to a free boundary tiling. All the statistics related to free boundary tilings will be collected in such central regions. Note that, as far as unstrained tilings are concerned, we restrict our study to polygonal boundaries of side p = 1 since we establish that such tilings behave like large p ones when D goes to infinity.

1. Monte Carlo algorithm

The algorithm is implemented as follows: at each Monte Carlo step, select a vertex at random, with uniform probability; if this vertex is flippable, then flip it; repeat sufficiently many times. The method is validated by the ergodicity of the space of configurations via flip dynamics: any configuration is reachable from any initial configuration.

A key point is to check that this algorithm samples configurations with uniform probability, since we are only interested in configurational entropy in which all configurations play the same role. Thus we need to establish that this algorithm defines a Markovian process which satisfies the detailed balance condition^{24}:

w(C_{1} → C_{2}) = w(C_{2} → C_{1}),  (9)

where w(C_{1} → C_{2}) denotes the transition probability from configuration C_{1} to configuration C_{2}. When two configurations differ by a single flip, there is only one way to go from one of them to the other one by a single flip. Therefore w = 1/N_{V} where N_{V} is the number of vertices. Since N_{V} is independent of the tiling, relation (9) is established.

Note that the flip acceptance rate equals the fraction (about 31% as reported in section IV 3) of vertices with threefold coordination. Vertices with more than threefold coordination cannot be flipped by an “elementary” phason flip. An alternative Monte Carlo dynamics, in which a threefold vertex is flipped in every step, does not obey detailed balance because the number of flippable vertices is not conserved, and hence eq. (9) is violated.

Another point to consider in Monte Carlo simulations is the correlation time τ (in units of Monte Carlo steps per vertex) which measures how many steps are necessary between samples to avoid excessive sample-to-sample correlations. Even though some recent work brought new results for small D (see^{26} and references therein), this question has no rigorous definitive answer for large D. However, our numerical calculations are all in agreement with the conjecture τ ≤ N_{V} /2.

2. Path counting and entropy estimate

To estimate the entropy σ_{D} of fixed boundary tilings, we shall use a path-counting algorithm based upon de Bruijn line enumeration^{12} directly derived from the general method described in section 2.2 of paper I^{1}. Figure 6 illustrates the method. At each step, we build a D + 1 → 2 tiling from a D → 2 one. Starting from any D → 2 tiling, we choose a path of length D along tile edges, going from bottom to top. We cut the tiling along this path, separating the two parts by length 1, then draw new bonds connecting previously identical vertices. Finally, we adjust all edge orientations to match the set of symmetry D + 1.

FIG. 6. Iterative construction of fixed boundary tilings (D → 2, p = 1).

Let P̄_{D} denote the mean number of bottom-to-top paths on the D → 2 tilings:

P̄_{D} = \frac{1}{B_{D}} ∑_{τ=1}^{B_{D}} P_{D}(τ),  (10)

where P_{D}(τ) is the exact number of paths on the tiling τ, and B_{D} is the number of such tilings. Then

B_{D+1} = P̄_{D}B_{D}.  (11)

Iterating relation (11), and taking the logarithm yields

log B_{D} = ∑_{k=1}^{D−1} log P̄_{k}.  (12)

For large D, the ratio \frac{1}{D} log P̄_{D} approaches a finite limiting value. Taking the limit of (12) as D → ∞, and noting eqs. (1) and (2), we find

lim_{D→∞} σ̄_{D} = lim_{D→∞} \frac{1}{D} log P̄_{D}.  (13)

Hence, by accumulating statistics on the number of paths allowed on tilings, we may evaluate the entropy.

We use Monte Carlo sampling to generate an ensemble of D → 2 tilings. On each tiling τ we can quickly evaluate P_{D}(τ) using a generalization of the Pascal Triangle construction: starting at the bottom of the tiling, assign each vertex an integer value equal to the number of paths that reach it from the bottom. The value at any vertex is iteratively the sum of the values at each prior vertex to which it is connected (see Fig. 7). At the end of the process, the value of the top vertex is P_{D}(τ).

FIG. 7. Path count construction, as described in the text.

We tested this algorithm on small D systems, where numerical values are in good agreement with the exactly known mean number of paths P̄_{D} = B_{D+1}/B_{D} (table I). For example, table II shows the convergence towards the exact value for D = 9 tilings. For large D systems, simulations are our only means of obtaining information about path count statistics, which are displayed in figure 8. In this figure, error bars are calculated as follows: the standard deviation Δ is measured numerically. The error is estimated by Δ/\sqrt{I}, where I is the number of independent measures in the simulation (that is^{24} its total length divided by twice the correlation time τ). Note that the actual errors calculated in table II are compatible with the previous error bar.

This graph presents a maximum near D = 20 which is related to finite D and finite size corrections. Indeed, it is demonstrated in paper I^{1} that for free boundary systems, finite D corrections to log P̄_{D}/D are of order 1/D^{2}. In addition to these finite D effects, we examine small systems with one de Bruijn line per family and therefore subject to strong boundary effects. For large D systems, this boundary effect should scale like 1/D. The combination of these corrections leads to the large D behavior

\frac{log P̄_{D}}{D} = σ̄_{∞} + \frac{Cst.}{D} − \frac{Cst.′}{D^{2}} + O(\frac{1}{D^{3}} ).  (14)

FIG. 8. Value of log P̄_{D}/D as a function of D up to D = 120; inset: large D behavior as a function of 1/D and quadratic fit. The numerical values as well as the fitted ones present a maximum near D = 20.

The numerical results in figure 8 reproduce this behavior; Fitting these numerical results with the above behavior, we observe that the fit is excellent. All numerical data beyond D = 15 coincide with the fitted ones up to error bars. In particular, the maximum observed near D = 20 is quite well reproduced by the fitted function.

Finally, the limiting value of log P̄_{D}/D can be estimated as σ̄_{∞} = 0.5676 ± 0.0001. The uncertainty estimate on this last limit comes from excluding the data with D = 120.

3. Vertex statistics

Other pertinent information available from Monte Carlo simulations concerns vertex statistics. They are indicators of the tiling’s local “microscopic” structure. They also are a fundamental ingredient of the mean-field approach of paper I^{1}. We focus on quantities which can be computed in the mean-field theory, especially coordination numbers and related quantities. We gather statistics both in the whole tiling (excluding vertices actually on the boundaries) and in a central region containing about 20% of the vertices. Presumably, the central region is only weakly influenced by the fixed boundaries, and we take that data as representative of free boundary tilings.

FIG. 9. Fractions of Z = 3 and Z = 5 vertices in D → 2 vertices as a function of 1/D, and quadratic fits. Diamonds concern the central region (20% vertices), whereas circles concern the whole tiling. The extrapolated large D values coincide.

We have run Monte Carlo simulations up to D = 100. After plotting the data as functions of 1/D, we extrapolate the limiting values via quadratic fits. Figure 9 illustrates two examples of coordination number statistics. We find the fractions of vertices (in the whole tiling and in the central region) whose coordination numbers Z = 3 and the fractions whose coordination numbers Z = 5. For finite D we see the central region has a relatively small fraction of tiles with Z = 3 in comparison with the whole tiling, and a relatively large fraction with Z = 5. The boundary regions are thus more likely to have Z = 3 and less likely to have Z = 5 than a free boundary tiling. However, the extrapolated values for D → ∞ agree to within the accuracy of the extrapolation. The same conclusion holds for Z = 4 (not represented here). For larger values of Z, the fractions of vertices are too small and we did not obtain relevant measures.

Hence, in the limit D → ∞ the whole tiling exhibits coordination statistics similar to a free boundary tiling. This point supports the existence of a thermodynamic limit - the local structure becomes uniform throughout the tiling, independent of proximity to the boundaries and the local strain they create.

Other quantities we examined relate directly to mean-field theory and were reported in paper I^{1}. Our mean-field theory was based upon the number of choices to be made while inscribing directed paths on tilings. The number of choices N_{c}(v) at a vertex v is the number p of “arms” emerging from a vertex. The number q of “legs” leading into a vertex together with the number of arms obey p + q = Z. Paper I reports the simulated probability distribution for p and mean value of pq, and shows these agree well with mean-field theory.

Here, we display in table IV the probability distribution for Z, and compare it with the predictions of mean-field theory. The agreement is also satisfactorily good.

4. Phason strained tilings

The data just reported in section IV 3 suggests that the tiling structure does not depend on the phason strain, since statistics are identical in the central region and in the boundary vicinity. We want to test this point directly by performing simulations on phason strained fixed boundary tilings.

We consider three different forms of phason strain in our numerical simulations: low frequency strain in which D/2 consecutive line families occur with k_{i} = 1 (i = 1, …, D/2), the remaining D/2 consecutive line families occur with k_{i} = 2 (i = D/2 + 1, …, D); the same low frequency with a larger amplitude so that k_{i} = 3 (i = D/2 + 1, …, D) replaces k_{i} = 2; high frequency strain with k_{i} = 1 (i odd) and k_{i} = 2 (i even). Figure 10 illustrates an equilibrated tiling with applied phason strain. The tiling has D = 30 and the strain is low frequency and large amplitude.

Computer simulations of path count statistics for D = 10, ..., 50 reveal no significant dependence of path count on phason strain amplitude or frequency. The small variations seen (see table V) cannot be separated from finite size effects. Thus entropy does not depend on phason strain, supporting a posteriori arguments developed for free boundary tilings^{1}. Recall that strain-independence of entropy in fixed boundary tilings relies upon strain-independence in free boundary tilings (see section III).

FIG. 10. Phason strained large D tiling. The added paths still go from bottom to top.

The same conclusion holds concerning vertex statistics on coordination number and choice distributions. No significant differences are found in comparison to unstrained fixed boundary tilings or free boundary tilings (see figure 11). These data support the existence of a universal local structure, independent of strain and boundary conditions in the limit D → ∞.

FIG. 11. Fractions of Z = 3 and Z = 5 vertices in strained tilings (p = 1, 1, …, 3, 3, …) (triangles), together with previous data on fixed boundary unstrained tilings (circles) and free boundary ones (diamonds). The limiting values coincide with the previous one, up to statistical errors. Lines are quadratic fits.

V. CONCLUSION

This paper tackles random tilings of high symmetry with fixed boundary conditions. In random tiling theory, boundary conditions are crucial because finite codimension fixed boundary tilings have a lower entropy than free boundary ones. However, we argue that boundaries become irrelevant when the codimension becomes large. We also demonstrate that p = 1 tilings with one de Bruijn line per family have the same entropy per tile as tilings with p large and even infinite, for sufficiently large D.

As a consequence, the numerical study of the entropy of large codimension random tilings can be concentrated on fixed boundary tilings filling a 2D-gon of side lengths set to 1. We perform exact enumeration for D ≤ 10, thanks to an analogy between fixed boundary tilings and some class of sorting algorithms. We use Monte Carlo simulations for larger tilings, up to D = 120. In both cases, the fact that tilings have a fixed polygonal boundary greatly simplifies their encoding and their manipulation in the memory of a computer. Our Monte Carlo data analysis is based on the same iterative process (i.e. concentrating on values of P̄_{D}) as our mean-field theory in paper I^{1}. We obtain a very accurate estimate of the entropy per tile, σ̄_{∞} = 0.5676 ± 0.0001, compared to the approximate mean-field value σ_{∞}^{MF} ≃ 0.598.

Insensitivity of the entropy per tile and vertex statistics to boundary conditions at large D restores a thermodynamic limit which does not exist at finite D. However, caution is required. Even if “topological” quantities such as entropy and vertex statistics become homogeneous in the tiling, its “metric” properties remain heterogeneous. For example, it seems that tiles become thinner near the boundary and that the mean area per tile decreases from the center of the tiling to its periphery (see figure 2). The entropy per unit area is not uniform in the tiling. This point will have to be clarified further.

Many of our methods and results can be generalized to higher dimensional tilings. However, the generalization may not be simple to implement in practice. Indeed the de Bruijn directed paths become directed surfaces built of faces of rhombohedra. The description of these surfaces in terms of successive choices must be generalized, which highly complicates their enumeration.

The restoration of the thermodynamic limit in the high symmetry limit should remain valid in three- or higherdimensions, because the vision of large regions with at most one de Bruijn surface per family still holds. The explicit proof should be similar to the appendix A, with only minor changes to account for the possible spontaneous decomposition into regions of lower effective codimension (as in the “arctic octahedron”, see Ref.^{14}). Presumably, the effective codimension still goes to infinity for almost all tiles.

One of the main motivations for the study of high codimension random tilings was the hope that mean-field theory might become exact in the high codimension limit. The principal conclusion of this numerical study, together with the mean-field results of paper I, is that the mean-field theory does not become exact, at least in the simplest version which neglects vertex correlations. Rather, the high codimension limit is quite nontrivial, and an exact analytic solution is yet to be found.

ACKNOWLEDGMENTS

We thank Chris Henley and Pavel Kalugin for useful discussions. This research is supported in part by the National Science Foundation under grants DMR-0111198 and INT-9603372 and by the CNRS.

APPENDIX A: THERMODYNAMIC LIMIT OF TWO-DIMENSIONAL FIXED BOUNDARY TILINGS

In this technical appendix, we demonstrate that the fixed boundary entropy σ̄_{D} attains a finite limit σ̄_{∞} when D → ∞, and that this limit coincides with the free boundary entropy σ_{∞}. Moreover, we argue that this limit is shared with finite side polygonal boundaries as well as non-diagonal tilings.

Our demonstration relies on the variational principle introduced in references [^{10,11}] to characterize the entropy gradient between the center and the boundary in typical tilings with given fixed boundaries. Typical tilings are those which maximize an entropy functional defined as the integral over the tiling of a local entropy per tile; This local entropy is the free-boundary entropy calculated with the local fractions of tiles. The entropy per tile of fixed boundary tilings is the maximum of this functional.

We show below that, given any domain in the tiling, when D goes to infinity, the tiling in this domain is a piece of (free) D′ → 2 tiling, such that D′ goes to infinity as D does. Hence the local entropy per tile tends to the large D free boundary entropy σ_{∞} (nearly) everywhere and fixed boundary tilings have the same entropy as free boundary ones.

We first consider the case where the sides of the polygonal boundary share the same length p (diagonal tilings), and take the large p limit before taking the large D limit. Secondly, we consider the case were D becomes large at fixed p. Finally, we discuss the non-diagonal case where these side lengths might be different. Our presentation is heuristic but can be made rigorous.

1. Diagonal tilings

The demonstration uses a particular tiling denoted as T_{0}. It is the dual of a particular de Bruijn grid^{1,3}, consisting of D families of p straight lines equally spaced out and forming a regular fan^{27} (see figure 12, left). The de Bruijn families are denoted by F_{1}, F_{2}, …, F_{D}, and are labeled counterclockwise so F_{k} makes an angle kπ/D with an arbitrary reference direction. Figure 12 displays such a grid and its dual tiling. For finite D, such a tiling is known^{11} not to maximize the entropy functional among the tilings with the same boundary (see below). By construction, the tiling T_{0} is divided into domains where only a fraction of the de Bruijn families intersect. In each domain, the tiling is homogeneous, with a well defined local entropy per tile at large p. Since only D′ < D de Bruijn families intersect in such a domain, this entropy a priori differs from the free boundary D → 2 one. However, by rotational symmetry, all the domains where exactly D′ families intersect have the same local entropy. In the following, the union of such domains where only D′ families meet will be denoted by A_{D′} and called the region of “effective dimensionality” D′.

FIG. 12. The de Bruijn straight lines (left) and the corresponding rhombus tiling in the 4 → 2 case with p = 20 (for clarity’s sake, we have only drawn 10 lines per family instead of 20 on the left figure). The regions where only 3 or even 2 families intersect clearly appear on the left figure. The corresponding regions are delimited by the two internal octagons on the right figure: in the outer region A_{2}, only 2 families meet; In the intermediate one A_{3}, 3 families meet; And in the central region A_{4}, all 4 families intersect.

Define S_{0}^{D→2} as the entropy functional evaluated at T_{0}. It satisfies S_{0}^{D→2} ≤ σ̄_{D}. Then

lim_{D→∞} S_{0}^{D→2} ≤ σ̄_{∞} ≤ σ_{∞}.  (A1)

Next we intend to show that for any real number h < 1,

lim_{D→∞} S_{0}^{D→2} ≥ h σ_{∞}.  (A2)

It then follows that

lim_{D→∞} S_{0}^{D→2} ≥ σ_{∞},  (A3)

and, owing to relation (A1), that

σ_{∞} ≤ lim_{D→∞} S_{0}^{D→2} ≤ σ̄_{∞} ≤ σ_{∞},  (A4)

and therefore that the 3 involved quantities are equal.

Let us prove the above statement (A2). We use the fact (see section IV 4 and also paper I^{1}), that the entropy σ_{∞} of free boundary tilings does not depend on strain in the large D limit. We show that the local entropy of T_{0} equals σ_{∞} in (nearly) all regions A_{D′}.

First we estimate how many tiles each region contains. Since a tile is defined as the intersection of two de Bruijn lines, we will count the number of such intersections in a given region, by calculating the number of intersections per unit area in each region on the one hand and the area of the regions on the other hand. However we work not in the “tiling metric” but in the “grid metric” instead. That is to say, in the figure where the de Bruijn lines were originally drawn straight (left-hand representation in fig. 12). Figure 13 (left) displays a schematic representation of a grid were the different regions under consideration clearly appear.

The number of intersections in a region A of area α is the sum over all pairs of de Bruijn families of the number of intersections of a pair. Consider a given pair of families F_{i} and F_{j}, i < j, which make an angle θ_{ij}. The number of such intersections is the area α divided by the area α_{0} of the unit cell of the lattice defined by these (only) two families (figure 13, right). If the distance which separates two lines of a family is set to 1, then α_{0} = 1/| sin θ_{ij}|. Hence the total number of intersections is

∑_{i<j} \frac{α}{1/| sin θ_{ij}|} = α ∑_{i<j} | sin θ_{ij}|,  (A5)

where the indices i and j run over the families present in A.

FIG. 13. Left: Schematical representations of the D families of de Bruijn lines and of their intersections, to be compared to fig. (12), left. When D tends to infinity, the polygonal boundary tends toward a circle of circumference 2Dp and the polygonal central region A_{D} tends toward a circle a diameter p. Right: Two families of de Bruijn lines making an angle θ_{ij}. The number of intersections of these two families in a region equals its area α divided by the area α_{0} of the grayed unit cell.

As displayed in figure 13 (left) a region A_{D′} is a crown made of 2D kites (or triangles for the inner crown) of equal areas. A kite (or triangle) of A_{D′} will be denoted by K_{D′} and its area by κ_{D′}. Henceforth the area of A_{D′} is equal to α_{D′} = 2Dκ_{D′}. In K_{D′}, there are D′ adjacent families of lines, for example F_{1}, …, F_{D′}. The angles that they make relatively to an arbitrary reference direction are therefore θ_{l} = θ_{0} + lπ/D, l = 1, …, D′ and the number of intersections in K_{D′} is

Q_{D′} = κ_{D′} ( ∑_{1≤i<j≤D′} sin \frac{(j − i)π}{D} ) .  (A6)

When D goes to infinity, x = D′/D becomes a continuous variable which represents a “fraction” of the D families of lines. For example, the previous sum in eq. (A6) can be estimated by an integral on the variables y = j/D and z = i/D. If E(x) denotes the integral part of x, then

ϕ_{D}(x) ≡ ∑_{1≤i<j≤E(xD)} sin \frac{(j − i)π}{D} ≃ D^{2} ∫_{0}^{x} dy ∫_{0}^{y} dz sin π(y − z) = \frac{D^{2}}{π^{2}} [πx − sin πx] .  (A7)

The central region A_{D} plays a particular role among the regions A_{D′}: when D is large, the fraction of tiles in A_{D} tends to 1/2. Indeed, as illustrated in figures 12 (left) and 13 (left), the central region tends toward a circle of diameter p. Its area is πp^{2}/4, and, owing to the above result applied to the case x = 1, the number of intersections that it contains grows like

\frac{D^{2}}{π} \frac{πp^{2}}{4} = \frac{D^{2}p^{2}}{4} .  (A8)

Since the total number of tiles grows like \frac{D^{2}p^{2}}{2} , this central region contains one half of the total number of tiles.

By comparison, the situation is completely different in the other regions A_{D′}, with D′ < D. In this case, when D goes to infinity, the fraction of tiles in such a region vanishes. Let us first explicitly compute the asymptotic behavior of this fraction and then discuss how to handle this vanishing character.

FIG. 14. (left) around the central region A_{D}, there is a first crown of triangles K_{D−1} (light gray), a second crown of kites K_{D−2} (medium gray), a third crown of kites K_{D−3} (dark gray) and so forth; in this figure, D = 6. Right: the three kinds of kites which define a kite-shaped region K_{D′} (white): a big kite minus two small kites (light gray) plus a smaller kite (dark gray).

Figure 14 (left) shows the geometry of the regions under consideration when D is finite. Around the central region A_{D}, there is a first crown triangular domains, then a second crown of kite-shaped regions, and as k increases, the regions A_{D′} lie in concentric crowns around the central one. Each crown contains 2D kite-shaped (or triangle-shaped) regions. As figure 14 (right) illustrates, such a kite can be seen as a greater rectangular kite^{28} minus two equal rectangular kites. These two latter kites have a non empty intersection, which is itself an even smaller kite. If the areas of these kites are respectively denoted by δ_{D′}, δ_{D′+1} and δ_{D′+2}, then the area of K_{D′} is

κ_{D′} = δ_{D′} − 2δ_{D′+1} + δ_{D′+2}.  (A9)

Now, for each of these 4 kites, its edges are two rays of the circular central region and two tangents to this central region, perpendicular to these two latter rays. Since the central region has diameter p, if the two rays make an angle θ, the area of the kite is

δ(θ) = \frac{p^{2}}{4} tan (\frac{θ}{2}) .  (A10)

For the domain K_{D′}, θ = π(D − D′)/D for the greater kite. Henceforth,

δ_{D′} = \frac{p^{2}}{4} tan (π \frac{D − D′}{2D} ) .  (A11)

The same formula holds for δ_{D′+1} and δ_{D′+2}. Thus if we set k = D − D′,

κ_{D′} = \frac{p^{2}}{4} [tan (\frac{kπ}{2D} ) − 2 tan (\frac{(k − 1)π}{2D} ) + tan (\frac{(k − 2)π}{2D} )] .  (A12)

Fixing the ratio x = D′/D, and approximating the above as a second derivative, we finally get

κ_{D′} = κ_{E(xD)} ≃ π^{2} \frac{p^{2}}{8} \frac{1}{D^{2}} tan (π(1 − x)/2) [1 + tan^{2} (π(1 − x)/2)] .  (A13)

Thus the number of tiles in A_{E(xD)} is

N_{E(xD)} = α_{E(xD)}ϕ_{D}(x) = 2Dκ_{E(xD)}ϕ_{D}(x).  (A14)

Finally, the fraction of tiles in the region A_{E(xD)} is

n_{E(xD)} = \frac{N_{E(xD)}}{N_{T}} = \frac{N_{E(xD)}}{D^{2}p^{2}/2} ≃ \frac{1}{2D} tan (π(1 − x)/2) [1 + tan^{2} (π(1 − x)/2)] [πx − sin πx]  (A15)

when D goes to infinity, since the number of tiles in the whole tiling is N_{T} = D^{2}p^{2}/2. As expected, these fractions vanish but they can be added to get a non vanishing number of intersections. More precisely, let us fix a real number 0 ≤ g ≤ 1 and let us compute the fraction of tiles which lie in regions where at least E(gD) families of lines intersect:

\frac{1}{2} + ∑_{l=E(gD)}^{D−1} n_{l} ≃ \frac{1}{2} + \frac{1}{2} ∫_{0}^{1−g} dx tan (πx/2) [1 + tan^{2} (πx/2)] [π(1 − x) − sin πx] ≡ \frac{1}{2} + \frac{1}{2}f(1 − g).  (A16)

The above integral is a continuous function f(g) for g ∈ [0, 1]. It fulfills the required condition f(0) = 1 (which means that the whole tiling contains a fraction 1 of tiles!). Moreover, f(g) < 1 when g > 0. Thus for any h < 1, there exists a real number g > 0 such that a fraction h of the tiles lie in regions where at least E(gD) families of lines meet. Since g > 0, when D goes to infinity, the number of families of lines also goes to infinity in such regions.

Moreover it is demonstrated in the following appendix B that the local tilings in such regions are true large codimension ones. Hence everywhere in the above region, the local entropy tends to the free boundary entropy σ_{∞}. A fraction h of the tiles lie in regions where the local entropy per tile tends to σ_{∞} when D → ∞. We get^{11} the expected relation (A2) and we conclude that free and fixed boundary entropies coincide. We also remark that the tiling T_{0} which was not supposed a priori to maximize the entropy functional indeed does so at the large D limit.

These arguments can also be extrapolated to finite p tilings provided D is large. Indeed, if p is finite, a domain K_{D′} can be very small and can contain very few tiles (and even no tile at all). Nevertheless in the previous demonstration, regions A_{E(xD)} with x_{0} ≤ x ≤ x_{0} + δx can be put together into larger regions, to which all the previous arguments can be applied.

2. Non-diagonal tilings

For large p non-diagonal tilings with side lengths k_{i} = x_{i}p where α ≤ x_{i} ≤ β, we use a particular tiling T′_{0} which is a variation of T_{0}. Its overall definition is the same except that the line spacing l_{i} depends on the de Bruijn family F_{i}. We follow the main steps of the prior demonstration, in particular the calculation of the areas κ_{D′} or α_{D′}.

We choose l_{i} = p/(x_{i}p − 1) ≃ 1/x_{i} so that the width l_{i}(k_{i} − 1) = p of the family F_{i} is independent of i in the grid representation (figure 13, left). Therefore the areas α_{D′} remain unchanged as compared to the previous section.

On the other hand, the density of intersections in each domain K_{D′} will vary because of varying line spacing. In particular the area α_{0} of the unit cell now depends on i and j:

α_{0}(i, j) = \frac{l_{i}l_{j}}{| sin θ_{ij}|} = \frac{1}{x_{i}x_{j}| sin θ_{ij}|}.  (A17)

The number of intersections (A6) becomes

Q̃_{D′} = κ_{D′} ( ∑_{1≤i<j≤D′} x_{i}x_{j} sin \frac{(j − i)π}{D} ) ≤ β^{2}Q_{D′}.  (A18)

Henceforth, the new numbers of tiles in the crown A_{D′} and in the whole tiling satisfies Ñ_{D′} ≤ β^{2}N_{D′} and Ñ_{T} ≥ α^{2}N_{T}. If we set again D′ = E(xD), the new fraction of tiles in A_{D′} satisfies

ñ_{E(xD)} = \frac{Ñ_{E(xD)}}{Ñ_{T}} ≤ \frac{β^{2}}{α^{2}} n_{E(xD)}.  (A19)

As a consequence, the outer crowns with small effective dimensionality D′, which had a vanishingly small contribution to the total number of tiles in the previous section, will again be negligible in the present case.

To finish the proof, we must check that the remainder of the tiling has a local entropy equal to σ_{∞}. The proof is explicited in appendix B in the diagonal case but can easily be adapted to the present case, leading to

\frac{n′_{ij}}{(n′)_{ij}^{∗}} ≥ \frac{α^{2}}{β^{2}} \frac{6}{π^{2}} ,  (A20)

where the notations are defined in appendix B. Even though the lower bound on tile fractions now depends on α and β because we are dealing with a nondiagonal case, it remains finite.

APPENDIX B: REGULARITY OF LARGE CODIMENSION FIXED BOUNDARY TILINGS

We check here that in each region A_{D′} of the diagonal tiling T_{0} the local entropy per tile tends to the free boundary one σ_{∞} as D → ∞. As discussed in paper I^{1}, a sufficient condition is that tile fractions should be “bounded” on A_{D′}. This means that all local tile fractions n_{ij} have the same order of magnitude as the strain-free ones n_{ij}^{∗}. In other words, there exist a finite constant a > 0 such that

n_{ij} ≥ a n_{ij}^{∗}  (B1)

for all i and j. We want to calculate the tile fractions n′_{ij} in A_{D′} and compare them to the strain-free ones. Recall that in a D → 2 tiling,

n_{ij}^{∗} = Cst. sin |π \frac{j − i}{D}| ,  (B2)

where the constant comes from the normalization relation ∑ n_{ij} = 1.

As in the previous appendix, we set D′ = E(xD) where D is large and x > 0 is finite. The main difficulty in this appendix comes from the fact that in the region A_{D′} under consideration, not all tile species occur. Without loss of generality, we assume that only the tiles with indices i, j = 1, …, D′ exist. These tiles appear in the region A_{E(xD)} with the fractions (B2), but with a different normalization constant because only some of them appear. The normalization relation now reads

∑_{1≤i<j≤D′} n′_{ij} = 1.  (B3)

At the large D limit, this sum can be replaced by an integral and we get the tile fractions in A_{E(xD)}:

n′_{ij} = \frac{1}{(D′)^{2}} \frac{π^{2}x^{2}}{πx − sin(πx)} sin |xπ \frac{j − i}{D′} | .  (B4)

Now we desire to check that this tiling has bounded local fractions of tiles when it is considered as a D′ → 2 tiling. We need to compare the above n′_{ij} to the corresponding quantities in an unstrained D′ → 2 tilings, in other words to

(n′)_{ij}^{∗} = \frac{π}{(D′)^{2}} sin |π \frac{j − i}{D′} | .  (B5)

By convexity of the sin function on the interval [0, π], we have: sin |xπ(j − i)/D′| ≥ x sin |π(j − i)/D′|, and

\frac{n′_{ij}}{(n′)_{ij}^{∗}} ≥ \frac{x}{π} \frac{π^{2}x^{2}}{πx − sin(πx)} ≥ \frac{6}{π^{2}} ,  (B6)

which achieves the proof: a = 6/π^{2} in condition (B1). From the results of paper I^{1}, we conclude that in such a region A_{D′}, the local entropy per tile equals the free boundary value σ_{∞} when D → ∞.

1 N. Destainville, M. Widom, R. Mosseri, F. Bailly, “Random tilings of high symmetry: I. Mean-field theory”, to be submitted to J. Stat. Phys. (2003). 2 V. Elser, Phys. Rev. Lett. 54, 1730 (1985). 3 N.G. de Bruijn, Ned. Akad. Wetensch. Proc. A84, 39 (1981); J. Phys. France 47, C3-9 (1986). 4 D. Shechtman, et al. Phys. Rev. Lett. 53,1951 (1984). 5 D. Levine, P.J. Steinhardt, Phys. Rev. Lett. 53, 2477 (1984). 6 C.L. Henley, in Quasicrystals, the State of the Art, Ed. D.P. Di Vincenzo, P.J. Steinhart (World Scientific, 1991), 429. 7 V. Elser, J. Phys. A 17, 1509 (1984). 8 G.H. Wannier, Phys. Rev. 79, 357 (1950); Phys. Rev. B 7 5017 (E, 1973).

9 D. Grensing, G. Grensing, J. Math. Phys. 24, 620 (1983); D. Grensing, I. Carlsen, H.Chr. Zapp, Phil. Mag. A 41, 777 (1980). 10 H. Cohn, M. Larsen, J. Propp, New York J. of Math. 4, 137 (1998); H. Cohn, R. Kenyon, J. Propp, J. Amer. Math. Soc.

14, 297 (2001). 11 N. Destainville, J. Phys. A: Math. Gen 31, 6123 (1998). 12 M. Widom, N. Destainville, R. Mosseri, F. Bailly, in Proceedings of the 6th International Conference on Quasicrystals (World

Scientific, Singapore, 1997). 13 N. Destainville, M. Widom, R. Mosseri, F. Bailly, Mat. Sci. Eng. A 294-296, 409 (2000). 14 M. Widom, R. Mosseri, N. Destainville, F. Bailly, J. Stat. Phys. 109, 945 (2002). 15 N. Destainville, R. Mosseri, F. Bailly, J. Stat. Phys. 102, 147 (2001). 16 R. Mosseri, F. Bailly, Int. J. Mod. Phys. B, Vol 7, 6&7, 1427 (1993). 17 N. Destainville, R. Mosseri, F. Bailly, J. Stat. Phys. 87, 697 (1997). 18 D.M. Knuth, Axioms and Hulls, in Lect. Notes in Computer Sci. 606, 35 (1992). 19 Oriented Matroids, A. Björner, M. Las Vergnas, B. Sturmfels, N. White, G.M. Ziegler (Cambridge University Press, 1993). 20 R. Kenyon, Algorithmica 9, 382 (1993). 21 S. Elnitsky, J. Combinatorial Theory A 77, 193–221 (1997). 22 Tilings of zonotopes: Discriminental arrangements, oriented matroids, and enumeration, G.D. Bailey, Ph. D. Thesis (Univ.

of Minnesota, 1997). 23 W. Li, H. Park and M. Widom, J. Stat. Phys. 66, 1 (1992). 24 Monte Carlo methods in statistical mechanics, M.E.J. Newman, G.T. Barkema (Clarendon Press, Oxford, 1999). 25 K.J. Strandburg, L.-H. Tang and M.V. Jaric, Phys. Rev. Lett. 63, 314 (1989); L.J. Shaw, V. Elser, C.L. Henley, Phys. Rev.

B 43, 3423 (1989); L.-H. Tang, Phys. Rev. Lett. 64, 2390 (1990); M. Oxborrow and C.L. Henley, Phys. Rev. B 48, 6966 (1993); F. Gahler, Proc. ICQ5 236 (1995); D. Joseph and M. Baake, J. Phys. A 29, 6709 (1996).

26 M. Luby, D. Randall and A. Sinclair, SIAM J. of Comp. 31, 167 (2001); N. Destainville, Phys. Rev. Lett. 88, 30601 (2002). 27 Strictly speaking, should this fan be exactly regular, there would be multiple intersections, for example at its very center. To

avoid this difficulty, each family must be slightly shifted by a random distance much smaller than the interline separation. Such shifts have been performed in the figure 12, even though it might not be clear because of resolution.

28 A rectangular kite has two right angles between its unequal sides.

TABLE I. Exact (except B_{11}/B_{10}) data for fixed boundary tilings, p = 1 and D = 1, 2, …, 10.

D  1  2  3  4  5
B_{D}(p = 1)  1  1  2  8  62
σ̄_{D}(p = 1)  0  0.231  0.347  0.413
B_{D+1}/B_{D}  1  2  4  7.75  14.65
D  6  7  8  9  10
B_{D}(p = 1)  908  24,698  1,232,944  112,018,190  18,410,581,880
σ̄_{D}(p = 1)  0.454  0.482  0.501  0.515  0.525
B_{D+1}/B_{D}  27.20  49.92  90.85  164.35  295.97 ± 0.04

TABLE II. Convergence of D = 9 path count data for increasing run length. Measured numerical errors on the latter quantity (Num. err.) are always smaller than estimated error bars (Estim. err.). The latter are calculated using numerically measured standard deviations and conjectured auto-correlation times (section IV 1; [^{24}]).

N_{MC}  10^{3}  10^{4}  10^{5}  10^{6}  10^{7}  exact
P̄_{9}  159.264  164.176  164.411  164.371  164.344  164.35
log P̄_{9}/9  0.563396  0.566771  0.566930  0.566903  0.566884  0.566891
Num. err.  3.5 10^{−3}  1.2 10^{−4}  3.9 10^{−5}  1.2 10^{−5}  6 10^{−6}
Estim. err.  8.5 10^{−3}  2.7 10^{−3}  8.5 10^{−4}  2.7 10^{−4}  8.5 10^{−5}

TABLE III. Exact data for fixed boundary tilings, p = 2 and D = 1, 2, …, 6.

D  1  2  3  4  5  6
B_{D}(p = 2)  1  1  20  5,383  16,832,230  327,674,047,931
σ̄_{D}(p = 2)  0  0.250  0.358  0.416  0.451

TABLE IV. The first values of the limiting distribution of coordination numbers Z, obtained both in the mean-field approximation and numerically, by Monte Carlo simulations.

Z  3  4  5  6  7
Mean-field  0.33  0.41  0.20  0.05  0.009
Numerical  0.31  0.43  0.21  0.04  0.005

TABLE V. Path count statistics for strained tilings.

D  p = 1, 1, ...  p = 1, 1, ..., 2, 2, ...  p = 1, 1, ..., 3, 3, ...  p = 1, 2, 1, 2, ...
10  0.569  0.572  0.572  0.573
20  0.573  0.572  0.571  0.574
30  0.572  0.571  0.570  0.572
40  0.572  0.570  0.570  0.571
50  0.571  0.571  0.570  0.570