# Thermodynamics of a Tiling Model

## 1 Introduction to the model

In this work we propose a model that shows a non-trivial thermodynamic behaviour due to its particular geometric structure.

This two-dimensional model is built using square tiles, called Wang’s tiles, on a square lattice. The edges of these tiles can be of six different ’colours’. But of all the possible types ($`6^4`$) that can be created changing the colours just a particular group of sixteen, found by Ammann , is taken into account (see figure 1).

This is one of the minimal set of Wang tiles such that the corresponding tilings exist and are non-periodic tilings. A tiling is a configuration of tiles placed edge-to-edge on the plane, where all contiguous edges have the same colour. If at least one tiling is allowed and none of the possible tilings shows a periodic pattern, then the set of tiles composing them is called aperiodic. Other aperiodic sets of tiles are often used as models for quasi-crystal materials , but this is not the present case.

The Wang tiles were the first aperiodic tiles to be discovered, in 1966 by Berger . They were initially important, and they still are, because of their use in problems of mathematical logic . The use that we make of them in this work is, anyway, far away from that point of view. Indeed we look at the behaviour of a system built by Wang tiles in a thermal bath, defining thermodynamic observables on these tilings. The main stimulus to study this model has been the high degeneracy of the perfectly matched configurations and their non-trivial aperiodic structure.

Putting the system in a thermal bath allows for translations of the tiles also in positions where there is no matching between edges of neighbours tiles, forming in this way a unmatched-tiling. We can assign an energy to each configuration which is equal to the number of links such that the colours of the facing edges are different. The energy of the exactly matched configurations, that from now on we will call ground states, is then zero.

If we label the type of tile in the position given by the coordinates $`(x,y)`$ in the plane by $`T_{x,y}`$ and the type of its four edges respectively towards South, East, North and West by $`T_{x,y}^{(S)}`$, $`T_{x,y}^{(E)}`$, $`T_{x,y}^{(N)}`$ and $`T_{x,y}^{(W)}`$ we can write an Hamiltonian for this system in the following way:

$$=\underset{(x,y)}{\overset{1,L1}{}}\left(1\delta \left(T_{x,y}^{(E)}T_{x+1,y}^{(W)}\right)\right)+\underset{(x,y)}{\overset{1,L1}{}}\left(1\delta \left(T_{x,y}^{(N)}T_{x,y+1}^{(S)}\right)\right).$$
(1)   
$`L`$ is the number of tiles in the lattice in one direction. The $`\delta (z)`$ is the delta of Kronecker.

## 2 Equilibrium analysis and critical behaviour

The numerical simulations done in order to study the equilibrium characteristics have been performed using the parallel tempering algorithm . We have chosen open boundary conditions on the two dimensional lattice. Indeed using periodic boundary conditions the fact that periodic tilings do not exist would have implied that the energy of the ground state would have been different from zero. With open boundary conditions the energy of the ground state is, by definition of (1), equal to zero. It is equivalent to say that all the tiles are edge by edge perfectly matched, forming some non-periodic structure.

### 2.1 Phase transition

At different temperatures, for every system we have computed the energy and the specific heath. Numerical simulations have been brought about on square lattices of different size: from a linear size of 8 to one of 32. Every equilibrium simulation has been carried out for a number of MC steps going from 10 millions to 100 millions, relatively to the size of the lattice. A range of temperature between 0 and 5 has been observed, then concentrating in small intervals in what comes out to be the critical region, around $`T=0.4`$. We controlled that at zero temperature the energy goes to zero and that the specific heat computed from the derivative of the energy coincides with the one computed from the energy fluctuations. These are checks of a correct thermalization.

The specific heat presents a smeared but clear change around the temperature $`0.4`$ (see figure 2), from a lower value for lesser temperatures to a higher value for greater. Increasing the size of the system we can observe that the crossing points between two curves of different size moves towards right and that the slope of the $`C_L(T)`$ line increases in this interval of temperatures as $`L`$ increases, but we also see that a peak grows around $`T=0.42`$.

