text:
# Monte Carlo Simulations of Spin Glasses at Low Temperatures

## I Introduction

There has recently been a renewed interest in the nature of the spin glass phase. Two principal theories have been investigated: the “droplet model” proposed by Fisher and Huse, (see also Refs. ), and the replica symmetry breaking (RSB) picture of Parisi. These scenarios and some others have also been considered by Newman and Stein. An important difference between these models concerns the number of large-scale, low energy excitations. RSB theory follows the exact solution of the infinite range Sherrington-Kirkpatrick (SK) model in predicting that there are excitations which involve turning over a finite fraction of the spins and which cost only a finite amount of energy in the thermodynamic limit. In addition, the surface of these large, finite energy excitations is expected to be space filling which means that the fractal dimension of their surface, $`d_s`$, is equal to the space dimension, $`d`$. The droplet theory argues that the lowest energy excitation involving a given spin and which has linear spatial extent $`L`$ typically costs an energy $`L^\theta `$, where $`\theta `$ is a (positive) exponent. Hence, in the thermodynamic limit, excitations which flip a finite fraction of the spins cost an infinite amount of energy. The droplet theory also predicts that $`d_s<d`$.

Recently, Krzakala and Martin (KM), and two of us (PY), have argued that a straightforward interpretation of their numerical results at zero temperature is intermediate between the droplet and RSB pictures in that there appear to be large scale excitations whose energy does not increase with size, but these have a surface with $`d_s<d`$. This interpretation of the results of KM and PY has, however, been recently challenged (though in opposite senses) by Marinari and Parisi and by Middleton. There has also been recent debate as to whether the $`\pm J`$ model has similar behavior to that of a model with a continuous distribution. In the scenario of KM and PY it is necessary to introduce two exponents which describe the growth of the energy of an excitation of scale $`L`$: (i) $`\theta (>0)`$ such that $`L^\theta `$ is the typical change in energy when the boundary conditions are changed, for example from periodic to anti-periodic, and (ii) $`\theta ^{}`$, which characterizes the energy of clusters excited within the system for a fixed set of boundary conditions. In this paper we test whether the picture proposed by KM and PY is compatible with finite temperature Monte Carlo simulations.

Several previous Monte Carlo simulations have found evidence for finite energy large scale excitations by looking at the order parameter function $`P(q)`$. In the thermodynamic limit this has delta functions at (plus or minus) the Edwards-Anderson order parameter, $`q_{EA}`$, corresponding to ordering within a single valley, and, according to RSB theory, a tail with a finite weight extending down to $`q=0`$. In the droplet theory, $`P(q)`$ is trivial, i.e. has only delta functions at $`\pm q_{EA}`$, though in a finite system there is a weight at the origin which vanishes with increasing $`L`$ like $`L^\theta `$.

These earlier Monte Carlo studies have found that the weight of $`P(q)`$ at the origin is independent of the systems size for temperatures down to $`T0.7T_c`$ in three dimensions and $`T0.6T_c`$ in four dimensions . However, these studies have been criticized as being too close to the critical point, so that the results are affected by critical fluctuations and very much larger sizes would be needed at these temperatures to see the asymptotic behavior of the low-temperature spin-glass state. Refs.  also argue, however, that clear evidence for droplet theory behavior could be seen even for quite small sizes at very low temperatures.

In this paper we check this prediction by performing Monte Carlo simulations in the low temperature region, though with an admittedly modest range of sizes, using the “parallel tempering” Monte Carlo method, also known as Exchange Monte Carlo. One difficulty with this approach is to ensure equilibration since the technique proposed earlier by one of us and Bhatt for conventional Monte Carlo does not work for parallel tempering. Here we use an alternative method, valid for the important case of a Gaussian distribution (which we use here), and which is closely related to the approach of Ref.  for the SK model.

Both in three and four dimensions, we find a tail in $`P(q)`$ which is independent of size (up to the sizes studied), for temperatures down to $`T0.2T_c`$ in 3D and $`T0.1T_c`$ in 4D, in contrast to the prediction of Refs. . We also find that data for the “link overlap”, defined below, fits well a description with $`d_s<d`$, though the extrapolation to the thermodynamic limit is quite large here. Thus our results are completely consistent with the earlier proposal of KM and PY.

We consider the short-range Ising spin glass in three and four dimensions, and, in addition, the Viana-Bray model. The latter is infinite range but with a finite average coordination number $`z`$, and is expected to show RSB behavior. All these models have a finite transition temperature. In the 3D case, the exponent $`\theta `$ obtained from the magnitude of the change of the ground state energy when the boundary conditions are changed from periodic to anti-periodic is about $`0.2`$, whereas in 4D it is much larger, about 0.7.

The Hamiltonian is given by

$$=\underset{i,j}{}J_{ij}S_iS_j,$$
(1)
where, for the short range case, the sites $`i`$ lie on a simple cubic lattice in dimension $`d=3`$ or 4 with $`N=L^d`$ sites ($`L8`$ in 3D, $`L5`$ in 4D), $`S_i=\pm 1`$, and the $`J_{ij}`$ are nearest-neighbor interactions chosen according to a Gaussian distribution with zero mean and standard deviation unity. Periodic boundary conditions are applied. For the Viana-Bray model each spin is connected with $`z=6`$ other spins on average chosen randomly (but with the constraint that the total number of bonds is exactly $`3N`$). We allowed the local coordination to fluctuate which is different from the more familiar Viana-Bray model in which each site has exactly the same coordination number, but we expect the properties of the two models to be very similar. The width of the Gaussian distribution is again unity, and the range of sizes is $`N700`$.

Our attention will focus primarily on two quantities: the spin overlap, $`q`$, defined by

$$q=\frac{1}{N}\underset{i=1}{\overset{N}{}}S_i^{(1)}S_i^{(2)},$$
(2)
where “$`(1)`$” and “$`(2)`$” refer to two copies (replicas) of the system with identical bonds, and the link overlap, $`q_l`$, defined by

$$q_l=\frac{1}{N_b}\underset{i,j}{}S_i^{(1)}S_j^{(1)}S_i^{(2)}S_j^{(2)}.$$
(3)
In the last equation, $`N_b`$ is the number of bonds ($`Nz/2`$ for the models considered here, where $`z`$ is the coordination number), and the sum is over all pairs of spins $`i`$ and $`j`$ which are connected by bonds. The advantage of calculating $`q_l`$ as well as $`q`$ is that if two spin configurations differ by flipping a large cluster then $`q`$ differs from unity by an amount proportional to the volume of the cluster while $`q_l`$ differs from unity by an amount proportional to the surface of the cluster.

## II Equilibration

Simulations of spin glasses at low temperatures are now possible, at least for modest sizes, using the parallel tempering Monte Carlo method. In this technique, one simulates several identical replicas of the system at different temperatures, and, in addition to the usual local moves, one performs global moves in which the temperatures of two replicas (with adjacent temperatures) are exchanged. It turns out to be straightforward to design an algorithm which satisfies the detailed balance condition, and it will also have a good acceptance ratio if the temperatures are fairly close together. In this way, the temperature of a given replica wanders up and down in a random manner, and each time the temperature goes low the system is likely to end up in a different valley of the energy landscape. Thus different valleys are sampled in much less time than it would take for the system to fluctuate between valleys if the temperature stayed fixed.

We choose a set of temperatures $`T_i,i=1,2,\mathrm{},N_T`$, in order that the acceptance ratio for the global moves is satisfactory, typically greater than about $`0.3`$. Since, at each temperature, we need two copies of the system to calculate $`q`$ and $`q_l`$ as shown in Eqs. (2) and (3), we actually run 2 sets of $`N_T`$ replicas and perform the global moves independently in each of these two sets.

To believe results of simulations carried out at low temperatures it is essential to have a sound criterion for equilibration. The technique pioneered by Ref.  does not work with parallel tempering Monte Carlo because the temperature does not stay constant. However, another method can be used for a Gaussian distribution of exchange interactions, which is a common and convenient choice. It depends on an identity first noted a long time ago by Bray and Moore  for the SK model. Here we give the corresponding result for the short range case. We start with the expression for the average energy per site,

$$U=\frac{1}{N}\underset{i,j}{}[J_{ij}S_iS_j__T]_{av},$$
(4)
where $`\mathrm{}_T`$ denotes the Monte Carlo average for a given set of bonds, and $`[\mathrm{}]_{av}`$ denotes an average over the (Gaussian) bonds $`J_{ij}`$. One can perform an integration by parts over the $`J_{ij}`$ to relate $`U`$ to the average link overlap defined in Eq. (3), i.e.

$$q_l\frac{1}{N_b}\underset{i,j}{}[S_iS_j__T^2]_{av}=1\frac{T|U|}{(z/2)J^2},$$
(5)
where the brackets $`\mathrm{}`$ indicate both a Monte Carlo average and and an average over disorder, $`J^2`$ is the variance of the interactions (set equal to unity in this paper), the sum is over sites $`i`$ and $`j`$ connected by bonds (each pair counted once), and the factor of $`z/2`$ arises because there are $`z/2`$ times as many bonds as sites. A very similar approach has been used to test equilibration of the parallel tempering method for the SK model, except that in that case the square of the spin glass order parameter appears rather than $`q_l`$.

We start the simulation by randomly choosing the spins in the $`2N_T`$ replicas to be uncorrelated with each other. This means that the two sides of Eq. (5) should approach the equilibrium value from opposite directions for the following reason. The data for $`q_l`$ will be too small if the system is not equilibrated because the random start means that the spins are initially further away from each other in configuration space than they will be in equilibrium, whereas initially the energy will not be as negative as in equilibrium so the right hand side of Eq. (5) will initially be too high. Hence we expect that if the two sides of Eq. (5) agree then the system is in equilibrium.

For illustration purposes we show in Fig. 1 how this works for three-dimensions with $`L=8`$ and $`T=0.5`$ (to be compared with $`T_c0.95`$). The data for $`q_l`$ increases as the length of the simulation increases while that determined from the energy decreases (by a lesser amount). Once the two agree they do not appear to change at longer times, indicating that they have reached equilibrium, as expected. Furthermore, the data for different moments of the spin overlap $`q`$ as well as the link overlap $`q_l`$ appear to saturate when $`q_l`$ has equilibrated so it does not appear that there are longer relaxation times for $`q`$ than for $`q_l`$. We also checked that the whole distribution $`P(q)`$ does not change with time once the two estimates for $`q_l`$ agree. By presenting results for an intermediate temperature, rather than the lowest temperature, in Fig. 1, we show that the results do not change for times longer than that needed for $`q_l`$ and $`12T|U|/z`$ to agree. The length of the simulation was chosen so that, at lowest temperature, the data for $`q_l`$ and $`12T|U|/z`$ just converged.

## III Results

### A Three Dimensions

In Table I, we show $`N_{\mathrm{samp}}`$, the number of samples, $`N_{\mathrm{sweep}}`$, the total number of sweeps performed by each set of spins (replicas), and $`N_T`$, the number of temperature values, used in the 3D simulations. For each size, the largest temperature is 2.0 and the lowest temperature is 0.1. However, for $`L=8`$, the data for the two lowest temperatures, $`T=0.10`$ and 0.15, are not fully equilibrated so data at these temperatures has been ignored and the lowest temperature used in the analysis is 0.20. This is to be compared with $`T_c0.95`$. The set of temperatures is determined by requiring that the acceptance ratio for global moves is satisfactory for the largest size, $`L=8`$, and for simplicity the same temperatures are also used for the smaller sizes. For $`L=3,4,5`$, and 6, the acceptance ratio for global moves is greater than about 0.6 in average and is always greater than 0.3 for each pair of temperatures. For $`L=8`$ the average acceptance ratio is 0.41, and the lowest value is 0.12.

