\documentstyle[aps,psfig]{revtex}

\def\sinf{\sigma_{\infty}}
\def\sinfbar{\bar{\sigma}_{\infty}}
\def\eh{\bf \hat{e}}
\def\be{\begin{equation}}
\def\ee{\end{equation}}
\def\ra{\rightarrow}

\begin{document}

\title{Entropy of three-dimensional integer partitions \\
and related three-dimensional rhombus tilings}

\author{M. Widom$^1$, R. Mosseri$^2$, N. Destainville$^3$}
\address{
$^1$Carnegie Mellon U. Department of Physics, Pittsburgh, PA  15213, USA\\
$^2$GPS, Universit\'e Paris 7 et 6, 2 place Jussieu, 
75251 Paris Cedex 05, France\\
$^3$Laboratoire de Physique Quantique, Universit\'e Paul Sabatier,
31062 Toulouse Cedex 04, France.}

\maketitle


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

\section{Introduction}
\label{intro}

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

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

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

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

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

\section{Partitions, tilings and boundary conditions}
\label{ptbc}

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

Solid partitions are defined as follows~\cite{destain97}: let us
consider a three-dimensional array of sides $k \times l \times m$.
Fix an integer $p>0$, called the {\em height} of the partition
problem. Put non-negative integers in that array, smaller than $p$,
with the following constraints: these integers are decreasing in each
of the 3 directions of space. Each integer partition is equivalent to
one tiling of the rhombic dodecahedron ``RD'' (outer frame in
Fig.~\ref{octa}) by rhombahedra. Because of the fixed boundaries we
call the entropy of this tiling $\sigma_{fixed}$.

\begin{figure}[tb]
\begin{center}
\ \psfig{figure=octarctic.eps,width=7.5cm} \
\end{center}
\label{octa}
\caption{The puckered octahedral boundary conditions. The octahedron 
is inscribed in the rhombic dodecahedron coming from the solid partition
method.}
\end{figure}

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

\begin{figure}[tb]
\begin{center}
\ \psfig{figure=bords1.eps,width=5cm} \

\ \psfig{figure=bords2.eps,width=5cm} \
\end{center}
\label{bords}
\caption{Boundary conditions for the partition problem associated 
with the strain free tiling problem with boundary conditions
of Fig.~\protect{\ref{octa}}. Two views of the truncated cube
are provided: two ``corners'' have been cut. The marked vertex
is the same in both views. On the other hand, we
have indicated the minimal and maximal values of parts.}
\end{figure}


\section{Monte Carlo}
\label{mc}

Conventional Monte Carlo simulations using Metropolis
sampling~\cite{qcmc} are useful for generating typical configurations.
They can reveal temperature driven phase transitions~\cite{PT} and
yield quantitative evaluations of phason elastic constants from
fluctuations. Evaluation of entropy by Monte Carlo simulation demands
specialized techniques.  Advanced methods using specialized dynamics
based on entropic sampling~\cite{entropic} and transfer
matrices~\cite{qctmmc} yield reasonably accurate free energies and
absolute entropies.

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

\subsection{Transition Matrix}
\label{tm}

For any legal partition $P=\{n_{ijk}\}$ we define its ``energy'' as
its total height
\be
\label{Ham}
E(P)=\sum_{ijk} n_{ijk}
\ee
and incorporate the constraint of legality by defining $E=\infty$ for
any partition $P$ that violates the ordering conditions inside the
partition or at the boundaries. Under single vertex flip dynamics a
single partition element increases or decreases by one unit resulting
in an energy change $\Delta E = \pm 1$.  The gound state of this model
is the lowest height legal partition. For the boundary conditions
employed here, the ground state is unique and has energy $E_{min}$.

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

For a given partition $P$, a number of partitions $n_{\pm}(P)$ can be
reached by single upwards or downwards steps. Additionally, a certain
number $n_0(P)$ of steps are forbidden due to the partition constraints.
The sum rule
\be
\label{sumrulen}
n_{-}(P)+n_{0}(P)+n_{+}(P)=N
\ee
holds for every partition $P$. In a single Monte Carlo step, the
probability for upwards or downwards transitions at infinite
temperature is $n_{\pm}(P)/N$, and the probabilty a move will be
rejected is $n_0(P)/N$. At finite temperatures the transition
probabilities must be multiplied by appropriate Boltzmann factors.

During a Monte Carlo simulation, many different partitions with a
given energy may be sampled. We can calculate exactly the infinite
temperature acceptance ratios of transitions by averaging over states
with energy $E$,
\be
\label{rates}
\omega_{\pm,0}(E)={{1}\over{W(E)}} \sum_{P(E)} n_{\pm,0}(P)/N
\ee
where the sum is over all partitions of energy $E$ and the
normalization $W(E)$ is the total number of partitions with energy
$E$. In general the set $P(E)$ and the value $W(E)$ are not
known. However, we may calculate the rates approximately by averaging
$n_{\pm,0}(P)/N$ over those partitions $P$ occuring during a Monte
Carlo simulation.  The acceptance ratios obey the sum rule
\be
\label{sumrulew}
\omega_{-}(E)+\omega_{0}(E)+\omega_{+}(E)=1
\ee
for any energy $E$.