Performing Finite Size Scaling (FSS) analysis of the crossing points of the specific heat curves we find that their abscissa tends to a $`T_{cross}(\mathrm{})=0.398\pm 0.007`$ as $`L\mathrm{}`$, with an exponent $`\nu =1.6\pm 0.5`$ and that the behaviour of the slope with increasing size is compatible with a diverging fit at that temperature. In this case we would have a jump at $`T_c=T_{cross}(\mathrm{})`$, which corresponds to a critical exponent $`\alpha `$ equal to zero.

But at the same time the fit of the peak height is consistent with a divergence at $`T_c=0.42\pm 0.01`$. Knowing that at criticality $`C|TT_c|^\alpha `$ and $`\xi |TT_c|^\nu L`$ we get $`CL^{\alpha /\nu }`$. Our data are both consistent with a power-law divergence ($`\alpha /\nu =0.35\pm 0.01`$), and with a logarithmic one ($`\alpha =0`$). A similar value of the critical temperature was found, in this last hypothesis, by Janowsky and Koch .

In any case there is evidence for a second order phase transition.

### 2.2 Order parameter

Wang’s tiles satisfy a particular property: if all the tiles of the two contiguous west and south sides (or equally of the north and east sides) of the square lattice are fixed, then at most one aperiodic tiling can be formed. This is due to the fact that each tile of the Wang set has different combinations of west-south (or north-east) colours. Then at most only one of them can be put, in the bottom-left (either top-right) corner. The same is true for the tiles to be put in the two corners created by putting the first tile. If any, there will be only one combination available also for them. Of course this is valid at $`T=0`$ in our system.

This condition implies that ground state degeneracy can not increase faster than $`16^{2L1}`$and therefore the entropy density is zero at zero temperature: $`s_L(0)a/L`$, where $`a<8\mathrm{log}2`$. Since just a few north-south or west-east combination are allowed the actual constant in the entropy is smaller than $`8\mathrm{log}2`$ (a stricter upper bound of $`a=\mathrm{log}12=3.585\mathrm{log}2`$ can be easily obtained).

In order to identify an order parameter we have to break this degeneracy.

We have then simulated the parallel evolution of two copies of the system, with the following procedure: one system reaches the equilibrium, then a copy of the equilibrium configuration is done and the evolution of this second copy is performed. One foundamental constraint is set: the boundary tiles on the south and east sides of the lattice stay unchanged in the dynamics towards equilibrium of the second copy.

For $`T>0`$ the equilibrium states are no more the exactly matched tilings. Because of the thermal noise some couple of neighbours sides of the Wang’s tiles can be unmatched. Furthermore there won’t be a unique tiling minimizing the energy. The second copy can then evolve to a different configuration.

We are interested in looking what happens when the temperature is increased from below to above the critical temperature, and how the behaviour is sensitive to the size of the system. With this aim we introduce an overlap that depends on the distance $`l`$ of the tiles from the two contiguous boundary sides that are fixed once that the first copy has reached equilibrium:

$$q(l)=\frac{1}{2(Ll)+1}\underset{yl,x=l}{}\underset{xl,y=l}{}\overline{\delta \left(T_{x,y}^{(1)}T_{x,y}^{(2)}\right)}$$
(2)   
Where the average $`\overline{(\mathrm{})}`$ has been performed over different realizations of the configurations, $`(\mathrm{})`$ is the time average at equilibrium and $`l=1,\mathrm{},L`$.

The overlap $`q(l)`$ is computed along the diagonal: starting from the vertex shared by the two fixed contiguous sides, $`q(1)`$, and ending at the opposite vertex of the lattice, $`q(L)`$. In order to gain more statistics the overlap is built averaging over the $`2(Ll)+1`$ elements in the row of ordinate $`l`$ with the coordinate $`xl`$ and in the column of abscissa $`l`$ with $`yl`$ and assigning this average value to the ’diagonal’ function $`q(l)`$.