In Table II we compare the average total energy $`E=NU`$ at $`T=0.20`$ with the average ground state energy obtained by finding the ground state of each sample with a hybrid genetic algorithm, as discussed elsewhere. The two energies are very close together, indicating that at this temperature our data are unlikely to be affected by the critical point.

Figs. 2 and 3 show (symmetrized) data for $`P(q)`$ at temperatures 0.20 and 0.50. There is clearly a peak for large $`q`$ and a tail down to $`q=0`$. At both temperatures one sees that the tail in the distribution is essentially independent of size. A more precise determination of the size dependence of $`P(0)`$ is shown in Fig. 4 where, to improve statistics, we average over the (discrete) $`q`$-values with $`|q|<q_{}`$, with $`q_{}=0.20`$. If general, we expect that $`P(0)L^\theta ^{}`$, where we allow $`\theta ^{}`$ to be different from $`\theta `$, the latter being obtained from boundary condition changes. In the droplet picture $`\theta ^{}=\theta `$. The dashed line in Fig. 4 has slope $`0.20`$ corresponding to the estimated value of $`\theta `$, so in the droplet picture the data is expected to follow a track parallel to this line. The actual size dependence is clearly much weaker than this, and consistent with a constant $`P(0)`$, which implies that the energy to create a large excitation does not increase with size, and therefore $`\theta ^{}0`$.

More precisely, a two-parameter fit of the data in Fig. 4 with the form $`aL^\theta ^{}`$, gives $`\theta ^{}=0.01\pm 0.05`$ for $`T=0.20`$, $`\theta ^{}=0.02\pm 0.02`$ for $`T=0.34`$, and $`\theta ^{}=0.01\pm 0.02`$, for $`T=0.50`$. Different values of $`q_{}`$ give similar results for the fits. We also tried a one-parameter fit in which $`\theta ^{}`$ is fixed. Assuming $`\theta ^{}=0`$ the goodness-of-fit parameter Q is $`0.34,0.77`$ and $`0.57`$ for $`T=0.20,0.34`$ and $`0.50`$ respectively, whereas assuming $`\theta ^{}=0.20`$ the goodness-of-fit parameters are very low, $`1.4\times 10^4,2.8\times 10^6`$ and $`7.6\times 10^{15}`$ respectively. Hence, just considering statistical errors for the sizes studied, the data is compatible with $`\theta ^{}=0`$ and not with $`\theta ^{}=0.20`$.

Refs.  studied $`P(0)`$ by the Migdal-Kadanoff approximation, which is known to yield the droplet picture asymptotically. They find that, although the behavior of $`P(0)`$ at higher temperatures is masked by critical point effects, data at low temperatures, such as those considered here, should show the droplet behavior. That we find quite different results indicates that the Migdal-Kadanoff approximation is not applicable to such small sizes. However, our data still do not rule out the possibility that the droplet theory, or some other theory, might be correct at larger sizes.

In Fig. 5 we show the data for $`P(0)`$ versus $`T`$. We see an approximately linear decrease of the data as $`T0`$. Note though, that there is some non-linearities as shown in the figure’s inset.

Fig. 6 shows the distribution of the link overlap $`q_l`$ at $`T=0.20`$. We see that there is a large peak at $`q_l`$ close to unity (with structure coming from the allowed discrete values of $`q_l`$) and a much weaker peak (note the logarithmic vertical scale) for smaller $`q_l`$ which grows slowly with increasing $`L`$ and moves to larger values of $`q_l`$. We will refer to this feature again below when we discuss the Viana-Bray model.

The variance of $`P(q_l)`$ is shown in Fig. 8 for several low temperatures. The data is consistent with a power law decrease to zero, i.e.

$$\mathrm{Var}(q_l)L^{\mu _l},$$
(6)
where the power $`\mu _l`$ seems to vary somewhat with $`T`$. This variation is probably due to corrections to scaling coming from the shift with $`T`$ of the complicated peak structure at large $`q`$ seen in Fig. 6.

The asymptotic value of $`\mu _l`$ is related to the exponents $`\theta ^{}`$ and $`d_s`$ that have been mentioned earlier. To see this, assume that the non-zero variance is largely due to the excitation of a single large cluster of size of order $`L`$. Its energy is of order $`L^\theta ^{}`$, and the probability that thermal fluctuations can create it is of order $`T/L^\theta ^{}`$, assuming a constant density of states for these excitations. One minus the link overlap between the two states is of order $`L^{(dd_s)}`$ because there is only a contribution to $`1q_l`$ from the surface of the cluster. Hence there is a probability $`TL^\theta ^{}`$ of getting a $`\delta q_l`$ of order $`L^{(dd_s)}`$, so the variance goes like $`TL^{\mu _l}`$ where

$$\mu _l=\theta ^{}+2(dd_s).$$
(7)
An extrapolation of our results for $`\mu _l`$ to $`T=0`$ gives $`0.76\pm 0.03`$. Assuming $`\theta ^{}=0`$ this implies

$$dd_s=0.38\pm 0.02$$
(8)
which is consistent with the $`T=0`$ result of PY, namely $`dd_s=0.42\pm 0.02`$.

We have also looked at more general fits of the form

$$\mathrm{Var}(q_l)=a+\frac{b}{L^c},$$
(9)
to see to what extent the data can rule out a non-zero value of $`a`$. This is of interest because a non-zero value for $`a`$ is required by standard replica symmetry breaking theory. We carried out fits of the form

$$\mathrm{ln}[\mathrm{Var}(q_l)a]=\mathrm{ln}bc\mathrm{ln}L,$$
(10)
in which $`a`$ is fixed, and $`\mathrm{ln}b`$ and $`c`$ are the fit parameters. The $`\chi ^2`$ of the fit is then determined as a function of $`a`$ and the results are shown in Fig. 9.

One sees that the minimum of $`\chi ^2`$ is for $`a=0`$ and the range of $`a`$ in which $`\chi ^2`$ has increased by less than unity relative to the $`a=0`$ value is $`a<5.3\times 10^4`$ for $`T=0.50`$ and $`a<1.3\times 10^3`$ for $`T=0.34`$. The width of the distribution of $`q_l`$ is $`\sqrt{a}`$ which has values $`0.023`$ and $`0.036`$ respectively. Thus while, our data cannot rule out a non-zero value for the width of the distribution of $`q_l`$ in the thermodynamic limit, it does suggest that this value, if non-zero, must be very small. We note, however, that a rather small value of $`a`$ is not unreasonable in RSB. For example, if $`P(q_l)`$ consists of two delta functions at a distance of 0.1, whose weights are 0.1 and 0.9 respectively, then the value of $`a`$ is 0.0009.

### B Four Dimensions

In four dimensions we present results down to a temperature of $`0.20`$, compared with $`T_c1.80`$. Parameters of the simulations are shown in Table III. For each size, the largest temperature is 2.80 and the lowest is 0.20. The acceptance ratio for global moves is always greater than about $`0.5`$ for $`L=3`$, about $`0.2`$ for $`L=4`$ and about $`0.3`$ for $`L=5`$. In Table IV we compare the average energy at $`T=0.2`$ with the average ground state energy obtained with the hybrid genetic algorithm of Ref. . As in 3D, the data indicate that the system is very close to the ground state.

Figs. 10 and 11 show data for $`P(q)`$ for temperatures 0.20 and 0.46. As in three-dimensions, the tail in the distribution is essentially independent of size. We display the full $`P(q)`$, rather than just the symmetric part as in 3D, in order to show that it has a symmetric form as expected (a symmetric form was also obtained in 3D).

The size dependence of $`P(0)`$, averaged over the range $`|q|<q_{}`$ with $`q_{}=0.2`$ is shown in Fig. 12. The dashed line has slope $`0.70`$ corresponding to the estimated value of $`\theta `$. In the droplet picture, the behavior should follow this form asymptotically. Clearly it does not for this small range of sizes. More precisely, performing a similar analysis as in 3D we find $`\theta ^{}=0.10\pm 0.12,0.08\pm 0.09`$ and $`0.17\pm 0.06`$ for $`T=0.20,0.32`$ and $`0.46`$ respectively. For the same temperatures, the goodness-of-fit parameter is $`0.67,0.65`$ and $`0.014`$, assuming $`\theta ^{}=0`$, which is acceptable, while assuming $`\theta ^{}=0.70`$ the goodness-of-fit parameters are tiny: $`2.1\times 10^{10},2.1\times 10^{12}`$ and $`4.2\times 10^{18}`$.

As in 3D we find different results from what is predicted by Refs.  at such low temperatures on the basis of the the Migdal-Kadanoff approximation. However, our data can not rule out the possibility that some other behavior may occur at larger sizes.

Note also, that our data for $`P(0)`$ decreases approximately linearly with temperature as $`T0`$, as shown in Figure 13.

Figs. 14 and 15 show the distribution of the link overlap $`q_l`$ at temperatures 0.20 and 0.46. As in 3D we see complicated structure at large $`q_l`$ and a subsidiary peak at smaller $`q_l`$ which grows with increasing $`L`$.

The variance of $`q_l`$ is shown in Fig. 16 at several low temperatures. The data is consistent with the power law decrease to zero shown in Eq. (6). The range of sizes is so small, and the values of $`\mu _l`$ also so small, that we are not able to rule out a non-zero value for $`L\mathrm{}`$ in 4D. However, the data is consistent with the asymptotic value being zero.

An extrapolation of our effective values of $`\mu _l`$ to $`T=0`$ gives $`0.35\pm 0.06`$, which, assuming $`\theta ^{}=0`$, gives

$$dd_s=0.17\pm 0.03.$$
(11)
This is just consistent with the $`T=0`$ results of PY who find $`dd_s=0.21\pm 0.01`$. However, the quoted error bars are from statistical errors only, so the difference may be partly due to systematic effects coming from the small range of sizes studied.

### C Viana-Bray Model

For the Viana-Bray model, $`T_c`$ is given by the solution of

$$\frac{1}{\sqrt{2\pi }}_{\mathrm{}}^{\mathrm{}}e^{x^2/2}\mathrm{tanh}^2\left(\frac{x}{T_c}\right)𝑑x=\frac{1}{z1},$$
(12)
where $`z(=6\mathrm{here})`$ is the coordination number. The solution is $`T_c=1.8075\mathrm{}`$, which is roughly twice the transition temperature of the 3D short range model considered here, which has the same coordination number.

Parameters of the simulations are shown in Table V. In each case, the largest temperature is 2.6 and the lowest temperature is 0.1. For $`N=700`$ the data is not equilibrated for temperatures lower than $`0.34`$, and is almost equilibrated at $`T=0.34`$. Except for $`N=700`$ the acceptance ratio for global moves is always greater than about 0.3. For $`N=700`$ the acceptance ratio is greater than 0.3 for most temperatures but there is one “bottleneck” where the acceptance ratio went down to 0.08.

First of all, in Fig. 17 we show that $`P(q)`$ has a weight at $`q=0`$ which appears to be independent of the system size, as expected.

