\documentstyle[aps,prb,psfig]{revtex} 
\input{macros.phys}

\def\sinf{\sigma_{\infty}}
\def\sinfbar{\bar{\sigma}_{\infty}}
\def\eh{{\bf \hat{e}}}
\def\fg{{\bf f}}

\renewcommand{\thesubsection}{\arabic{subsection}}

\begin{document}

\title{Random tilings of high symmetry: II. Boundary conditions and numerical
studies}

\author{M. Widom}
\address{Department of Physics, Carnegie Mellon University, 
Pittsburgh, PA  15213, USA.}

\author{N. Destainville}
\address{Laboratoire de Physique Th\'eorique, UMR CNRS-UPS 5152,
Universit\'e Paul Sabatier,
31062 Toulouse Cedex 04, France.}

\author{R. Mosseri} \address{Groupe de Physique des Solides,
CNRS et Universit\'e 6, 
Campus Boucicaut, 140 Rue de Lourmel, 75015 Paris, 
France.}
 
\author{F. Bailly}
\address{LPSB-CNRS, 92195 Meudon Cedex, France.}


\date{\today}
\maketitle


\begin{abstract}
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.
\end{abstract}


\section{Introduction}
\label{intro}

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~\cite{paperI}, 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~\cite{Elser}, the membrane being the union of
hypercube faces (see paper I~\cite{paperI}, section II.1).