Consider the infinite temperature detailed balance condition.  At
infinite temperature the probability that the energy of partition has
value $E$ is $W(E)/W$, where $W$ is the total number of legal
partitions. Multiplying this probability by the
acceptance rate of upwards transitions $\omega_{+}(E)$ yields the
total forwards transition rate. The backwards rate is obtained in
similar fashion. Detailed balance requires that the total rate of
forward transitions from energy $E$ to energy $E+1$ must equal the
total rate of backwards transitions, hence
\be
\label{balance}
\omega_{+}(E) W(E)/W=\omega_{-}(E+1) W(E+1)/W.
\ee
It is useful to rearrange the detailed balance eq.~(\ref{balance}) to
find
\be
\label{iterate}
W(E+1)={{\omega_{-}(E+1)}\over{\omega_{+}(E)}} W(E)
\ee
which allows us to extract the full density of states function $W(E)$
using uniqueness of the ground state, $W(E_{min})=1$.

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

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

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

\subsection{Numerical Data}
\label{data}

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

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

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

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

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

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

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

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

\vspace{2cm}

\begin{thebibliography}{99}

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

\bibitem{elser84} {\sc V. Elser}.~-- Solution of the Dimer Problem on
an Hexagonal Lattice with Boundary, {\em J. Phys. A: Math. Gen.} {\bf
17}, 1509 (1984).

\bibitem{proppandco} H. Cohn, M. Larsen, J. Propp, The Shape of a Typical
Boxed Plane Partition, {\em New York J. of Math.} {\bf 4}, 137 (1998);
H. Cohn, R. Kenyon, J. Propp, A Variational Principle for Domino
Tilings, {\em Journal of the AMS} {\bf 14}, 297 (2001).

\bibitem{destain98} N. Destainville, Entropy and boundary conditions
in random rhombus tilings, {\em J. Phys. A: Math. Gen} {\bf 31}, 6123
(1998).

\bibitem{destain97} N. Destainville, R. Mosseri, F. Bailly,
Configurational Entropy of Codimension-One Tilings and Directed Membranes, 
{\em J. Stat. Phys.} {\bf 87}, Nos 3/4, 697 (1997).

\bibitem{destain01} N. Destainville, R. Mosseri, F. Bailly,
Fixed-boundary octagonal random tilings: a combinatorial 
approach, {\em J. Stat. Phys.} {\bf 102}, 147 (2001).

\bibitem{TMMC} J.S. Wang, T.K. Tay and R.H. Swendsen, Phys. Rev. Lett. 
{\bf 82}, 476 (1999); R.H. Swendsen, B. Diggs, J.-S. Wang, S.-T. Li
and J.B. Kadane, Int. J. Mod. Phys. C{\bf 10}, 1563 (1999);

\bibitem{exact} M. Widom {\em Phys. Rev. Lett.} {\bf 70}, 2094 (1993);
P. Kalugin {\em J. Phys. A} {\bf 27}, 3599 (1994);
J. de Gier and B. Nienhuis, {\em Phys. Rev. Lett.} {\bf 76}, 2918 (1996);
J. de Gier and B. Nienhuis, {\em J. Phys. A} {\bf 31}, 2141 (1998);

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

\bibitem{qctm} M. Widom, D.P. Deng and C.L. Henley
{\em Phys. Rev. Lett.} {\bf 63}, 310 (1989); W. Li, H. Park and
M. Widom, {\em J. Stat. Phys.} {\bf 66}, 1 (1992); M.E.J. Newman and
C.L. Henley, {\em Phys. Rev. B} {\bf 52}, 6386 (1995)

\bibitem{qctmmc} 
L.J. Shaw and C.L. Henley, {\em J. Phys. A} {\bf 24}, 4129 (1991);

\bibitem{entropic}
F. Gahler, Proc. ICQ5, 236 (1995); 

\bibitem{PT}
H.-C. Jeong and P.J. Steinhardt, Phys. Rev. B{\bf 48}, 9394 (1993);


\newpage

\end{thebibliography}


\begin{table}[htb]
\caption{Convergence of entropy data for increasing run length. The
residual $R=\log{W(E_{max})}$ measures the cumulative error in $W(E)$.}
\label{converge}
\begin{tabular}{||l|r|r|r||}
$N$              &   $10^4$    &  $ 10^5$  &  $10^6$  \\
$\sigma_{fixed}$ &   0.147385  &  0.147400 &  0.147349 \\
$R_{rixed}$      &   0.180     &  0.429    &  0.002    \\
$\sigma_{free}$  &   0.197300  &  0.197234 &  0.197273 \\
$R_{free}$       &   0.143     & -0.202    & -0.036  \\
\end{tabular}
\end{table}

\begin{table}[htb]
\caption{Size-dependent entropies. Values in parentheses are uncertainties
in final digit. Values without uncertainties are exact.}
\label{sigmas}
\begin{tabular}{||r|l|l|l||}
$L$   & $\sigma_{free}$ & $\sigma_{fixed}$ & $\sigma_{free}/\sigma_{fixed}$ \\
\hline
1     & 0.1732868       &  0.1732868       & 1.000 \\
2     & 0.1732868       &  0.1601239       & 1.080 \\
3     & 0.17947(2)      &  0.1545769       & 1.161 \\
4     & 0.18455(6)      &  0.1516677       & 1.217 \\
5     & 0.18829(3)      &  0.1502(1)       & 1.254 \\
6     & 0.19108(4)      &  0.1492(1)       & 1.281 \\
7     & 0.193(1)        &  0.1486(1)       & 1.299 \\
8     & 0.19486(2)      &  0.147997(3)     & 1.317 \\
9     & 0.196183(1)     &  0.14762(2)      & 1.329 \\
10    & 0.19727(4)      &  0.14735(5)      & 1.339 \\
\hline
$\infty$& 0.215(5)      &  0.145(1)        & 1.48(3) \\
\end{tabular}
\end{table}

\end{document}