text:
# Mechanistic approach to generalized technical analysis of share prices and stock market indices

## I Introduction

Technical indicators as moving average and $`momentum`$ are part of the classical technical analysis and are much used in efforts to predict market movements. One question is whether these techniques provide adequate ways to read the trends, and later on allow for an investment strategy development. It has been shown that moving average trading rules can be utilized  for USA and UK markets. In both futures and spot foreign currency markets significant profits can be earned along these lines. Parisi and Vasquez recently used the moving average technique on emerging markets  to show that buy signals consistently generate higher returns than sell signals. Moreover, returns following sell signals are shown to be negative, which is not easily explained by any of the currently existing equilibrium models. Related studies by Gunasekarage and Power  also showed that technical trading rules have predictive ability in South Asian stock markets.

One surprise to a physicist is the neglect of the volume of transactions in the classical way of predicting the evolution of a share price or a market index. Yet there was a sort of thermodynamic relationship between share price and exchanged volume of shares in estimating the value of some company from an investor point of view. Such an addition to measure and if possible predict the evolution of stocks is introduced here. Indeed we present a generalization of the classical technical analysis concepts taking into account the volume of transactions. The latter becomes a time dependent or generalized mass. The notions of pressure, acceleration and force are deduced. In that spirit, a generalized (kinetic) energy is easily defined. It is pointed out that the momentum correlations take into account the sign of the fluctuations, while the energy is geared toward the absolute value of the fluctuations. They have different patterns which are checked by searching for the crossing points of their respective moving averages. The evolution of IBM share price and volume between Jan 01, 1990 and Dec 31, 2000 is used for illustrations.

## II Technical Analysis

### A Moving Average

Consider a time series $`x(t)`$ given at N discrete times $`t`$. The series (or signal) moving average $`M_\tau (t)`$ over a time interval $`\tau `$ is defined as

$$M_\tau (t)=\frac{1}{\tau }\underset{i=t}{\overset{t+\tau 1}{}}x(i\tau )t=\tau +1,\mathrm{},N$$
(1)
i.e. the average of $`x`$ over the last $`\tau `$ data points. For simplicity we suppose that the ticking times are equally spaced. One can easily show that if the signal $`x(t)`$ increases (decreases) with time, $`M_\tau (t)<x(t)`$ ($`M_\tau (t)>x(t)`$). Thus, a moving average captures the past trend of the signal over a given time interval $`\tau `$. The IBM daily closing value price signal between Jan 01, 1990 and Dec 31, 2000 is shown in Fig. 1 (top figure) together with the $`\tau =50`$ day moving average taken from Yahoo. The bottom figure shows the daily $`volume`$ of transactions given in millions.

The moving average notion has already been discussed. It is obvious that like any other statistical $`mean`$, a moving average $`M_\tau (t)`$, depends on the number of data points taken into account. There can be as many moving averages as $`\tau `$ intervals. The shorter the $`\tau `$ interval the closer to a signal is the moving average. However, a too short moving average may give false messages about the long time trend of the signal. In Fig. 2 two moving averages of the IBM signal for $`\tau `$=5 days (i.e. 1 week) and 21 days (i.e. 1 month) are compared for illustration.

The intersections of the $`price`$ signal with a moving average can define so-called lines of $`support`$ or $`resistance`$. A support (resistance) line occurs when the local minimum (maximum) of $`x(t)`$ bounces on $`M_\tau (t)`$. The lines are supposed to indicate the price level at which $`\tau `$-investors believe that prices will move higher or lower respectively. In Fig. 2, IBM lines of resistance happen e.g. around June 1993 and lines of support around mid Oct. 1993 for the $`\tau `$ = 5 days and 21 days cases respectively. Support and resistance levels depend on $`\tau `$ and are based in principle on investment horizon strategy, but are in fact containing much psychological fancy.

Other features of the moving average prone investor framework are the intersections between $`two`$ moving averages $`M_{\tau _1}`$ and $`M_{\tau _2}`$. They might occur or not at drastic changes in the trends of $`x(t)`$. Consider again the two moving averages of IBM price signal for $`\tau _1=5`$ days and $`\tau _2=21`$ days (Fig. 2). If $`x(t)`$ increases for a long period of time before decreasing rapidly, $`M_{\tau _1}`$ will cross $`M_{\tau _2}`$ from above. This event is called a death cross in empirical finance. In contrast, when $`M_{\tau _1}`$ crosses $`M_{\tau _2}`$ from below, the crossing point is a gold cross. They appear more drastic when the respective slopes have different signs. In Fig.2 a death cross and a gold cross occur near March 93 and Oct. 93 respectively. The density of such crossing points between two moving averages as a function of the difference in the characteristic $`\tau `$’s defining the moving averages has been discussed elsewhere and is shown in Fig.3 for the case of interest here. Based on this idea, a new and efficient approach has been suggested in order to estimate an exponent that characterizes the roughness of a signal. From Fig.3 and Ref. the IBM roughness exponent has been found equal to $`0.44\pm 0.02`$ for the time interval considered.

### B Momentum Indicator

The so called momentum is another instrument of the technical analysts. We will refer to it here as the classical momentum for reasons to become obvious later. The classical momentum of a stock is defined over a time interval $`\tau `$ as

$$R_\tau (t)=\frac{x(t)x(t\tau )}{\tau }=\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}t=\tau +1,\mathrm{},N$$
(2)
For $`\mathrm{\Delta }t`$ =$`\tau `$ = 1, the momentum is nothing else than the volatility. The momentum $`R_\tau `$ for three time intervals, $`\tau `$=5, 21 and 250 days, i.e. one week, one month and one year, are shown in Fig. 4 for IBM between 1990 and 2000, together with a blow up for the years 1999-2000. The longer the period, the smoother the momentum signal. Relevant information on the price trend turnovers is usually considered to be found in a few moving averages of the $`momentum`$, or momentum indicator, i.e. in

$$R_\tau ^\mathrm{\Sigma }(t)=\underset{i=t}{\overset{t+\tau 1}{}}\frac{x(i)x(i\tau )}{\tau }t=\tau +1,\mathrm{},N$$
(3)
like those for $`\tau `$ = 1 week, 1 month and 1 year. Notice that those are calculated over the same time intervals over which the momentum is calculated. For IBM, these $`R_\tau ^\mathrm{\Sigma }`$ are shown in Fig. 5. The density of intersections between moving averages of these momenta could be calculated. The result is displayed in Fig.3.

### C Classical Strategy

In Fig. 6 the IBM signal and its weekly (short-term), monthly (medium-term) and yearly (long-term) moving averages are compared to the weekly (short-term), monthly (medium-term) and yearly (long-term) momentum indicators to show the 1999 bullish and also beginning of bearish trends. The strategic message that is coming out of reading the combination of these six indicators (Fig. 6) is that one could start buying at some momentum bottom and sell at a maximum. A buy position is found for both monthly and weekly momentum indicators around (1) Feb. 99. Another buy position occurs around mid-April 99 when the price surge confirms the momentum uptrend, at a gold cross (2). Selling signals are given at the maximum of the monthly momentum indicator in mid-Jan. 99 (6), near the second half of mid-May (7) and mid-July (8). A death cross (3) occurs between the short and medium term moving averages in Sept. 99 and is subsequently followed by the maximum (4) of the monthly momentum, even turning negative. In Oct. 99, occurs the maximum of the long-term momentum, just before the Oct. 99 crash (5). A good technical analyst would have strongly recommended to sell the position since the price is also falling down below the weekly moving average. Hope (or faith?) for prospect reoccurs after Nov. 99 during one month.

## III Generalized technical analysis

Stock markets do have another component beside prices or volatilities. This is the volume of transactions. It is introduced here as the physical mass of stocks. Remember that the number of shares is constant over rather long time intervals, i.e. usually between splits, like the mass of an object.

Consider $`V(t)`$ to be the volume of transactions of a stock with price $`x(t)`$ at time $`t`$, - Fig. 1 (bottom). A generalized momentum $`\stackrel{~}{R}_\tau `$ over a time interval $`\tau `$ can be defined (see Fig. 7) as in physics through