A plot of $`P(q_l)`$ is shown in Figs. 18 and 19 for $`T=0.34`$ and $`0.70`$. Note that the data for $`N=700`$, $`T=0.34`$ show a dip around $`q0.5`$, due to imperfect equilibration of the Monte Carlo runs for some samples (this explains also the fluctuations of the data for $`P(q)`$ in Fig. 17). As in the 3D and 4D data discussed above we see a 2-peak structure develop as the size increases. This is also clearly visible in the earlier work of Ciria et al. in 4D. For the Viana-Bray model, the position of the smaller peak shifts neither with $`L`$ or $`T`$ whereas for the 3D data, see Fig. 6, it clearly shifts to larger values of $`q_l`$ with increasing size. In 4D, see Figs. 14 and 15, the range of sizes is sufficiently small that it is difficult to tell whether there is a shift in the position of the peak or not, but if it is present, it appears to be a smaller effect than in 3D. It would be helpful to understand the physics behind the two peak structure. For the hierarchical lattice used in the Migdal-Kadanoff approximation, Bokil et al. have given an explanation, but it is not clear to us how this goes over to the models discussed here.

A plot of $`\mathrm{Var}(q_l)`$ against $`N`$ is shown in Fig. 20 for different temperatures. In contrast to the data for 3D shown in Fig. 8, (which is for a similar number of spins) the data is clearly tending to a constant at large $`N`$.

This is confirmed by the $`\chi ^2`$ analysis of the fits corresponding to Eq. (10) shown in Fig. 21. Clearly the asymptotic value of $`a`$ is large and finite. Compare this figure with Fig. 9, which shows the corresponding results in 3D, and where $`a=0`$ gives the best fit.

We conclude this section by pointing out the sizes studied (which covers a similar range to that in 3D and 4D) are sufficient to determine the correct asymptotic behavior for the Viana-Bray model.

## IV Conclusions

To conclude, Monte Carlo simulations at low (but finite) temperatures agree with earlier $`T=0`$ studies of KM and PY that there appear to be large-scale low energy excitations which cost a finite energy, and whose surface has fractal dimension less than $`d`$. However, since the sizes that we study are quite small, there could be a crossover at larger sizes to different behavior, such as the droplet theory (with $`\theta ^{}=\theta (>0)`$) or an RSB picture (where $`\theta ^{}=0,dd_s=0`$). We note, however, that our results for short range models are quite different from those of the mean-field like Viana-Bray model for samples with a similar number of spins, and, furthermore, our results for the Viana-Bray model do predict the correct asymptotic behavior for that model.

###### Acknowledgements.

We would like to thank D. S. Fisher, G. Parisi, E. Marinari, O. Martin, M. Mézard, M. A. Moore, H. Bokil and A. J. Bray for helpful discussions and correspondence. This work was supported by the National Science Foundation under grants DMR 9713977 and 0086287. The numerical calculations were made possible by a grant of time from the National Partnership for Advanced Computational Infrastructure, and by use of the UCSC Physics graduate computing cluster funded by the Department of Education Graduate Assistance in Areas of National Need program.


format: text

text:
# Accretion column disruption in GX 1+4

## 1 Introduction

The X-ray pulsar GX 1+4 is unusual in several respects. It is the only confirmed neutron star in a symbiotic binary Belczynski et al. (2000) which appears to be accreting from the giant companions’ stellar wind. The neutron star spin period evolution exhibits long periods of almost constant spin-up or -down, punctuated by ‘torque reversals’ lasting tens of days (e.g. Nelson et al., 1997). Indirect estimates indicate a surface magnetic field strength of 2–3$`\times 10^{13}`$ G Beurle et al. (1984); Dotani et al. (1989); Mony et al. (1991); Greenhill et al. (1993); Cui (1997). These values are significantly greater than the $`10^{12}`$ G inferred for other accreting X-ray pulsars, generally from measurements of cyclotron scattering resonance features (CSRFs; e.g. Heindl et al., 2000) in the X-ray spectra. As the presence of such a strong field in GX 1+4 awaits confirmation by such methods, the possibility that the field strength may in fact be one or more orders of magnitude lower remains.

Historically X-ray pulsar spectra have generally been fitted with variations on a power law, which is well known to be a signature of unsaturated Comptonization (e.g Pozdnyakov et al., 1983). Model fitting to the available archival RXTE data show that the spectra of GX 1+4 are broadly consistent with an analytic approximation to a Comptonization continuum component combined with a gaussian to represent iron line emission, both attenuated by neutral absorbtion with variable column density Galloway (2000).

The archival RXTE data, which span two orders of magnitude in flux for the source, has also revealed that the X-ray lightcurves and mean pulse profiles of all but one of the observations exhibit sharp dips spanning $`\mathrm{\Delta }\varphi 0.1`$–0.2 in phase. During 1996 July the countrate dropped to $`10`$ $`\mathrm{count}\mathrm{s}^1`$ over $`6`$ h; sharp dips were observed in the mean pulse profiles even then Giles et al. (2000). Similar dips are observed in only two other X-ray pulsars, RX J0812.4-3114 Reig & Roche (1999) and A 0535+262 Čemeljić & Bulik (1998). A recent analytic and modelling study suggests that the sharp dips are signatures of eclipses of the emitting region by the accretion column (Galloway et al. 2000, in preparation). Pulse-phase spectroscopy using a continuum spectral model simulating thermal Comptonization of both GX 1+4 and RX J0812.4$``$3114 shows a significant increase in the fitted optical depth $`\tau `$ coincident with the dip phase. If the inclination angle $`i`$ of these systems are roughly equal to their magnetic colatitude $`\beta `$, such an increase may result when one of the accretion columns approaches alignment with the line of sight once each rotation period. At the phase of closest alignment, photons emitted from the polar cap must propagate a greater distance through the column to escape and ultimately reach the observer. Numerical modeling using a simplified column geometry confirms that the resulting increase in optical depth will have just the effect observed in GX 1+4 and RX J0812.4$``$3114 Galloway & Wu (2000).

## 2 Observations and Data Analysis

RXTE consists of three instruments, the proportional counter array (PCA) sensitive to photons in the energy range 2–60 keV, the high-energy X-ray timing experiment (HEXTE) covering 16–250 keV, and the all-sky monitor (ASM) which spans 2–10 keV Giles et al. (1995). Analysis of RXTE data presented in this letter was carried out using lheasoft 5.0, released 23 February 2000 by the RXTE Guest Observer Facility (GOF). The data were first screened to ensure that the pointing offset was $`<0.02\mathrm{°}`$ and the source was $`>10\mathrm{°}`$ from the sun. This introduces additional gaps to the data. Instrumental background from cosmic ray interactions and as a result of satellite passages close to the SAA are estimated using the pcabackest software which is included in the lheasoft package. The ‘bright’ source background models from the 23 February 2000 release were used, which are appropriate when the net source countrate is $`60`$ $`\mathrm{count}\mathrm{s}^1`$. Spectral fitting was undertaken using the xspec spectral fitting package version 11 Arnaud (1996). For more details of the spectral analysis see Galloway (2000).

Between 1997 May 16–20 GX 1+4 was observed daily for $`12`$ min. The phase-averaged background subtracted countrate varied from day to day, from $`100`$ $`\mathrm{count}\mathrm{s}^1`$ (with all 5 proportional counter units operating) on 17 May to $`500`$ $`\mathrm{count}\mathrm{s}^1`$ on 18 and 20 May (Figure 1). For four of the five observations the sharp dips which form the primary minimum in the pulse profile were clearly observed in the background subtracted PCA countrate. The pulse fraction is defined as

$`f_\mathrm{p}={\displaystyle \frac{F_{\mathrm{max}}F_{\mathrm{min}}}{F_{\mathrm{max}}}}\mathrm{where}`$
$`F_{\mathrm{min}}<F(\varphi _i)<F_{\mathrm{max}}i1\mathrm{}n`$
and where $`F(\varphi _i)`$ is the mean pulse profile calculated over $`n`$ bins. On May 16 and 18–20 $`f_\mathrm{p}0.8`$, as is typical for the source. On 17 May however, the sharp dips in the pulse profile were weaker or absent, resulting in a pulse fraction $`f_\mathrm{p}0.5`$ (Figure 1, second panel from top). The minima are instead broad and irregular, lacking the consistency over successive cyles observed on other days.

Spectral model fitting of the Standard-2 mode spectra averaged over each daily observation was undertaken with the Comptonization continuum model (‘compTT’ in xspec; Titarchuk, 1994) following Galloway (2000). Note that it was not possible to obtain a background spectrum for the May 20 observation due to problems with the filter file; the background spectrum from the previous days’ observation was used instead. Several of the spectral fit parameters show significant variability (Table 1). The column density $`n_\mathrm{H}`$ varied between $`10^{23}\mathrm{cm}^2`$ and $`3.2\times 10^{23}\mathrm{cm}^2`$. The source spectrum temperature $`T_0`$ was somewhat higher than is typical on 16 May but decreased to a more typical 1.2–1.3 keV by 18 May. For two of the spectra the scattering plasma temperature $`T_\mathrm{e}`$ similarly could not be constrained by the fitting algorithm, and instead was fixed at 10 keV. For the other three days $`T_\mathrm{e}`$ varied between $`6.5`$–7 and $`10`$ keV. The most dramatic variation was in the fitted optical depth $`\tau `$, which was 3.6–4 on May 16 and 18–20 but $`19`$ on May 17. The equivalent width of the Fe gaussian component was also unusually high for this spectrum, at more than 1.8 keV. These variations occurred as the luminosity varied between 1.5–$`3.6\times 10^{37}\mathrm{erg}\mathrm{s}^1`$, assuming a source distance of 10 kpc. The actual distance of the source is thought to be in the range 3–15 kpc Chakrabarty & Roche (1997).

The $`\chi _\nu ^2`$ values indicate good fits for each spectrum except on May 17. The spectrum obtained on this day does not resemble any of the other $`50`$ spectra extracted over the course of the RXTE observations between 1996 and 1997. The model fit was problematic in several respects. Initial fits with all parameters free to vary resulted in a fit value for the $`\tau _{\mathrm{max}}`$ parameter (associated with the Xe absorption edge at 4.83 keV) an order of magnitude greater than the maximum for all the other spectra. Thus the value for this parameter was instead fixed (‘frozen’) at the fit value for the spectrum extracted over the whole of the 1997 May observation. The fitted $`T_0`$ value, when free to vary, also became much larger than was consistent than for the other spectra and so was fixed at 1.3 keV. This value is typical for the source at other times Galloway (2000). While the resulting fit statistic $`\chi _\nu ^22.60`$ indicates an unacceptable fit, no improvement was found using alternative continuum models, including power law, power law with exponential cutoff and broken powerlaw. A comparable fit was obtained with a blackbody component resulting in $`\chi _\nu ^2=2.74`$ with fitted $`T_{\mathrm{bb}}=6.65_{6.54}^{6.76}`$ keV in rough agreement with $`T_\mathrm{e}`$ for the compTT fit of $`6.51_{6.40}^{6.60}`$ keV.