The statistic sample is formed repeating the simulation several times with different initial configurations, in order to obtain different equilibrium configurations (we have always checked if the same equilibrium configuration could possibly appear more than once starting from different initial configurations to avoid annoying biases of the statistic sample, but this never happened). Every size has been simulated for at least 100 different initial configurations.

At zero temperature the overlap is always one.

At higher temperature but below the phase transition it goes to a constant lesser than one (far enough from the boundary $`L`$, where finite size effects are overwhelming).

Above the transition point $`q(l)`$ decays very rapidly, compatible with an exponential, to the lowest possible value, corresponding to completely uncorrelated copies (hot or disordered phase). Since the tiles are of sixteen different types and the probability distribution of the tile-types is uniform, this lowest value is equal to $`1/16`$.

In figure 3 we show the behaviour of $`q(l)`$ at different temperatures, both above and below the critical one, for different sizes. From this figures we can see that the $`q(l)`$ seems to approach some plateau for $`lL/2`$ at temperatures below $`T0.42`$. This is more evident as we go to bigger sizes.

We can compute a $`q(T)`$ for every size, averaging over the plateau values of $`q(l)`$. The plateau is every time chosen taking a small window around $`L/2`$ and then enlarging it as long as the plateau value stays constant. As $`l`$ grows further finite size effects destroy this plateau. We observe that this parameter shows a change approaching $`T_c`$: when $`T<T_c`$, it increases from the value of $`1/16`$ that it has at high temperature to the value of 1 that it reaches at zero temperature. The profiles of $`q_L(T)`$ are plotted in figures 4.

The probability distribution of the overlap’s values at the plateau has a non-trivial shape as soon as $`T<T_c`$ (figure 5). This is not so astonishing since fixing the boundary conditions we have broken the degeneracy of the equilibrium states. (We stress that due do the overlap’s definition the distribution doesn’t show any symmetry $`qq`$). The form of the $`P(q)`$ averaged over different realizations of the tile’s configurations at two continuous edges of the lattice is not strongly dependent from the size of the samples, at least for the simulated cases.

The probability distributions of the overlap for different realizations of the equilibrium boundary conditions in the bottom and left edges of the tiling are represented in figure 6.

## 3 Off-equilibrium analysis

In order to study systems with very slow relaxation times, much greater than the experimental times, it is very important to look at the off-equilibrium dynamics. This approach describes glasses, spin-glasses and, in general, any system with many equilibrium states divided by high free energy barriers, in a more realistic way in comparison with an equilibrium approach, since equilibrium is, in practice, never reached during experiments.

To keep the sample out of equilibrium during the run of a numerical simulation we need the typical distance $`\xi (t)`$ over which the system has reached equilibrium at a certain time to be always smaller than the linear size $`L`$ of the lattice. This distance is also called dynamical correlation distance. Practically, in order to be sure of avoiding thermalization, we choose sizes bigger than those used in the static analysis. The greater is the size the longer is the time needed to the correlation distance to reach the size of the system. More over we used in this case a standard Monte Carlo algorithm instead of the parallel tempering, so that the thermalization times increased sensitively under the transition point.

In order to study the response of the system to an external field we can embed the system into a perturbative field. Namely a field directed in one of the sixteen possible directions in the ’tile-types space’, where a particular tile of the set of sixteen can be favoured (positive field) or disfavoured (negative field). This can be a uniform or non-uniform field. To avoid any preference towards one particular type of tile we have chosen a non-uniform random field, whose value at each site is independent from the other sites.

Thus we add to the Hamiltonian (1) the perturbative term:

$$\underset{(x,y)}{}h_{x,y}^{sign}\left(1\delta \left(T_{x,y}h_{x,y}^{type}\right)\right).$$
(3)   
where $`h_{x,y}=h_{x,y}^{type}h_{x,y}^{sign}`$ is the random external field pointing along one of the tile-types or opposite to it.

$`h_{x,y}^{type}`$ has a uniform distribution of the sixteen possible choices of tile-type, and $`h_{x,y}^{sign}`$ gives the magnitude of the field and its sign (it can be randomly positive or negative, always according to a uniform distribution).

The probability distribution can then be written as:

$$𝒫(h_{x,y})=\frac{1}{16}\underset{n}{\overset{0,15}{}}\delta \left(h_{x,y}^{type}e^{i\pi \frac{n}{16}}\right)\times \frac{1}{2}\left(\delta (h_{x,y}^{sign}h_o)+\delta (h_{x,y}^{sign}+h_o)\right)$$
(4)   
with $`\overline{h_{x,y}}=0`$ and $`\overline{h_{x,y}^2}=h_o^2`$. $`n=0,\mathrm{},15`$ gives the ’direction’ of the field in a representation in the complex plane (see figure 7)

Once the system is cooled from high temperature, it is left evolving until a certain time $`t_w`$, usually called waiting time. At $`t=t_w`$ the field is switched on and we begin recording the values of the temporal correlation function $`C(t,t_w)`$ and of the integrated response function $`m(t,t_w)`$.

For our model they are defined as follows:

$$C(t,t_w)=\frac{1}{L^2}\underset{x,y}{}\delta (T_{x,y}(t)T_{x,y}(t_w)),$$
(5)   
$$m(t,t_w;h)=\frac{1}{L^2}\underset{x,y}{}\frac{\overline{h_{x,y}^{sign}(t_w)\delta (T_{x,y}(t)h_{x,y}^{type}(t_w))}}{h_o}$$
(6)   
and the susceptibility is

$$\chi (t,t_w)=\underset{h_o0}{lim}\frac{m(t,t_w;h)}{h_o}.$$
(7)   
For numerical computation it becomes:

$$\chi (t,t_w)\frac{m(t,t_w;h)}{h_o}.$$
(8)   
Here $`(\mathrm{})`$ is the average over different dynamical processes and $`\overline{(\mathrm{})}`$ the average over the random realizations of the perturbative external field. In a system at equilibrium the Fluctuation-Dissipation Theorem (FDT) holds:

$$\chi (tt_w)=\frac{1C(tt_w)}{T}.$$
(9)   
where we have made explicit use of the fact that the correlation function is defined in such a way that $`C(t_w,t_w)=C(0)=1`$.

Out of equilibrium this theorem is no more valid. Anyway a generalization is possible , at least in the early times of the dynamics. This generalization is made introducing a multiplicative factor $`X(t,t_w)`$ depending on two times such that :

$$\chi (t,t_w)=\frac{X(t,t_w)}{T}\left(1C(t,t_w)\right)$$
(10)   
In a certain regime, called aging regime, $`X(t,t_w)<1`$ and the FDT is violated. An important assumption is that this modified coefficient $`X(t,t_w)/T`$ depends on $`t`$ and $`t_w`$ only through $`C(t,t_w)`$. Its inverse is also called effective temperature $`T_eT/X\left[C(t,t_w)\right]`$ since the system in this regime, for a given time-scale, seems to behave like a system in equilibrium at a temperature different from the heat-bath temperature.

In terms of susceptibility and correlation functions we can write a general functional dependence:

$$\chi (t,t_w)=\frac{1}{T}S\left[C(t,t_w)\right]$$
(11)   
The $`S[C]`$ here defined would be $`1C`$ if we were at equilibrium.

From our probe we find that we first have a regime where the relation is linear with a coefficient equal to the heat-bath temperature. This early regime is sometimes called stationary, since the observables computed at very short times (compared with $`t_w`$) do not depend on the age of the system. During this regime the system goes fast towards a local minimum.

After this first time the $`S[C]`$ bends and the coefficient of $`X(t,t_w)/T`$ is no more the inverse heat-bath temperature. $`X(t,t_w)`$ is now less than one, as if the system would be at an effective temperature bigger than the one of the heat-bath. It also seems to change continuously as the system evolves, until it reaches zero (figures 8 \- 10).

The value of the autocorrelation function at which $`S[C]`$ leaves the $`1C`$ line, equal to the average overlap value $`q`$ in the static, increases continuously with temperature, as we already observed in the static analysis.