$$\stackrel{~}{R}_\tau (t)=\frac{V(t)}{<V(t)>_\tau \tau }\frac{x(t)x(t\tau )}{\tau }=m(t)\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t},t=\tau +1,\mathrm{},N$$
(4)
where the total volume of transactions over the interval $`\tau `$ is $`<V>_\tau \tau =_{i=1}^\tau V(i)`$. In so doing, we introduce some financial analogy to a generalized time dependent mass $`m(t)`$ of a diffusing object. The total volume in the denominator is introduced for a normalization purpose. Notice that we could choose other normalizations or definitions of the mass or the generalized momentum : one could take the $`log(V(t))`$, but this introduces a nonlinear transformation. One could use $`log[V(t+\tau )/V(t)]`$ but that can be negative. One could normalize with respect to $`<V>_\tau `$ only, but $`\stackrel{~}{R}_\tau (t)`$ values would be larger, ….

In order to search for more definite indications on the stock trend changes, represented as the complex influence between price and volume of transactions, we further consider a moving average of the generalized momentum, i.e. a generalized momentum indicator (Fig.8) which is

$$\stackrel{~}{R}_\tau ^\mathrm{\Sigma }(t)=\underset{i=t}{\overset{t+\tau 1}{}}\frac{V(t)}{<V>_\tau \tau }\frac{x(i)x(i\tau )}{\tau }t=\tau +1,\mathrm{},N$$
(5)
The IBM generalized momentum indicator $`\stackrel{~}{R}_\tau `$ for three time intervals $`\tau `$, i.e. 1 week, 1 month and 1 year and their corresponding moving average $`\stackrel{~}{R}_\tau ^\mathrm{\Sigma }`$ are shown in Fig. 8. It is observed that the multiplicative factor results only in enhancing features, i.e. the mass being here positive. However the slopes in the generalized momentum indicator are much enhanced due to volume variations. A comparison of moving averages of the IBM signal and its IBM generalized momentum for three different horizons together with the volume of transactions are shown in Fig. 9. The volume of transactions is plotted in tenths of millions on the generalized momentum base line. The long-term (yearly) generalized momentum signal is marked by dot-dash curve and has small positive values, with a steady increase up to the maximum of the IBM price around mid Sept. 1999 marked by the index (3). See the well marked crash of Oct. 99; also displayed in Fig.9, just after (5). The crash coincides with the maximum in the volume of transactions, and a negative momentum.

A closer look of the short (weekly) and medium-term (monthly) momentum in terms of sell/buy message point of views again suggests to buy around Feb. 1999 (see Fig. 6 as well) at the minimum of the momentum curve. The new feature is that the peaks of the generalized medium term (monthly) and long-term (yearly) momentum enhances several changes in the price trend represented by its monthly or yearly moving average. Peaks are observed at mid Jan. 99 (6), mid May (7) (after the gold cross (2) ) and mid July 1999 (8). Notice that (3) and (8) coincide with the monthly momentum turning to negative values after reaching a weak maximum. Notice the differences between Fig. 6 and Fig. 9, i.e. the enhanced, and even modified structures, between (6) and (1), and between (4) and the Oct. crash, due to the volume effect.

### A Pressure, Acceleration and Force Indicators

Introducing a ”mass” $`m(t)=V(t)/(<V>_\tau \tau )`$ and the ”velocity” $`v(t)=\mathrm{\Delta }x/\tau `$ through the standard notion of physical momentum

$$\stackrel{~}{R}_\tau ^\mathrm{\Sigma }=p(t)=\underset{i=t}{\overset{t+\tau 1}{}}m(i)v(i)$$
(6)
enables us to take into account both sides of the market coin and their impact on one another.

Note that $`<V>_\tau \tau `$ ensures to include a $`pressure`$*** A global kinetic theory for prices has been derived considering an equilibrium market (with actors having all identical relaxation times). In closing the set of equations, an equation of state, with a pressure and a temperature, between the price, as the order parameter of a stock, and the volume of exchanged shares were introduced. contribution from the variation of the accumulated volume of transactions during the time period $`\tau `$ on the change of the momentum of the stock. From a mathematical point of view at the time of the peaks of generalized momentum curve $`\stackrel{~}{R}_\tau ^\mathrm{\Sigma }`$, the first derivative of the momentum are obviously equal to zero and the second derivative is negative or positive depending whether there is a peak or a dip. Because of the time dependent ”mass” this implies that

$$p^{}(t)=\underset{i=t}{\overset{t+\tau 1}{}}m^{}(i)v(i)+\underset{i=t}{\overset{t+\tau 1}{}}m(i)v^{}(i)=0$$
(7)
or

$$p^{}(t)=\underset{i=t}{\overset{t+\tau 1}{}}m^{}(i)v(i)+\underset{i=t}{\overset{t+\tau 1}{}}m(i)a(i)=0,$$
(8)
where $`a(i)`$ can be thought as an $`acceleration`$.Price velocity and price acceleration are two fundamental indicators which of course already exist in the economy literature and econophysics, e.g. they were recently used to construct a general classification of market indices possible patterns deviating from the random walk. Therefore the second term in Eq. 8 represents a classical $`force`$ that acts upon the object (stock or market index). It originates from the speed of change of stock price whose sign is either positive or negative. At times of extrema this force balances the first term in Eq. 8, determining the rate in the volume of transactions. The sign of this first term depends on the accumulated sum of $`m^{}(i)v(i)`$, since both $`m^{}(i)`$ and $`v(i)`$ can be negative or positive.

Therefore, the fact that the first derivative of the momentum $`p^{}(t)`$ vanishes at time $`t_1`$ can lead to a deeper reading of the price/volume interactions. The maximum of the (generalized or not) momentum indicator at $`t_1`$, mid Jan. 1999, mark (6) in Fig. 6 and Fig.9, is related to a sharp increase of the volume of transactions, as the rescaled volume of transactions is of the order of $`\stackrel{~}{R}_\tau ^\mathrm{\Sigma }`$. For Eq. 8 to be satisfied, the second term in the equation should be negative. Thus the derivative of the $`velocity`$, i.e. the acceleration $`\mathrm{\Delta }v(t)/\mathrm{\Delta }t`$ should be negative, which means a change in the derivative of the price trend. This is easily seen at the resistance point, that also coincides with a death cross between the IBM signal and its moving averages for the one week and one month time averaging at (1). The same discussion pertains to the maximum of the momentum curve in mid July 1999, mark (8) in Fig. 6 and Fig. 9, when the IBM share price breaks the resistance line, just before a death cross when the price drops below its monthly moving average, while the weekly moving average becomes negative.

However, the cause of the peak in the generalized momentum indicator at mid-May, 1999, mark (7) in Fig. 6 and Fig. 9, is more complex than for mid-Jan. 99, mark (6), or mid-July 99, mark (8). Indeed the acceleration during most of the two prior months is positive and so is the second term in Eq. 8. At the time of the price maximum on May 13, 1999 the classical force has to be (and is) balanced by a (negative) contribution from the first term. Because the velocity $`v(i)`$ is positive during that period, the derivative of the volume of transactions (the mass) has to be (and is) most of the time negative.

Thus the increase of momentum is not entirely due to the price increase, but does have another component, hidden if one does not consider the above generalization. After May 13, 1999, the price drops to the monthly moving average and rebounds on a support line. Thus this can be interpreted as a continuous price increase with corrections to the moving average price due to the influence of the volume of transactions. The observed evolution of the generalized momentum of the IBM stock implies that some generalized force can be considered as a cause of these changes.

### B Energy

Mechanically speaking it can be thought that some $`energy`$ is also accumulated through the interplay between the price and the volume of transactions and causes a generalized force to act. Therefore, a kinetic ”energy” can be introduced as

