- Email: [email protected]
Finite-size scaling of the three-dimensional Ising model
Volume 248, number 3,4
PHYSICS LETTERS B
4 October 1990
Finite-size scaling of the three-dimensional Ising model R.A. W e s t o n Centrefor Particle Theory, Universityof Durham, Durham DH1 3LE, UK Received 12 March 1990; revised manuscript received 17 July 1990
The finite-size scaling of the spin and energy density correlation lengths is examined for an Ising model on a 128 × L 2 lattice, with L=4, 6, 8 and 10. The results are found to be consistent with a relationship of the form ~/L=A/x (where ~is the correlation length associated with the field with scaling dimension x) for antiperiodic but not for periodic boundary conditions on the L 2 lattice. These Monte Carlo results confirm the findings of Henkel, who examined the model in its extreme anisotropic limit on smaller lattices. In addition, the value of A is determined to be 0.117 (3).
1. Introduction The conformal invariance o f certain types o f statistical-mechanical systems at their critical point may be exploited in order to relate their behaviour in infinite euclidean space to that in pseudo-finite geometries (i.e. finite in one or more but not all o f the dimensions). In particular, Cardy has shown [ 1-3 ] how correlation functions o f scaling fields transform covariantly under the conformal map ~N--,SN-1 × ~. I f such a correlation function has the power law behaviour, 1 ( O ( r l ) O ( r 2 ) ) oc irl_r2lEx,
(1.1)
in 0~N, then in SN- 1× ~ it will decay along the R direction with a correlation length
~/r= 1 / x ,
(1.2)
where r is the radius o f the SN- 1 hypersphere. Strictly, this prediction only applies as the separation along the R direction tends to infinity. This prediction has been confirmed for many statistical models in two dimensions, either from the exact solution [4], from numerical evaluation o f the transfer matrix eigenvalues [ 5 ], or by direct numerical simulation of the action [ 6 ] (i.e. using a Monte Carlo algorithm to generate a set of lattice configurations with the correct probability distribution). The relationship (1.2) has not yet been confirmed in three 340
dimensions, primarily because o f the difficulty of obtaining transfer matrix solutions on lattices large enough to adequately cover the sphere. Alcaraz and H e r r m a n n [7] have examined the transfer matrix limit o f the three-dimensional Ising model on different platonic solids (i.e. the regular polyhedra which m a y be inscribed in a sphere, the largest o f which is the dodecahedron with only 20 points), and have had to conclude that no useful information about the critical behaviour of the system can be obtained using these small lattices. The way ahead may be to try and simulate the action on a lattice consisting of a large n u m b e r of repetitions of a r a n d o m lattice covering $2; an approach currently being investigated by the author. While no numerical data currently exists to either confirm or deny (1.2), in 1987 Henkel obtained some interesting and unexpected results for the finite-size scaling behaviour of the three-dimensional Ising model on another similar geometry, that o f a tor u s × R [ 8 - 1 0 ] . Henkel determined the relevant eigenvalues of the two-dimensional hamiltonian of the model (the log of the transfer matrix o f the extreme anisotropic limit o f the model). He thus obtained the correlation lengths associated with the exponential decay of the correlation functions o f both the spin field a and the local energy density e (~ being the product o f a spin with the sum o f its nearest neighhours). He found that the correlation lengths scaled with the edge size o f the lattice L, and that their ratio
0370-2693/90/$ 03.50 © 1990 - Elsevier Science Publishers B.V. (North-Holland)
Volume 248, number 3,4
PHYSICS LETTERSB
(extrapolated to infinite lattices) was [ 10 ] ~/~ =3.62(7) for periodic boundary conditions,
( 1.3 )
~,/~ =2.76(4) for antiperiodic boundary conditions.
(1.4)
The values of the spin and energy density scaling dimensions of the three-dimensional Ising model may be obtained from a finite-size scaling evaluation of the critical exponents and use of the hyperscaling relationships. They may also be evaluated by a more direct finite-size scaling analysis [9]. The values quoted by Henkel are [9] xo=0.515(9) and x,=1.42(0.02). The ratio of these values is x~/ xo= 2.76 (2). A relationship of the form
~/L=A/x,
(1.5)
is thus consistent with (1.4) for antiperiodic boundary conditions on the toms, but not with (1.3) for periodic boundary conditions. This was Henkel's result. Unlike $2× E, the geometry of a torus × ~ is not related to R 3 by a conformal transformation. And so it is not possible to map the correlation functions covariantly onto this geometry. No prediction exists for a relationship of the form of (1.5). The fact that Henkel's observations are inconsistent with such a relationship for periodic boundary conditions but consistent for antiperiodic boundary conditions is surprising. For the Ising model on the two-dimensional geometry of a cylinder (that is on S × ~ , a geometry of the form discussed above and which is conformally related to E 2) the effect of boundary conditions is well understood. While ( 1.2 ) is true for the spin and energy correlation lengths for periodic boundary conditions around the cylinder, it is no longer true for antiperiodic boundary conditions. The reason is that the spin and energy density no longer correspond to the primary fields present in the antiperiodic theory (the effect of boundary conditions on the operator content of two-dimensional theories is discussed fully by Cardy [ 11 ] ). Conversely Henkel's results are consistent with (1.5) for antiperiodic boundary conditions on his geometry, but inconsistent for periodic boundary conditions - a curious situation.
4 October 1990
The hamiltonian used by Henkel was only defined up to an arbitrary normalisation. So while he was able to determine whether the ratio of the spin and energy correlation lengths was inversely proportional to the ratio of the scaling dimensions, he was unable to find the constant of proportionality A in (1.5) (if such a relationship held). The purpose of the simulation described in this paper was thus two-fold: firstly, to test Henkel's findings on larger lattices by simulating the action, and secondly, if Henkel's observations were correct, to determine the constant of proportionality A in (1.5) for antiperiodic boundary conditions on the torus. I have obtained results on 128×42, 128×62, 128×82 and 128>( l02 lattices. The mean values of the ratios of the correlation lengths were found to be 3.7( 1 ) for the periodic lattice and 2.6(2) for the antiperiodic lattice. These values are consistent with Henkel's results (1.3) and (1.4). Using the value xo= 0.515 (9), and assuming ( 1.5 ) was true for antiperiodic boundary conditions, it was then possible to use the measured values of the spin correlation length on the antiperiodic lattice in order to determine a value for A. The value determined in this way was A = 0.117 (3). In section 2, there is a description of the numerical method used to carry out this simulation. The results are presented in section 3. Finally, some conclusions are drawn in section 4.
2. Numerical method The action for the Ising model is
S=-K E
(ij)
aitTj,
(2.1)
where tri takes the value ( + 1 ) at each lattice vertex, and (ij) indicates a sum over all nearest-neighbour pairs. This action was simulated using a heat bath algorithm [ 12 ] in which the probability of the ith spin turning over, if all the others remained fixed, was 1
Pi = 1 +exp[2Kai~s (aj) ] "
(2.2)
The sum in the denominator is over all nearest neighbours of the ith spin. A random number was chosen from a uniform distribution between 0 and 1. I f this number was less than P~ then the ith spin was flipped. 341
Volume 248, number 3,4
PHYSICS LETTERS B
If it was greater than P~ then the spin was left alone. Since the sum in (2.2) involves only nearest-neighbour interactions, it was correct to fix the spins on one non-nearest neighbour sublattice while updating all the spins on the other sublattice. Thus in practice the algorithm was implemented by calculating the nearest-neighbour sum for every point on the lattice and then using (2.2) to update all the spins on one complete non-nearest neighbour sublattice. The nearest-neighbour sums were then re-calculated and used to update the other sublattice. This complete process constituted one iteration of the lattice. The simulation was carried out at K=Kc=0.22168, the critical point of the infinite three-dimensional Ising model [ 13 ]. The spins were stored in bit form on one long vector. This wound first around one torus and then over successive toil to finally link up with itself (the lattice was taken to be periodic in the R direction, on a scale much larger than that of the torus). So the nearest-neighbour shifts could be carried out be cyclically shifting the entire spin vector by different amounts. Periodic (and antiperiodic) boundary conditions were implemented by identifying five subregions within the lattice, corresponding to sites in the bulk and those lying on the four boundaries of the tori. The cyclical shifts involved in calculating the nearest-neighbour sums were different for each of these subregions of sites. This is the fastest and most general way of implementing periodic boundary conditions on the Cambridge AMT DAP machine used for the simulation. This is a parallel machine, made up of a 64 × 64 array of microprocessors which simultaneously perform the same processing tasks on different data (it is a SIMD "Single Instruction Multiple Data" machine). It can achieve an effective rate of approximately 30 Megaflops. Correlation lengths were extracted by performing chi-squared linear fits to the approximately linear portions of the log of the correlation functions versus separation. In all cases the statistical error on the correlation function dominated the deviation from linearity expected at large separations due to the periodicity of the lattice. The connected correlation functions were calculated as averages over the set of sampled configurations in two ways. Firstly and most simply, for each sampled configuration the spins and energy densities were averaged over each torus to 342
4 October 1990
form two new vectors of the same length, N, as the number of tori along the R direction. In terms of such a vector S, the relevant connected correlation function was evaluated as
(s(n)s(O))c=(~i~=l(SiSli+nlu)) -
Si
i 1
,
(2.3)
where the brackets on the right-hand side refer to an average over the sampled configurations. Summing over the individual tori had the effect of removing non-zero momentum modes, which decay along the R direction more rapidly than the state characterised by the correlation length at infinite separation. The problem with this method of measuring the connected correlation functions became apparent when the statistical errors were estimated. The error on the average Q of N measurements Q~ was estimated by distributing the Q~ into 10 bins. Measurement Qi was placed in bin number I i l to. The error on Q was then the standard deviation of the 10 bin averages. This binning technique was designed to take into account the possible correlation of successive sampled configurations. A straight standard deviation of all the data would have underestimated the statistical error if such a correlation existed. The assumption was made that such a correlation was negligible after 10 measurements. For several test runs a more sophisticated technique was also used to estimate the errors. This technique varied the bin size, and is fully described in ref. [ 14 ]. The estimates of both techniques were similar (neither was systematically larger, and they generally differed by less than a factor of two). The simpler technique, which required much less storage space and was more suited to analysing the errors on differing numbers of measurements, was therefore used. The method described above calculated the connected correlation functions by obtaining the total and disconnected parts from averages over a set of configurations and then subtracting these two parts. The two averages were obtained by averaging over the same set of configurations. As such the displacements of the estimated values of the two parts from their true values would have been correlated. So although the errors obtained with the rebinning technique should
Volume 248, number 3,4
PHYSICSLETTERSB
have reflected the uncertainty in each of the pans, they overestimated the error in their difference. The second method of obtaining the connected correlation function and its associated error was designed in order to produce a more realistic value for the error. The estimate of the disconnected part of the correlation function may be written 2
N----i
Si
= ~-~ .. ( S , . ) ( S j . ) ,
(2.4)
where again Zi represents a sum over the N tori, and ( ) represents an average over all measurements. If N' is the number of separations A for which there is negligible correlation between the field averages on different tori, then (2.4) is
So the connected energy correlation function was evaluated by first calculating the quantity
1 ~ (SiSli+jl~,) '
(2.6)
for separations j = 0 ~ m (for some integer m). The mean of this quantity was evaluated over the range of separations j = J + 1~ m. This mean value was then subtracted from the value of (2.6) for separations 0~A. Averaging the resultant A+ 1 values over all measurements provided a direct estimate of the connected correlation function for these separations. The errors on these quantities were estimated from the standard deviation of ten bins as above. The advantage of this method is that the error in the total Green's function estimator and the error in the estimator of its asymptotic value partly cancel on each configuration. This method should produce a value for the uncertainty in the final connected correlation function which takes the correlation of the errors in its two parts into account.
3. Results The details and results of this simulation are summarised in tables 1, 2 and 3. Table 1 displays all the information about the number of equilibration iterations, number of iterations, number of measure-
4 October 1990
ments and iteration times for the different sized lattices. Tables 2 and 3 show the ratio of the correlation lengths for periodic and antiperiodic boundary conditions on these lattices. The mean values of the ratios are Rp=3.7(1 ) and RAp=2.6(2), for periodic and antiperiodic boundary conditions respectively. These ratios are consistent with the extrapolated values ofRp = 3.62 (7) and RAp = 2.76 (2) determined by Henkel from lattices with up to 52 points [ 10 ]. The errors on the correlation lengths shown in tables 2 and 3 were obtained in several ways. First, all the connected spin-spin correlation functions were evaluated using the first of the techniques described above. It was not possible to use the second more realistic technique because the spin-spin correlation function did not reach an asymptotic value (i.e. d in (2.5) was greater than the length of the lattice). The total and disconnected parts of the correlation function were obtained as averages over all the configurations. The error on the connected correlation function was simply the sum of the errors in its two parts. A chi-squared fit of the log of the connected correlation function versus separation was carried out using these errors. The correlation length was the inverse gradient obtained from this fit. Initially this same approach was taken in order to work out the energy correlation length and its associated error. However, typically all correlation of the energy density disappeared after 8 or 9 lattice spacings, leaving the total function to oscillate about its asymptotic value. The log of the correlation function only achieved linearity after 2 or 3 lattice spacings (see fig. 2). For the intermediate separation the log plots were approximately linear. However, the errors on the two parts of the correlation function obtained using the first technique were of the same order of magnitude or larger than this difference. These errors were obviously completely unrealistic in this instance. After this problem was realised, the second technique was used to evaluate the energy-energy correlation function and its associated errors for subsequent runs. The correlation lengths measured in this way are shown in the starred entries in tables 2 and 3. Both methods were used to obtain the energy correlation functions for the 128 × 6 × 6 antiperiodic lattice. The correlation functions were the same (to well within the statistical errors evaluated using the more sophisticated technique) out to a lattice spacing of at least ten. 343
Volume 248, number 3,4
PHYSICS LETTERS B
4 October 1990
Table 1 Simulation details: Periodic and antiperiodic. Edge l e n g t h L
Equilibrium iterations ( X 103)
Further iterations ( X 1 0 3 )
Iterations per second
4
100
600
4
100
6
100 100 100 100 100 100
600 900 1100 900 300 300
4 3 11 9 1 1
51
8 10
Table 2 Periodic boundary conditions.
L
~spin
(-~energy
Ratio
~,pi./L
4 6 8 10
3.07(1) 4.72(1) 6.28(1) 7.45(4)
0.85(2) 1.25(3) 1.64(3) 2.09(6)*
3.61(7) 3.8(1) 3.82(7) 3.6(1)
0.768(3) 0.787(2) 0.785(1) 0.745(4)
Table 3 Antiperiodic boundary conditions. L
~spin
(-~¢n©rg'y
Ratio
~pin/L
6 8 10
1.375(5) 1.812(4) 2.24(2)
0.50(1)* 0.68(3) 0.94(15)*
2.76(6) 2.7(1) 2.4(4)
0.229(1) 0.2265(5) 0.224(2)
In order to obtain energy correlation lengths from the data already collected using the first naive approach, chi-squared fits were carried out on the linear portions o f the log plots. All points in such a fit were given the same weight. This weight was the standard deviation of the points from the best fit line. The error on the gradient obtained from this fit was simply an estimate o f the deviation from linearity. The errors on the correlation lengths were also estimated from the starred data in tables 2 and 3 by assuming that the errors would decrease with the square root o f the number o f iterations, and increase linearly with the number o f lattice points. The errors evaluated with this latter method were generally about twice as big as those evaluated with the former method. These are the values which are quoted in the unstarred entries in tables 2 and 3. (All errors quoted in this simulation are intended to provide only rough indica344
Number of measurements ( X 103)
40 27
tions o f the statistical uncertainties involved. ) Plots o f the log o f spin and energy correlation functions for the periodic 102 lattice are shown in figs. 1 and 2. The errors on the energy correlation function were obtained with the more sophisticated method. The linear portion selected from fig. 2 was between lattice separations o f 3 and 8, that is, the first two points were excluded. The spin and energy correlation lengths evaluated by a chi-squared fit to these graphs were ~ = 7.45 (4) and 2.09 (6). In order to analyse the possible effects of the lengthwise periodicity of the lattice, a simulation was also carried out using a periodic 76 X 6 × 6 lattice and the results were compared with those shown in table 2 for the 1 2 8 X 6 X 6 periodic lattice. The smallest value o f the ratio of the edge length of the torus to the length of the lattice, in the simulations described so far, was 128/10. This is approximately the same as 76/6. So any finite length errors present with the 1 2 8 × 10X 10 lattice should also have been present with the 76 X 6 X 6 lattice. A periodic lattice was used because periodic lattices had the longer correlation lengths and so presumably also the larger finite length errors. With the same number of iterations and measurements as shown in table 1 for the 128 × 6 × 6 periodic lattice, the spin and energy density correlation lengths measured on the 76 X 6 X 6 lattice were measured to be ~s=4.71 (2) and ~ = 1.25(3). The errors on the energy correlation function were obtained using the second, more sophisticated technique described above. The correlation lengths are the same, to within the quoted statistical errors, as those for the 128 × 6 × 6 lattice (see table 2). Finite length correlations should thus be negligible for all the correlation lengths given in tables 2 and 3.
4 October 1990
PHYSICS LETTERS B
Volume 248, number 3,4 -1 -1.5 -2 -2.5 -3 -3.5
-4 -4.5 -5 -5.5 "60
a
*
h
i
t
i
2
3
4
5
6
7
separation Fig. 1. Log of spin correlation function versus lattice separation: 192 periodic lattice.
-1.5[ -2 -2.5 -3 -3.5 __0
-4 -4.5
-5.5
-60
;
i
10
15
~
separation Fig. 2. Log of energy correlation function versus lattice separation: 102 periodic lattice. 345
Volume 248, number 3,4
PHYSICS LETTERS B
Scaling has been checked, for both periodic and antiperiodic boundary conditions, by measuring ~/L for the spin correlation length. The values are shown in tables 2 and 3. There is a spread in the scaled correlation lengths o f 6% for the periodic case and 2% for the antiperiodic case. Finally if a relationship o f the form (1.5) is assumed to be correct for the antiperiodic case (which is consistent with the observations o f Henkel and myself), then the mean values o f ~/L given in table 3 may be used in order to obtain a value o f A~ x = 0.227 (1) (recall that the spin correlation length could always be obtained much more accurately than the energy correlation length - hence the small errors). Using Henkel's value of x a = 0 . 5 1 5 ( 9 ) [9], gives a value for the constant of proportionality A--0.117(3).
4. Conclusions On lattices with up to 102 points on the toms, I have measured the ratio o f the spin and energy correlation lengths. The mean values obtained for both periodic and antiperiodic boundary conditions on the torus are consistent with the results obtained by Henkel on smaller lattices of up to 52 points [ l 0 ]. For periodic and antiperiodic boundary conditions, the spin correlation length was found to scale with the edge length o f the t o m s to within accuracies o f 6% and 2% respectively. Assuming the scaling law to be o f the form (1.5) for the antiperiodic case (which is consistent with the observations o f Henkel and myself) and using x , = 0 . 5 1 5 (9), produced a value A = 0.117(3) for the constant in ( 1.5 ). There are several obvious theoretical questions that need to be addressed. Why does (1.5) not apply to the t o m s with periodic boundary conditions, as one might naively suppose? Why does it apply for the to-
346
4 October 1990
rus with antiperiodic conditions? And finally, where does this value of 0.117 (3) for the constant A come from? It is not the 1/2r~ factor that appears for the geometry S z × ~. In order to approach these theoretical problems, it will be necessary to have a fuller understanding of the operator content o f the three-dimensional model, and in particular the effect of boundary conditions on this operator content. This was the knowledge that was required in order to make sense o f the effect o f the various boundary conditions on the two-dimensional Ising model [ 11 ].
Acknowledgement I would like to thank R o n Horgan, Ian D r u m m o n d and Simon Catterall for their helpful advice, and to acknowledge the SERC for funding. This work was carried out in The Department o f Applied Mathematics and Theoretical Physics in Cambridge.
References [ 1] J.L. Cardy, in: Phase transition and critical phenomena, Vol. 11 eds. C. Domb and J.L. Lebowitz (Academic Press, London, 1986). [2] J.L. Cardy, J. Phys. A 17 (1984) L385. [3] J.L. Cardy, J. Phys. A 18 (1985) L757 [4] P. Nightingale and H. Blote, J. Phys. A 16 (1983) L657. [ 5] K.A. Penson and M. Kolb, Phys. Rev. B 29 (1984) 2854. [6] R.A. Weston, PhD Thesis, Cambridge (1989 ), unpublished; S.M. Catterall, DAMTP preprint (1989). [ 7 ] F.C. Alcarazand H.J. Herrmann, J. Phys. A 20 ( 1987) 5735. [8] M. Henkel, J. Phys. A 19 (1986) L247. [9] M. Henkel, J. Phys. A 20 (1987) 995. [ 10 ] M. Henkel, J. Phys. A 20 ( 1987) L769. [ 11 ] J.L. Cardy, Nucl. Phys. B 275 [FS 17 ] (1986 ) 200. [ 12] M. Creutz, Nucl. Phys. B (Proc. Suppl.) 4 (1988) 547. [ 13 ] G.S. Pawley, R.H. Swendsen, D.J. Wallace and K.G. Wilson, Phys. Rev. B 29 (1984) 4039. [ 14] R.A. Weston, Phys. Len. B 219 (1989) 315.
PHYSICS LETTERS B
4 October 1990
Finite-size scaling of the three-dimensional Ising model R.A. W e s t o n Centrefor Particle Theory, Universityof Durham, Durham DH1 3LE, UK Received 12 March 1990; revised manuscript received 17 July 1990
The finite-size scaling of the spin and energy density correlation lengths is examined for an Ising model on a 128 × L 2 lattice, with L=4, 6, 8 and 10. The results are found to be consistent with a relationship of the form ~/L=A/x (where ~is the correlation length associated with the field with scaling dimension x) for antiperiodic but not for periodic boundary conditions on the L 2 lattice. These Monte Carlo results confirm the findings of Henkel, who examined the model in its extreme anisotropic limit on smaller lattices. In addition, the value of A is determined to be 0.117 (3).
1. Introduction The conformal invariance o f certain types o f statistical-mechanical systems at their critical point may be exploited in order to relate their behaviour in infinite euclidean space to that in pseudo-finite geometries (i.e. finite in one or more but not all o f the dimensions). In particular, Cardy has shown [ 1-3 ] how correlation functions o f scaling fields transform covariantly under the conformal map ~N--,SN-1 × ~. I f such a correlation function has the power law behaviour, 1 ( O ( r l ) O ( r 2 ) ) oc irl_r2lEx,
(1.1)
in 0~N, then in SN- 1× ~ it will decay along the R direction with a correlation length
~/r= 1 / x ,
(1.2)
where r is the radius o f the SN- 1 hypersphere. Strictly, this prediction only applies as the separation along the R direction tends to infinity. This prediction has been confirmed for many statistical models in two dimensions, either from the exact solution [4], from numerical evaluation o f the transfer matrix eigenvalues [ 5 ], or by direct numerical simulation of the action [ 6 ] (i.e. using a Monte Carlo algorithm to generate a set of lattice configurations with the correct probability distribution). The relationship (1.2) has not yet been confirmed in three 340
dimensions, primarily because o f the difficulty of obtaining transfer matrix solutions on lattices large enough to adequately cover the sphere. Alcaraz and H e r r m a n n [7] have examined the transfer matrix limit o f the three-dimensional Ising model on different platonic solids (i.e. the regular polyhedra which m a y be inscribed in a sphere, the largest o f which is the dodecahedron with only 20 points), and have had to conclude that no useful information about the critical behaviour of the system can be obtained using these small lattices. The way ahead may be to try and simulate the action on a lattice consisting of a large n u m b e r of repetitions of a r a n d o m lattice covering $2; an approach currently being investigated by the author. While no numerical data currently exists to either confirm or deny (1.2), in 1987 Henkel obtained some interesting and unexpected results for the finite-size scaling behaviour of the three-dimensional Ising model on another similar geometry, that o f a tor u s × R [ 8 - 1 0 ] . Henkel determined the relevant eigenvalues of the two-dimensional hamiltonian of the model (the log of the transfer matrix o f the extreme anisotropic limit o f the model). He thus obtained the correlation lengths associated with the exponential decay of the correlation functions o f both the spin field a and the local energy density e (~ being the product o f a spin with the sum o f its nearest neighhours). He found that the correlation lengths scaled with the edge size o f the lattice L, and that their ratio
0370-2693/90/$ 03.50 © 1990 - Elsevier Science Publishers B.V. (North-Holland)
Volume 248, number 3,4
PHYSICS LETTERSB
(extrapolated to infinite lattices) was [ 10 ] ~/~ =3.62(7) for periodic boundary conditions,
( 1.3 )
~,/~ =2.76(4) for antiperiodic boundary conditions.
(1.4)
The values of the spin and energy density scaling dimensions of the three-dimensional Ising model may be obtained from a finite-size scaling evaluation of the critical exponents and use of the hyperscaling relationships. They may also be evaluated by a more direct finite-size scaling analysis [9]. The values quoted by Henkel are [9] xo=0.515(9) and x,=1.42(0.02). The ratio of these values is x~/ xo= 2.76 (2). A relationship of the form
~/L=A/x,
(1.5)
is thus consistent with (1.4) for antiperiodic boundary conditions on the toms, but not with (1.3) for periodic boundary conditions. This was Henkel's result. Unlike $2× E, the geometry of a torus × ~ is not related to R 3 by a conformal transformation. And so it is not possible to map the correlation functions covariantly onto this geometry. No prediction exists for a relationship of the form of (1.5). The fact that Henkel's observations are inconsistent with such a relationship for periodic boundary conditions but consistent for antiperiodic boundary conditions is surprising. For the Ising model on the two-dimensional geometry of a cylinder (that is on S × ~ , a geometry of the form discussed above and which is conformally related to E 2) the effect of boundary conditions is well understood. While ( 1.2 ) is true for the spin and energy correlation lengths for periodic boundary conditions around the cylinder, it is no longer true for antiperiodic boundary conditions. The reason is that the spin and energy density no longer correspond to the primary fields present in the antiperiodic theory (the effect of boundary conditions on the operator content of two-dimensional theories is discussed fully by Cardy [ 11 ] ). Conversely Henkel's results are consistent with (1.5) for antiperiodic boundary conditions on his geometry, but inconsistent for periodic boundary conditions - a curious situation.
4 October 1990
The hamiltonian used by Henkel was only defined up to an arbitrary normalisation. So while he was able to determine whether the ratio of the spin and energy correlation lengths was inversely proportional to the ratio of the scaling dimensions, he was unable to find the constant of proportionality A in (1.5) (if such a relationship held). The purpose of the simulation described in this paper was thus two-fold: firstly, to test Henkel's findings on larger lattices by simulating the action, and secondly, if Henkel's observations were correct, to determine the constant of proportionality A in (1.5) for antiperiodic boundary conditions on the torus. I have obtained results on 128×42, 128×62, 128×82 and 128>( l02 lattices. The mean values of the ratios of the correlation lengths were found to be 3.7( 1 ) for the periodic lattice and 2.6(2) for the antiperiodic lattice. These values are consistent with Henkel's results (1.3) and (1.4). Using the value xo= 0.515 (9), and assuming ( 1.5 ) was true for antiperiodic boundary conditions, it was then possible to use the measured values of the spin correlation length on the antiperiodic lattice in order to determine a value for A. The value determined in this way was A = 0.117 (3). In section 2, there is a description of the numerical method used to carry out this simulation. The results are presented in section 3. Finally, some conclusions are drawn in section 4.
2. Numerical method The action for the Ising model is
S=-K E
(ij)
aitTj,
(2.1)
where tri takes the value ( + 1 ) at each lattice vertex, and (ij) indicates a sum over all nearest-neighbour pairs. This action was simulated using a heat bath algorithm [ 12 ] in which the probability of the ith spin turning over, if all the others remained fixed, was 1
Pi = 1 +exp[2Kai~s (aj) ] "
(2.2)
The sum in the denominator is over all nearest neighbours of the ith spin. A random number was chosen from a uniform distribution between 0 and 1. I f this number was less than P~ then the ith spin was flipped. 341
Volume 248, number 3,4
PHYSICS LETTERS B
If it was greater than P~ then the spin was left alone. Since the sum in (2.2) involves only nearest-neighbour interactions, it was correct to fix the spins on one non-nearest neighbour sublattice while updating all the spins on the other sublattice. Thus in practice the algorithm was implemented by calculating the nearest-neighbour sum for every point on the lattice and then using (2.2) to update all the spins on one complete non-nearest neighbour sublattice. The nearest-neighbour sums were then re-calculated and used to update the other sublattice. This complete process constituted one iteration of the lattice. The simulation was carried out at K=Kc=0.22168, the critical point of the infinite three-dimensional Ising model [ 13 ]. The spins were stored in bit form on one long vector. This wound first around one torus and then over successive toil to finally link up with itself (the lattice was taken to be periodic in the R direction, on a scale much larger than that of the torus). So the nearest-neighbour shifts could be carried out be cyclically shifting the entire spin vector by different amounts. Periodic (and antiperiodic) boundary conditions were implemented by identifying five subregions within the lattice, corresponding to sites in the bulk and those lying on the four boundaries of the tori. The cyclical shifts involved in calculating the nearest-neighbour sums were different for each of these subregions of sites. This is the fastest and most general way of implementing periodic boundary conditions on the Cambridge AMT DAP machine used for the simulation. This is a parallel machine, made up of a 64 × 64 array of microprocessors which simultaneously perform the same processing tasks on different data (it is a SIMD "Single Instruction Multiple Data" machine). It can achieve an effective rate of approximately 30 Megaflops. Correlation lengths were extracted by performing chi-squared linear fits to the approximately linear portions of the log of the correlation functions versus separation. In all cases the statistical error on the correlation function dominated the deviation from linearity expected at large separations due to the periodicity of the lattice. The connected correlation functions were calculated as averages over the set of sampled configurations in two ways. Firstly and most simply, for each sampled configuration the spins and energy densities were averaged over each torus to 342
4 October 1990
form two new vectors of the same length, N, as the number of tori along the R direction. In terms of such a vector S, the relevant connected correlation function was evaluated as
(s(n)s(O))c=(~i~=l(SiSli+nlu)) -
Si
i 1
,
(2.3)
where the brackets on the right-hand side refer to an average over the sampled configurations. Summing over the individual tori had the effect of removing non-zero momentum modes, which decay along the R direction more rapidly than the state characterised by the correlation length at infinite separation. The problem with this method of measuring the connected correlation functions became apparent when the statistical errors were estimated. The error on the average Q of N measurements Q~ was estimated by distributing the Q~ into 10 bins. Measurement Qi was placed in bin number I i l to. The error on Q was then the standard deviation of the 10 bin averages. This binning technique was designed to take into account the possible correlation of successive sampled configurations. A straight standard deviation of all the data would have underestimated the statistical error if such a correlation existed. The assumption was made that such a correlation was negligible after 10 measurements. For several test runs a more sophisticated technique was also used to estimate the errors. This technique varied the bin size, and is fully described in ref. [ 14 ]. The estimates of both techniques were similar (neither was systematically larger, and they generally differed by less than a factor of two). The simpler technique, which required much less storage space and was more suited to analysing the errors on differing numbers of measurements, was therefore used. The method described above calculated the connected correlation functions by obtaining the total and disconnected parts from averages over a set of configurations and then subtracting these two parts. The two averages were obtained by averaging over the same set of configurations. As such the displacements of the estimated values of the two parts from their true values would have been correlated. So although the errors obtained with the rebinning technique should
Volume 248, number 3,4
PHYSICSLETTERSB
have reflected the uncertainty in each of the pans, they overestimated the error in their difference. The second method of obtaining the connected correlation function and its associated error was designed in order to produce a more realistic value for the error. The estimate of the disconnected part of the correlation function may be written 2
N----i
Si
= ~-~ .. ( S , . ) ( S j . ) ,
(2.4)
where again Zi represents a sum over the N tori, and ( ) represents an average over all measurements. If N' is the number of separations A for which there is negligible correlation between the field averages on different tori, then (2.4) is
So the connected energy correlation function was evaluated by first calculating the quantity
1 ~ (SiSli+jl~,) '
(2.6)
for separations j = 0 ~ m (for some integer m). The mean of this quantity was evaluated over the range of separations j = J + 1~ m. This mean value was then subtracted from the value of (2.6) for separations 0~A. Averaging the resultant A+ 1 values over all measurements provided a direct estimate of the connected correlation function for these separations. The errors on these quantities were estimated from the standard deviation of ten bins as above. The advantage of this method is that the error in the total Green's function estimator and the error in the estimator of its asymptotic value partly cancel on each configuration. This method should produce a value for the uncertainty in the final connected correlation function which takes the correlation of the errors in its two parts into account.
3. Results The details and results of this simulation are summarised in tables 1, 2 and 3. Table 1 displays all the information about the number of equilibration iterations, number of iterations, number of measure-
4 October 1990
ments and iteration times for the different sized lattices. Tables 2 and 3 show the ratio of the correlation lengths for periodic and antiperiodic boundary conditions on these lattices. The mean values of the ratios are Rp=3.7(1 ) and RAp=2.6(2), for periodic and antiperiodic boundary conditions respectively. These ratios are consistent with the extrapolated values ofRp = 3.62 (7) and RAp = 2.76 (2) determined by Henkel from lattices with up to 52 points [ 10 ]. The errors on the correlation lengths shown in tables 2 and 3 were obtained in several ways. First, all the connected spin-spin correlation functions were evaluated using the first of the techniques described above. It was not possible to use the second more realistic technique because the spin-spin correlation function did not reach an asymptotic value (i.e. d in (2.5) was greater than the length of the lattice). The total and disconnected parts of the correlation function were obtained as averages over all the configurations. The error on the connected correlation function was simply the sum of the errors in its two parts. A chi-squared fit of the log of the connected correlation function versus separation was carried out using these errors. The correlation length was the inverse gradient obtained from this fit. Initially this same approach was taken in order to work out the energy correlation length and its associated error. However, typically all correlation of the energy density disappeared after 8 or 9 lattice spacings, leaving the total function to oscillate about its asymptotic value. The log of the correlation function only achieved linearity after 2 or 3 lattice spacings (see fig. 2). For the intermediate separation the log plots were approximately linear. However, the errors on the two parts of the correlation function obtained using the first technique were of the same order of magnitude or larger than this difference. These errors were obviously completely unrealistic in this instance. After this problem was realised, the second technique was used to evaluate the energy-energy correlation function and its associated errors for subsequent runs. The correlation lengths measured in this way are shown in the starred entries in tables 2 and 3. Both methods were used to obtain the energy correlation functions for the 128 × 6 × 6 antiperiodic lattice. The correlation functions were the same (to well within the statistical errors evaluated using the more sophisticated technique) out to a lattice spacing of at least ten. 343
Volume 248, number 3,4
PHYSICS LETTERS B
4 October 1990
Table 1 Simulation details: Periodic and antiperiodic. Edge l e n g t h L
Equilibrium iterations ( X 103)
Further iterations ( X 1 0 3 )
Iterations per second
4
100
600
4
100
6
100 100 100 100 100 100
600 900 1100 900 300 300
4 3 11 9 1 1
51
8 10
Table 2 Periodic boundary conditions.
L
~spin
(-~energy
Ratio
~,pi./L
4 6 8 10
3.07(1) 4.72(1) 6.28(1) 7.45(4)
0.85(2) 1.25(3) 1.64(3) 2.09(6)*
3.61(7) 3.8(1) 3.82(7) 3.6(1)
0.768(3) 0.787(2) 0.785(1) 0.745(4)
Table 3 Antiperiodic boundary conditions. L
~spin
(-~¢n©rg'y
Ratio
~pin/L
6 8 10
1.375(5) 1.812(4) 2.24(2)
0.50(1)* 0.68(3) 0.94(15)*
2.76(6) 2.7(1) 2.4(4)
0.229(1) 0.2265(5) 0.224(2)
In order to obtain energy correlation lengths from the data already collected using the first naive approach, chi-squared fits were carried out on the linear portions o f the log plots. All points in such a fit were given the same weight. This weight was the standard deviation of the points from the best fit line. The error on the gradient obtained from this fit was simply an estimate o f the deviation from linearity. The errors on the correlation lengths were also estimated from the starred data in tables 2 and 3 by assuming that the errors would decrease with the square root o f the number o f iterations, and increase linearly with the number o f lattice points. The errors evaluated with this latter method were generally about twice as big as those evaluated with the former method. These are the values which are quoted in the unstarred entries in tables 2 and 3. (All errors quoted in this simulation are intended to provide only rough indica344
Number of measurements ( X 103)
40 27
tions o f the statistical uncertainties involved. ) Plots o f the log o f spin and energy correlation functions for the periodic 102 lattice are shown in figs. 1 and 2. The errors on the energy correlation function were obtained with the more sophisticated method. The linear portion selected from fig. 2 was between lattice separations o f 3 and 8, that is, the first two points were excluded. The spin and energy correlation lengths evaluated by a chi-squared fit to these graphs were ~ = 7.45 (4) and 2.09 (6). In order to analyse the possible effects of the lengthwise periodicity of the lattice, a simulation was also carried out using a periodic 76 X 6 × 6 lattice and the results were compared with those shown in table 2 for the 1 2 8 X 6 X 6 periodic lattice. The smallest value o f the ratio of the edge length of the torus to the length of the lattice, in the simulations described so far, was 128/10. This is approximately the same as 76/6. So any finite length errors present with the 1 2 8 × 10X 10 lattice should also have been present with the 76 X 6 X 6 lattice. A periodic lattice was used because periodic lattices had the longer correlation lengths and so presumably also the larger finite length errors. With the same number of iterations and measurements as shown in table 1 for the 128 × 6 × 6 periodic lattice, the spin and energy density correlation lengths measured on the 76 X 6 X 6 lattice were measured to be ~s=4.71 (2) and ~ = 1.25(3). The errors on the energy correlation function were obtained using the second, more sophisticated technique described above. The correlation lengths are the same, to within the quoted statistical errors, as those for the 128 × 6 × 6 lattice (see table 2). Finite length correlations should thus be negligible for all the correlation lengths given in tables 2 and 3.
4 October 1990
PHYSICS LETTERS B
Volume 248, number 3,4 -1 -1.5 -2 -2.5 -3 -3.5
-4 -4.5 -5 -5.5 "60
a
*
h
i
t
i
2
3
4
5
6
7
separation Fig. 1. Log of spin correlation function versus lattice separation: 192 periodic lattice.
-1.5[ -2 -2.5 -3 -3.5 __0
-4 -4.5
-5.5
-60
;
i
10
15
~
separation Fig. 2. Log of energy correlation function versus lattice separation: 102 periodic lattice. 345
Volume 248, number 3,4
PHYSICS LETTERS B
Scaling has been checked, for both periodic and antiperiodic boundary conditions, by measuring ~/L for the spin correlation length. The values are shown in tables 2 and 3. There is a spread in the scaled correlation lengths o f 6% for the periodic case and 2% for the antiperiodic case. Finally if a relationship o f the form (1.5) is assumed to be correct for the antiperiodic case (which is consistent with the observations o f Henkel and myself), then the mean values o f ~/L given in table 3 may be used in order to obtain a value o f A~ x = 0.227 (1) (recall that the spin correlation length could always be obtained much more accurately than the energy correlation length - hence the small errors). Using Henkel's value of x a = 0 . 5 1 5 ( 9 ) [9], gives a value for the constant of proportionality A--0.117(3).
4. Conclusions On lattices with up to 102 points on the toms, I have measured the ratio o f the spin and energy correlation lengths. The mean values obtained for both periodic and antiperiodic boundary conditions on the torus are consistent with the results obtained by Henkel on smaller lattices of up to 52 points [ l 0 ]. For periodic and antiperiodic boundary conditions, the spin correlation length was found to scale with the edge length o f the t o m s to within accuracies o f 6% and 2% respectively. Assuming the scaling law to be o f the form (1.5) for the antiperiodic case (which is consistent with the observations o f Henkel and myself) and using x , = 0 . 5 1 5 (9), produced a value A = 0.117(3) for the constant in ( 1.5 ). There are several obvious theoretical questions that need to be addressed. Why does (1.5) not apply to the t o m s with periodic boundary conditions, as one might naively suppose? Why does it apply for the to-
346
4 October 1990
rus with antiperiodic conditions? And finally, where does this value of 0.117 (3) for the constant A come from? It is not the 1/2r~ factor that appears for the geometry S z × ~. In order to approach these theoretical problems, it will be necessary to have a fuller understanding of the operator content o f the three-dimensional model, and in particular the effect of boundary conditions on this operator content. This was the knowledge that was required in order to make sense o f the effect o f the various boundary conditions on the two-dimensional Ising model [ 11 ].
Acknowledgement I would like to thank R o n Horgan, Ian D r u m m o n d and Simon Catterall for their helpful advice, and to acknowledge the SERC for funding. This work was carried out in The Department o f Applied Mathematics and Theoretical Physics in Cambridge.
References [ 1] J.L. Cardy, in: Phase transition and critical phenomena, Vol. 11 eds. C. Domb and J.L. Lebowitz (Academic Press, London, 1986). [2] J.L. Cardy, J. Phys. A 17 (1984) L385. [3] J.L. Cardy, J. Phys. A 18 (1985) L757 [4] P. Nightingale and H. Blote, J. Phys. A 16 (1983) L657. [ 5] K.A. Penson and M. Kolb, Phys. Rev. B 29 (1984) 2854. [6] R.A. Weston, PhD Thesis, Cambridge (1989 ), unpublished; S.M. Catterall, DAMTP preprint (1989). [ 7 ] F.C. Alcarazand H.J. Herrmann, J. Phys. A 20 ( 1987) 5735. [8] M. Henkel, J. Phys. A 19 (1986) L247. [9] M. Henkel, J. Phys. A 20 (1987) 995. [ 10 ] M. Henkel, J. Phys. A 20 ( 1987) L769. [ 11 ] J.L. Cardy, Nucl. Phys. B 275 [FS 17 ] (1986 ) 200. [ 12] M. Creutz, Nucl. Phys. B (Proc. Suppl.) 4 (1988) 547. [ 13 ] G.S. Pawley, R.H. Swendsen, D.J. Wallace and K.G. Wilson, Phys. Rev. B 29 (1984) 4039. [ 14] R.A. Weston, Phys. Len. B 219 (1989) 315.