The data from each of the observations on May 16–20 were divided into 10 phase bins based on the pulse period $`P=126.4319(0)\pm 0.0005(5)`$ s estimated from regular BATSE monitoring of the source. The ephemeris for each day was chosen arbitrarily so that the first phase bin was centered on the primary minimum, defined as phase 0.0. Each of the resulting phase-selected spectra were fitted with the usual spectral model with at most $`T_\mathrm{e}`$, $`\tau `$ and the compTT and gaussian component normalisations free to vary. The other parameters were fixed at the fit values for the phase-averaged spectrum on the same day. This approach is essentially identical to that of Galloway et al. (2000). For the May 16 and 19 observations $`T_\mathrm{e}`$ was also fixed at the same value as for the phase-averaged spectrum (see Table 1). Variation of $`T_\mathrm{e}`$ and $`\tau `$ with phase is small except close to the primary minimum. The fitted $`\tau `$ around phase 0.0 on May 16 and 18–20 was significantly greater than the value for the phase-averaged spectrum, by 20–60 per cent. For the May 17 observation, at the much shallower primary minimum $`\tau `$ is instead marginally lower than the phase-averaged fit values, while $`T_\mathrm{e}`$ is somewhat higher (Figure 2).

## 3 Discussion

For sources with $`i\beta `$, assuming that $`i`$, $`\beta `$ and the footpoints of the accretion column on the neutron star do not change, dips due to column eclipses must always be observed while accretion continues through the columns. The apparent disappearance of the dip during the 1997 May 17 observations clearly challenges this picture. The flux at this time is more than an order of magnitude larger than the lowest level measured by RXTE; complete cessation of accretion seems unlikely. Another possibility is that the geometry of the system is changing such that the accretion column drops out of alignment briefly. The inclination angle $`i`$ may vary due to geodetic precession of the pulsar, as is suggested for some radio pulsars (e.g. Kramer, 1998). Alternatively, the magnetic colatitude $`\beta `$ may change due to evolution of the neutron star magnetic field. Variations due to the former effect, however, are only expected on timescales of years. And while the latter possibility cannot be discounted due to the lack of knowledge regarding such evolution, it seems much more likely that each of these mechanisms would cause a much longer-lasting, if not permanent, cessation of eclipses rather than what is observed in GX 1+4.

Variations in the position of the column itself are much more plausible. As the accretion rate $`\dot{M}`$ (and hence the luminosity $`L_\mathrm{X}`$) decreases, the balance of magnetic and gas pressures at the inner edge of the accretion disc means that the disk will be truncated further out from the neutron star Ghosh & Lamb (1979a, b). Field lines threading the boundary region where disk material decelerates and becomes bound to the accretion column will thus meet the neutron star at progressively higher latitudes on the star. In principle at least this mechanism may cause the column to move in and out of alignment, in which case eclipses (and dips) would be expected only over a finite range of $`L_\mathrm{X}`$. Instead what appears to be the case for GX 1+4 is that dips are normally observed at the full range of luminosities at which the source was observed by RXTE, and are perhaps absent at a particular flux level only within that range. Observations at flux levels comparable to that of 1997 May 17 however do generally exhibit sharp dips. If the accretion column has indeed been disrupted rather than moved, accretion must be occurring over a much larger region of the neutron stars’ magnetosphere to give the observed X-ray flux. That material can enter the magnetosphere of the star without being channeled along the field lines suggests either that the magnetic field strength is in actual fact low for an X-ray pulsar, or that some global magnetohydrodynamic instability is operating which can completely disrupt the flow.

The unique phase-averaged RXTE spectrum observed coincident with the unusually low pulse fraction underlines the significance of this event. That the optical depth has increased fivefold compared to values obtained before and after suggests that the source has been ‘smothered’ by hot accreting material, which presumably also contributes to the unusually strong Fe line emission. Such a large $`\tau `$ suggests that the photon spectrum will approach thermalization with the scattering material; the close agreement of the fitted blackbody temperature and the scattering plasma temperature provides corroboration for the compTT component fits during this interval. The pulse-phase spectral variation was also different to that normally observed in the source. Typically $`\tau `$ is significantly greater than the mean value around the primary minimum; this is the case for previously published results (e.g. Galloway et al., 2000) as well as the observations on May 16 and 18–20. On May 17 $`\tau `$ was instead somewhat less than the phase-averaged value around the primary minimum. These results are strongly suggestive of a substantially different distribution of matter around the source. The increase in the equivalent width also observed on May 17 was not accompanied by the usual increase in the neutral column density $`n_\mathrm{H}`$ Kotani et al. (1999); Galloway (2000), further indicating varying conditions close to the star rather than in the outer circumstellar matter.

At least one similar event has been observed previously from this source. On 1993 December 11 a pulse fraction $`35`$ per cent (20–100 keV) was measured using balloon-borne instruments, along with an unusually hard spectrum with power law energy index $`0.54\pm 0.18`$ Rao et al. (1994). For comparison, a power law fit was made to RXTE spectra from 17 May in the energy band 20–40 keV only. The resulting photon index $`1.15_{1.03}^{1.28}`$ is consistent with the fit value for the full PCA energy range, and represents a significantly harder spectrum than on 1993 December 11. Flux measurements by the Burst and Transient Source Experiment (BATSE) aboard the Compton Gamma-ray Observatory satellite indicate that following that observation the source entered a hard X-ray low state with $`L_\mathrm{X}0.1\mathrm{keV}\mathrm{cm}^2\mathrm{s}^1`$ (20–60 keV) which persisted (with brief interruptions) for the next $`250`$ d. In contrast, the May 17 observation was made during a flare lasting $`20`$ d where the flux peaked at $`0.2\mathrm{keV}\mathrm{cm}^2\mathrm{s}^1`$. The other spectra from 1997 May give fitted $`\tau `$ values at intermediate levels between the two distinct spectral states suggested by analysis of the archival RXTE observations Galloway (2000). It is possible that the disruption of the accretion column is also related to transitions between distinct spectral and luminosity states inferred from analysis of the full set of RXTE archival data Galloway (2000).

This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center, and also the BATSE pulsar group home page at http://www.batse.msfc.masa.gov.


text:
# SHORT-TIME CRITICAL DYNAMICS

## Abstract

An introductory review to short-time critical dynamics is given. From the scaling relation valid already in the early stage of the evolution of a system at or near the critical point, one derives power law behaviour for various quantities. By a numerical simulation of the system one can measure the critical exponents and, by searching for the best power law behaviour, one can determine the critical point. Critical slowing down as well as finite size corrections are nearly absent, since the correlation length is still small for times far before equilibrium is reached. By measuring the (pseudo) critical points it is also possible to distinguish (weak) first-order from second-order phase transitions.

In this talk I like to give an introductory report about the main features of short-time critical dynamics. The topic exists since about ten years ago. It will be shown that it is possible to measure the critical temperature as well as all the static and dynamic critical exponents already in the short-time regime of the evolution of statistical systems, far before the equilibrium is reached. This can be done by starting either from an initial system of very high temperature, suddenly quenched to the critical temperature, or by starting with a completely ordered state, leading to the same results. For a comprehensive review we refer to references and. It will also be shown that short-time critical dynamics can in addition be used as a tool to distinguish first-order from second order phase transitions by comparing the critical temperature obtained from the different starting conditions.

It has for long been known that statistical system for particular values of their coupling show critical behaviour. This behaviour is mainly characterised by the fact that the correlation length $`\xi _x`$ becomes infinite as well as the corresponding correlation time $`\xi _t`$. If we define the difference between the coupling and the critical one by $`\tau `$, we have

$$\xi _x\tau ^\nu ,\xi _t\tau ^{\nu z};\tau =\frac{KK_c}{K_c}.$$
(1)
Other quantities approach zero or infinity when $`\tau 0`$, e.g., the magnetisation of a spin system behaves as

$$M(\tau )(\tau )^\beta .$$
(2)
The exponents $`\nu `$, $`\beta `$ and $`z`$ are static rsp. dynamic critical exponents. Mainly due to the large correlation length and -time there exists a dynamical scaling form. We do not explicitly quote the scaling form here, but refer to it in the subsequent discussion.

In 1989 Janssen, Schaub and Schmittmann have shown, that this scenario is not only valid in equilibrium, but a scaling relation is already valid in the short-time regime of the evaluation of a critical system. The system should be prepared at a very high temperature, with a small magnetisation $`m_0`$ remaining. It is then suddenly quenched to the critical temperature and released to the dynamical evolution of model A. The authors show, using renormalisation group analysis, that after a macroscopic small time $`t_{mic}`$ a scaling form is valid. For the $`kth`$ moment of the magnetisation it reads

$$M^{(k)}(t,\tau ,L,m_0)=b^{k\beta /\nu }M^{(k)}(t/b^z,b^{1/\nu }\tau ,L/b,b^{x_0}m_0).$$
(3)
In this equation $`b`$ is an arbitrary scaling factor, $`t`$ is the time, $`\tau `$ is defined in Eq.(1), L is the linear lattice size, and $`m_0\genfrac{}{}{0pt}{}{>}{}0`$ the initial magnetisation. Except of this additional last argument the scaling form looks exactly like that in equilibrium. This argument $`m_0`$ gives rise to a new, independent critical exponent $`x_0`$, the scaling dimension of the initial magnetisation.

The questions arises how to exploit the scaling relation (3) in order to get the desired information about the critical exponents. We investigate the system numerically with Monte-Carlo simulation at a temperature at or near the critical temperature $`T_c`$, measuring the magnetisation and its second moment, $`M(t)`$ and $`M^{(2)}(t)`$, as well as the autocorrelation $`A(t)=_iS_i(t)S_i(0)`$, where a time unit is defined as a complete update of the whole lattice. Since the system is in the early stage of the evolution the correlation length is still small and finite size problems are nearly absent. Therefore we generally consider $`L`$ large enough and skip this argument. An occasional check is made to prove this assumption. We now choose the scaling factor $`b=t^{1/z}`$ so that the main $`t`$-dependence on the right is cancelled. Expanding the scaling form (3) for $`k=1`$ with respect to the small quantity $`t^{x_0/z}m_0`$, one obtains

$$M(t)m_0t^\theta F(t^{1/\nu z}\tau )=m_0t^\theta \left(1+at^{1/\nu z}\tau +O(\tau ^2)\right),$$
(4)
where $`\theta =(x_0\beta /\nu )/z`$ has been introduced. In most cases $`\theta >0`$, i.e., for $`\tau =0`$ the small initial magnetisation increases in the short-time region. Figure 1 shows the initial increase for the 3-dimensional Ising model as an example. For small $`\tau 0`$ we have expanded the function $`F(t^{1/\nu z}\tau )`$. There appear corrections to the simple power law dependent on the sign of $`\tau `$. Therefore simulating the system for values of the coupling $`K`$ in the neighbourhood of the critical point one obtains $`M_K(t)`$ with non-perfect power law behaviour, and the critical point $`K_c`$ can be obtained by interpolation. Figure 2 shows a plot of the magnetisation for the 2-dimensional 3-state Potts model for three different values of the coupling. The curve with the best power law behaviour is found for $`K_c=1.0055(8)`$, while the exact value is $`1.00505`$.

| Figure 1. Three-dimensional Ising model: Initial increase of the magnetisation for different initial magnetisations $`m_0`$, plotted vs. time on log-log scale. $`L=128`$. Result is $`\theta =0.108(2)`$ | Figure 2. Two-dimensional 3-state Potts model: Magnetisation for three different values of the coupling near the critical point. From above $`J=1.00750`$, $`1.00505`$, and $`1.0025`$. Dotted curve $`J_c=1.0055(8)`$ |
| --- | --- |

