A numerical study of finite-size scaling for first-order phase transitions

NUCLEAR P H VS I C S B

Nuclear Physics B 370 (1992) 773—796 North-Holland

________________

A numerical study of finite-size scaling for first-order phase transitions A. Billoire, R. Lacaze and A. Morel Service de Physique Théorique de Saclay * 91191 Gif-sur-Yvette Cedex, France Received 3 October 1991 Accepted for publication 1 November 1991

We have performed a very high statistics numerical simulation of the q = 10 and q = 7 two-dimensional Potts models, in the transition region. Results are compared to the rigorous asymptotic predictions of Borgs, Koteck~’and Miracle-Sole. The leading behavior as a function of the lattice size does agree with predictions. The expected next to leading terms are not always seen, and we observe sizable non-scaling corrections, which theory predicts to become exponentially small on large lattices.

1. Introduction With the advent of powerful high-speed computers, the Monte Carlo technique [1] has become widely used to investigate phase transitions of statistical systems. It allows us to locate the transitions, decide about their order and extract the critical exponents or the latent heat. However, the very determination of the order of the transition become hazardous if the transition is weakly first order. In such a case, a safe method to determine the order of the transition is to look at the finite-size behavior of a suitably chosen cumulant like the specific heat CV the susceptibility x~or Binder’s fourth-order cumulant BL. If the transition is first order, C V/La x/L” and BL will go, in the vicinity of the transition point, to a non zero limit as the lattice size L goes to infinity, with finite-size corrections behaving as 1 IL”. If the transition is second order, they all go to zero with non-trivial exponents. Numerical data that show that for instance the minimum of BL is approaching a non-zero limiting value, with a 1/La finite-size behavior, give a strong case for a first-order transition. If the 1 IL” behavior is not observed, however, the possibility remains of a cross-over at a larger lattice size towards a zero limiting value typical of a non first-order behavior. ,

*Laboratoire de la Direction des Sciences de la Matière du CEA, SDI 6177 CNRS. 0550-3213/92/$05.0O

© 1992 Elsevier Science Publishers B.V. All rights reserved —

A. Billoire et a!. / Finite-size scaling

774