$$E_\tau ^\mathrm{\Sigma }(t)=\underset{i=t}{\overset{t+\tau 1}{}}\left(\frac{x(i)x(i\tau )}{\tau }\right)^2t=\tau +1,\mathrm{},N.$$
(9)
The generalized mass introduced above leads also to generalize the ”energy” of the stock signal to be like in mechanics (we drop the factor $`1/2`$)

$$\stackrel{~}{E}_\tau ^\mathrm{\Sigma }(t)=\frac{1}{\tau }\underset{i=t}{\overset{t+\tau 1}{}}\frac{V(t)}{<V>_\tau }\left(\frac{x(i)x(i\tau )}{\tau }\right)^2t=\tau +1,\mathrm{},N$$
(10)
or

$$\stackrel{~}{E}_\tau ^\mathrm{\Sigma }(t)=m(t)\left(\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}\right)^2=mv^2.$$
(11)
Theoretical and generalized $`energy`$ of the IBM stock are shown in Figs. 10-11 for the same three usual time interval averages $`\tau `$, i.e. one week, one month, one year. Note (insert of Fig. 10) two large peaks in the monthly theoretical energy : the high one in mid-May 1999 corresponding to (7) in Fig. 9, and one at the crash time in October. A small peak in mid-July, mark (8) in Fig. 9, corresponds to the maximum of the price and the generalized momentum. The appearance of a huge increase in energy in Oct. 99 together with the negative slope of the generalized momentum indicator should be interpreted now as a serious warning for an incoming crash. A mimimum in kinetic energy can be seen to provide information on strong continuing price increase or decrease depending on the momentum rate of change. The local maximum of a kinetic energy as in classical physics indicates some accumulation of energy at that time for the evolution of the stock, - accumulation which must be dissipated ! It is interesting to note that the monthly generalized energy (Fig. 11, inset) emphasizes even further the huge accumulation of energy in Oct. 1999.

## IV Conclusions

A thermodynamic-like phase diagram between the daily closing price ($`P`$) and the daily transaction volume ($`V`$) resulted in a fundamental phase diagram for companies quoted on stock exchanges. This pointed out to search for correlations between volume and price. In the same spirit of econophysics connecting physics and financial data for share prices or stock market indices we have developed a generalization of the classical method by technical analysts, i.e starting from the notion of moving averages and momentum indicators. One can take into account not only the price of a stock but also the volume of transactions which is similar to a time dependent generalized mass in mechanics. This might interestingly serve for generalizing ideas on anomalous diffusion model(s) for price evolutions as well. A notion of pressure, acceleration and force concepts are deduced and identified to their classical mechanics counterparts. The terms have sometimes been used in the literature but outside their classical meaning. A generalized (kinetic) energy is also easily defined, in the same lines of an analogy. A correspondence between mechanical terms, classical technical analysis terms and to be seen generalized concepts of technical analysis is found in Table I. This might help in devising a free energy like approach, as also introduced elsewhere.

Introducing the product of price and mass of transactions in indicators renders the new signals very rough. The number of intersections, mimima and maxima is highly modified. Whether the above generalization of technical analysis concepts can be used to increase the number of buy/sell orders, and develop new strategies through nine rather than six indicators (shown as a summary in Figs. 12-13) is left for further work. The interest of such considerations will be obvious if cumulated (good !) predictions from the above stand over usual findings.

Finally it seems that the (generalized) momentum and energy concepts so introduced concern other effects studied in economy, i.e. the correlation existence of signs and amplitudes of share price variations. It is easily understood that the momentum (and momentum correlations) take into account the sign of the fluctuations (and their correlations), while the energy (and energy correlations) is geared toward the absolute value of the fluctuations (and their correlations). In thermodynamics of non-equilibrium processes, these correlations are described by the $`viscosity`$ and thermal conductivity coefficient respectively. Thus the above concepts might also serve in a dynamic equation framework.

Acknowledgements

MA thanks Etienne Labie and Denis Vanderborght (Leleux Ass.) for practical comments and references.


format: text

text:
# There might be superluminal particles in nature The paper has been published in Journal of Shaanxi Normal University (Natural Science Edition),Vol.29,No.3

## 1 INTRODUCTION

Nothing can travel faster than light. Is this statement as true now as it ever was? In 2000 there were two experiments showing superluminal propagation of microwave  or laser pulse . Yet physicists believe that the law of physics have remained intact . However, the experimental discovery of negative mass square of neutrino in recent years ,

$$E^2=c^2p^2+m^2c^4,$$
$`(1)`$
$$m^2(\nu _e)=2.5\pm 3.3eV^2,$$
$`(2)`$
though far from accurate, does strongly hint that neutrino might be a particle moving faster than light. Actually, rewriting (1) as

$$E^2=c^2p^2m_s^2c^4$$
$`(3)`$
and using the quantum relations $`E=\mathrm{}\omega `$ and $`p=\mathrm{}k`$, one easily derives the kinematic relation for superluminal particle ( also named as tachyon in literature ) as follows ($`m_s`$ is called ”proper mass”):

$$p=\stackrel{~}{m}u=\frac{m_su}{\sqrt{\frac{u^2}{c^2}1}},E=\stackrel{~}{m}c^2=\frac{m_sc^2}{\sqrt{\frac{u^2}{c^2}1}}.$$
$`(4)`$
Here the velocity of particle $`u`$ is identified with the group velocity $`u_g=\frac{d\omega }{dk}`$ of wave. To derive (3), a quantum equation for neutrino was proposed in Ref.  (see also ). In this paper, we will compare three kinds of equation—Dirac equation and two ”superluminal” equations and discuss relevant problems.

## 2 Dirac equation versus superluminal equation

The Dirac equation for fermion with rest mass $`m_0`$ reads

$$i\mathrm{}\frac{}{t}\psi =ic\mathrm{}\stackrel{}{\alpha }\psi +\beta m_0c^2\psi $$
$`(5)`$
where $`\psi `$ is a four-component spinor wave function and $`\alpha _i`$ and $`\beta `$ are 4$`\times `$4 matrixes. In Dirac representation, they are expressed as

$$\psi =\left(\begin{array}{c}\phi \\ \chi \end{array}\right),\alpha _i=\left(\begin{array}{cc}0& \sigma _i\\ \sigma _i& 0\end{array}\right),\beta =\left(\begin{array}{cc}I& 0\\ 0& I\end{array}\right)$$
$`(6)`$
where $`\phi `$ and $`\chi `$ are two-component spinors, $`\sigma _i`$ are Pauli matrixes. Now we perform a unitary transformation as

$$\psi U\psi =\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)\psi =\psi ^{}=\left(\begin{array}{c}\xi \\ \eta \end{array}\right).$$
$`(7)`$
Let the matrix $`U`$ act on Eq.(5) from the left, due to noncommutativity between $`U`$ and $`\alpha _i`$ ( or $`\beta `$ ), we find Dirac equation in ”Weyl representation” with

$$\psi =\left(\begin{array}{c}\xi \\ \eta \end{array}\right),\alpha _i^{(W)}=\left(\begin{array}{cc}\sigma _i& o\\ 0& \sigma _i\end{array}\right),\beta ^{(W)}=\left(\begin{array}{cc}0& I\\ I& 0\end{array}\right)$$
$`(8)`$
Note that, however, the representation transformation leads to important change in physical interpretation. While $`\phi `$ and $`\chi `$ in (6) represent the ”hidden particle and antiparticle fields” in a particle state, $`\xi `$ and $`\eta `$ in (8) characterize the ”hidden left-handed and right-handed rotating fields” in a particle with 100% left-handed helicity explicitly.

Now let’s perform a unitary transformation on Dirac equation (5) in Weyl representation:

$$\psi =\left(\begin{array}{c}\xi \\ \eta \end{array}\right)U_s\psi =\left(\begin{array}{cc}i& 0\\ 0& 1\end{array}\right)\psi =\psi _s=\left(\begin{array}{c}\xi _s\\ \eta _s\end{array}\right).$$
$`(9)`$
After setting

$$m_s=im_0,$$
$`(10)`$
we obtain an equation with real proper mass $`m_s`$:

$$\begin{array}{c}i\mathrm{}\frac{}{t}\xi _s=ic\mathrm{}\stackrel{}{\sigma }\xi _sm_sc^2\eta _s,\hfill \\ i\mathrm{}\frac{}{t}\eta _s=ic\mathrm{}\stackrel{}{\sigma }\eta _s+m_sc^2\xi _s.\hfill \end{array}$$
$`(11)`$
A substitution of plane-wave solution $`\xi _s\eta _sexp[i(pxEt)/\mathrm{}]`$ leads to Eq. (3) immediately. Why so simple transformation like (9) can change Dirac equation into a ”superluminal” equation? It is because the $`U_s`$ in (9) is a nonhermitian matrix. While a unitary but hermitian matrix like that in (7) merely amounts to a representation transformation, a unitary but nonhermitian matrix is capable of changing the physical essence of equation. Eq.(11) can be expressed in Dirac representation as

$$\begin{array}{c}i\mathrm{}\frac{}{t}\psi _s=ic\mathrm{}\stackrel{}{\alpha }\psi _s+\beta _sm_sc^2\psi _s\hfill \end{array}$$
$`(12)`$
$$\psi _s=\left(\begin{array}{c}\phi _s\\ \chi _s\end{array}\right)=U\left(\begin{array}{c}\xi _s\\ \eta _s\end{array}\right),\beta _s=\left(\begin{array}{cc}0& I\\ I& 0\end{array}\right).$$
$`(13)`$
Note that the matrix $`\beta _s`$ and its counterpart in Eq. (11) are antihermitian.

On the other hand, if we directly set $`m_0=im_s`$ in Eq. (5) to get the ”Dirac equation with imaginary mass” as

$$\begin{array}{c}i\mathrm{}\frac{}{t}\psi ^{(i)}=ic\mathrm{}\stackrel{}{\alpha }\psi ^{(i)}+\beta _s^{(i)}m_sc^2\psi ^{(i)}\hfill \end{array},$$
$`(14)`$
then we have

$$\psi ^{(i)}=\left(\begin{array}{c}\phi ^{(i)}\\ \chi ^{(i)}\end{array}\right),\beta _s^{(i)}=\left(\begin{array}{cc}iI& 0\\ 0& iI\end{array}\right).$$
$`(15)`$
The fact that equation (14) is wrong can easily be seen by its plane-wave solution when $`p0`$, yielding a decoupling solution:

$$\begin{array}{c}\phi ^{(i)}e^{i(im_st)}e^{m_st}\hfill \\ \chi ^{(i)}e^{i(im_st)}e^{m_st}\hfill \end{array}$$
$`(16)`$
which show the violation of unitarity. By contrast, such kind of meaningless solution is definitely excluded in Eq. (12). The sharp difference is stemming from the following fact. Three kinds of Eqs. (5), (12) and (14), all respect the basic symmetry: under the space-time inversion ($`\stackrel{}{x}\stackrel{}{x}`$, $`tt`$)

$$\phi (\stackrel{}{x},t)\chi (\stackrel{}{x},t),\chi (\stackrel{}{x},t)\phi (\stackrel{}{x},t),$$
$`(17)`$
the theory remains invariant (while a concrete solution of particle transforms into that of antiparticle). However, we should consider a smaller symmetry—the space-inversion ($`\stackrel{}{x}\stackrel{}{x}`$, $`tt`$):

$$\xi (\stackrel{}{x},t)\eta (\stackrel{}{x},t),\eta (\stackrel{}{x},t)\xi (\stackrel{}{x},t).$$
$`(18)`$
Then we see that Eqs. (5) and (14) remain invariant whereas Eq. (12) (or (11)) fails to do so. In other words, our new superluminal equation reflects the maximum parity violation, a property exactly explaining the permanent helicity of neutrino—while $`\nu __L`$ and $`\overline{\nu }__R`$ are physically realized, $`\nu __R`$ and $`\overline{\nu }__L`$ are forbidden strictly—as verified by experiments .

## 3 Virtual particle as a superluminal particle

It is no surprise to see solution like (16) when Eq. (15) contains an imaginary mass $`im_s`$ explicitly. Actually, in quantum mechanics, we often endow an imaginary part of mass $`m`$ to wave function for describing an unstable particle:

$$mmi\frac{\mathrm{\Gamma }}{2},\mathrm{\Gamma }=\mathrm{}/\tau ,\psi ^2e^{t/\tau }$$
$`(19)`$
with $`\tau `$ being the decay lifetime. Of course, the unitarity of wave function is destroyed.

Hence, it is hopeless to set $`m_0=im_s`$ directly in Dirac equation for describing a stable superluminal particle. However, we often discuss the ”virtual particle” in covariant perturbation theory in the sense that its momentum $`p`$ and energy $`E`$ can vary independently and so are not subjected to the constraint $`E^2p^2=m_0^2`$. Both $`E^2>p^2`$ and $`E^2<p^2`$ cases can occur. We might as well say that a virtual particle could be superluminal but it is unstable too. A stable superluminal particle can only be realized as $`\nu __L`$ or $`\overline{\nu }__R`$ ensured by the violation of parity symmetry.

## 4 Maximum parity violation displays the beauty of nature

Usually, a symmetry, i.e., an invariance under some transformation yields an affirmative guarantee for a theory and its violation seems to spoil the validity of the theory. Now we see that for a meaningful discrete symmetry like space-inversion, what happens in nature is either keeping intact ( for subluminal particle ) or being violated to maximum ( for superluminal neutrino ). In the latter case the violation must reach its maximum because the constraint of larger symmetry—the space-time invariance—must be held at the same time. \[In the first line of Eq. (11), the coefficient of $`m_s`$ must be fixed as $`(1)`$ to ensure the invariance of transformation (17) together with the violation of transformation (18).\] We might as well look at the maximum violation of parity as some ”antisymmetry” ( rather than ”asymmetry” ) which also provides affirmative guarantee for the validity of superluminal theory for neutrino. In fact, the normal consequence of violation of hermitian property being the violation of unitarity ( as shown by (16)) is now recast into a strange realization—of four would be unstable solutions for a same momentum, two of them ($`\nu __R`$ and $`\overline{\nu }__L`$) are eliminated whereas other two ($`\nu __L`$ and $`\overline{\nu }__R`$) are stabilized.

## 5 Summary and discussion

We can have three kinds of Dirac-type equation:

$$\begin{array}{c}i\mathrm{}\frac{}{t}\psi =ic\mathrm{}\stackrel{}{\alpha }\psi +\beta mc^2\psi \hfill \end{array}$$
$`(20)`$
as listed in the following table.

Finally, two remarks are in order:

(a) According to our present understanding, no boson but fermion can be superluminal as long as the parity symmetry is violated to maximum. Most likely, one kind of known particles, the neutrino, is just a superluminal particle, a tachyon with spin 1/2.

(b)The physical implication of nonhermitian transformation $`U_s`$ amounts to an extra phase difference $`\pi `$/2 between $`\xi _s`$ and $`\eta _s`$ in comparison with that in subluminal case. Hence we see onceagain that it is the phase which plays a dominant role in quantum mechanics.

(c)The origin ( or mechanism ) of mass generation has been remained in physics as a mystery for hundred years. It is time to say something now. According to detailed analysis in Ref. , the finite and changeable mass of a fermion is not due to one kind of excitation which varies continuously inside but a manifestation of coherent cancellation effect between two fields rotating with opposite helicities implicitly, either one of them can be excited infinitely in essence. It is precisely this fantastic effect explains the amazing result that a tachyon’s energy (mass) approaches zero when its velocity $`u`$ increases to infinity (as shown in Eq. (4)).

Acknowledgements: The author wishes to thank Dr. Tsao Chang for bringing the superluminal problem of neutrino to his attention and relevant discussions.


text:
# Correlated defects, metal-insulator transition, and magnetic order in ferromagnetic semiconductors

(January 22, 2002)