For the second moment of the magnetisation one can deduce $`M^{(2)}(t)L^d`$ ($`d`$ = dimension of the system), because the correlation length is still small in the early time region even at the critical point. Combining this with the result of the scaling form for $`\tau =0`$ and $`b=t^{1/z}`$, $`M^{(2)}(t)t^{2\beta /\nu z}M^{(2)}(1,t^{1/z}L)`$, one obtains the power law

$$M^{(2)}(t)t^{c_2}c_2=\left(d2\frac{\beta }{\nu }\right)\frac{1}{z}.$$
(5)
An example for the power law behaviour of this quantity is given in Fig. 3.

From the scaling form (3), setting again $`b=t^{1/z}`$, one derives for $`_\tau \mathrm{ln}M(t,\tau )|_{\tau =0}`$ the power law

$$_\tau \mathrm{ln}M(t,\tau )|_{\tau =0}t^{\frac{1}{\nu z}}.$$
(6)
Approximating numerically the derivative of $`M(t)`$ involves the difference of small quantities and is therefore affected more by uncertainties. However, Fig. 4 shows for the 2-dimensional 3-state Potts model as example that the data here also, for $`t>20`$, show perfect power law behaviour.

Another interesting quantity is the autocorrelation

$$A(t)=\frac{1}{L^d}\underset{i}{}S_i(t)S_i(0).$$
(7)
An analysis (for $`m_0=0`$) leads to

$$A(t)t^{c_a}c_a=\frac{d}{z}\theta .$$
(8)
The plot in Fig. 5 shows a nearly perfect power law. Again the 3-dimensional Ising model has been taken as an example.

| Figure 3. Three-dimensional Ising model: Second moment of the magnetisation, plotted vs. time on log-log scale. $`L=128`$. From the slope one gets $`c_2=0.970(11)`$ | Figure 4. Two-dimensional 3-state Potts model: $`_\tau \mathrm{ln}M(t,\tau )`$, plotted vs. time on log-log scale. Dots represent the straight line fitted to the curve within the time interval $`[20\mathrm{}100]`$ |
| --- | --- |

Care has to be taken to exclude the time region $`t<t_{mic}`$, as well as the large time region where deviations of the the power law behaviour occur due to finite size effects or due to too large fluctuations. An example for a more detailed analysis of the straight line (on log-log plot) in Fig. 5 is given in Fig. 6. The individual slope in bins between $`t`$ and $`1.5t`$ is plotted versus $`t`$. It is seen that the slope for $`t\genfrac{}{}{0pt}{}{<}{}10`$ is not yet constant and therefore the region $`t<10`$ should be excluded. It is also seen that the fluctuations within the higher bins increase although much more points enter.

Till now a completely disordered initial state has been considered as starting point, i.e., a state of very high temperature. The question arises how a completely ordered initial state evolves, when heated up suddenly to the critical temperature. In the scaling form (3) one can skip besides $`L`$, also the argument $`m_0=1`$:

$$M^{(k)}(t,\tau )=b^{k\beta /\nu }M^{(k)}(t/b^z,b^{1/\nu }\tau )(m_0=1).$$
(9)
The system is simulated numerically by starting with a completely ordered state, whose evaluation is measured at or near the critical temperature. The quantities measured are $`M(t)`$ and $`M^{(2)}(t)`$. With $`b=t^{1/z}`$ one avoids the main $`t`$-dependence on the right of eq.(9), and for $`k=1`$ one has

$$M(t,\tau )=t^{\beta /\nu z}M(\mathrm{\hspace{0.33em}1},t^{1/\nu z}\tau )=t^{\beta /\nu z}\left(1+at^{1/\nu z}\tau +O(\tau ^2)\right).$$
(10)
| Figure 5. Three-dimensional Ising model: Autocorrelation plotted vs. time on log-log scale. $`L=128`$. From the slope one gets $`c_a=1.36(1)`$ | Figure 6. Three-dimensional Ising model: slope $`c_a`$ measured between $`t`$ and $`1.5t`$ plotted vs. time on log-log scale. |
| --- | --- |

Similarly as in case of the a completely disordered start, see eq.(4), one can choose couplings in the neighbourhood of the critical point and fix from these measurements the critical point by looking for the best power law behaviour. We would like to mention that measurements starting from $`m_0=1`$ are much less affected by fluctuations, because the quantities measured are rather big in contrast to those from a random start. Figure 7 is an example of the 3-dimensional Ising model. Both the critical point as well as the slope at that point can be determined from the measurement. Figure 8 shows the decay of the magnetisation at the critical point for different lattice sizes from $`L=32`$ to $`256`$. The data points for all lattice sizes agree completely up to $`t=1000`$, only the data for $`L=32`$ show some deviation for $`t\genfrac{}{}{0pt}{}{>}{}300`$. This justifies our statement that finite size effects are unimportant in short-time critical dynamics.

Also for the cold start one derives from (9)

$$_\tau \mathrm{ln}M(t,\tau )|_{\tau =0}t^{\mathrm{\hspace{0.33em}1}/\nu z}$$
(11)
which allows to measure the ratio $`1/\nu z`$. With the magnetisation and its second moment the time dependent Binder cumulant

$$U(t)=\frac{M^{(2)}}{(M)^2}1t^{d/z}$$
(12)
is defined. From its slope one can directly measure the dynamic exponent $`z`$. For the last two quantities we give examples in Figs. 9 and 10, using the 3-dimensional Ising model.

| Figure 7. Three-dimensional Ising model: Decay of the magnetisation for initial magnetisations $`m_0=1`$, for three different values of the coupling near the critical point. | Figure 8. Three-dimensional Ising model: Decay of the magnetisation for different lattice sizes at the critical point. |
| --- | --- |

In summarising the topics discussed up to now, we would like to make the following remarks:

* From an investigation of the system from a high-temperature initial state, a new independent critical exponent $`\theta `$ can be determined.
* A determination of the critical point and of all the critical exponents is possible starting from a high-temperature or from a zero-temperature initial state. The results are the same for second-order transitions.
* In contrast to investigations in equilibrium, the correlation length is still finite in the short-time regime, therefore finite size effects are strongly reduced.
* Also the correlation time is still finite, therefore critical slowing down is absent.
* We work directly with sample averages, not with time averages.

As is well known, critical slowing down can be overcome in some cases by the successful non-local cluster algorithm. However, the cluster algorithm is not applicable, e.g., for systems with randomness or frustration. With short-time dynamic simulations we have, e.g., successfully investigated the chiral degree of freedom in the 2-dimensional fully frustrated XY model, where critical slowing down makes measurements in equilibrium extremely difficult.

| Figure 9. Three-dimensional Ising model: Logarithmic derivative of the magnetisation with respect to $`\tau `$, obtained from an ordered start. | Figure 10. Three-dimensional Ising model: Time evolution of the Binder cumulant $`U(t)`$ for different lattice sizes. |
| --- | --- |

Till now we have investigated critical systems which undergo a second order phase transition. The relaxation of initial states with very high temperature suddenly quenched to the critical temperature, or of initial states with zero temperature, suddenly heated up to the critical temperature, is measured in the short-time regime. Various quantities exhibit a power law behaviour. From measurements with couplings in the neighbourhood of the critical point we can find by interpolation the optimal power law behaviour and determine in this way the critical point. The result is the same for both initial conditions within the statistical error.

If the phase transition is first order, at least if it is weak, the behaviour should be similar to second order. By cooling down a very hot initial state, one should find a pseudo critical point $`K^{}`$, and by heating up a cold initial state another pseudo critical point $`K^{}`$ is found. We expect

$$K^{}<K_c<K^{}.$$

A difference between $`K^{}`$ and $`K^{}`$ is a clear signal for a first-order transition.

The p-state Potts model undergoes a second order phase transition for $`p4`$ at the critical point $`K_c=\mathrm{ln}(1+\sqrt{p})`$. For $`p5`$ the transition is of first order, but still weak for small $`p`$. We have investigated this model for $`p=5`$ and $`p=7`$. For the initial condition $`m_0=0`$ the second moment $`M^{(2)}(t)`$ of the magnetisation has been measured, and for initial condition $`m_0=1`$ the quantities $`M(t)`$ and $`M^{(2)}(t)`$. In case of $`p=7`$ the time interval extended up to $`t=\mathrm{6\hspace{0.33em}000}`$ for $`L=280`$ and $`560`$, while in case of $`p=5`$ a more careful analysis is necessary. Here we have chosen time intervals up to $`\mathrm{40\hspace{0.33em}000}`$ for $`L=560`$. In both cases we find a clear difference between $`K^{}`$ and $`K^{}`$. In case of $`K^{}`$, both the magnetisation and the Binder cumulant (12) have been used. Figure 11 shows the result for $`p=7`$. The results for $`p=5`$ are similar. This example shows that short-time critical dynamics can be successfully used as a tool to distinguish between first- and second order phase transitions.

Figure 11. 7-state Potts model: Determination of the pseudo-critical point $`K^{}`$ determined from $`M^{(2)}(t)`$ and of $`K^{}`$ determined from $`M(t)`$ and $`U(t)`$. For the time interval $`[t\mathrm{}t_{max}]`$ used for the fit various starting points between $`t=100`$ and $`800`$ and the maximum time have been used.


text:
# Glassiness and constrained dynamics of a short-range non-disordered spin model

## I Introduction

Understanding the nature of the low-temperature dynamics of glasses and other strongly interacting many body systems remains one of the outstanding open problems in condensed matter and statistical physics .

Consider, for example, the case of glass-forming supercooled liquids. One of the reasons for glassiness in such systems is that the rearrangements of atoms necessary for their relaxation involves activation over energy barriers. It is natural to plot the logarithm of the viscosity or equilibration time against the inverse temperature, giving a so-called Arrhenius plot, which should take a straight-line form if the barrier heights remain constant with varying temperature. For “strong” liquids , such as SiO2, the Arrhenius plot is indeed a straight line, but for “fragile” liquids it is not, following instead the empirical Vogel–Fulcher law $`\mathrm{exp}[\text{const.}/(TT_0)]`$, although other forms not displaying a finite $`T`$ singularity can be fit too . If the dynamics of fragile liquids is due to activation (which is not the only possible explanation), this implies that the energy barriers grow with decreasing temperature, presumably because they increase with the increasing size of correlated regions in the system.

Another set of interesting questions, most relevant for fragile liquids, revolves around the possible existence of an ideal (continuous) phase transition to a true thermodynamic glass state at some temperature $`T_K`$ (the Kauzmann temperature) lying below the glass temperature . Although the existence of such a transition would resolve Kauzmann’s paradox, in which the extrapolation of the configurational entropy of the supercooled liquid appears to pass below that of the crystal , and is also supported by the analogy between fragile glasses and discontinuous mean field spin-glasses , there is some evidence against it. For example, it has been shown numerically  that a thermodynamic phase transition is absent for a system of polydisperse hard-disks, a typical model glass former. It has also been argued that the extrapolations which yield a positive $`T_K`$ in fragile liquids composed of molecules of finite size are flawed .