Rhombus tilings are dual to de Bruijn grids~\cite{paperI,debruijn}. 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 {\em 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~\cite{Shechtman,Levine}, metal alloys that exhibit
rotational symmetry forbidden by conventional crystallography. One of
the key properties of random tilings~\cite{Henley91} 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 \ra \infty$
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.

\begin{figure}[htb]
\begin{center}
\ \psfig{figure=fixed_D2.eps,width=12.5cm} \
\end{center}
\caption{Fixed boundary $D \ra 2$ tilings with $D=3,5$ and $7$. The
  side lengths of the polygons are all $p=8$.}
\label{pavagesD2}
\end{figure}

As usual in statistical mechanics, tiling systems can have different
boundary conditions, such as free, fixed or periodic. For example, the
tilings in figures~\ref{pavagesD2}, \ref{ex40.2} or
\ref{DeBruijn4.2}, 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 \ra
2$ tilings (see Fig.~\ref{pavagesD2}, left) equals
0.261~(\cite{elsershape}), while the free boundary entropy equals
0.323~(\cite{Wannier}). This phenomenon can be
understood~\cite{elsershape,Grensing,matheux,Bibi97B} by inspection of
Fig.~\ref{pavagesD2}. 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~\cite{paperI} developed a mean-field theory for free boundary
tilings, applicable to the limit of high rotational symmetry.  Two
earlier papers~\cite{Widom97,Bibi00} presented initial studies of this
problem. The first paper~\cite{Widom97} 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~\cite{Bibi00} 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 \ra \infty$. 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~\ref{fixed} 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~\ref{thermo} where we argue that free and fixed boundary
tilings attain the same thermodynamic limit as $D \ra \infty$. Then,
section~\ref{simul} 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.

\section{Fixed boundary tilings}
\label{fixed}

Fixed boundary tilings such as those of figure~\ref{pavagesD2} 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 {\em 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~\cite{elsershape,Grensing,matheux,Bibi97B,4to3}. 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~\ref{thermo} that when $D$
becomes large (figure~\ref{ex40.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 $\bar{\sigma}_D$ and their limiting value by
$\sinfbar$. 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~\cite{octo01} (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
\begin{equation}
N_T = \sum_{1\leq i<j \leq D} k_i k_j.
\label{NT}
\end{equation}
If all $k_i$ are set to 1, $N_T = D(D-1)/2$.


\begin{figure}[ht]
\begin{center}
\ \psfig{figure=rt_40_1_12E6.eps,width=5cm} \
\end{center}
\caption{A $40 \ra 2$ fixed boundary tiling with one de Bruijn line in
  each family.}
\label{ex40.2}
\end{figure}

The relationship between tilings and partitions has been widely
explored in references~\cite{elsershape,octo01,Mosseri93B,Bibi97}. The
idea is to code a random rhombus tiling by an array of integers
satisfying certain ordering constraints.  Figure~\ref{ex.part}
illustrates this point in the simple hexagonal case~\cite{elsershape}:
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 \times 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.

\begin{figure}[ht]
\begin{center}
\ \psfig{figure=mapping32.eps,height=4.5cm} \
\end{center}
\caption{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.}
\label{ex.part}
\end{figure}

This point of view can be extended to any $D \ra d$
problem~\cite{octo01,Mosseri93B,Bibi97}: any $D \ra d$ tiling can be
coded with a generalized partition on an array related to a $D-1 \ra
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~\ref{pavagesD2}).

We now present an analogy between fixed boundary tilings and
algorithms for sorting lists~\cite{octo01}. Although Computer Science is
not our motivation, we are interested in the same enumeration problem
\cite{Knuth92,Bjorner93}. 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\cite{matheux,Kenyon93,Elnitsky97,Bailey97}. We shall use results from
these fields throughout this paper.

In the sorting language, a {\em comparator} $[i;j]$ acts on a list
$(x_1, x_2,\ldots, 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 \cite{Knuth92}, we call a {\em complete} sorting
algorithm a sequence of such comparators which sorts in the increasing
order any list of real numbers $(x_1,x_2,\ldots,x_D)$. This sorting
algorithm will be called {\em 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 \leq x_j$. Knuth shows that a sequence of comparators is a
sorting algorithm if it correctly sorts the completely reversed list
$(D,D-1,\ldots,1)$. This means that a complete primitive sorting
algorithm is a sequence of comparators $[i;i+1]$ that transforms the
list $(D,D-1,\ldots,1)$ into the list $(1,2,\ldots,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~\ref{sort} 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 $\simp{D}{2}$ crossings.

\begin{figure}[ht]
\begin{center}
\ \psfig{figure=sort.eps,height=2.5cm} \
\end{center}
\caption{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]$.}
\label{sort}
\end{figure}

We now establish the link between sorting algorithms and the de Bruijn
representations~\cite{paperI,debruijn} of $D \ra 2$ tilings. Each line
of the diagram~\ref{sort} 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~\ref{sort}, the fourth
and the fifth comparator ({\em 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
\cite{Knuth92,Elnitsky97}. 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
\geq B_D$). Since each pattern of crossings (equivalence class of
algorithms) defines a tiling, the tiling entropy density per tile is
\begin{equation}
\label{sigmaD}
\bar{\sigma}_D={{1}\over{N_T}} \log{B_D}.
\end{equation}

\begin{figure}[ht]
\begin{center}
\ \psfig{figure=SdeD.eps,width=7cm} \
\end{center}
\vspace{-4mm}
\caption{Entropy per tile $\bar{\sigma}_{D}$ of fixed boundary
tilings ($p=1$). We indicate the limiting values $0.568$ (full line)
and the close upper bound $\log 2 \simeq 0.693$ (dashed line).}
\label{plot.pavages.finis}
\end{figure}

Sorting algorithms are easy to manipulate and provide efficient
enumeration numerical tools.  This analogy provides us with the number
$B_D$ of $D \ra 2$ tilings for small values of $D \le 10$ from the
work of Knuth \cite{Knuth92} as extended by us~\cite{Widom97} (see
Table~\ref{pavages.finis}). Unfortunately it is not possible to reach
the next value ($D=11$) for the foreseeable future using our current
algorithms.  Figure~\ref{plot.pavages.finis} suggests that the entropy
$\sinfbar(p=1)$ has a well defined limit when $D$ tends to infinity.
Indeed, Knuth~\cite{Knuth92} gives lower and upper bounds on $B_D$
from which we get
\begin{equation}
{1 \over 3} \log 2 \leq \lim_{D \ra \infty} \bar{\sigma}_{D}
\leq 2 \log 2. 
\end{equation}
Moreover, Bj\"orner (see reference \cite{Bjorner93}, p.~270) 
derives a better upper bound: 
\begin{equation}
\lim_{D \ra \infty}  \bar{\sigma}_{D} \leq 1.44 \log 2.
\end{equation}
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~\cite{paperI} that for
free boundary entropies $\lim_{D \ra \infty} \sigma_D \leq
\log{2}$, we finally get:
\begin{equation}
\label{Ubound}
\lim_{D \ra \infty} \bar\sigma_D  \leq \log{2}.
\end{equation}

This point of view generalizes to systems with more than one line per
de Bruijn family~\cite{octo01}, leading to the definition of {\em
partial} sorting algorithms for merging pre-ordered lists of
numbers. Suppose that we have $D$ families of $k_i$ numbers each
($i=1,\ldots,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,\ldots,k_D$. In
reference~(\cite{octo01}), we derive an analytic expression for
$A_D(k_1,\ldots,k_D)$. However, it proves impossible to calculate
$B_D$ analytically.

The diagonal case, with $k_1=k_2=\cdots=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~\ref{pavages.finis2}).

\section{Thermodynamic limit}
\label{thermo}

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
\ref{fixe} and \ref{encadre}. 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~\cite{Wannier} and
0.261~\cite{elsershape}. For octagonal tilings with $D=4$ the values
are 0.434~\cite{lpw} and 0.36(1)~\cite{octo01}. Therefore the relative
difference $(\sigma-\bar{\sigma})/\sigma$ decreases between $D=3$ and
$D=4$, consistent with a vanishing difference at large $D$.
Data on $p=1$ and $p=2$ tilings
(Tables~\ref{pavages.finis} and~\ref{pavages.finis2}) suggest that
when $D$ becomes large, the entropy becomes independent of $p$.

Furthermore, a na\"{\i}ve observation of figures~\ref{pavagesD2} 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 \cite{Bibi97}, 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 $\sinf$.

More precisely, let $(x_i)$ be any set of numbers $0 < \alpha
\leq x_i \leq \beta$, where $\alpha$ and $\beta$ are fixed positive
real numbers. Consider fixed boundary tilings of a $2D$-gon with edge
lengths $k_i=x_i p$ as $p \ra \infty$. 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~\cite{paperI} (section II.3)).
The tile fractions are
\begin{equation}
\bar{n}_{ij} = \frac{k_i k_j}{N_T} = \frac{k_i k_j}{k_1 k_2 + \ldots +
k_{D-1} k_D},
\end{equation}
owing to eq.~(\ref{NT}). In an unstrained tiling, where all $k_i$ are
equal,
\begin{equation}
\bar{n}_{ij}^* = \frac{2}{D(D-1)}
\end{equation}
whereas the condition $0 < \alpha \leq x_i \leq \beta$ imposes
\begin{equation}
\bar{n}_{ij} \geq \frac{\alpha^2}{\beta^2} \; \bar{n}_{ij}^*
= \frac{\alpha^2}{\beta^2} \; \frac{2}{D(D-1)}.
\end{equation}
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~\ref{fixe} that the
sequence $\bar{\sigma}_D$ reaches a thermodynamic limit
$\bar{\sigma}_{\infty}$ 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 $\sinf$, and $\lim_{D \ra \infty} \bar{\sigma}_D =
\sinf$. 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 \ra \infty$ and $p \ra \infty$ commute with each
other. The proofs of all these results are given in
appendix~\ref{fixe}.

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~\cite{paperI}. 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.

\section{Simulations}
\label{simul}

Monte Carlo numerical simulations~\cite{Newman} are widespread in the
fields of random tilings and quasicrystals where there exist many
difficult unsolved theoretical questions~\cite{mc}. 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 {\em via}
mean-field arguments~\cite{paperI}, we perform Monte Carlo
calculations on large $D$ tilings.  The $D \ra 2$ configuration space
is sampled for large $D$ systems {\em via} single-vertex flip
dynamics. This method is validated by the connectivity of the space of
configurations in two dimensions~\cite{Kenyon93,octo01}.

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.


\subsection{Monte Carlo algorithm}
\label{MCalgo}

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
{\em 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 {\em detailed balance condition}~\cite{Newman}:
\begin{equation}
w(\CC_1 \ra \CC_2) = w(\CC_2 \ra \CC_1),
\label{det:bal}
\end{equation}
where $w(\CC_1 \ra \CC_2)$ denotes the transition probability from
configuration $\CC_1$ to configuration $\CC_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~(\ref{det:bal}) is established.

Note that the flip acceptance rate equals the fraction (about 31\% as
reported in section~\ref{vertex:stats}) 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.~(\ref{det:bal}) is violated.

Another point to consider in Monte Carlo simulations is the
correlation time $\tau$ (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~\cite{Mixing} and references
therein), this question has no rigorous definitive answer for large
$D$. However, our numerical calculations are all in agreement with the
conjecture $\tau \leq N_V/2$.

\subsection{Path counting and entropy estimate}

To estimate the entropy $\sigma_D$ of fixed boundary tilings, we shall
use a path-counting algorithm based upon de Bruijn line
enumeration~\cite{Widom97} directly derived from the general method
described in section 2.2 of paper I~\cite{paperI}.
Figure~\ref{iter.tilings} illustrates the method.  At each step, we
build a $D+1 \ra 2$ tiling from a $D \ra 2$ one.  Starting from any $D
\ra 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$.

\begin{figure}[ht]
\begin{center}
\ \psfig{figure=worm.eps} \
\end{center}
\caption{Iterative construction of fixed boundary tilings ($D \ra 2$,
  $p=1$).}
\label{iter.tilings}
\end{figure}

Let $\bar{P}_D$ denote the mean number of bottom-to-top paths on the $D \ra
2$ tilings:
\begin{equation}
\bar{P}_D = {1 \over B_D} \sum_{\tau=1}^{B_D} P_D(\tau),
\end{equation}
where $P_D(\tau)$ is the exact number of paths on the tiling $\tau$,
and $B_D$ is the number of such tilings. Then
\begin{equation}
B_{D+1} = \bar{P}_D B_D.
\label{iter.rel}
\end{equation}
Iterating relation~(\ref{iter.rel}), and taking the logarithm yields
\begin{equation}
\label{logBD}
\log{B_D} = \sum_{k=1}^{D-1} \log{\bar{P}_k}.
\end{equation}
For large $D$, the ratio ${{1}\over{D}} \log{\bar{P}_D}$ approaches a finite
limiting value.  Taking the limit of~(\ref{logBD}) as $D \ra \infty$,
and noting eqs.~(\ref{NT}) and~(\ref{sigmaD}), we find
\begin{equation}
\label{sigmalimit}
\lim_{D \ra \infty} \bar{\sigma}_{D}
= \lim_{D \ra \infty} {{1}\over{D}} \log{\bar{P}_D}.
\end{equation}
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 \ra 2$
tilings. On each tiling $\tau$ we can quickly evaluate $P_D(\tau)$
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.~\ref{pathcount}). At the end of the
process, the value of the top vertex is $P_D(\tau)$.

\begin{figure}[ht]
\begin{center}
\ \psfig{figure=pathcount.eps,width=4cm} \
\end{center}
\caption{Path count construction, as described in the text.}
\label{pathcount}
\end{figure}

We tested this algorithm on small $D$ systems, where numerical values
are in good agreement with the exactly known mean number of paths
$\bar{P}_D=B_{D+1}/B_D$ (table~\ref{pavages.finis}). For example,
table~\ref{test.accuracy} 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~\ref{Mean}. In this figure, error bars
are calculated as follows: the standard deviation $\Delta$ is measured
numerically. The error is estimated by $\Delta/\sqrt{I}$,
where $I$ is the number of {\em independent} measures in the
simulation (that is~\cite{Newman} its total length divided by twice
the correlation time $\tau$). Note that the actual errors calculated
in table~\ref{test.accuracy} are compatible with the previous error
bar.

This graph presents a maximum near $D=20$ which is related to finite
$D$ {\em and} finite size corrections. Indeed, it is demonstrated in
paper I~\cite{paperI} that for free boundary systems, finite $D$
corrections to $\log \bar{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
\begin{equation}
{\log \bar{P}_D \over D} = \bar{\sigma}_{\infty} + {Cst. \over D} - {Cst.'
\over D^2} + \OC({1 \over D^3}).
\end{equation}

\begin{figure}[ht]
\begin{center}
\ \psfig{figure=Snew.eps,width=10cm} \
\end{center}
\caption{Value of $\displaystyle{\log \bar{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$.}
\label{Mean}
\end{figure}

The numerical results in figure~\ref{Mean} 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 $\displaystyle{\log \bar{P}_D /D}$ can be
estimated as $\bar{\sigma}_{\infty} = 0.5676 \pm 0.0001 $. The uncertainty
estimate on this last limit comes from excluding the data with
$D=120$.



\subsection{Vertex statistics}
\label{vertex:stats}

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~\cite{paperI}.  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.

\medskip

\begin{figure}[ht]
\begin{center}
\ \psfig{figure=Z_3.eps,height=6.3cm}~\hspace{1cm}~\psfig{figure=Z_5.eps,height=6.3cm} \
\caption{Fractions of $Z=3$ and $Z=5$ vertices in $D \ra 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.}
\label{Z_fractions}
\end{center}
\end{figure}

We have run Monte Carlo simulations up to $D=100$. After plotting the
data as functions of $1/D$, we extrapolate the limiting values {\em
via} quadratic fits.  Figure~\ref{Z_fractions} 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 \ra \infty$ 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 \ra \infty$ 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~\cite{paperI}. 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~\ref{Z_num} the probability distribution
for $Z$, and compare it with the predictions of mean-field theory.
The agreement is also satisfactorily good.



\subsection{Phason strained tilings}
\label{phason:strain}

The data just reported in section~\ref{vertex:stats} 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,\ldots,D/2)$, the remaining $D/2$
consecutive line families occur with $k_i=2$ $(i=D/2+1,\ldots,D)$; the
same low frequency with a larger amplitude so that $k_i=3$
$(i=D/2+1,\ldots,D)$ replaces $k_i=2$; high frequency strain with
$k_i=1$ ($i$ odd) and $k_i=2$ ($i$ even).  Figure~\ref{fig:strain}
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~\ref{tab:strain})
cannot be separated from finite size effects. Thus entropy does not
depend on phason strain, supporting {\em a posteriori} arguments
developed for free boundary tilings~\cite{paperI}. Recall that
strain-independence of entropy in fixed boundary tilings relies upon
strain-independence in free boundary tilings (see
section~\ref{thermo}).

\begin{figure}[ht]
\begin{center}
\ \psfig{figure=strain.eps,width=6.5cm}
\end{center}
\caption{Phason strained large $D$ tiling. The added paths still go from
bottom to top.}
\label{fig:strain}
\end{figure}

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~\ref{coord:strained}).  These data
support the existence of a universal local structure, independent of
strain and boundary conditions in the limit $D \ra \infty$.


\begin{figure}[ht]
\begin{center}
\
\psfig{figure=Z_3str.eps,height=6.3cm}~\hspace{1cm}~\psfig{figure=Z_5str.eps,height=6.3cm}
\
\end{center}
\caption{Fractions of $Z=3$ and $Z=5$ vertices in strained tilings 
($p=1,1,\ldots,3,3,\ldots$) (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.}
\label{coord:strained}
\end{figure}


\section{Conclusion}
\label{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\leq 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 $\bar{P}_D$) as our
mean-field theory in paper I~\cite{paperI}.  We obtain a very
accurate estimate of the entropy per tile, $\bar{\sigma}_{\infty} =
0.5676 \pm 0.0001 $, compared to the approximate mean-field value
$\sigma_{\infty}^{MF} \simeq 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~\ref{ex40.2}). The entropy
{\em 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 higher-dimensions, because the vision
of large regions with at most one de Bruijn {\em surface} per family
still holds. The explicit proof should be similar to the
appendix~\ref{fixe}, with only minor changes to account for the
possible spontaneous decomposition into regions of lower effective
codimension (as in the ``arctic octahedron'', see Ref.~\cite{4to3}).
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 {\em 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.

\section*{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

\section{Thermodynamic limit of two-dimensional fixed boundary
  tilings}
\label{fixe}

In this technical appendix, we demonstrate that the fixed boundary
entropy $\bar{\sigma}_D$ attains a finite limit $\sinfbar$ when $D
\ra \infty$, and that this limit coincides with the free boundary entropy
$\sinf$. 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~[\cite{matheux,Bibi97B}] 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 {\em 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' \ra 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 $\sigma_{\infty}$ (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.

\subsection{Diagonal tilings}
\label{diagtilings}

The demonstration uses a particular tiling denoted as $\TC_0$.  It is
the dual of a particular de Bruijn grid~\cite{paperI,debruijn},
consisting of $D$ families of $p$ {\em straight} lines equally spaced
out and forming a regular fan~\cite{fan} (see
figure~\ref{DeBruijn4.2}, left). The de Bruijn families are denoted by
$F_1,F_2,\ldots,F_{D}$, and are labeled counterclockwise so $F_k$
makes an angle $k \pi / D$ with an arbitrary reference
direction. Figure~\ref{DeBruijn4.2} displays such a grid and its dual
tiling. For finite $D$, such a tiling is known~\cite{Bibi97B} not to
maximize the entropy functional among the tilings with the same boundary
(see below).  By construction, the tiling $\TC_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 {\em a
priori} differs from the free boundary $D \ra 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 ${\cal
A}_{D'}$ and called the region of ``effective dimensionality'' $D'$.

\begin{figure}[ht]
\begin{center}
\ \psfig{figure=octolines.eps,height=6cm} \hspace{1.5cm} \psfig{figure=octo.eps,height=6cm} \
\end{center}
\caption{The de Bruijn straight lines (left) and the
corresponding rhombus tiling in the $4 \ra 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 ${\cal A}_2$, only 2 families meet; In the intermediate
one ${\cal A}_3$, 3 families meet; And in the central region ${\cal
A}_4$, all 4 families intersect.
\label{DeBruijn4.2}}
\end{figure}

Define $S_0^{D \ra 2}$ as the entropy functional evaluated at $\TC_0$.
It satisfies $S_0^{D \ra 2} \leq \bar{\sigma}_{D}$.  Then
\begin{equation}
\lim_{D \ra \infty} S_0^{D \ra 2} \leq \bar{\sigma}_{\infty}
\leq \sigma_{\infty}.
\label{1st.ineq}
\end{equation}
Next we intend to show that for any real number $h < 1$,
\begin{equation}
\lim_{D \ra \infty} S_0^{D \ra 2} \geq h \; \sigma_{\infty}.
\label{alpha.ineq}
\end{equation}
It then follows that 
\begin{equation}
\lim_{D \ra \infty} S_0^{D \ra 2} \geq \sigma_{\infty},
\end{equation}
and, owing to relation~(\ref{1st.ineq}), that 
\begin{equation}
\sigma_{\infty} \leq 
\lim_{D \ra \infty} S_0^{D \ra 2} \leq \bar{\sigma}_{\infty}
\leq \sigma_{\infty},
\end{equation}
and therefore that the 3 involved quantities are equal.

Let us prove the above statement~(\ref{alpha.ineq}). We use the fact
(see section~\ref{phason:strain} and also paper I\cite{paperI}), that
the entropy $\sigma_{\infty}$ of free boundary tilings does not depend
on strain in the large $D$ limit. We show that the local entropy of
$\TC_0$ equals $\sigma_{\infty}$ in (nearly) all regions $\AC_{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.~\ref{DeBruijn4.2}).
Figure~\ref{2familles} (left) displays a schematic representation of a
grid were the different regions under consideration clearly appear.

The number of intersections in a region $\AC$ of area $\alpha$ 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 $\theta_{ij}$. The number of such
intersections is the area $\alpha$ divided by the area $\alpha_0$ of
the unit cell of the lattice defined by these (only) two families
(figure \ref{2familles}, right). If the distance which separates two
lines of a family is set to 1, then $\alpha_0 = 1/|\sin
\theta_{ij}|$. Hence the total number of intersections is
\begin{equation}
\sum_{i<j} {\alpha \over 1/|\sin \theta_{ij}|} = \alpha \sum_{i<j} |\sin
\theta_{ij}|,
\end{equation}
where the indices $i$ and $j$ run over the families present in $\AC$.

\begin{figure}[ht]
\begin{center}
\ \psfig{figure=familles2.eps,height=4cm} \hspace{1cm} 
\psfig{figure=familles1.eps,height=4cm} \
\end{center}
\caption{Left: Schematical representations of the $D$ families of de
Bruijn lines and of their intersections, to be compared to
fig.~(\ref{DeBruijn4.2}), left.  When $D$ tends to infinity, the
polygonal boundary tends toward a circle of circumference $2 D p$ and
the polygonal central region ${\cal A}_D$ tends toward a circle a diameter
$p$. Right: Two families of de Bruijn lines making an angle
$\theta_{ij}$.  The number of intersections of these two families in a
region equals its area $\alpha$ divided by the area $\alpha_0$ of
the grayed unit cell.
\label{2familles}}
\end{figure}

As displayed in figure~\ref{2familles} (left) a region ${\cal A}_{D'}$
is a crown made of $2D$ kites (or triangles for the inner crown) of
equal areas. A kite (or triangle) of ${\cal A}_{D'}$ will be denoted
by $\KC_{D'}$ and its area by $\kappa_{D'}$. Henceforth the area of
${\cal A}_{D'}$ is equal to $\alpha_{D'} = 2D \kappa_{D'}$. In
$\KC_{D'}$, there are $D'$ {\em adjacent} families of lines, for
example $F_1,\ldots,F_{D'}$. The angles that they make relatively to
an arbitrary reference direction are therefore $\theta_l = \theta_0 +
{l \pi / D}$, $l=1,\ldots,D'$ and the number of intersections in
$\KC_{D'}$ is
\begin{equation}
Q_{D'} = \kappa_{D'} \left( \sum_{1 \leq i<j \leq D'} \sin {(j-i) \pi
\over D} \right) .
\label{Q}
\end{equation}
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.~(\ref{Q}) 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
\begin{equation}
\phi_D(x) \equiv
\sum_{1 \leq i<j \leq E(x D)}  \sin {(j-i) \pi \over D} 
\simeq D^2 \int_0^{x} dy \int_0^y dz \sin \pi (y-z)  \\
= {D^2 \over \pi^2 } \left[ \pi x -\sin \pi x \right].
\end{equation}			

The central region ${\cal A}_D$ plays a particular role
among the regions ${\cal A}_{D'}$: when $D$ is large, the fraction of
tiles in ${\cal A}_D$ tends to $1/2$. Indeed, as illustrated in figures
\ref{DeBruijn4.2} (left) and \ref{2familles} (left), the
central region tends toward a circle of diameter $p$. Its area is $\pi
p^2/4$, and, owing to the above result applied to the case $x=1$, the
number of intersections that it contains grows like
\begin{equation}
{D^2 \over \pi} {\pi p^2 \over 4} = {D^2 p^2 \over 4}.
\end{equation}
Since the total number of tiles grows like $\displaystyle{D^2 p^2
\over 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.

\begin{figure}[ht]
\deuxfigs{regions.eps}{kites.eps}
\caption{(left) around the central region ${\cal A}_D$, there is a
first crown of triangles $\KC_{D-1}$ (light gray), a second crown of
kites $\KC_{D-2}$ (medium gray), a third crown of kites $\KC_{D-3}$
(dark gray) and so forth; in this figure, $D=6$. Right: the three
kinds of kites which define a kite-shaped region $\KC_{D'}$
(white): a big kite {\em minus} two small kites (light gray) {\em
plus} a smaller kite (dark gray).
\label{regions}}
\end{figure}

Figure~\ref{regions} (left) shows the geometry of the regions under
consideration when $D$ is finite. Around the central region ${\cal
A}_D$, there is a first crown triangular domains, then a second crown
of kite-shaped regions, and as $k$ increases, the regions ${\cal
A}_{D'}$ lie in concentric crowns around the central one. Each crown
contains $2D$ kite-shaped (or triangle-shaped) regions. As
figure~\ref{regions} (right) illustrates, such a kite can be seen
as a greater rectangular kite\cite{rectkite} {\em 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 $\delta_{D'}$, $\delta_{D'+1}$
and $\delta_{D'+2}$, then the area of ${\cal K}_{D'}$ is
\begin{equation}
\kappa_{D'}=\delta_{D'} - 2 \delta_{D'+1} + \delta_{D'+2}.
\end{equation}

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 $\theta$, the area of the kite is
\begin{equation}
\delta(\theta)={p^2 \over 4} \tan\left( {\theta \over 2} \right).
\end{equation}
For the domain ${\cal K}_{D'}$, $\theta=\pi (D-D')/ D$ for the
greater kite.  Henceforth,
\begin{equation}
\delta_{D'}={ p^2 \over 4} \tan \left( \pi {D-D'\over 2 D} \right).
\end{equation}
The same formula holds for $\delta_{D'+1}$ and $\delta_{D'+2}$. Thus
if we set $k=D-D'$,
\begin{equation}
\kappa_{D'}=
{p^2 \over 4} \left[
\tan \left( {k \pi \over 2 D} \right)
-2 \tan \left( {(k-1) \pi \over 2 D} \right)
+\tan \left( {(k-2) \pi \over 2 D} \right) \right].
\end{equation}
Fixing the ratio $x=D'/D$, and approximating the above as a second
derivative, we finally get 
\begin{equation}
\kappa_{D'} = \kappa_{E(xD)} \simeq \pi^2 {p^2 \over 8} {1 \over D^2} 
\tan \left( \pi (1-x)/ 2  \right) \left[ 1 +
\tan^2 \left( \pi (1-x) / 2 \right) \right].
\label{alpha.k} 
\end{equation}
Thus the number of tiles in $\AC_{E(xD)}$ is
\begin{equation}
N_{E(xD)} = \alpha_{E(xD)} \phi_D(x) = 2D \kappa_{E(xD)} \phi_D(x).
\end{equation}
Finally, the fraction of tiles in the region
${\cal A}_{E(xD)}$ is
\begin{equation}
n_{E(xD)} = \frac{N_{E(xD)}}{N_T}= \frac{N_{E(xD)}}{D^2 p^2 /2} \simeq
{1 \over 2D } \tan \left( \pi (1-x)/ 2  \right) \left[ 1 +
\tan^2 \left( \pi (1-x) / 2 \right) \right]
\left[\pi x - \sin \pi x
\right]
\end{equation}
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 \leq g \leq 1 $ and let us
compute the fraction of tiles which lie in regions where at least
$E(g D)$ families of lines intersect:
\begin{equation}
{1 \over 2} + \sum_{l=E(g D)}^{D-1} n_{l} \simeq 
{1 \over 2} + {1 \over 2} \int_{0}^{1-g} dx 
\tan \left( \pi x / 2 \right)
\left[ 1 + \tan^2 \left( \pi x / 2 \right) \right] \left[ \pi (1-x) - 
\sin \pi x \right] 
\equiv {1 \over 2} + {1 \over 2} f(1-g).
\end{equation}
The above integral is a continuous function $f(g)$ for $g \in [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(g D)$
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~\ref{encadre}
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 $\sigma_{\infty}$. A fraction $h$
of the tiles lie in regions where the local entropy per tile tends to
$\sigma_{\infty}$ when $D \ra \infty$.  We get~\cite{Bibi97B} the
expected relation~(\ref{alpha.ineq}) and we conclude that free and
fixed boundary entropies coincide. We also remark that the tiling
$\TC_0$ which was not supposed {\em 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 $\KC_{D'}$
can be very small and can contain very few tiles (and even no tile at
all).  Nevertheless in the previous demonstration, regions
$\AC_{E(xD)}$ with $x_0 \leq x \leq x_0 + \delta x$ can be put
together into larger regions, to which all the previous arguments can
be applied.

\subsection{Non-diagonal tilings}
\label{nondiagtilings}

For large $p$ non-diagonal tilings with side lengths $k_i = x_i p$
where $\alpha \leq x_i \leq \beta$, we use a particular tiling
$\TC'_0$ which is a variation of $\TC_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 $\kappa_{D'}$ or
$\alpha_{D'}$.

We choose $l_i = p/(x_i p -1) \simeq 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~\ref{2familles}, left). Therefore the areas $\alpha_{D'}$
remain unchanged as compared to the previous section.

On the other hand, the density of intersections in each domain
$\KC_{D'}$ will vary because of varying line spacing. In particular
the area $\alpha_0$ of the unit cell now depends on $i$ and $j$:
\begin{equation}
\alpha_0(i,j) = \frac{l_i l_j}{| \sin \theta_{ij} |} = \frac{1}
{x_i x_j | \sin \theta_{ij} |}.
\end{equation}
The number of intersections~(\ref{Q}) becomes
\begin{equation}
\tilde{Q}_{D'} = \kappa_{D'} \left( \sum_{1 \leq i<j \leq D'} x_i x_j
\sin {(j-i) \pi \over D} \right) \leq \beta^2 Q_{D'} .
\end{equation}
Henceforth, the new numbers of tiles in the crown $\AC_{D'}$ and in
the whole tiling satisfies $\tilde{N}_{D'} \leq \beta^2 N_{D'}$ and
$\tilde{N}_{T} \geq \alpha^2 N_{T}$. If we set again $D'=E(xD)$, the
new fraction of tiles in $\AC_{D'}$ satisfies
\begin{equation}
\tilde{n}_{E(xD)} = \frac{\tilde{N}_{E(xD)}}{\tilde{N}_{T}} \leq 
\frac{\beta^2}{\alpha^2} \; n_{E(xD)}.
\end{equation}
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 $\sinf$. The proof is explicited in
appendix~\ref{encadre} in the diagonal case but can easily be
adapted to the present case, leading to
\begin{equation}
\frac{n'_{ij}}{(n')^*_{ij}} \geq \frac{\alpha^2}{\beta^2} \; \frac{6}{\pi^2},
\end{equation}
where the notations are defined in appendix~\ref{encadre}. Even though
the lower bound on tile fractions now depends on $\alpha$ and $\beta$
because we are dealing with a nondiagonal case, it remains finite.


\section{Regularity of large codimension fixed boundary tilings}
\label{encadre}

We check here that in each region $\AC_{D'}$ of the diagonal tiling
$\TC_0$ the local entropy per tile tends to the free boundary one
$\sigma_{\infty}$ as $D\ra\infty$.  As discussed in paper
I~\cite{paperI}, a sufficient condition is that tile fractions should
be ``bounded'' on $\AC_{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
\begin{equation}
n_{ij} \geq a \; n_{ij}^*
\label{encadreappend}
\end{equation} 
for all $i$ and $j$. 
We want to calculate the tile fractions $n'_{ij}$ in $\AC_{D'}$ 
and compare them
to the strain-free ones. Recall that in a $D \ra 2$ tiling,
\begin{equation}
n^*_{ij} = Cst. \; \sin\left| \pi \; \frac{j-i}{D} \right|,
\label{prop}
\end{equation}
where the constant comes from the normalization relation $\sum 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 $\AC_{D'}$ under consideration,
not all tile species occur. Without loss of generality, we assume that
only the tiles with indices $i,j=1,\ldots,D'$ exist. These tiles
appear in the region ${\cal A}_{E(xD)}$ with the
fractions~(\ref{prop}), but with a different normalization
constant because only some of them appear. The normalization relation now reads
\begin{equation}
\sum_{1\leq i<j \leq D'} n'_{ij}=1.  
\end{equation}
At the large $D$ limit, this sum can be replaced by an integral and we
get the tile fractions in ${\cal A}_{E(xD)}$:
\begin{equation}
n'_{ij} = \frac{1}{(D')^2} \; \frac{\pi^2 x^2}{\pi x - \sin (\pi x)} \;
\sin \left| x \pi \; \frac{j-i}{D'} \right|.
\end{equation}

Now we desire to check that this tiling has bounded local fractions of
tiles {\em when it is considered as a $D' \ra 2$ tiling}. We need to
compare the above $n'_{ij}$ to the corresponding quantities in an
unstrained $D' \ra 2$ tilings, in other words to
\begin{equation}
(n')^*_{ij} = \frac{\pi}{(D')^2} \; \sin \left| \pi \; \frac{j-i}{D'} \right|.
\end{equation}
By convexity of the sin function on the interval $[0,\pi]$, we have: 
$\sin | x \pi (j-i)/D' | \geq x \sin | \pi (j-i)/D' |$, and
\begin{equation}
\frac{n'_{ij}}{(n')^*_{ij}} \geq \frac{x}{\pi} \; \frac{\pi^2 x^2}{\pi
x - \sin (\pi x)} \geq \frac{6}{\pi^2},
\end{equation}
which achieves the proof: $a=6/\pi^2$ in
condition~(\ref{encadreappend}). From the results of paper
I~\cite{paperI}, we conclude that in such a region $\AC_{D'}$, the
local entropy per tile equals the free boundary value $\sigma_{\infty}$
when $D \ra \infty$.

\bigskip

\begin{thebibliography}{99}

\bibitem{paperI} N. Destainville, M. Widom, R. Mosseri, F. Bailly,
``Random tilings of high symmetry: I. Mean-field theory'', to be
submitted to {\em J. Stat. Phys.} (2003).

\bibitem{Elser} V. Elser, {\em Phys. Rev. Lett.} {\bf 54}, 1730 (1985).

\bibitem{debruijn} N.G. de Bruijn, {\em Ned. Akad. Wetensch. Proc.}
{\bf A84}, 39 (1981); {\em J. Phys. France} {\bf 47}, C3-9 (1986).

\bibitem{Shechtman} D. Shechtman, et al.  {\em Phys. Rev. Lett.} {\bf
53},1951 (1984).

\bibitem{Levine} D. Levine, P.J. Steinhardt, {\em Phys. Rev. Lett.}
{\bf 53}, 2477 (1984).

\bibitem{Henley91} C.L. Henley, {\em in} {\em Quasicrystals, the State
of the Art}, Ed. D.P.  Di Vincenzo, P.J.  Steinhart (World Scientific,
1991), 429.

\bibitem{elsershape} V. Elser, {\em J. Phys. A} {\bf 17}, 1509 (1984).

\bibitem{Wannier} G.H. Wannier, {\em Phys. Rev.} {\bf 79}, 357 (1950);
{\em Phys. Rev.} B {\bf 7} 5017 (E, 1973).

\bibitem{Grensing} D. Grensing, G. Grensing, {\em J. Math. Phys.}
{\bf 24}, 620 (1983); D. Grensing, I. Carlsen, H.Chr. Zapp, {\em
Phil. Mag. A} {\bf 41}, 777 (1980).

\bibitem{matheux} H. Cohn, M. Larsen, J. Propp, {\em New York J. of
Math.} {\bf 4}, 137 (1998); H. Cohn, R. Kenyon, J. Propp, {\em
J. Amer. Math. Soc.} {\bf 14}, 297 (2001).

\bibitem{Bibi97B} N. Destainville, {\em J. Phys. A: Math. Gen} {\bf
31}, 6123 (1998).

\bibitem{Widom97} M. Widom, N. Destainville, R. Mosseri, F. Bailly,
{\em in} {\em Proceedings of the 6th International Conference on
Quasicrystals} (World Scientific, Singapore, 1997).

\bibitem{Bibi00} N. Destainville, M. Widom, R. Mosseri, F. Bailly,
{\em Mat. Sci. Eng.~A} {\bf 294-296}, 409 (2000).

\bibitem{4to3} M. Widom, R. Mosseri, N. Destainville, F. Bailly,
{\em J. Stat. Phys.} {\bf 109}, 945 (2002).

\bibitem{octo01} N. Destainville, R. Mosseri, F. Bailly, 
{\em J. Stat. Phys.} {\bf 102}, 147 (2001).

\bibitem{Mosseri93B} R. Mosseri, F. Bailly, {\em
Int. J. Mod. Phys. B}, Vol 7, {\bf 6}\&{\bf 7}, 1427 (1993).

\bibitem{Bibi97} N. Destainville, R. Mosseri, F. Bailly,
{\em J. Stat. Phys.} {\bf 87}, 697 (1997).


\bibitem{Knuth92} D.M. Knuth, Axioms and Hulls, {\em in} {\em Lect. Notes 
in Computer Sci.} {\bf 606}, 35 (1992).

\bibitem{Bjorner93} {\em Oriented Matroids}, 
A. Bj\"orner, M. Las Vergnas, B. Sturmfels, N. White,
G.M. Ziegler (Cambridge University Press, 1993).

\bibitem{Kenyon93} R. Kenyon, {\em Algorithmica} {\bf 9}, 382 (1993).

\bibitem{Elnitsky97} S. Elnitsky, {\em J. Combinatorial
  Theory} A {\bf 77}, 193--221 (1997).

\bibitem{Bailey97} {\em Tilings of zonotopes: Discriminental
arrangements, oriented matroids, and enumeration}, G.D. Bailey, Ph. D. Thesis
(Univ. of Minnesota, 1997).

\bibitem{lpw} W. Li, H. Park and M. Widom,
{\em J. Stat. Phys.} {\bf 66}, 1 (1992).

\bibitem{Newman} {\em Monte Carlo methods in statistical mechanics},
M.E.J. Newman, G.T. Barkema (Clarendon Press, Oxford, 1999).

\bibitem{mc} K.J. Strandburg, L.-H. Tang and M.V. Jaric, {\em
Phys. Rev. Lett.} {\bf 63}, 314 (1989); L.J. Shaw, V. Elser,
C.L. Henley, {\em Phys. Rev. B} {\bf 43}, 3423 (1989); L.-H. Tang,
{\em Phys. Rev. Lett.} {\bf 64}, 2390 (1990); M. Oxborrow and
C.L. Henley, {\em Phys. Rev. B} {\bf 48}, 6966 (1993); F. Gahler,
Proc. ICQ5 236 (1995); D. Joseph and M. Baake, {\em J. Phys. A} {\bf
29}, 6709 (1996).

\bibitem{Mixing} M. Luby, D. Randall and A. Sinclair, {\em SIAM J. of
Comp.}  {\bf 31}, 167 (2001); N. Destainville, {\em Phys. Rev. Lett.}
{\bf 88}, 30601 (2002).

\bibitem{fan} 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~\ref{DeBruijn4.2}, even though it might
not be clear because of resolution.

\bibitem{rectkite} A rectangular kite has two right angles between its
unequal sides.



\end{thebibliography}

\newpage

\begin{table}[ht]
\caption{Exact (except $B_{11}/B_{10}$) data for fixed boundary
tilings, $p=1$ and $D=1, 2, \ldots, 10$.}
\vspace{0.3cm}
\begin{center}
\begin{tabular}{|l|r|r|r|r|r|}
$D$ & 1 & 2 & 3 & 4 & 5 \\
\hline
$B_D(p=1)$ & 1 & 1 & 2 & 8 & 62 \\
$\bar{\sigma}_{D}(p=1)$ & & 0 & 0.231 & 0.347 & 0.413 \\
$B_{D+1}/B_D$ & 1 & 2 & 4 & 7.75 & 14.65 \\
\hline
\hline
$D$ & 6 & 7 & 8 & 9 & 10 \\
\hline
$B_D(p=1)$ & 908 & 24,698 & 1,232,944 & 112,018,190 & 18,410,581,880 \\
$\bar{\sigma}_{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 \pm 0.04$\\
\end{tabular}
\end{center}
\label{pavages.finis}
\end{table}

\begin{table}[ht]
\caption{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~\ref{MCalgo};
[\cite{Newman}]).
}
\vspace{0.3cm}
\begin{center}
\begin{tabular}{|l|r|r|r|r|r|r|}
$N_{MC}$ & $10^3$ & $10^4$ & $10^5$ & $10^6$ & $10^7$ & exact \\
\hline
$\bar{P}_9 $ & 159.264& 164.176& 164.411 & 164.371 & 164.344 & 164.35 \\
$\log \bar{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}$ & \\
\end{tabular}
\end{center}
\label{test.accuracy}
\end{table}

\begin{table}[ht]
\caption{Exact data for fixed boundary tilings, $p=2$ and 
$D=1, 2, \ldots, 6$. }
\vspace{0.3cm}
\begin{center}
\begin{tabular}{|l|r|r|r|r|r|r|}
$D$ & 1 & 2 & 3 & 4 & 5 & 6 \\
\hline
$B_D(p=2)$ & 1 & 1 & 20 & 5,383 & 16,832,230 & 327,674,047,931 \\
$\bar{\sigma}_{D}(p=2)$ & & 0 & 0.250 & 0.358 & 0.416  & 0.451 \\
\end{tabular}
\end{center}
\label{pavages.finis2}
\end{table}

\begin{table}[ht]
\caption{The first values of the limiting distribution of coordination
numbers $Z$, obtained both in the mean-field approximation and
numerically, by Monte Carlo simulations.}
\begin{center}
\begin{tabular}{|l|r|r|r|r|r|}
$Z$ & 3 & 4 & 5 & 6 & 7 \\
\hline
Mean-field & 0.33 & 0.41 & 0.20 & 0.05 & 0.009\\
\hline 
Numerical & 0.31 & 0.43 & 0.21 & 0.04 & 0.005 \\
\end{tabular} 
\end{center}
\label{Z_num}
\end{table}

\begin{table}[ht]
\caption{Path count statistics for strained tilings.}
\label{tab:strain}
\begin{center}
\begin{tabular}{|l|r|r|r|r|}
$D$ & $p=1,1,...$ & $p=1,1,...,2,2,...$ & $p=1,1,...,3,3,...$ 
& $p=1,2,1,2,...$ \\
\hline
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 \\
\end{tabular}
\end{center}
\end{table}

\end{document}