## Abstract

The effect of disorder on transport and magnetization in ferromagnetic III–V semiconductors, in particular (Ga,Mn)As, is studied theoretically. We show that Coulomb-induced correlations of the defect positions are crucial for the transport and magnetic properties of these highly compensated materials. We employ Monte Carlo simulations to obtain the correlated defect distributions. Exact diagonalization gives reasonable results for the spectrum of valence-band holes and the metal-insulator transition only for correlated disorder. Finally, we show that the mean-field magnetization also depends crucially on defect correlations.

*Introduction*.—Recently, there has been substantial interest in diluted ferromagnetic III–V semiconductors due to observations of relatively high Curie temperatures Ohno ; Ohno1 ; Reed . This makes these materials promising for applications as well as interesting from the physics point of view. They could allow the incorporation of ferromagnetic elements into semiconductor devices, and thus the integration of processing and magnetic storage on a single chip. To be specific, we consider here manganese-doped GaAs. The properties of this material rely on the dual role played by the Mn impurities: They carry a local spin due to their half-filled d-shell and dope the system with holes, which mediate a ferromagnetic indirect exchange interaction between the spins. Furthermore, the materials are highly compensated presumably due to antisite defects (As substituted for Ga) Potashnik , which drives them towards a metal-insulator transition (MIT). It is clearly important to understand the interplay between transport properties, magnetic ordering, and the defect configuration. In this letter we show that *correlated* defects are required for a description consistent with experiments.

There is a growing body of theoretical work on (Ga,Mn)As Koenigrev . Most approaches either start from the heavily doped regime Dietl97 ; Jungwirth ; Koenig ; Abolfath ; Dietl ; CDM ; Schlie  or from the weak-doping limit BB . In the case of heavy doping Dietl97 ; Jungwirth ; Koenig ; Abolfath ; Dietl ; CDM ; Schlie  one considers valence-band holes moving in the disorder potential of the defects. For large hole concentration the Fermi energy $`E_F`$ lies deep in the valence band compared to the Coulomb-potential fluctuations and the latter are neglected. On the other hand, for light doping BB  the local impurity states overlap only weakly, forming an impurity band. This model is problematic since for the typical doping range the impurity band would be much broader than the energy gap comment . Both approaches have in common that the effect of antisites on the disorder and correlations of defects are neglected.

The outline of the remainder of this paper is as follows. First, we discuss the configuration of Mn impurities and antisite defects with the help of Monte Carlo simulations and obtain the disorder potential and its spatial correlations. Then, we derive the spectrum and localization properties of holes in this potential and discuss the MIT. Finally, we calculate the polarizations of Mn and hole spins within a selfconsistent mean-field theory.

*Defect configurations*.—In the theory of doped semiconductors one often assumes the dopants to be randomly distributed SE ; Koenigrev . However, the present system is highly compensated. This leads to defects of either charge being present (relative to Ga, Mn impurities have charge $`e`$ while antisites carry $`+2e`$) and to a low density of carriers, which only weakly screen the Coulomb interaction between the charged defects SE . In the resulting *Coulomb plasma* the defect positions are expected to be highly correlated. We study here the limiting case where the defects come into thermal equilibrium during growth Keldysh . This is motivated by experiments Potashnik  suggesting that defects diffuse rapidly at $`250^{}`$C, which is a typical growth and annealing temperature for (Ga,Mn)As grown by molecular beam epitaxy. The resulting configuration is assumed to be quenched at low temperatures.

To find typical defect configurations close to equilibrium we perform Monte Carlo simulations for the classical Coulomb system of Mn impurities and antisites on the Ga sublattice with the Hamiltonian $`H=1/2_{i,j}q_iq_j/(ϵr_{ij})e^{r_{ij}/r_{\mathrm{scr}}}`$, where $`q_i`$ are the defect charges, $`r_{ij}`$ is their separation, and $`ϵ11`$. The screening length $`r_{\mathrm{scr}}`$ is obtained from nonlinear screening theory SE . $`r_{\mathrm{scr}}`$ is much larger than the nearest-neighbor separation of Ga sites for realistic parameters so that it hardly affects the small-scale defect correlations relevant here. We employ the Metropolis algorithm at $`250^{}`$C for systems of $`20\times 20\times 20`$ conventional cubic unit cells with periodic boundary conditions, unless stated otherwise.

A typical result for Mn concentration $`x=0.05`$ and $`p=0.3`$ holes per Mn, as suggested by experiments Ohno , is shown in Fig. 1. The concentration of antisites is obtained from $`x`$ and $`p`$ under the assumption of charge neutrality. Most of the defects have formed clusters. This leads to screening of the disorder potential by the defects, as can be seen from the potential correlation function $`D(r)V(𝐫)V(𝐫^{})_{\left|𝐫𝐫^{}\right|=r}V^2`$ plotted in Fig. 2 for $`x=0.05`$ and $`p=0.3`$. Note that $`\mathrm{\Delta }V\sqrt{D(0)}`$ is the width of the distribution of $`V(𝐫)`$. For correlated defects, $`D(r)`$ and thus $`\mathrm{\Delta }V`$ are strongly reduced and the Coulomb potential is screened on the much shorter length scale $`r_{\mathrm{ion}}`$, which is of the order of the nearest-neighbor Ga site separation rem.rscr . (This ionic screening is only possible due to the presence of defects of either charge. Thus it does not take place, *e.g.*, in ferromagnetic II–VI semiconductors since Mn is isovalent in these materials.) Importantly, however, $`\mathrm{\Delta }V0.284`$eV is still *not* small compared to the Fermi energy $`|E_F|0.329`$eV.

Recent high resolution x-ray diffraction Schott  as well as resistivity and magnetization measurements Potashnik  support the formation of small clusters. The more direct determination of the local structure around Mn dopants by extended x-ray-absorption fine structure Shioda  does not give direct information on the abundance of specific elements around the Mn. However, the observed loss of order around Mn is certainly consistent with cluster formation. Finally, we show in the following that experimental results for the band gap, transport, and the magnetization Ohno ; Beschoten  strongly support the clustering.

*Hole states*.—To find the spectrum and localization properties of valence-band holes moving in the disorder potential $`V(𝐫)`$, we start from the Hamiltonian $`H=_i(\mathrm{}^2/2m^{})_i^2+V(𝐫_i)`$, where we use the envelope function and parabolic-band approximations rem.KL . The calculations are done for spin-less holes, which is justified since the additional disorder introduced by the exchange interaction with the Mn spins is much smaller than $`\mathrm{\Delta }V`$.

The Hamiltonian is solved by exact diagonalization in a plane-wave basis of $`13^3=2197`$ states rem.number , giving the energy spectrum and normalized eigenfunctions. The Fermi energy $`E_F`$ is obtained from the hole concentration $`px`$. We also obtain the inverse participation ratio $`\mathrm{IPR}(n)=1/_𝐫|\psi _n(𝐫)|^4`$ of the states $`\psi _n(𝐫)`$. The IPR allows to estimate the position of the mobility edge, since it is of the order of the volume for extended states but is essentially independent of system size for localized states. Figure 3 shows a comparison of the IPR as a function of energy for uncorrelated defects and for the equilibrium configuration obtained above, for 5% Mn and $`0.3`$ holes per Mn. The valence-band edge is strongly smeared out by uncorrelated defects so that the energy gap is completely filled, which is not observed experimentally, whereas correlated defects only lead to a small tail in the gap. We have also studied the dependence of the IPR on system size $`L^3`$ over the limited feasible range in $`L`$, finding that $`\mathrm{IPR}(n)L^3`$ for the flat part of the curve, signifying extended states, whereas the IPR is independent of $`L`$ for states sufficiently far in the tail, which are thus localized. Hence, the mobility edge lies on the slope of the curve. The states are much more localized for uncorrelated defects, due to the larger disorder potential, whereas for correlated defects the states close to $`E_F`$ are clearly extended. Thus the spectrum and the localization properties agree with experiments Ohno  only for correlated defects.