The above expectations are based on the heuristic theory of finite-size scaling close to a first-order transition of refs. [2—51which unfortunately does not allow to compute corrections to the leading behavior. A rigorous theory of finite-size effects close to a first-order transition point has been developed in refs. [6,7], when it is possible to write the partition function in terms of a contour model with small activities. This has been applied to the q-state Potts model in ref. [81. One does not know, however, if the not so large systems that can be simulated on existing computers are in the asymptotic regime described in refs. [6—81.We have performed a very high statistics numerical simulation [9—11],in order to investigate this point, using a model which has a firstorder phase transition and for which the transition inverse temperature fit, and the internal energies of the coexisting phases are exactly known, namely the two-dimensional Potts models (with q = 10, and q = 7), on a square lattice [5,12,13]. In sect. 2 we review the predictions of finite-size scaling for the first-order transition of the two-dimensional Potts model. In sect. 3 we give details about the simulation and the data analysis. In sect. 4 we discuss the q = 10, and in sect. 5 the q = 7 results. Conclusions are drawn in sect. 6. 2. Theory The q-state Potts model [14,151 is a generalization of the Ising spin model. The spin variables a, associated with the sites of a d-dimensional lattice, can take q different values. The hamiltonian of the model is

H=—

~

ö~,.

(1)

nearest neighbor

The two-dimensional Potts model has a first-order phase transition [18] for q > 4 at the temperature

fit=ln(~’~+l).

(2)

As shown in ref. [81,the partition function Z (/3, L) of the Potts model on a square lattice of size L” with periodic boundary conditions can be written, for large enough q, as Z(fl,L)

=

exp{—/3L”f~(/3)}+ qexp{—/3L”f0(/3)} +O(e_b~~) exp{—flL”min{f~(/3),fo(fl)}},

(3)

A. Billoire et al. / Finite-size scaling

775

where b is some strictly positive number, f0 (/3) and fa (/1) are smooth functions such that for /3 < fit, fa(13) is identical with the free energy and fo (/3) > fd(/.J). For /3 > /J~,f~(/3)is identical with the free energy and fd (/3) > f0 (/3). Since f0 and fd are many times differentiable functions, one can expand them around the transition point

fifi(fl)

=

2( ~

thf~(Pt)+

+ ...

(4)

(fl-fi~E~+ (fl_flt)

It follows that the specific heat

CV

=

/32L”((E2)

—

(E)2)

(5)

has a maximum at (2)

/3(CVmax)

3”).

/3~ Ed_E

(6)

0Ld+ -~+O(l/L

The height of this maximum increases linearly with L’~’,

2~ +O(1/Ld),

CVmax _Ld~(E

(7)

0_E~)2 + CV~

whereas for fixed /3 ~ fit, CV(/3) goes to a constant, as L goes to infinity. Expressions for the coefficients and CV (2) as functions of the E 1 and C~ can be found in ref. [131. Another very interesting observable is the ratio [5] fl~T~

1l- (E4) 3~ (E2)2 BL—~

(8)

It is non-positive, and vanishes in the infinite-volume limit, for all values of /3, even at a second-order critical point, since fluctuations disappear in this limit ((E”) (E)’1). Close to a first-order phase transition, however, fluctuations between the coexisting phases prevent BL from vanishing. One finds that BL reaches a minimum equal to [7,13,19]

BLmin

=

-

(E~

+ BL~2~ + O(l/L2d)

(9)

A. Billoire et a!. / Finite-size scaling

776

at the point (10)

~

as functions of the E, and C, Expressions of the coefficients fi~and BL~2~ can be found in ref. [13]. Although the use of BLmin as an indicator of the order of phase transitions has been much publicized, C Vmax/Ld is as good an indicator indeed. The value of BL depends on the choice made of the arbitrary constant one can add to the definition of the energy. This leads to introduce the quantity *

U4

— —

((E—(E))4) ((E (E))2)2’

11

-

which is independent of such a constant. U4 is strictly larger than one, but when the energy probability distribution is

P(E)

=

~(ó(E—E 1) +ó(E—E2)),

(12)

where E1 and E2 are some constants, in which case U4 = 1. It means that U4 can reach its limiting value only at a first-order transition point, in the infinite-volume limit. For large but finite volumes, U4 reaches a minimum

2(EQ—Ed)2 + O(1/L2”).

U4min

=

(13)

1 + Ld fl~

at the point fl(U4min)

=

fit— ~(q)

~

+

~2~~d~38)

+

(or CVmax/L’~) and U4mjn are dual since BLmin (or CVmax/Ld) going to zero means second (or higher) order, whereas U4mjn going to one means first order. The above formulae for the extrema C Vmax, BLmin, U4mjn and the corresponding effective fl’s have higher power-law corrections that may hide the asymptotic behavior on lattices that can be simulated today. In contrast, the BLmin

*This ratio was suggested to us independently by S. Gupta and A. Irbäck, and by A. Sokal.

A. Billoire et a!. / Finite-size scaling

777

TABLE 1 Internal energies (E = (H) p, /L’~) and specific heats of the two coexisting phases for the two-dimensional q-state Potts model for q = 10 and 7. The “physical correlation length” ~ (fit) is also given

q 7 10

E

0 —1.55460 —1.66425

Ed —1.20133 —0.96820

E~— E0 0.35328 0.69605

Cd — C0 0.22343 0.44763

=

Cd 50±10 12.7±0.3

‘~(flt) 30 -.~ 6

expressions for bulk averages evaluated at the (infinite-volume limit) transition point fi = fit do not have power-law corrections, as a consequence of eq. (3). The average energy is given by EL(flt)

=

E~+qE0 +O(e~L),

(15)

and the value of the specific heat is CVL(flt)

=

Co1~~+qCd + (1~q)2(EoEd)flt+O(e).

(16)

The energy at fit does not depend on the lattice size, up to exponentially small corrections. This provides an efficient estimate of the transition temperature [8] by the following “two-lattice method”. One simulates lattices of increasing sizes L1 < L2 < L~< and consider fleff(Li,Li÷i), the solution of the fixed-point equation ...,

EL1(fl)

=

EL,~1(fl)

(17)

for i = 1,2 The estimate fleff(LE,Li-i-i) converges towards fit exponentially. In order to implement the above procedure, one needs to know the value of the function EL, (/3) for arbitrary fl fit. We use the spectral density method (see sect. 3) to reconstruct this function for each L from a single run performed at flMc fi (C Vm~3. To the order considered, all quantities are expressed in terms of fit and of the energies and specific heats of the two coexisting phases. The transition temperature, E0, Ed and the difference C0 Cd are known exactly [18]. Numerical values for q = 7 and 10 can be found in table 1, together with our estimate of Cd and the estimate of the “physical correlation length” ~(q), defined as the pure phase, infinite-volume limit, correlation length. It has the same value on each side of the transition, from duality, and is measured [20] ‘~

—

A. Bi!!oire et a!. / Finite-size scaling

778

to be equal to 5.9 ±0.7 at fit for the q = 10 model. This value is extrapolated to smaller q-values with the formula [211

f

~2

‘~

(18)

The physical correlation length should not be confused with the exponentially large correlation lengths, due to tunneling between vacua, that appear on finite lattices [2,161. Before going to the discussion of our results, we would like to comment on the original heuristic theory of Challa et al. [5]. In this “two gaussian peaks” model, the shape of the energy probability density is assumed to be PL(E)

2 a0 ,~ ~ L”fl?(E—E0—C0oT) 7~exPl 2C0

(19)

2l

+\/~_exP 2Cd I’ ad J. L’~fl~(E_Ed—CdöT) 1

where öT = fl_1 _fl~i•The value for the ratio ao/ad is not predicted. One can obtain the two gaussian model, in the large-L limit, from the rigorous theory of refs. [6—8],and derive the value for aO/ad: by inverse Laplace transform of eq. (3), one obtains eq. (19) with ao/ad = qexp[flL”(f~(fl) —fo(fl))}. However, it does not mean that the two gaussian peaks model is correct. It is not. As already remarked in ref. [13], the observed positions of the maxima of the two peaks of PL (E) at /3 = fit have a strong L-dependence, in disagreement with the two gaussian peaks model. Terms that are negligible in eq. (3) are responsible for this, after the inverse Laplace transform. 3. The simulation Configurations are generated with the one hit Metropolis algorithm, the trial value for a spin variable a being chosen randomly among [1,. q 1. The lattice is divided into four sub-lattices which are updated successively. We use a multi-spin coding technique in order to gain speed, sixteen independent copies of the system being simulated in parallel. The sixteen copies of a given site are updated using different pseudo-random numbers in order to avoid possible correlations between the sixteen histories. On a 502 lattice for /3 fit, q = 10 or 7, our program updates 3.0 x 106 spins per second on a CRAY X-MP EA 28. We use the CRAY supplied congruential pseudo-random number generator RANF. After a sub-lattice update, which uses 16(N2/4) successive pseudorandom numbers, the next M 4N2 pseudo-random numbers are generated, . . ,

“~

—

A. Bi!!oire et a!. / Finite-size scaling

779

where M is the first prime number larger that 4N2. These extra numbers are not used. We do this extra work in order to reduce correlations between the successive random numbers used in order to update a given lattice site. Observables are computed for each of the sixteen systems. Statistical error estimates are obtained from a jackknife analysis [221 of these sixteen data. This method provides reliable error estimates and bias reduction. We use the so called “spectral density method” [231 in order to determine the thermodynamical averages for any fi, using data taken at one or a few fl-values. This method is invaluable for our purpose. Let ~ be the lattice average of some quantity for a given value of /3, at Monte Carlo time t. An estimate of (~) for any value of /3’ is ~,

(~)

>exp[—(fl’ =

fl)L”Et] >exp[_(fih_fl)LdE,] —

~Z

(20)

,

where E 1 is the energy of the configuration #t at coupling /3. The expectation value ~ is computed separately for each of the sixteen independent systems. This enables a straightforward (jackknifed) error estimate. The spectral density method has to be used with some care. It actually works for Ifl fl’IL” not too large, otherwise only a tail of the ~ distribution contributes in the t average in eq. (20) and the results become unreliable. As an illustration let us consider the maximum of the specific heat CVmax. Each run, performed at a fi value flMc = p1, /32..., gives an estimate of C Vma~namely C Vn~ax,C V~ax... and an estimate of the statistical error on C Vmax. These estimates are combined using weighted averages to yield our final results. 2 test. We We monitor managedtheto consistency have a x2 of the estimates CV~ax... by the of x freedom in all cases (but for the roughly equal to CVr~ax, the number of degrees L = 50 lattices for q = 10) by increasing the statistics and/or restricting the interval IflMc filL” (we had to exclude some data from the weighted averages, as shown in table 2). For q = 10, L = 50, the dynamics is so slow that several 108 sweeps are not enough to have full proof results, and the two simulations made at fit and fi = 1.425 give inconsistent results. The results from these two fl-values are plotted separately in what follows. We chose not to increase the statistics reached in ref. [11] at fit and to concentrate on the point /3 = 1.425 for reasons explained below. For q = 10, we took data on lattices of sizes L = 12 to L = 50. Lattice sizes, fl-values, and the number of iterations performed can be found in table 2. We have found more efficient to concentrate our efforts, for each lattice size, on a single point close to fi (C Vmax), i.e. a point where the system spends roughly half of the time in an ordered phase and the other half in the disordered phase. —

—

780

A. Billoire et a!. / Finite-size scaling

TABLE 2 The q = 10 simulation. Data with an “E” in the fourth column are excluded from the spectral density analysis. Here, one sweep means a sweep of the 16 systems that are updated in parallel

L 12 12 12 12 12 16 20 20 20 20 20 24 36 36 44 50

The q

=

/J Nsweeps 1.3500 1000000 1.3800 1000000 1.4000 1000000 1.4100 1000000 1.4200 1000000 1.4130 42000000 1.3800 100000 1.4000 2900000 1.4100 3000000 1.4200 48000000 1.4250 2000000 1.4200 41000000 1.4240 54000000 1.4250 5000000 1.4245 65000000 1.4250 44 000 000

J E E

E E

L 12 12 12 12 12 16 20 20 20 20

fi 1.4250 fit 1.4270 1.4300 1.4400 fit fit 1.4270 1.4300 1.4400

24 fit 36 fit 36 1.4270 50

fit

Nsweeps 1000000 2000000 1000000 1000000 1000000 2000000 3000000 3000000 3000000 100000 E 4000000 10000000 5000000 20 400 000

TABLE 3 7 simulation. Here, one sweep means a sweep of the 16 systems that are updated in parallel

fi

L 16 16 24 24 32 32 48 64

Nsweeps 1000000 1000000 2000000 2000000 4000000 7000000 6000000 15000000

1.2800 fit 1.2870 fit 1.2890 fit 1.2910 1.2925

For q = 7, we took data on lattices of sizes L = 16 to L = 64. Lattice sizes, fl-values, and the number of iterations performed can be found in table 3. We have evaluated the auto-correlation functions [24] of the total spin B, and of the order parameter M defined as

S

=

~ q—i

N~~2ITh0I~,

(21)

where N,~,is the fraction of the spins of the configuration having the value a,

A. Billoire et a!. / Finite-size scaling

781

TABLE 4 2NS’ and average The q = 10 simulation. Singlet and non-singlet dynamical masses ~uSand / numbers of order—disorder flips seen during one history. Numbers between parenthesis are the estimated statistical errors on the last digit of the number on their left

L 16 20 20 24 24 36 44 50

The q

=

fi

PS

PNS

1.4130 1.4200 1.4250 1.4200 fit 1.4240 1.4245 1.4250

0.000265(5) 0.000108(1) 0.000120(8) 0.0000538(8) 0.000065(4) 0.0000078(2) 0.00000260(9) 0.00000110(9)

0.000129(1) 0.0000366(6) 0.0000133(12) 0.00000286(6) 0.0000058(6) 0.0000030(1) 0.00000117(7) 0.00000040(6)

Flips 3348(12) 1677(13) 40(2) 773(8) 32(2) 158(3) 66(2) 19(1)

TABLE 5 7 simulation. Singlet and non-singlet dynamical masses

9us and PNS, and average number seen during one history. Numbers between parenthesis are the estimated statistical errors on the last digit of the number on their left

of order—disorder flips L 16 16 24 24 32 32 48 64

P

Ps 0.00100(2) 0.00110(4) 0.000318(8) 0.000355(15) 0.000138(4) 0.000166(4) 0.0000425(15) 0.0000152(4)

1.28000 1.29356 1.28700 1.29356 1.28900 1.29356 1.29100 1.29250

N0

— —

~1

z_~ ~

va,,o

PNS

0.00033(1) 0.000078(4) 0.000112(4) 0.0000295(15) 0.000059(3) 0.0000134(6) 0.0000223(15) 0.0000058(3)

~ ,

IVI

—

Flips 276(3) 113(3) 184(3) 131(4) 169(2) 131(4) 78(2) 76(2)

qmaxa{Na}—l —

The large Monte Carlo time behavior of the auto-correlation functions of M and B defines the singlet (rs) and non-singlet (rNS) exponential auto-

correlation times [11]respectively. The dynamical masses ~t, = l/r, can be found in table 4 for q = 10 (table 4 contains only those results that changed since ref. [11]), and table 5 for q = 7. Extreme slowing down is observed on the largest lattices for q = 10. We only see a few order—disorder flips [11,17] on the L = 50 lattices, and this explains why we got unreliable results from those lattices. The dynamics is much faster for q = 7, as expected since the order—disorder surface tension and thus the barrier between ordered and disordered states are smaller. This is the reason why we put much less computational efforts on q = 7 than on q = 10.

782

A. Billoire et a!. / Finite-size scaling

4. Results for the q

=

10 model

Our results for CVmax, BLmin, U4mjn and the corresponding fl’s can be found in table 6. They have very small statistical errors, and the improvement is spectacular when compared to previously published results [5,13]. This is due to the use of the spectral density method together with much larger statistics. Our results for the maximum of the specific heat C Vmax as a function of are summarized in fig. 1. The data points are perfectly aligned within the precision of the drawing. This way to present data for CVmax, although frequently used, is somewhat deceiving as shown in fig. 2 where CVmax/L2 is presented as a function of l/L2, showing some curvature that fig. 1 did not show. An enlarged view of the small 1 /L2 part of fig. 2 can be found in fig. 3. The three curves drawn are the asymptotic predictions of eq. (7) for Cd = 12.7 ±0.3 which is our estimate of Cd for the q = 10 model. This estimate comes from various fits to eq. (7) with a l/L” or a e~L term added. Figs. 2 and 3 are in nice agreement with the theoretical expectations; CVmax/L2 does go to a non-zero constant as L grows, and the leading finite-size corrections do behave like 1 /L2. The two L = 50 data are plotted separately. They are mildly inconsistent. Since the data altogether approach the theoretical curve from below, we are not that happy with the L = 50, /3 = 1.425 point (which is supposed to be more reliable than the fl~point) that overshoots the theoretical curve. Taken alone, it would imply the value Cd = 13.5 + 0.5.

Results of the q L

50 50 44 36 24 20 16 12

=

TABLE 6 10 simulation. Values of the extrema of CV(fi), BL(fi) and U4(fi). The first L = 50 line is for fi = 1.425, the second for fi~

CVmax

628.89(54) 625.94(163) 489.26(42) 331.26(34) 153.24(12) 109.62(08) 73.64(05) 44.95(06)

BLmin

U4min

—0.11250(11)

1.09619(59) 1.09881(127) 1.11802(51) 1.15865(71) 1.26316(77) 1.31719(66) 1.38473(72) 1.47147(122)

—0.11225(37)

—0.11368(12) —0.11594(11) —0.12443(11) —0.13071(12) —0.14159(11) —0.16225(25)

fi(CVmax)

fi(BLmin)

50 50

1.42479(06) 1.42454(12)

1.42414(06) 1.42389(12)

1.42478(06) 1.42454(12)

44 36

1.42434(04) 1.42361(04)

1.42350(04) 1.42233(04)

1.42434(04) 1.42360(04)

24 20

1.42068(04) 1.41846(03)

1.41767(04) 1.41401 (04)

1.42065(04) 1.41840(03)

16 12

1.41457(05) 1.40683(11)

1.40734(05) 1.39321(12)

1.41444(05) 1.40646(12)

fi(U4min)

A. Billoire et a!. / Finite-size scaling

783

800 CVmax 600

—

—

400

—

—

200—

—

0 I

0

1000

I

I

I

2000

I

I

I

3000

2, for the q = 10 two-dimensional Fig. 1. The maximum of the specific heat, as a function of L Potts model. The straight line is the theoretical prediction. The intercept is fitted.

This would spoil the agreement with all other points. Our results for the Binder parameter as function of 1 /L2 can be found in fig. 4 together with the prediction of eq. (9). The figure is very similar to fig. 2. An enlarged view of the small l/L2 part is shown in fig. 5. The /3 = 1.425, L = 50 data point spoils the picture somewhat again. Alone, it would imply Cd = 13.4 + 0.3. Results for the minimum of U4 can be found in fig. 6. In marked contrast with the situation for CVmax/L2 and BLmin, the deviations from the theoretical expectation are huge at small L’s, and the large L data over-shoot substantially the theoretical curve (the L = 50, fi = 1.425 point would be on the theoretical curve for Cd = 14.58 ±0.09). Taken altogether the data are consistent with a value one for the infinite-volume limit of U4min, but are far from providing a convincing evidence that the limit is indeed one. The effective fl’s, fl(CVmax) and fl(BLmin) as functions of l/L2, can be found in fig. 7, together with the prediction of eqs. (6) and (10). We did not plot our data for fl(U4min), that are nearly identical with the data for fl(CVmax), in agreement with eqs. (6) and (14), but an amazing fact, given the difference between the behaviors of U4min and CVmaX/L2. Note that fi~ depends on C 0 and Cd only through their exactly known difference Cd C0. The agreement with the theoretical expectations is good, up to corrections that fade away rapidly as L grows. In order to take a closer look, we have plotted —

A. Billoire et a!. / Finite-size scaling

784

0.32

2

—

CVmax/L

—

0.30—

—

0.28—

—

0.26—

—

0.24 I

0.000

I

I

I

I I

0.002

I

I

I I

0.004

I

I

I

I

0.006

I

I

I

0.008

Fig. 2. The maximum of the specific heat divided by L2, as a function of i/L2, for the q = 10 two-dimensional Potts model. The three curves are drawn using the central value and the one standard deviation estimates for Cd.

0.0000

0.0002

0.0004

0.0006

0.0008

Fig. 3. The maximum of the specific heat divided by L2, as function of i/L2, for the q = 10 two-dimensional Potts model (enlarged view of the small i/L2 part of fig. 2). The three curves are drawn using the central value and the one standard deviation estimates for Cd.

Bi!loire et a!. / Finite-size scaling

A.

—0.10

785

I

BLmin —0.12

—

—

x

—0.14

—

—

—0.16

—

—

—0.18

I

I

I

I I

0.000

I

I

I

I

0.002

I

I

I

I

0.004

I

I

I

0.006

0.008

2, for the q = 10 Fig. 4. The minimum of the Binder’s fourth-order cumulant, as a function of 1/L two-dimensional Potts model. The three curves are drawn using the central value and the one standard deviation estimates for Cd.

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

BLmin

—0.110

—

—0.115

—

—

—

a

0.0000

0.0002

0.0004

0.0006

0.0008

Fig. 5. The minimum of the Binder’s fourth-order cumulant, as a function of 1 /L2, for the q = 10 two-dimensional Potts model (enlarged view of the small l/L2 part of fig. 4). The three curves are drawn using the central value and the one standard deviation estimates for Cd. The L = 50, fit point has been slightly shifted to the right, to improve readability.

786

A. Billoire et a!. / Finite-size scaling

1.5

I

1.4

—

I

I

4) / I

1111111

I I

I

~1

I

—

/

1.3-

x -

x

1.2—

—

1.1

—

1.0 I 0.000 I

I

I I

I

I

0.002

I I

I

I

I

I

0.004

I

I

I

0.006

I

0.008

2, for the q = 10 two-dimensional Potts model. Fig. 6. The minimum of U4, as a function of 1 /L The three curves are drawn using the central value and the one standard deviation estimates for Cd.

1.43

I

I

I

I I

I

I

I

I

I

)< ~(CV~~)

1.42

I

I

I

I

I

I

I

I

+ l~(BLmin)

—

—

+

>1

1.41—

1.40

—

1.39 I 0.000 Fig.

—

I

I

I

I I

0.002

I

I

I

I

0.004

I

I

I

I

0.006

I

I

I

0.008

and fi(BLmin), as a function of ilL2, for the q = 10 two-dimensional Potts model. The i/L4 contributions to the theoretical curves are barely visible.

7. fi(CVmax),

A. Billoire et a!. / Finite-size scaling

1.43ii~_

I

I

I-

I

x

I I

I

I

~

1y(CVm~x)

+

I

I

787

I

I

~y(BLmin)

1.429

—

—

1.428

—

— a

1—

1.427— I

+

I

::0.00000 NIlIlI,,

0.00002

0.00004

4, for the q = 10 two-dimensional Potts model. Fig. 8. The y(CVmax), nearly horizontal and y(BLmin), line isas the function prediction of l/L for y (CVmax), the other line the prediction for y (BLmin). The effect on y (BLmjn) of the uncertainty on the value for Cd is barely visible.

the data for

y(CVmax)

=

fl(CVmax) + Ed—EOL2’

y(BLmin)

=

fl(BLmin) + ln(q(E

(23)

0~~~)i

(24)

4 in fig. 8 together with the 1/L4 term computed in as function ref.a[13]. Our of datal/Lare in blatant disagreement with the latter. They are not precise enough, however, to determine the functional form of the correction to the leading 1 /L2 behavior, and would be consistent with e.g. either a 1 /L4 or l/L3 behavior. The effect is anyway quite small (see the scale of fig. 8). Our results for E(flt), CV(fl 1) and 2, fleff can be found Fig. 9 shows together with a in fit table to the7.simple ansatz our data for E(flt) as function of 1/L EL(flt)

with E~

=

=

E~+ A

(E~+ qE 0)/(l + q), A

=

(25)

e~~/L0

=

—0.187 + 0.007, L0

=

5.28 + 0.07,2

8.1 for six degrees of freedom, with a notable contribution to this

x

A. Billoire et a!. / Finite-size scaling

788

I

—1.59

I

I

I

I

I I

I

I

I

I

I

I

I

I

E(~~)

I 0.000

II

II

0.002

0.004

0.006

II

I

0.008

2, for the q = 10 two-dimensional Potts model. The curve drawn Fig. 9. E (Pt), as function of 1 /L is a fit to the data of the form EL (fit) = E~ + A e1-1~° with parameters as given in the text.

coming from the suspicious point L = 50, fit despite its large errors. Two comments are in order, first the fitted value for L 0 is close to ~(fit), next the fact that corrections are “exponentially small” does not mean that they are not visible. They are as large as 1.2% for L = 12. The small statistical errors obtained allow for a stringent test of the twolattice method to estimate the transition temperature (using eq. (17)). This is shown in fig. 10, where one finds fl~a(L~,L~+i) as function of 1/L~for {L,} = [12, 16,20,24,36,44]. The mean energy EL (fl) on the lattice of size L, is computed with the spectral density method using only the simulation performed at the value flMc the closest to fi (C Vmax); in all cases, but for L = 12, this is the simulation with the largest statistics. The other simulations performed on the same lattice size (if any) are not used. The data of fig. 10 are well represented by the shape

fiL

= floo

+

fib

e/L0,

(26)

with the parameters fl,~= 1.42596 ±0.00007, in remarkable agreement with what we know indeed (fit = 1.42606), fib = 0.025 ±0.001, L0 = 5.02 ±0.04. By excluding the point L = 12, one would obtain the value fl = 1.42596 ± 0.00010. In conclusion, flecr from the two-lattice method is by far closer to

A. Billoire et a!. / Finite-size scaling

1.429

1.426

789

~IiIIIIIIIII

—

1.425 I 0.000

—

I

I

I

I I

0.002

I

I

I I

I

I

I

I

0.004

I I

0.006

I

I

I

0.008

Fig. 10. fieff(Li,Li+1), as a function of l/L~,for the q = 10 two-dimensional Potts model. The curve drawn is a fit to the data of the form fieff(Lj,Lj+i) = fit + A e~i/’-owith parameters as given in the text.

the infinite-volume limit fit than traditional estimators like fi (C Vmax). The theoretical prediction that the deviation of flea from fit is not a power law in 1 /L is essential in getting such a nice result for fl~. A glance at fig. 10 teaches 2 would perfectly reproduce the data, leading us that a linear behavior in 1 /L however to a wrong estimate (fi~~ 1.4255). We finally present in fig. 11 our results for CV(fit), as a function of i/L2. They are in disagreement with our expectations. Although eq. (16) for CV (fit) is believed to be less sensitive to corrections than eq. (7) for C Vm~,we see much larger deviations from the theoretical curve. In order to reproduce the data, we would need quite a larger value of Cd, e.g. the ansatz

CV(fi~)= C + (1±q)2(E

2fi~+ A L2 eL~

0

—

(27)

Ed)

gives a stable fit with C 18, A —0.85, L 2 7. In such a fit, 0 as4.5 and as x 46% of CV(flt) at the O(e~’-~)terms coming from eq. (3) make much L = 12 ! A value of C~ 18 would ruin the agreement observed for L 36 between our data for CVmax (BLmin) and eq. (7) (eq. (9)). It would mean that the asymptotic regime sets in for larger lattice.

790

A. Billoire et a!. / Finite-size scaling

~IIIIIIIiiIIIIII

111111111

0.000

0.002

I

I

I

I

0.004

I I

I

I

0.006

2, as a function of i/L2, for the q

Fig. 11. CV(fit)/L

=

I

I

0.008

10 two-dimensional Potts model.

TABLE 7

Results of the q

=

10 simulation. Values of E(fit), CV(fit)/L2 and fi~ff(Lj,L~+i).The first L = 50 line is for fi = 1.425, the second for fit

L 50 50 44 36 24

E(fit) —1.5974(67) —1.6198(83) —1.6035(32) —1.5993(20) —1.6022(09)

CV(fit)/L2 0.0928(76) 0.0660(102) 0.0877(37) 0.0971(23) 0.1090(11)

1.42643(49) 1.42595(19) 1.42615(08)

20 16

—1.6050(07) —1.6103(06)

0.1167(07) 0.1286(07)

1.42650(19) 1.42697(16)

12

—1.6200(09)

0.1491(11)

1.42822(50)

5. Results for the q

=

fieff

7 model

The q = 7 model has a weaker transition than the q = 10 model, and we expect the asymptotic regime to be pushed to higher values of L. We cannot tell exactly how far it will be pushed, before doing the experiment, since there are two length scales in the problem: the inverse latent heat, which is nearly twice larger for q = 7, and the correlation length, which is roughly five times larger according to eq. (18). The inverse surface tension, a third relevant length scale, presumably goes like the correlation length.

A. Billoire et al. / Finite-size scaling

0.14

~I

I

I

I

I

I

I

I

I I

I

I

11111

791

I I

2

CVmax/L

0.04 ~IIIIIII 0.000

III

0.00 1

0.002

1111

0.003

0.004

Fig. 12. The maximum of the specific heat, as a function of l/L2, for the q Potts model.

=

7

two-dimensional

The data for the maximum of the specific heat C Vmax/L2 as a function of 1 /L2 as found if fig. 12 indeed look less first-order like (results for BLmin are very similar). The statistical errors are extremely small again. We have drawn the theoretical estimates using the educated guess C~= 50±10, obtained from various fits, like in the q = 10 case. CVmax/L2 does go to a non-zero constant as L grows, but we do not see the expected asymptotic 1 /L2 behavior, and if we were to estimate the latent heat from these data, we would overestimate its value. Results for the minimum of U4 can be found in fig. 13. The energy probability distribution on lattices of size ~ 64 is clearly far from the asymptotic limit. Although U4 is a neat quantity, from a theoretical point of view, it is, alas, not very useful to decide about the order of the transition. Recall that one has to decide whether U4min reaches one, or a value larger than one, in the infinite-volume limit. The effective fl’s, fi(CVmax), and fi(BLmjn) as a function of l/L2, are to be found in fig. 14. We did not plot our data for fi(U4min), which are again very close to the data for /3 (C Vmax). The theoretical curves drawn do not include the 1 /L4 terms from eqs. (6) and (10) which would made things worse. Fig. 15 shows our data for E(fit) as function of l/L, together with a simple fit,

792

A. Billoire et a!. / Finite-size scaling

1.8 I

I

I

11111111

I

I

I

I

I

I

I

I

I

I

U4min 1.6

*

—

x

*

1.4

— *

1.2

-

1.0 I 0.000 I

—

I

I

I I

0.00 1

I

I

I

I

I

I

I

0.002

I

I

0.003

2, for the q

Fig. 13. The minimum of U4, as a function of ilL

x 1.29—

fi(~Vmax)

+

I

I

=

7

two-dimensional Potts model.

<

1.27 11111111111 0.000 0.00 1 0.002

I

0.004

fi(BLmtn) —

1.28—

Fig. 14.

I

*

III

I

0.003

II

—

III

0.004

0.005

fi(CVmax), and fi(BLmjn), as a function of l/L2, for the q = 7 two-dimensional Potts model. The upper straight line is the prediction for fi(CVmax), the lower for fi(BLmin).

A. Billoire et a!. / Finite-size scaling

~

~:: —1.56

I

—

0.00

793

I

— 0.02

0.04

0.06

0.08

Fig. 15. E(fit), as a function of 1/L, for the q = 7 two-dimensional Potts model. The curve 1’--~with parameters as given in the drawn is a fit to the data of the form EL (fit) =text. E~+ A e_’—

EL(fit)

=

E~+ A e~~/L0

(28)

with A = —0.067 ±0.002, L 2 = 10 (the L = 16 data = 27.3 ±0.4 and x being excluded from the fit).0 The length scale Lb is close to ‘~(flt)as in the q = 10 case. In fig. 16 we finally present our results for CV(fi~),as function of l/L2. In this case (q = 7) they are not worse than for CVmax when compared to theory. 6. Conclusions We have performed a high-statistics numerical simulation of the q = 10 and q = 7 two-dimensional Potts models, on lattices of various sizes L. These models are known to present a first-order phase transition at a known value fit of fi. By high statistics, we mean that during the MC run the systems did flip many times between ordered and disordered phases. However, in the L = 50 case, the number of flips was too low for the data to be fully reliable.

A. Billoire et a!. / Finite-size scaling

794

0.08

///

—

0.06—

CV(fi~)

—

*

—

a I

0.04—

0.02 I 0.000

—

I

I

I

I

I

0.00 1

I

I

I

I

I

0.002

I

I

I

I

2, for the q

Fig. 16. C V (fit), as a function of 1 /L

I

0.003 =

I

I

I

I I

0.004

I

I

I

0.005

7 two-dimensional Potts model.

We have made a comparison of our data with theoretical predictions on the L-dependence of (a) extrema of moments of the energy distribution (CV/L2, BL, U4) and of the fi values where these extrema are reached; (b) moments evaluated at fi = fit. The theoretical predictions of interest here are the following. Given the values at fit of the energies and specific heats of the pure phases in the thermodynamical limit (three out of these four parameters are exactly known), theory provides the first terms of an expansion in terms of x = 1 /L2 for the above quantities. The estimates so obtained are accurate either to the next order in x (case (a), eqs. (6)—(14), or up to O(e’~)terms (case (b), eqs. (15) and (16). The result of our comparison between theory and data is the following: Although the familiar plot of the maximum of the specific heat versus the volume L2 seems to exhibit perfect linearity (fig. 1), a closer look, allowed by the very high precision reached, shows significant deviations at small L. (i) In case (a), our data do behave linearly in x for small x. The x 0 limiting value agree well with the prediction. The observed slope is reproduced for the value Cd = 12.7 ±0.3 of the only parameter left free, although the data for U4min would prefer a slightly larger value (see figs. 2—8 and 12—14). (ii) However, whenever the x2 term is known, as it is for the /3 values —~

A. Billoire et al. / Finite-size scaling

795

where energy moments are extremum, this term does not explain at all the observed departure from linearity. This finding is clearly exemplified in fig. 8. (iii) In principle, the theory is especially powerful at fit [case (b)] since the power expansion in x contains only a few terms with a O(e’~) rest. We found that the L-behavior of the average energy is in beautiful agreement with the expected behavior: the exponential ansatz of eq. (26) works very well (figs. 9 and 15). Unfortunately, the situation is not so nice for CV(fit)/L2: the expected small x linear behavior is not seen. The latter findings can be interpreted as evidence for large “exponentially small” corrections. Ifsuch is the case, the true value of Cd must be substantially larger than thought before (rather 18 than 13), and similar corrections are present in case (a). Although this contradicts the optimistic conclusions drawn in (i) above, it may account at the same time for the deviations discussed in (ii) . It would mean that the lattice sizes which can be studied on presently available computers are not large enough to actually observe the approach to asymptotics predicted by theory for (a class of) first-order transitions. The leading power-law behavior is probably obscured by 0 (e~L) terms. This finding is in fact quite interesting since there is more physics in the latter terms, which reflect the existence of coexisting phases close to the transition, than in power-law expansions in x related to the expansion of the pure phase free energies in powers of (fl fir). If the problem at hand is to learn by numerical means about a transition whose characteristics are not already known, the conclusions of our study would be: it is possible to convince oneself that the transition is first order when quantities like C Vmax/L2 or BLmin clearly have non-zero limits when L cx, but one has to relax on the request that the l/L2 behavior must be seen. In such a case, the actual extrapolations to the infinite-volume limit require hints about the nature and behavior of sub-leading contributions. We have shown how theoretical predictions at fit can be used to measure the fit value with a precision which dwarf conventional methods. ~

—

—*

We would like to thank R. Balian, C. Borgs, S. Gupta, A. Irbäck, W. Janke, R. Koteck~’,A. Messager, S. Miracle-Sole, B. Petersson, J. Ruiz, A. Sokal and J. Zinn-Justin for conversations. References [1] K. Binder and D.W. Heermann, Monte Carlo simulations in statistical physics, An introduction, Springer Series in Solid State Physics 80 (1988) [2] V. Privman and M.E. Fisher, J. Stat. Phys. 33 (1983) 385 [3] Y. Imry, Phys. Rev. B21 (1980) 2042; M.E. Fischer and A. Berker, Phys. Rev. B26 (1982) 2507; E. Brézin and J. Zinn-Justin, Nucl. Phys. B257 (1985) 867;

A. Billoire et a!. / Finite-size scaling

796

K. Binder, Rep. Prog. Phys. 50 (1987) 783; Finite-size scaling and numerical simulations of statistical systems, ed. V. Privman (World Scientific, Singapore, 1990) [4] K. Binder and D.P. Landau, Phys. Rev. B30 (1984) 1477 [5] M.S. Challa, D.P. Landau and K. Binder, Phys. Rev. B34 (1986) 1841 [6] C. Borgs and J. Imbrie, Commun. Math. Phys. 123 (1989) 305 [7] C. Borgs and R. Koteck~,preprint HUTMP 90/B258; J. Stat. Phys. 61(1990) 79 [8] C. Borgs, R. Koteck~’and S. Miracle-Sole, J. Stat. Phys. 62 (1991) 529; C. Borgs and W. Janke, preprint FUB-HEP 6/91 [9]A. Billoire, R. Lacaze and A. Morel, Nucl. Phys. B (Proc. Suppl.)20 (1991) 661 [10] A. Billoire, Talk given at the Workshop on Monte Carlo Methods in Theoretical Physics, June 27 - July 6 1990, Elba International Physics Center, Isola d’Elba, Italy, to be published [11] A. Billoire, S. Gupta, A. Irbãck, R. Lacaze, A. Morel and B. Petersson, NucI. Phys. B358 (1991) 231

[12]M. Fukugita, H. Mino, M. Okawa and A. Ukawa, J. Phys. A23 (1990) L561 [13] J. Lee and J.M. Kosterlitz, Phys. Rev. B43 (1991) 3265 [14] R.B. Potts, Proc. Camb. Phil. Soc. 48 (1952) 106 [15] F.Y. Wu, Rev. Mod. Phys. 54 (1982) 235; [Erratum: 55 (1983)] [16] C. Borgs and J. Imbrie, preprint FUB 91 [17] R. Horsley, preprint HLRZ 91-12 [18] R.J. Baxter, J. Phys. Al 5 (1982) 3329; Exactly solved models in statistical mechanics (Academic Press, New York, 1982) [19] A. Billoire, S. Gupta, A. Irbäck, R. Lacaze, A. Morel and B. Petersson, Phys. Rev. B42 (1990) 6743 [20] P. Peczak and D.P. Landau, Phys. Rev. B39 (1989) 11932. [21] J.L. Black and V.J. Emery, Phys. Rev. B23 (1981) 429 [22] R.G. Miller, Biometrika 61 (1974) 1; B. Efron, Biometrika 68 (1981) 589 [23] A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 61 (1988) 2635; Phys. Rev. Lett. 63 (1989) 1195 [24] A.D. Sokal, Cours de Troisième cycle de ia physique en Suisse Romande, Lausanne, June 1989