Given the many interesting open questions in this field, it is important to find simple and if possible solvable models for the various features displayed by glasses. The discontinuous mean-field spin-glasses, such as the spherical $`p`$-spin model , provide a good example of such a model system. However, while some dynamical features of these models  are also observed in simulations of Lennard-Jones glasses , the presence of disorder in the models and, crucially, of long-range interactions limits their usefulness as models of systems which are intrinsically short-ranged and disorder-free, like hard spheres in two and three dimensions.

An alternative approach to modeling these latter systems is the “hierarchically constrained dynamics” of Palmer et al. , in which it is hypothesized that, for a strongly interacting system displaying glassy behaviour, it should be possible to describe the dynamics in terms of hierarchies of degrees of freedom, from fast to slow, independent of the presence of disorder or even frustration. These hierarchies would be weakly interacting in the energetic sense, but their dynamics would be constrained, the faster modes constraining the slower ones. It is known that the presence of kinetic constraints in the dynamics can directly induce glassiness. A good example is the kinetically constrained Ising chain , for which the Hamiltonian is trivial, but the transition rates for the flipping of individual Ising spins depend on the states of neighbouring spins. Kinetic constraints can also arise as a result of entropic barriers, as in the Backgammon model, for instance . Fragile glass behaviour has also been observed in kinetic lattice-gas models . Although interesting in their own right, these models have somewhat ad hoc dynamics, and represent only one side of the scenario discussed in Ref. . A recent comparative study of one-dimensional constrained kinetic models has been given by Crisanti et al. .

In this paper, we study in detail a model introduced by Newman and Moore , which is a form of two-dimensional Ising model with uniform short-range ferromagnetic interactions. Despite the absence of either disorder or frustration, this model displays glassy behaviour at low temperatures. The cause of this behaviour is the presence of energy barriers which grow logarithmically with the size of correlated regions. The model has a dual description in terms either of strongly interacting spins subject to simple single-spin-flip dynamics, or of free “defects” subject to a constrained dynamics. It is thus an explicit realization of the constrained dynamics scenario of Palmer et al. . Other non-disordered short-range spin systems displaying glassy features are three-dimensional Ising models with competing nearest and next-nearest neighbour interactions , or with ferromagnetic four-spin plaquette interactions .

The paper is organized as follows. In Section II we describe the model and extend the exact solution of the statics given in Ref.  by calculating all equilibrium spin correlation functions, and showing that the phase-transition to an ordered state occurs only at zero temperature. In Section III, we give our principal results, which concern the out-of-equilibrium dynamics of the system following a quench from a random configuration to low temperatures. We show that the equilibration time diverges with an exponential inverse temperature squared law, similar to that found for the asymmetrically constrained Ising chain. We study the behaviour of one-time quantities, like energy density and spin correlations and two-time quantities, like autocorrelation and response functions. We also study the out-of-equilibrium fluctuation-dissipation relations and relate them to the configurational entropy, which can be calculated exactly. Our conclusions are given in Section IV.

## II Model and Static Solution

We consider the model introduced in Ref. , which consists of Ising spins $`\sigma =\pm 1`$ on a triangular lattice with uniform short-range three-spin ferromagnetic interactions: each spin interacts only with its nearest-neighbours, and only in groups of three lying at the vertices of a downward-pointing triangle on the lattice. Note that this is distinct from the model of Baxter and Wu , which has interactions on upward-pointing triangles also (and which is not glassy). The Hamiltonian for the model is

$$H=\frac{1}{2}J\underset{mn}{}\sigma _{mn}\sigma _{m,n+1}\sigma _{m1,n+1}+\frac{1}{2}NJ,$$
(1)
where the indices $`m`$ and $`n`$ run along the unit vectors of the lattice $`\stackrel{}{a}_1\widehat{x}`$ and $`\stackrel{}{a}_2\frac{1}{2}(\widehat{x}+\sqrt{3}\widehat{y})`$. The constant in the Hamiltonian is added for convenience to make the minimum possible energy equal to zero. The model can also be formulated using the defect variables

$$\tau _{mn}\sigma _{mn}\sigma _{m,n+1}\sigma _{m1,n+1},$$
(2)
in terms of which the Hamiltonian is

$$H=\frac{1}{2}J\underset{mn}{}\tau _{mn}+\frac{1}{2}NJ.$$
(3)
On lattices of size a power of two in at least one direction, with periodic boundary conditions, there is a one-to-one correspondence between spin and defect configurations, and hence the partition function is given by

$$Z=(2e^{\frac{1}{2}\beta }\mathrm{cosh}\frac{1}{2}\beta )^N,$$
(4)
where $`N`$ is the number of spins, and we set $`J=1`$ from here on. The equilibrium energy density is then

$$\epsilon _{\mathrm{eq}}=\frac{1}{2}\left(1+\tau \right)=\frac{1}{2}(1\mathrm{tanh}\frac{1}{2}\beta ).$$
(5)
We now wish to invert Eq. (2) and write the spins as functions of the defects. It is known  that if the spins are represented by the variables $`s_{mn}\frac{1}{2}(\sigma _{mn}+1)\{0,1\}`$ the spins below an isolated defect on the lattice form a Pascal’s triangle mod 2 (i.e., a triangular array with the binomial coefficients $`\left(\genfrac{}{}{0pt}{}{n}{r}\right)`$ mod 2 as entries). This means that the value of the spin $`s_{mn}`$ is given by the superposition mod 2 of the Pascal’s triangles of all defects $`\tau _{kl}=+1`$ with $`mlkm`$ and $`ln`$. The inverse of the transformation (2) then reads

$$\sigma _{mn}=\underset{\genfrac{}{}{0pt}{}{nl}{mlkm}}{}(\tau _{kl})^{\left(\genfrac{}{}{0pt}{}{ln}{mk}\right)}.$$
(6)
With this result we can compute the equilibrium correlation functions of the spins. First, the magnetization is given by averaging (6):

$$\sigma _{mn}=(\mathrm{tanh}\frac{1}{2}\beta )^{𝒩_{mn}}$$
(7)
where

$$𝒩_{mn}=\underset{\genfrac{}{}{0pt}{}{nl}{mlkm}}{}\left(\genfrac{}{}{0pt}{}{ln}{mk}\right)\text{ mod 2}$$
(8)
is the total number of ones in the inverted Pascal’s triangle with its tip at site $`(m,n)`$. This number diverges faster than the linear size of the triangle, so in the thermodynamic limit we obtain