In Fig. 4 we show the IPR as a function of energy for correlated defects at various Mn dopings and hole concentrations. The hole concentration for given Mn doping was chosen in accordance with experiments Ohno ; Matsukura . The spectrum and the IPR change little with Mn doping, since the potential screened by the defects on the length scale $`r_{\mathrm{ion}}`$ is much less affected by a change of concentration than the bare Coulomb potential.

The dominant effect comes from the hole concentration. Our result for the scaling of the IPR suggests that for 5% and 3% Mn the states at $`E_F`$ are extended and we expect metallic behavior. For 1% Mn $`E_F`$ is close to the mobility edge, while for 0.5% Mn the states at $`E_F`$ are localized. From this we estimate the MIT to take place at a Mn concentration of the order of 1%. This result is in reasonable agreement with experiments Ohno . For comparison, uncorrelated defects would lead to an MIT at about 5% Mn, see Fig. 3.

*Magnetization*.—Finally, we discuss the spontaneous magnetization and its dependence of the defect configuration. The exchange interaction between hole and Mn spins is described by the Hamiltonian $`H=_i(\mathrm{}^2/2m^{})_i^2+V(𝐫_i)J_{\mathrm{pd}}_{i,l}𝐬_i𝐒_l\delta (𝐫_i𝐑_l)`$, where $`V(𝐫)`$ is the disorder potential obtained above, $`J_{\mathrm{pd}}45`$eV Å3 is the exchange integral Okaba ; Bhatta ; Abolfath , $`𝐬_i`$ and $`𝐒_l`$ are the hole and Mn spin operators, respectively, and $`𝐑_l`$ are the Mn impurity positions. In second quantized form, the Hamiltonian reads

$`H`$ $`=`$ $`{\displaystyle \underset{n\sigma }{}}c_{n\sigma }^{}ϵ_nc_{n\sigma }`$
$`J_{\mathrm{pd}}{\displaystyle \underset{n\sigma ,n^{}\sigma ^{}}{}}{\displaystyle \underset{l}{}}c_{n\sigma }^{}\psi _n^{}(𝐑_l){\displaystyle \frac{𝝉_{\sigma \sigma ^{}}}{2}}𝐒_l\psi _n^{}(𝐑_l)c_{n^{}\sigma ^{}},`$
where $`ϵ_n`$ are the eigenenergies obtained in the absence of exchange as discussed above, $`\psi _n(𝐫)`$ are the corresponding eigenfunctions, and $`𝝉`$ is the vector of Pauli matrices. The exchange interaction is decoupled at the mean-field level, introducing the averaged Mn spin polarizations $`𝐌_l𝐒_l`$ and hole spin polarizations at the Mn sites, $`𝐦_l_{n\sigma ,n^{}\sigma ^{}}c_{n\sigma }^{}\psi _n^{}(𝐑_l)(𝝉_{\sigma \sigma ^{}}/2)\psi _n^{}(𝐑_l)c_{n^{}\sigma ^{}}`$. We do *not* perform a spatial average, i.e., the disorder is retained. The hole and Mn sectors can now be diagonalized separately, $`𝐌_l`$ and $`𝐦_l`$ are calculated, and the procedure is iterated. This selfconsistent mean-field theory is similar to the one employed by Berciu and Bhatt BB . The main difference is that we start from realistic hole states for the disorder potential.

Figure 5(a) shows the Mn and hole spin polarizations for $`x=0.05`$ and $`p=0.3`$. Also shown are the polarizations obtained neglecting disorder by using plane waves for the hole states and the virtual crystal approximation (VCA) for the Mn spins Jungwirth ; Koenigrev ; rem.KL . We see that disorder strongly enhances $`T_c`$ BB . This increase is easily understood: $`T_c`$ is determined by the indirect exchange interaction of nearest-neighbor Mn pairs, which is inversely proportional to the separation Dietl97 ; Jungwirth ; Koenig . Without disorder, the typical distance is $`r_{\mathrm{vca}}n_{\mathrm{Mn}}^{1/3}`$, where $`n_{\mathrm{Mn}}`$ is the density of Mn impurites. For $`x=0.05`$, $`r_{\mathrm{vca}}9.67`$Å. On the other hand, with disorder the nearest-neighbor separation on the Ga sublattice is $`r_{\mathrm{dis}}=3.99`$Å. Therefore, the indirect exchange, and thus $`T_c`$, is larger by a factor of order $`r_{\mathrm{vca}}/r_{\mathrm{dis}}2.4`$ in the disordered case.

For *correlated* defects the polarization curve looks more “mean-field-like,” in qualitative agreement with experiments Ohno ; Dietl ; Beschoten . This can be explained by the extended nature of the relevant hole states for correlated defects, which leads to a long-range indirect exchange interaction. Thus the effective field seen by a Mn spin is averaged over many other Mn spins and the spatial fluctuations of this effective field are small. On the other hand, for *uncorrelated* defects the holes are partly localized, the indirect exchange is of shorter range, and the fluctuations are larger. We believe that this leads to the smeared-out magnetization curve. It is clearly important to study the effect of fluctuations on the magnetization. We expect fluctuations to reduce $`T_c`$ in particular for the uncorrelated case. Figure 5(b) shows that $`T_c`$ is reduced for smaller Mn and hole concentrations and that the typical temperature scale is approximately proportional to $`x`$, in accordance with experiments Ohno . Note also that in all cases the holes are only partially polarized.

To conclude, we have obtained typical defect configurations in (Ga,Mn)As and studied the effect of the resulting disorder potential on the valence-band holes. The strong Coulomb interactions of charged defects leads to the formation of defect clusters during growth. The previously neglected antisite defects are crucial for this effect. We have shown that such a correlated defect distribution is required to understand the hole spectrum, the metal-insulator transition, and the shape of the magnetization curves. Our results should also help to clarify why the properties of ferromagnetic semiconductors depend so strongly on details of the growth process.—We profited from discussions with P. J. Jensen, J. König, J. Schliemann, and in particular M. E. Raikh, who emphasized to us the importance of compensation.

*Note added*.—Recently, Yang and MacDonald Yang  have studied the effect of the hole-hole interaction, assuming a random distribution of Mn dopants and antisites.


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

## I Introduction

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

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

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

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

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

This quantitative result has never been generalized to higher dimension or codimension tilings, because its generalization requires the knowledge of the free boundary entropy if one wants to use the same proof as in the hexagonal case. However, thanks to numerical simulations, it has recently been conjectured by Linde, Moore, Nordhal  that in dimensions higher than 2, the corresponding arctic region should be a polytope itself, with flat boundaries. It was further conjectured (by Destainville and Mosseri  and independently by Propp ) that in this case the entropy density should be spatially uniform and maximal in the unfrozen region. These conjectures renew the interest for the partition method since the relation between both entropies becomes amazingly simple in this case.

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

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

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

## II Partitions, tilings and boundary conditions

### A Generalities

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

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

We now define three-dimensional solid partitions. Consider a three-dimensional array of sides $`k_1\times k_2\times k_3`$. Fix an integer $`p>0`$, called the height of the partition problem. Put non-negative integers in the array, no larger than $`p`$, with the constraint that these integers decrease in each of the three directions of space. More precisely, if $`i_1`$, $`i_2`$ and $`i_3`$ are indices attached to the boxes of the array ($`1i_\alpha k_\alpha `$), we denote by $`n_{i_1,i_2,i_3}`$ the integral variables attached to these boxes (the parts). Our partition array contains $`N_p=k_1k_2k_3`$ parts. The partition constraint is

$$0n_{i_1,i_2,i_3}p$$
(1)
and

$$n_{i_1,i_2,i_3}n_{j_1,j_2,j_3}$$
(2)
whenever $`i_1j_1`$, $`i_2j_2`$ and $`i_3j_3`$.

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