The $`S[C]`$ shows no strong dependence from the magnitude of the field, at least for the values that we have used to perturbate the system. Furthermore the system shows only a slight dependence from the size (measures on $`L=64,128,256`$) (fig. 9).

From the figures 10 we can see that $`S[C]`$ moves towards some asymptotic line increasing $`t_w`$.

The initial difference between $`1C`$ and $`S[C]`$ is due to lack of statistics: if $`N_{fr}`$ is the number of field realization performed in the simulation there is a difference of order $`N_{fr}/L^2`$ between the value of the integrated response function in the thermodynamic limit and the value at finite size.

We can also look at the link with the static analysis. Like in the case of mean field spin glasses we can suppose, following , that for $`t,t_w\mathrm{}`$, $`C(t,t_w)q`$ and $`X[C(t,t_w)]x(q)`$, where $`x(q)`$ is the cumulative distribution of $`P(q)`$:

$$x(q)=_0^q𝑑q^{}P(q^{})$$
(12)   
If this link is valid we can connect the susceptibility multiplied by the heat-bath temperature at a certain correlation value ($`S[C]`$) to the integral $`_1^C𝑑qx(q)`$.

In our case there is an agreement between dynamic data and the values of this last integral computed on the static data. The two approaches are consistent if the external field is not too large, in order to avoid the non-linear effects that are neglected in the derivation of formula 6. The field cannot even be too small, though, in order to make the system move from the meta-stable state in which it has gone during $`t_w`$. The lower is the temperature at which we cool the system, the bigger is the external field necessary to make it explore other parts of the space of states. A probe at too low temperatures then is not really working because non-linear effects interfere heavily or because the system does not reach the aging regime.

Two examples of the agreement between dynamic and static data are shown in figure 10.

## 4 Conclusions

We have studied the thermodynamics of a tiling model built by Wang tiles. For this kind of system we have found evidence of a phase transition from a completely disordered phase, in which the tiles on the plane are completely uncorrelated between each other, to a phase in which they begin to present an organized, also if very complicated, structure. For $`T0`$ this structure becomes an exactly matched tiling.

In order to characterize the phase of the system we have determined an order parameter with a non-zero mean value below $`T_c`$. It is the overlap between two tilings built with the same boundary conditions on the lower and the left sides of the lattices. Our data hint that it could have a non-trivial probability distribution under the phase transition.

Under the critical point the tiling system shows the aging phenomenon: the answer of the system to an external perturbation and the time autocorrelation function depend on the history of the system. This brings to a violation of the fluctuation-dissipation theorem for $`T<T_c`$. The integrated response function times the temperature ($`S[C]`$) shows a progressive bending when the system leaves the stationary regime until it reaches a constant value ($`X=0`$) and the dynamically and statically determined $`q`$ values increase continuously with decreasing temperature.

From this kind of behaviour it is not clear wether the model belongs to the class of systems showing domain growth or it is rather more similar to a spin glass in magnetic field. Indeed for very long times (small values of the correlation function) the fluctuation-dissipation ratio goes eventually to zero and it can not be excluded that the dynamics evolves through domains growth , even though in our case the nature of the domains of tiles should still be theoretically understood. Nevertheless for a very large interval of time, i.e. of values of $`C`$ in the plot of figure 10, the $`S(C)`$ is continuosly bending ($`0<X<1`$) like in a spin glass model , especially in the regions $`0.3<C<0.6`$, for $`T=0.35`$, and $`0.4<C<0.7`$, for $`T=0.3`$, as shown in figure 10. For values of $`C`$ smaller than these the $`S(C)`$ curves flatten but the predictions obtained from the static behaviour (the $`P(q)`$ are shown in figure 5) are also nearly flat, making quite difficult the distinction between a domain growth dynamics, where the response function is constant in the aging regime, and a more complicated behaviour where even the response function shows aging.

Eventually we show that our data are consistent with the equivalence between the equilibrium function $`_C^1𝑑qx(q)`$ and the dynamically determined $`S[C]`$ as $`t,t_w\mathrm{}`$.

Acknowledgments We warmly thank J.Kurchan for having stimulated this work and for his advises.