$$\sigma _{mn}=\{\begin{array}{cc}1& \text{for }T=0\hfill \\ 0& \text{for }T>0\text{,}\hfill \end{array}$$
(9)
for all $`(m,n)`$, implying that the system has a $`T=0`$ static phase transition.

Arbitrary correlation functions can be computed in a similar manner:

$$\sigma _{mn}\mathrm{}\sigma _{m^{}n^{}}=(\mathrm{tanh}\frac{1}{2}\beta )^{𝒩_{mn\mathrm{}m^{}n^{}}}$$
(10)
where now $`𝒩_{mn\mathrm{}m^{}n^{}}`$ is the total number of ones in the superposition mod 2 of the inverted Pascal’s triangles starting in positions $`(m,n)`$ through $`(m^{},n^{})`$. Since $`\left(\genfrac{}{}{0pt}{}{l}{0}\right)=\left(\genfrac{}{}{0pt}{}{l}{l}\right)=1`$, the left and right edges of a Pascal’s triangle contain only ones, so any superposition of two triangles always has an infinite number of ones coming from the edges. This implies that all 2-spin correlations vanish at $`T>0`$. The first non-zero correlations are those for three spins at the vertices of inverted equilateral triangles of side $`2^k`$. In this case $`𝒩=3^k`$, and the correlation is

$$C_k^{(3)}=\sigma _{mn}\sigma _{m,n+2^k}\sigma _{m2^k,n+2^k}=(\mathrm{tanh}\frac{1}{2}\beta )^{3^k}.$$
(11)
Notice that

$$C_{k+1}^{(3)}(\mathrm{tanh}\frac{1}{2}\beta )=C_k^{(3)}(\mathrm{tanh}^3\frac{1}{2}\beta ).$$
(12)
A similar scaling relation is seen in the one-dimensional Ising model.

At low temperatures the system is in a scaling region, and using Eq. (11) we find a correlation length of

$$\xi =(\mathrm{ln}\mathrm{coth}\frac{1}{2}\beta )^{\mathrm{ln}2/\mathrm{ln}3}.$$
(13)
## III Dynamics

Newman and Moore  found that the model studied here shows glassy behaviour under a single-spin-flip dynamics. They showed that following a quench from $`T=\mathrm{}`$ the system was unable to equilibrate in finite time at low enough temperature ($`T0.2`$ for $`t_{\mathrm{max}}10^9`$ in their simulations). The system also fell out of equilibrium for exponential cooling with a variety of cooling rates. The reason for this glassiness is to be found in the details of the system’s dynamics. Single-spin flips correspond to flips of the defect variables on triples of sites forming upward-pointing triangles on the lattice, and sets of such flips can be combined to flip upward-pointing triangles of side $`2^k`$ for any integer $`k`$. Any isolated such triangle with $`k>0`$ is locally stable; in order to remove it we have to cross an energy barrier of height $`k`$. Thus as the system relaxes it has to cross energy barriers which grow logarithmically with the size of the equilibrated regions.

Another way of looking at this is to observe that at low temperatures the flipping of an isolated defect is heavily suppressed, since such a flip requires the creation of two other defects and so incurs a net energy penalty. Thus, a defect requires the presence of another neighbouring defect to be able to flip, a situation highly reminiscent of the facilitated kinetic Ising model of Ref. ; the defects are non-interacting in the Hamiltonian, but their low temperature dynamics is effectively constrained.

At low temperatures excitations of linear size larger than 1 can only be annihilated via activation. From the observation that barriers grow logarithmically with linear size it is straightforward to estimate equilibration time. The rate of relaxation of an excitation of linear size $`d`$ is given by the Arrhenius formula $`\mathrm{\Gamma }(d)\mathrm{exp}(\mathrm{ln}d/T\mathrm{ln}2)`$. Thus after time $`t`$ the average linear distance between defects is $`dt^{T\mathrm{ln}2}`$. The equilibrium value of this distance at low $`T`$ is $`d_{\mathrm{eq}}\mathrm{exp}(\beta /2)`$ \[see Eq. (5)\], and hence the equilibration time is

$$t_{\mathrm{eq}}\mathrm{exp}\left(\frac{1}{2T^2\mathrm{ln}2}\right).$$
(14)
This exponential inverse temperature squared (EITS) form is similar to that obtained in Ref.  for the asymmetrically constrained Ising chain (ACIC) , except for the 2 in the denominator which is due to the two-dimensional nature of our model. The EITS form has no finite-temperature singularity, which is consistent with the fact that our model has no finite-temperature phase transition. The behaviour of the relaxation, however, is “fragile”. In fact, the EITS form can be approximated by a Vogel–Fulcher form with $`T_0\frac{2}{3}T_g`$. A system which has an analogous but “strong” behaviour is the Frenkel-Kontorova model : at low temperatures above its $`T=0`$ phase transition, the relaxation of defects is hindered by the presence of constant energy barriers, and the equilibration time diverges in an Arrhenius way.

Equation (14) clarifies the observations made in Ref. . In the case of relaxation following a quench from $`T_0=\mathrm{}`$ to $`T`$, the highest temperature at which the system was found to equilibrate within time $`t_{\mathrm{max}}=10^9`$ was $`T=0.2`$, while for $`T=0.18`$ or smaller the system was unable to reach equilibrium. This is explained by the fact that $`t_{\mathrm{eq}}(T=0.2)10^8`$, while $`t_{\mathrm{eq}}(T=0.18)10^{10}`$. In the case of exponential cooling, our expression for $`t_{\mathrm{eq}}`$ implies a maximum allowable cooling rate of $`\gamma _{\mathrm{max}}(T)1/t_{\mathrm{eq}}(T)`$ if we wish to remain in equilibrium at temperatures $`T`$ and greater. For the different annealing simulations of , the temperatures at which the energy ceases to follow the equilibrium curve are given by the solution of $`\gamma _{\mathrm{max}}(T)=\gamma `$ for the cooling rate $`\gamma `$ used.

### A One time quantities

We now consider the behaviour of one-time quantities for the model, such as energy density and spin correlations, following a quench from $`T=\mathrm{}`$ to a low temperature.

The $`T=\mathrm{}`$ configuration is random, and therefore has zero correlation length. As the system relaxes following the quench it tries to increase its correlation length to the equilibrium value, Eq. (13), which is finite but large at low temperature. In terms of the defects, the system is trying to decrease the defect density or equivalently the internal energy per spin $`\epsilon `$ from its initial value of $`\epsilon _0=\frac{1}{2}`$. In Fig. 1 we show numerical results for the evolution in time of the energy density after a quench from $`T=\mathrm{}`$ to a variety of final temperatures. The simulations were performed on a $`256\times 256`$ rhombic lattice using a Bortz–Kalos–Lebowitz continuous time algorithm . As the figure shows, after an initial temperature-independent exponential relaxation corresponding to the (barrierless) removal of pairs of neighbouring defects, the energy displays “plateaus”, which are more pronounced the lower the temperature. These plateaus are the result of the system becoming trapped in locally stable configurations. As we can see, the time taken by the system to escape these plateaus becomes larger with decreasing temperature.

As shown in the previous section, the typical length-scale in the system (i.e., the typical distance separating defects) grows with time as $`t^T`$, and in the absence of other length-scales we might thus expect the data in Fig. 1 approximately to collapse when plotted against $`t^T`$, or equivalently against $`\nu =T\mathrm{ln}t`$ on the logarithmic scales used in the figure. (A collapse of this kind was found for the ACIC in Ref. .) In the inset of Fig. 1 we show that this is indeed the case for our model. Save for the initial exponential transient, the collapse of the curves for different temperatures is very good. Moreover, the first plateau seems to extend from $`\nu =0`$ to $`1`$, the second from $`\nu =1`$ to $`2`$, and so on (a behaviour also seen in the ACIC).

At low temperatures, the relaxation of the system can also be regarded as an anomalous coarsening process. For the case considered here of periodic boundary conditions, there exists a unique minimum energy configuration of the model in which all the spins point down and there are no defects. For free boundary conditions, however, there are $`2^{2L1}`$ degenerate minima. (The spins along the bottom and right hand side of our rhombic system, for instance, may be chosen arbitrarily, with the rest being uniquely fixed by the requirement that there be no defects.) Following a quench, spins start to rearrange locally to eliminate defects, regardless of the boundary conditions, and thus form domains of the various free-boundary minima. The imperfect matching of these domains will be marked by the presence of defects.

The coarsening dynamics for $`T`$ close to zero can be studied approximately using the generating function method of Sollich and Evans . Consider a description of a configuration of the system in terms of domains containing no defects ($`\tau =1`$) bounded by sites containing defects ($`\tau =1`$). In contrast with the one-dimensional model studied in , such a description cannot be made exact for our model. Nevertheless, as far as average properties go, such as distance between defects, we can think of low energy configurations as an approximate tiling of defect-free parallelograms each delimited by a defect in, say, its top right-hand corner. This effectively maps the problem into one dimension. As we now show, this crude approximation, which allows us to apply the method of Ref. , gives reasonable results.

When $`T0`$, we have $`\epsilon _{\mathrm{eq}}0`$ and the plateaus of Fig. 1 become distinct stages in the dynamics. During stage $`n`$, all domains of linear size $`d`$ in the interval $`2^{k1}<d2^k`$ are annealed away. The time-scale for this process is $`\mathrm{O}(e^{\beta k})`$, and differs by a factor of $`e^\beta `$ from the time-scale associated with the following stage. In this limit, the master equation for the coarsening process can be solved using generating functions. If $`P_k(d)`$ is the probability distribution of lengths $`d`$ at the beginning of $`k`$-th stage, $`G_k(z)=_{2^{k1}<d}P_k(d)z^d`$ is the probability generating function for this distribution, and $`H_k(z)_{2^{k1}<d2^k}P_k(d)z^d`$ is the equivalent function for the active domains—those which will be annealed away during the current stage—the following recursive relations are obtained

$$G_{k+1}(z)=1+[G_k(z)1]\mathrm{exp}\left(H_k(z)\right).$$
(15)
Given an initial distribution $`P_0(d)`$, these equations can be solved numerically.

In Fig. 2 we show the average linear length of the domains $`d(t)=1/\sqrt{\epsilon (t)}`$ as a function of $`\nu =T\mathrm{ln}t`$ from our simulations. As $`T0`$, we expect $`d(t)`$ to become a “staircase” function, with the steps extending between integer values of $`\nu `$. The approximate heights of these steps are given by the derivative of the generating function at $`z=1`$: $`d_k=_ddP_k(d)=G_k^{}(1)`$. These steps are also plotted in the figure (dot-dashed line) and, as we can see, the agreement is quite good.

We now turn to the behaviour of the spin correlation functions as we approach equilibrium. In equilibrium, as shown above, the magnetization and all two-spin correlations vanish at any finite temperature. Moreover, at all times following a quench, the average magnetization and all equal-time 2-spin correlations vanish as well, as shown in Fig. 3. This is a consequence of the three-spin interactions: regardless of the value of a spin $`i`$, the Hamiltonian favours a neighbouring spin $`j`$ equally to be up or down, since the only interaction between the two spins also includes a third spin of unknown orientation.

In equilibrium, the first non-zero correlations are those for triplets of spins at the vertices of downward pointing triangles of linear size $`2^k`$. This is also the case in the out-of-equilibrium regime following a quench. In Fig. 4 we show the absolute value of the $`3`$-spin correlations $`|C_k^{(3)}|`$ as a function of time after a quench, for various values of $`k`$. The behaviour of these curves illustrates the nature of the stages by which relaxation takes place. Initially, correlations at all length-scales are zero. The system has no barriers to relaxation of excitations of length-scale 1, so $`|C_0^{(3)}|`$ starts to grow immediately after the quench. At $`\nu =T\mathrm{ln}t=0`$, $`|C_1^{(3)}|`$ grows exponentially fast to a first quasi-stationary value, while all $`C_{k>1}^{(3)}`$ remain zero: this first plateau corresponds to partial equilibration up to length-scales of 2. At $`\nu =1`$, the system starts to relax up to length-scales of 4, with $`C_2^{(3)}`$ becoming nonzero, while $`C_{k>2}^{(3)}=0`$. And so forth.

Within the plateaus all one-time quantities have roughly stationary values, and so it is natural to ask to what extent the system is in a quasi-equilibrium state. One answer to this question comes from examining the agreement between our non-equilibrium correlation functions and the exact relations (12) for the equilibrium correlation function. In Fig. 5 we compare $`|C_0^{(3)}|`$ with $`|C_k^{(3)}|^{(1/3^k)}`$, and find that the out-of-equilibrium correlation functions still collapse approximately when scaled appropriately, the collapse getting better as time progresses. This may indicate that the later stages can be described by an approximate equilibrium at an appropriate effective temperature. We discuss this point further in Section III C.

### B Two-time quantities

We now turn to the behaviour of two-time quantities in the out-of-equilibrium regime following a quench from high temperature. We first consider the local two-time spin autocorrelation function

$$C(t,t_w)=\frac{1}{N}\underset{mn}{}\sigma _{mn}(t)\sigma _{mn}(t_w).$$
(16)
In Fig. 6 we show $`C(t,t_w)`$ for quenches from $`T=\mathrm{}`$ to $`T=0.12`$ for three different values of the weighting time $`t_w`$, as a function of the rescaled time difference $`\nu =T\mathrm{ln}\tau `$, with $`\tau =tt_w`$. The three values of $`t_w`$ used correspond to $`\nu =0`$, $`1`$, $`2`$, i.e., to the starting points of the first three of the plateaus discussed in Section III A. In Fig. 7 we present the same correlations as a function of $`t/t_w`$.

Since the system does not reach equilibrium at $`T=0.12`$ on the timescales simulated, we do not in general expect the correlation functions to be functions of the time difference $`\tau =tt_w`$ only. In fact, as Fig. 6 shows, the behaviour of the correlation functions has a clear dependence on the waiting time. On the other hand, as we can see from Fig. 7, neither does it obey a simple aging form, scaling with $`t/t_w`$, as observed in other cases . This is because the length-scale associated with the relaxation of the spins does not grow as a simple power of time (see Figs. 2 and 4). In general, if equilibration is associated with the growth of a length-scale $`l(t)`$, then the two-time correlation functions should scale with $`l(t)/l(t_w)`$ . Behaviour of this kind was found in the autocorrelation functions for the ACIC , where the appropriate length-scale is the average distance between upward-pointing spins. In our model, however, rescaling the spin correlation functions by the distance between defects fails, and we have not been able to find a suitable length-scale such that the scaling with $`l(t)/l(t_w)`$ holds.

Another observable of interest is the response function $`\chi (t,t_w)`$, which measures the response of the spins at time $`t`$ to a small field applied at time $`t_w`$. Comparison of such a response function with the two-time correlations studied above can reveal violations of the fluctuation–dissipation theorem (FDT) which are expected in systems displaying aging. In general, in order to obtain the response function corresponding to the local correlation function calculated above, we would need to apply a field to only a single spin on the lattice, or equivalently we could apply a random field and measure the staggered response . Both of these approaches present significant numerical challenges, and require a substantial investment of CPU time in simulation to extract clean results. For our model, however, this turns out to be unnecessary because, as mentioned in Section III A, all off-diagonal spin-spin correlations vanish, implying that two-time autocorrelations of the magnetization $`mN^1_{mn}\sigma _{mn}`$ are proportional to the local spin correlations thus:

$$m(t)m(t^{})=\frac{1}{N}C(t,t^{}).$$
(17)
This means that we need only to measure the response to a uniform field.

We have performed simulations in which the Hamiltonian was perturbed with a small uniform magnetic field applied at a time $`t_w`$ following a quench to finite temperature: $`\mathrm{\Delta }H(t)=h_0\theta (tt_w)_{mn}\sigma _{mn}`$. We measured the integrated response (or relaxation function)

$$\chi (t,t_w)=\frac{1}{Nh_0}\underset{nm}{}\sigma _{nm}(t)_h,$$
(18)
which is equivalent to measuring the change in magnetization per spin of the system, divided by the strength of the applied field. We have experimented with a variety of different values for the field strength $`h_0`$; the results presented here are for $`h_0=0.02`$, which we find to be well within the linear regime. In Fig. 8 we show the measured response for $`T=0.12`$ and $`\nu _w=T\mathrm{ln}t_w=0`$ and 1. Notice that the response is not in general a monotonically increasing function of time, but instead shows cusps at times corresponding to the jumps between plateaus. For higher temperatures the cusps are less pronounced and eventually disappear (see inset).

A similar non-monotone behaviour is seen in the response function for the defects in our model, and can be understood physically by the following argument. The energy density $`\epsilon `$ plays a role for the defects equivalent to that played for the spins by the magnetization $`m`$, and the equivalent of a small applied magnetic field $`\delta h`$ is a small change in temperature $`\delta T=T^2\delta h`$. Consider then the behaviour of the energy density shown in Fig. 1. At low temperatures we expect $`\epsilon `$ to follow a staircase function, with steps extending between integer values of $`\nu `$. When the field is applied (i.e., when the temperature is changed), the perturbed energy density is another staircase function, with each step having a slightly higher (lower) value for negative (positive) field, and lasting a longer (shorter) time, the steps now extending between integer values of $`\nu (1T\delta h)`$. The corresponding integrated response is given by the difference between the perturbed and unperturbed energy densities, divided by the field. Within each plateau this difference is mostly small and roughly constant, but a sharp change is seen at the boundaries between steps. The unperturbed ($`T0`$) density jumps to the next plateau when $`\nu `$ is an integer, but the perturbed one only does so when $`\nu (1T\delta h)`$ is an integer, and thus the absolute value of the response is large during the short interval between these two events. Once the perturbed density jumps to the next plateau, the response becomes small again. This increase and decrease at integer $`\nu `$ is responsible for the humps seen in the integrated response. The argument generalizes to perturbations with a staggered random field.

The response function for the ACIC was studied in Ref. , where it was found to be a monotonically increasing function of time. However, as mentioned above, the defect representation of our model behaves very similarly to the ACIC, and we therefore conjecture that non-monotone behaviour will be seen in the ACIC also at sufficient low temperatures, for precisely the reasons given above.

### C Configurational entropy and fluctuation–dissipation relations

A configuration of our model is a local energy minimum if and only if no two defects occupy adjacent sites on the lattice. In the context of glasses, such states are known as inherent structures . A nice feature of the model is that the set of inherent structures is isomorphic to the set of allowed states of Baxter’s hard-hexagon model , allowing the distribution of inherent structures to be calculated exactly.

The grand-partition function of the hard-hexagon model on a lattice of $`N`$ sites is

$$Z_N(z)=\underset{M=0}{\overset{N}{}}g(M,N)z^M,$$
(19)
where $`z`$ is the fugacity and $`g(M,N)`$ is the number of ways of placing $`M`$ non-adjacent particles on the triangular lattice. In our case, $`M=N\epsilon `$ corresponds to the total energy of our model, so $`g_N(\epsilon )g(N\epsilon ,N)`$ gives the density of states for inherent structures with energy $`\epsilon `$. For large $`N`$, $`g(n,N)`$ is exponential in $`N`$, so that $`g_N(\epsilon )=e^{NS_c(\epsilon )}`$, where $`S_c(\epsilon )`$ is the configurational entropy of the model, i.e., the entropy density of metastable states. In this case, Eq. (19) becomes

$$Z_N(z)=𝑑\epsilon \mathrm{exp}\left[N\left(\epsilon \mathrm{ln}z+S_c(\epsilon )\right)\right].$$
(20)
For large $`N`$ the integral is dominated by the saddle point in the exponent, and we obtain

$`\epsilon (z)`$ $`=`$ $`{\displaystyle \frac{\mathrm{ln}\kappa }{\mathrm{ln}z}},`$ (21)
$`S_c[\epsilon (z)]`$ $`=`$ $`\mathrm{ln}\kappa (z)\epsilon (z)\mathrm{ln}z,`$ (22)
where $`\kappa (z)lim_N\mathrm{}[Z_N(z)]^{1/N}`$ is the partition function per site of the hard-hexagon model, which is known exactly in the thermodynamic limit . Between them, Eqs. (21) and (22) determine $`S_c(\epsilon )`$ parametrically.

At low defect densities, $`S_c(\epsilon )`$ reduces to

$$S_c(\epsilon )=\epsilon \mathrm{ln}\epsilon +\mathrm{O}(\epsilon ),$$
(23)
which is the general form for a low concentration of non-interacting point defects in an ordered structure . The configurational entropies of the disordered Ising chain, constrained Ising chains, and the Backgammon model all have this asymptotic form . Note that $`S_c`$ has infinite slope with respect to $`\epsilon `$ at $`\epsilon =0`$ where it vanishes, which is consistent with the fact that none of these models has a finite temperature phase transition. However, since $`S_c^{}`$ increases logarithmically as $`\epsilon `$ is decreased, an extrapolation of $`S_c^{}(0)`$ from the behaviour at finite $`\epsilon `$ would wrongly suggest a finite slope and therefore a finite Kauzmann temperature. It has been argued that this mechanism would rule out an ideal-glass phase transition in materials composed of limited size molecules and conventional molecular interactions .

Having calculated correlation and response functions for the spins in our model (Section III B), we can use our results to study the relation between fluctuations and responses, which in some other systems is related to the configurational entropy. The natural way to do this is by means of a parametric plot of $`\chi (t_w+\tau ,t_w)`$ versus $`C(t_w+\tau ,t_w)`$ for fixed $`t_w`$ . For a system in equilibrium at temperature $`T`$, such a plot would be linear with slope $`1/T`$, in accordance with the fluctuation–dissipation theorem. In a glassy system on the other hand, the FDT is normally violated in the out-of-equilibrium low-temperature regime as $`t_w\mathrm{}`$, and this violation encodes important information about the system’s dynamical behaviour (see Ref.  for a review). The classic example is the mean-field $`p`$-spin spin-glass , for which the FDT plot is piecewise linear: for large values of $`C`$, corresponding to relaxation of fast degrees of freedom, the plot has a slope of $`1/T`$ as it would in equilibrium, but for smaller $`C`$, corresponding to the relaxation of slow modes, it has slope $`1/T_{\mathrm{eff}}`$, where $`T_{\mathrm{eff}}>T`$ is interpreted as a effective temperature for these modes . Furthermore, in the $`p`$-spin model $`T_{\mathrm{eff}}`$ is numerically equal to the inverse of the slope of the configurational entropy at the asymptotic energy density . Similar behaviour has also been observed in more realistic model glass formers .

At sufficiently low temperatures our model is clearly out of equilibrium, and although it does not reach the long-time asymptotic regime corresponding to approach to equilibrium in our simulations, it is possible that the FDT plots can still provide information on the relaxation process, and maybe that there exist effective temperatures associated with the different time-scales in the problem . In Fig. 9 we show the FDT plot for the spin response and correlation functions. The unusual non-monotonic shape is a consequence of the non-monotonicity of the response function.

The curves should be “read” from right to left in the plot, and are composed of a sequence of segments, each associated with one of the plateaus. The starting point of each segment and the part of the curve in which the response is increasing correspond to the time the system spends within the relevant plateau. The maximum and the downward portion of the segment correspond to the transition to the next plateau. The first part of each segment has a shape similar to that found in Refs.  and  for the one-dimensional Ising model following a quench into the temperature scaling region: the FDT curve there has slope $`1/T`$ at $`C=1`$, and bends continuously with decreasing $`C`$ to reach a slope of $`1/2T`$ at $`C=0`$. (Note that this behaviour is very different from that seen in domain growth in higher dimensions .) In our model, each of the segments of the FDT plot starts with slope approximately $`1/T`$ (indicated by the dotted line in the figure), but the subsequent shape of the curve varies from one segment to another. The similarity within the plateaus to the behaviour of the one-dimensional Ising model is consistent with what was found for the statics of the model in Section II.

An important open question is whether the configurational entropy of Eqs. (21) and (22) plays any role in the out-of-equilibrium dynamics. A plausible explanation for the shape of the FDT curves of Fig. 9 is that the initial part of each segment corresponds to quasi-stationary thermal excitations of fast modes, while the latter part corresponds to slower large-scale rearrangements arising from occasional jumps between local minima before the transition to the next plateau takes place. If the slope of the latter part of each segment corresponds to an effective temperature for that segment, we might expect these temperatures to be related to the rate of change of $`S_c(\epsilon )`$ at the energy density of the associated plateau. It is noteworthy that the final slopes of each segment in Fig. 9 are roughly equal to $`1/S_c^{}(\epsilon _k)`$, with values approximately $`\frac{1}{2}`$, $`\frac{1}{3}`$, and $`\frac{1}{4}`$ for the first three plateaus. This observation is however speculative and furthermore contrasts with the conclusions of Ref.  for the ACIC, so it deserves further investigation.

## IV Conclusions

In this paper we have studied the low temperature behaviour of the glassy two-dimensional spin model with uniform ferromagnetic short-range three-spin interactions introduced by Newman and Moore. The model has a dual description in terms either of interacting spins or of free defects, the mapping between the spin and defect representations being one-to-one. This allows us to compute exactly all equilibrium correlation functions, both for spins and defects. We have also shown that the model displays no static phase transition at finite temperature.

Despite the simplicity of its statics, the model’s low-temperature single-spin-flip dynamics is highly non-trivial. A spin flip in the model corresponds to the flipping of three neighbouring defects, which implies that at low temperatures the dynamics of the defects is constrained: a defect can only be flipped if another neighbouring defect is present. This in turn implies that the relaxation of isolated defects is an activated process, and the size of the corresponding energy barriers are found to grow logarithmically with the distance between defects. The dual representation of the model in terms of strongly interacting spins with a simple dynamics, and of free defects subject to kinetic constraints, is an explicit realization of the hierarchical constrained dynamics scenario of Palmer et al. .

We have studied in detail the dynamics of the model after a quench to low temperatures. The presence of logarithmically growing barriers leads to non-Arrhenius relaxation, the equilibration time being of the exponential inverse temperature squared form $`t_{\mathrm{eq}}\mathrm{exp}(1/[2T^2\mathrm{ln}2])`$. One-time quantities such as internal energy display “plateaus” in their equilibration profiles, which correspond to the trapping of the system in local energy minima. Each plateau is associated with a specific stage in the dynamics, the $`k^{\mathrm{th}}`$ plateau corresponding to the partial equilibration of length-scales up to $`2^k`$. The behaviour of observables related to the defects is strikingly similar to that seen in the asymmetrically constrained Ising chain. For example, in the $`T0`$ limit the average distance between defects can be well approximated using the analytic methods of Sollich and Evans  which yield exact results for the asymmetrically constrained model.

We have also studied two-time quantities for the model, like spin autocorrelations and response functions. At low temperatures, the response functions have the unusual property of being non-monotonic: they display humps at exactly those times at which the system jumps between plateaus. This behaviour has also been observed in other models at times well within the activated regime, such as the constrained Ising chains and the Backgammon model in one dimension , models for two-dimensional froths , and vibrated granular media . An important open question is whether this is a generic feature of the out-of-equilibrium dynamics of activated processes.

## Acknowledgements

The authors would like to thank Jean-Philippe Bouchaud, Paul Goldbart, Jorge Kurchan, Felix Ritort, David Sherrington, and Peter Sollich for useful discussions. JPG wishes to thank the Santa Fe Institute and MEJN the Subdepartment of Theoretical Physics, University of Oxford, for their kind hospitality while this work was carried out. This work was funded in part by EC Grant No. ARG/B7-3011/94/27 and EPSRC Grant No. GR/M04426.