Following the same construction in three dimensions, four-dimensional hypercubes are stacked above the three-dimensional partition array, with the heights of the stacks equal to the corresponding parts. Then project into three dimensions along the (1,1,1,1) direction of the hypercubic lattice. Like in the hexagonal case, the so-obtained tilings fill a polyhedron, a “rhombic dodecahedron” ($`RD`$) of integral sides $`k_1`$, $`k_2`$, $`k_3`$ and $`p`$ (see the outer frame in figure 3). The total number of tiles,

$$N_t=k_1k_2k_3+k_1k_2p+k_1k_3p+k_2k_3p.$$
(3)
We call tilings with rhombic dodecahedron boundaries “$`RDB`$-tilings” and denote their configurational entropy per tile  by $`\sigma _{fixed}`$.

The source of configurational entropy can be easily understood in either the tiling or the partition representation. Figure 2 illustrates the basic “tile flip” in both the $`32`$ and the $`43`$ tilings. In each case interior tiles are rearranged without disturbing the surface of a region. In terms of the equivalent partition, the height of one partition element increases or decreases by one unit. This elementary local move is ergodic, and every legal tiling can be reached by a succession of such flips.

An alternative description of the same tilings and partitions is in terms of de Bruijn lines  and membranes. Figure 1 (right) shows the de Bruijn lines of the tiling (center). These are formed by connecting centers of parallel edges on every rhombus. There are three families of de Bruijn lines, each family with an average orientation perpendicular to the tile edges. Although de Bruijn lines within a family never cross, de Bruijn line crossings among different families occur at the center of each rhombus. For three-dimensional tilings, we have instead de Bruijn membranes. Elementary flips such as those illustrated in Fig. 2 bend the membrane locally, and membrane fluctuations become the source of configurational entropy.

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

### B Boundary conditions

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

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

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

To exploit the calculational advantages of a partition representation, while achieving the physical free-boundary entropy in the thermodynamic limit, we adapt the partition method so that the corresponding tilings exhibit no frozen regions. The new boundary, even though fixed, has no macroscopic effect on tiling entropy in the thermodynamic limit. The tilings become homogeneous, displaying the free-boundary entropy density throughout.

The idea is to consider tilings (we focus on “diagonal” tilings with $`k_1=k_2=k_3=p`$) that fill a regular octahedron $`O`$ instead of the rhombic dodecahedron $`RD`$. Eight vertices of the $`RD`$ must be truncated to produce the $`O`$. We call tilings with octahedron boundaries $`OB`$-tilings. Such an octahedron is displayed in figure 3. It is inscribed in an $`RD`$ and has puckered boundaries instead of flat ones. Despite this puckering, the boundaries are effectively flat in the thermodynamic limit. The same kind of idea is developed in reference  in the hexagonal case where a puckered two-dimensional hexagonal boundary is introduced so that the tilings are homogeneous and exhibit no frozen regions. It is demonstrated (and numerically checked) that the random tilings filling this puckered hexagon display a free-boundary entropy in the thermodynamic limit.

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

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

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

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

Using the methods discussed in section II A, such partitions generate tilings filling the octahedron $`O`$. For an original partition cube of sides $`k_1=k_2=k_3=p`$, the full $`RD`$ contains $`N_t=4p^3`$ tiles (see equation (3)). The octahedron $`O`$ contains $`N_t=4p^34(p1)p(p+1)/3`$ tiles (these are the $`4p^3`$ tiles in the $`RD`$ minus the tiles in the truncated corners). These partitions are used in section III to estimate numerically $`\sigma _{free}`$.

### C Arctic octahedron conjecture (Linde, Moore and Nordahl)

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

In three dimensions, Linde, Moore and Nordahl have recently explored numerically the typical shape of a $`RDB`$-tiling and have conjectured  that the two-dimensional circle then becomes not a sphere but a regular octahedron, inscribed in $`RD`$ like $`O`$ in figure 3. More precisely, the unfrozen region is not exactly $`O`$ but tends towards $`O`$ at the large size limit.

This conjecture has a crucial consequence : the statistical ensemble of $`RDB`$-tilings is dominated by tilings periodic outside $`O`$ and random inside $`O`$, equivalently $`OB`$-tilings like in the previous section completed by eight periodically tiled pyramids to fill $`RD`$. Note that James Propp  has also conjectured independently of the present work that $`RDB`$-tilings should be homogeneous inside $`O`$. Taking into account the tiles in the frozen regions to calculate an entropy per tile, one finally gets

$$\sigma _{fixed}=\frac{N_O}{N_{RD}}\sigma _{free}=\frac{2}{3}\sigma _{free},$$
(4)
since the ratio of the numbers of tiles in $`RD`$ and $`O`$ is 3/2. The next section is dedicated to the calculation of this ratio via Monte Carlo simulations.

## III Monte Carlo

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

We employ a variant of the transition matrix method  that couples a conventional Metropolis Monte Carlo simulation with a novel data collection and analysis scheme to construct a numerical approximation to the transition matrix, described below in section III A. The density of states is an eigenvector of the exact transition matrix, and the sum of the density of states yields the total number of states. This method yields highly accurate absolute entropies with impressive efficiency.

### A Transition Matrix

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

$$E(P)=\underset{ijk}{}n_{ijk},$$
(5)
and incorporate the constraint of legality by defining $`E=\mathrm{}`$ for any partition $`P`$ that violates the ordering conditions inside the partition or at the boundaries. Under single vertex flip dynamics (see Fig. 2) a single partition element increases or decreases by one unit resulting in an energy change $`\mathrm{\Delta }E=\pm 1`$. The ground state of this model is the lowest height legal partition. For the boundary conditions employed here, the ground state is unique and we denote its energy as $`E_{min}`$. There is also a unique maximum energy state of energy $`E_{max}`$.

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 $`\mathrm{exp}(E/T)`$. Note that illegal partitions are never accepted. Low temperatures bias the partitions towards low heights, while high temperatures remove any bias according to height. As $`T`$ becomes large, the Metropolis simulation faithfully reproduces the random tiling ensemble in which all configurations have equal weight. If we choose, we may take negative values of temperature $`T`$. Large negative temperatures again reproduce the random tiling ensemble, while small negative temperatures favor partitions of maximal height.

For a given partition $`P`$, a number of partitions $`n_\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

$$n_{}(P)+n_0(P)+n_+(P)=2N_p$$
(6)
holds for every partition $`P`$. This value is twice the number of parts because each part can be raised or lowered. For any given partition $`P`$, the calculation of $`n_{\pm ,0}`$ is easy (though not fast) and exact.

We define the transition matrix $`\omega _{\pm ,0}(E)`$ as a matrix with $`E_{max}E_{min}+1`$ rows (one for each allowed energy $`E`$) and three columns (labeled “+” for upwards transitions, 0 for forbidden transitions and “-” for downwards transitions. Alternatively, we can think of these three columns as the diagonal and two off-diagonals of a square matrix $`\omega (E,E^{})`$ of dimension $`E_{max}E_{min}+1`$. Formally, we define

$$\omega _{\pm ,0}(E)\frac{1}{W(E)}\underset{P(E)}{}n_{\pm ,0}(P)/2N_p$$
(7)
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, preventing us from actually calculating the transition matrix. However, we may calculate the matrix numerically with high precision by averaging $`n_{\pm ,0}(P)/2N_p`$ over those partitions $`P`$ occurring during a Monte Carlo simulation. By virtue of eq. (6), the matrix elements obey the sum rule

$$\omega _{}(E)+\omega _0(E)+\omega _+(E)=1$$
(8)
for any energy $`E`$, so each row of the transition matrix can be normalized independently without knowledge of $`W(E)`$.

The transition matrix can be interpreted in terms of the rates at which Monte Carlo moves would be accepted in a hypothetical simulation at infinite temperature (even if our actual simulation is performed at finite temperature). In a single Monte Carlo step, the probability for upwards or downwards transitions at infinite temperature is $`n_\pm (P)/2N_p`$, and the probability a move will be rejected is $`n_0(P)/2N_p`$. We can predict finite temperature transition probabilities by multiplying $`\omega _{\pm ,0}(E)`$ with appropriate Boltzmann factors, and in fact these agree well with actual observed acceptance rates.

To extract the density of states and the entropy, consider the infinite temperature detailed balance condition. At infinite temperature the probability a randomly chosen partition has energy $`E`$ is $`W(E)/Z`$, where

$$Z=\underset{E}{}W(E)$$
(9)
is the total number of legal partitions (the “partition function”). Multiplying the probability $`W(E)/Z`$ 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

$$\omega _+(E)W(E)/Z=\omega _{}(E+1)W(E+1)/Z.$$
(10)
It is useful to rearrange the detailed balance eq. (10) to find

$$W(E+1)=\frac{\omega _{}(E+1)}{\omega _+(E)}W(E)$$
(11)
which allows us to iteratively extract the full density of states function $`W(E)`$ using uniqueness of the ground state, $`W(E_{min})=1`$. Finally, the total entropy

$$S=\mathrm{ln}Z,$$
(12)
and the entropy density $`\sigma =S/N_t`$.

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

The accuracy of our result is controlled by the accuracy with which we determine $`\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. We store the matrix as the integer valued sum of the integers $`n_{\pm ,0}(P)`$, and impose the normalization (8) 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, ordinary 4 byte integers can not hold the data. We implemented special procedures to handle storage and algebraic manipulations of large integers.

### B Numerical Data

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

For both fixed and free boundary tilings, symmetries of the model dictate that the density of states is symmetric about the midpoint energy $`E_{mid}=(E_{min}+E_{max})/2`$. The density of states $`W(E)`$ has a strong maximum at the midpoint. By including both positive and negative temperature sweeps we sample both low and high energy states. By reversing our sweeps we mitigate possible systematic sampling errors associated with the direction of the sweep.

At each temperature during a sweep, we accumulate data on $`N_{MC}`$ sample configurations. For a given configuration (partition $`P`$), the time required to calculate $`n_{\pm ,0}(P)`$ is proportional to the partition size $`N_p`$, and hence to the time required to attempt $`N_p`$ flips. Because the data accumulation is time consuming, we can perform many single vertex flips between sample without significantly slowing down the simulation. We take $`N_{FL}=N_p/5`$ single vertex flips, which yields rough equality between time spent flipping and collecting data. Prior to data collection at any temperature, we anneal at fixed temperature for $`N_{MC}\times N_{FL}/100`$ steps. The total number of attempted flips in an $`N_{MC}=10^6`$ run on a $`10\times 10\times 10`$ partition is thus in excess of $`4\times 201\times 10^6\times 10^3/5=1.6\times 10^{11}`$ attempted flips. A run of this length takes 5 days on a 1.7GHz Pentium 4 processor.

This protocol was chosen to ensure approximately uniform coverage of energies from $`E_{min}`$ to $`E_{max}`$. The energy distribution of partitions visited at each temperature overlaps strongly the energy distribution of partitions visited at the previous and following temperatures. The value of $`T_{min}`$ is chosen sufficiently low to guarantee some coverage of the extreme energy states. Because there are few states at the energy extremes, there is an entropic barrier to reaching these extremes. This causes a spike in the coverage close to the extreme energies that cannot be avoided using Metropolis sampling. The value of $`T_{max}`$ is chosen sufficiently high that the midpoint energy is a local maximum in the coverage. Fig. 5 plots the number of partitions sampled as a function of energy during the longest runs (length $`10^6`$) for the $`p=4`$ partition. Note that the number of hits exceeds $`10^6`$ uniformly for each energy in the range $`0E256`$.

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

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

Table I shows the convergence of entropy data as run length grows. We show the convergence only for the largest system size, which represents our worst case. The residual $`R=\mathrm{log}W(E_{max})`$ measures the failure of the calculated microcanonical entropy to return to zero at high energy. It reflects the cumulative error in $`W(E)`$, and thus, assuming statistical independence of the errors at each energy, it provides an upper bound on the error in $`W(E)`$ for any energy. Converting the error in $`W(E)`$ into an error in the entropy, we estimate an uncertainty

$$|\mathrm{\Delta }\sigma |R/N_t.$$
(13)
Our numerical simulation results in Table I, and also simulations of smaller systems for which the entropy is known exactly, are in good numerical agreement with this estimate.

Table II summarizes our data for all system sizes studied. We report the entropy resulting from the run of length $`N_{MC}=10^6`$, and use the difference between that run and the length $`N_{MC}=10^5`$ run for the quoted uncertainty.

To extrapolate our data to infinity, we fit to functional forms. For the fixed boundary entropy, we expect finite-size corrections of order $`\mathrm{log}(p)/p`$ and $`1/p`$ because of boundary effects. Indeed, these two terms are the first correction terms in the exactly solvable one- and two-dimensional random tilings. Note that they have the same order of magnitude for the values of $`p`$ considered here. This logarithmic form fits the data much better than a simple correction of order $`1/p`$. For fixed boundary tilings the data fit well to

$$\sigma _{fixed}(p)0.1450.0049\frac{\mathrm{log}(p)}{p}+\frac{0.034}{p},$$
(14)
from which we conclude that $`\sigma _{fixed}=0.145(3)`$. The value is obtained from a fit excluding the data point with $`p=1`$, while the uncertainty estimate comes from excluding instead the data point with $`p=10`$. This value is very close to a conjectured limit  of 0.139, but given our small uncertainty, we believe the conjectured value is not exact.

For the free boundary data we find

$$\sigma _{free}(p)0.2140.052\frac{\mathrm{log}(p)}{p}\frac{0.046}{p},$$
(15)
from which we conclude that $`\sigma _{free}=0.214(2)`$. As in the fixed boundary case, the logarithmic finite size correction fits the data better than a simple correction of order $`1/p`$. Should the octahedral fixed boundary have a trivial short-range effect, it would create only a $`1/p`$ leading term (from surface vs bulk contributions), as indeed occurs in the corresponding two-dimensional hexagonal case  (see II B). Therefore, it appears that the strain-free octahedral boundary has a nontrivial long-range effect on the local entropy density, presumably by limiting the fluctations of the height function within a “penetration depth” that depends on the side length $`p`$. A penetration depth growing like $`\mathrm{log}(p)`$ is consistent with our leading finite-size correction.

For the ratio, we fit to

$$\frac{\sigma _{free}(p)}{\sigma _{fixed}(p)}1.480.37\frac{\mathrm{log}(p)}{p}\frac{0.54}{p},$$
(16)
and we conclude the ratio $`\sigma _{free}/\sigma _{fixed}=1.48(3)`$. This value equals 3/2 within our uncertainty.

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

## IV Discussion

Returning to the arctic octahedron conjecture, we recall (see equation (4)) we expected a ratio $`\sigma _{free}/\sigma _{fixed}=3/2`$. In fact, with few extra hypotheses, we now demonstrate that our numerical results provide a strong support to the arctic octahedron conjecture. Suppose that, as in the two-dimensional case , the function $`\varphi _{\text{max}}`$ maximizing the entropy functional $`s[\varphi ]`$ subject to a given boundary condition is unique . If $`\varphi _0`$ is the piecewise linear function which vanishes identically in the central octahedral region $`O`$ and has maximum strain (and 0 entropy) in the eight pyramidal regions comprising $`RDO`$, then $`s[\varphi _0]=2/3\sigma _{free}`$, and thus $`s[\varphi _0]=\sigma _{fixed}`$. Therefore $`\varphi _0=\varphi _{\text{max}}`$, by uniqueness of $`\varphi _{\text{max}}`$.

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

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

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

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

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

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

To finish with, we mention that this transition matrix Monte Carlo technique can easily be adapted to the numerical calculation of the entropy of two-dimensional rhombus tilings. Indeed, the structure of the configuration space of such tiling problems in terms of flips has been characterized in reference . Note however that no such simple result as the “arctic octahedron phenomenon” is expected in these two-dimensional classes of tilings, but the calculation of their configurational entropy is a challenge in itself. This work is in progress.

###### Acknowledgements.

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