Critical behaviour of two Ising models with multispin interactions

Physica 139A (1986) 395-411 North-Holland, Amsterdam

CRITICAL BEHAVIOUR OF TWO ISING MODELS WITH MULTISPIN INTERACTIONS H.W.J. BL(3TE, A. COMPAGNER, P.A.M. CORNELISSEN, A. HOOGLAND Laboratorium voor Technische Natuurkunde, Postbus 5046, 2600 GA Delft, The Netherlands F. MALLEZIE and C. VANDERZANDE* Instituut voor Theoretische Fysica, Universiteit Leuven, B-3030 Leuven, Belgium

Received 9 June 1986

We have performed extensive Monte Carlo simulations on the two-dimensional Ising models with n-spin interactions described recently by Debierre and Turban. Results for n = 3 models with sizes up to 128 x 128 are analyzed by means of finite-size scaling. This yields a value of the magnetic exponent Yh close to ~. Direct estimates of the temperature exponent YTdo not converge convincingly. However, assuming the presence of logarithmic corrections such as in the 4-state Potts model, we obtain an improved estimate of the temperature exponent close to a2, in agreement with 4-state Potts universality. This result is further supported by an exact mapping between the n = 3 model and the 4-state Potts model in an anisotropic limit. For the n = 4 model, we confirm that the phase transition is first order, and we estimate the discontinuities in the energy and the magnetization.

I. Introduction W e c o n s i d e r a t w o - d i m e n s i o n a l Ising m o d e l with two-spin i n t e r a c t i o n s in the y d i r e c t i o n a n d n - s p i n i n t e r a c t i o n s in the x d i r e c t i o n with a r e d u c e d H a m i l tonian

H =~ x,y

t

n--1 i=0

K2Sx,ySx,y+l + K n 1-I Sx+i,y ) ,

(1)

i.e. the c o u p l i n g s d e n o t e d by K are e q u a l to i n t e r a c t i o n s J d i v i d e d b y - k T . This m o d e l is self-dual1'2). F o r fixed ratio K 2 / K n the self-dual p o i n t is determined by * Present address: Limburgs Universitair Centrum, Departement WNF, Universitaire Campus, B-3610 Diepenbeek, Belgium. 0 3 7 8 - 4 3 7 1 / 8 6 / $ 0 3 . 5 0 (~) E l s e v i e r Science P u b l i s h e r s B.V. ( N o r t h - H o l l a n d Physics P u b l i s h i n g Division)

396

H.W.J. BLOTE et al. sinh[2K2] sinh[2K.] = 1

independent of n. If a single phase transition exists, it must be located at the self-dual point. Only for n = 2 (the square lattice Ising model) and for n = oo (see Maritan et al.3)) proofs exist that such a single phase transition indeed is present. Numerical work and approximations indicate that this is also the case for n = 3 and n = 4. Debierre and Turban 2) have applied finite-size scaling and transfer matrix techniques to the n = 3 model with finite sizes up to 9. They suggested that the n = 3 model belongs to the same universality class as the 4-state Potts 4) and B a x t e r - W u 5) models. Furthermore, the phase transition of the n = 4 model was concluded to be probably first order. Using finite-size scaling and other methods, Penson et al. 6) and Iglri et al. 7) found similar results for the one-dimensional quantum version of this model. However, the numerical results for the temperature exponent of the n = 3 model were low in comparison with the 4-state Potts model, namely YT ~ 4, while one has YT = for the 4-state Potts model. Recently, Vanderzande 8) has confirmed part of these results for the one-dimensional quantum version of the n = 3 model, using finite-size scaling and exact numerical diagonalization of the Hamiltonian for systems with linear sizes up to 15, and also other techniques. Standard finite-size analysis was found to produce estimates of the temperature exponent close to those obtained earlier6'7). A difficult problem in the analysis of models with 4-state Potts universality is the existence of a marginal operator 9) which may lead to logarithmic corrections, and thus prevent rapid finite-size convergence. For the 4-state Ports model the consequences were found to b e serious1°'11): finite-size scaling of small systems produces values of the temperature exponent that are too small ( y T ~ 4). It may be difficult to detect the presence of the logarithmic corrections if the finite-size data are restricted to small systems, in particular if they are subject to statistical errors. In order to check for the presence of logarithmic corrections, it is necessary to take into account much larger systems than those which are manageable with the transfer matrix and Hamiltonian methods. For this reason, we have performed Monte Carlo simulations on L × L models with linear sizes up to L = 128. In order to obtain results with sufficient statistical accuracy, rather long simulations were necessary. These were performed on the DISP (Delft Ising System Processor) 12) which has a speed over 10 6 spin updates/second. Particular attention has been given to the randomnumber generator 12) (in order to reduce unwanted correlations). Moreover, the spins were visited randomly, which tends to suppress any remaining effects of a possibly less than ideal random-number generator. The results for the n = 3 model are given in section 2. We have also

CRITICAL BEHAVIOUR OF TWO ISING MODELS

397

investigated the n = 4 model: section 3 gives the results. A discussion concerning the character of the transitions of the two models can be found in section 4.

2. R e s u l t s f o r t h e n = 3 m o d e l

Simulations on the D I S P are subject to certain restrictions: the boundary conditions are periodic, and the lattice sizes have to be powers of 2. This m a y influence the finite-size convergence of our data, as becomes clear when one considers the possible ground states of the n = 3 model. These are labelled + + + , + - - , - + - and - - + ; the ground states are periodic in the x and y directions. Three of the ordered states do not fit on lattices with sizes that are powers of 2. As a consequence, the + + + phase is favoured, leading to a positive expectation value of the magnetization. F u r t h e r m o r e , we may expect alternation in the finite-size results because odd powers of 2 are equal to a multiple of 3 minus one, whereas even powers of 2 equal a multiple of 3 plus one. H o w e v e r , it seems reasonable to assume that these effects b e c o m e unimportant for sufficiently large system sizes. The simulations of the n = 3 models with periodic boundaries on the D I S P t o o k place at the expected critical point K 2 --- K 3 --- K D on the basis of duality: KD = 1 ln(1 + V~). The simulations had a total length of 2 × 108 sweeps for most system sizes. D a t a for the magnetization mL, energy uL, etc. were taken at intervals of 100 sweeps. The susceptibility XL was calculated from the haean square magnetization (the magnetization of the infinite system is expected to vanish). The specific heat C L was obtained from the fluctuations in the energy. T h e r m a l and magnetic data are shown in tables I and II respectively, together with exact results for 2 × 2 and 4 × 4 systems. The results for the critical energy are, within their numerical accuracies, equal to - 2 - 1 / 2 1 n ( l + V ~ ) = - 0 . 6 2 3 2 . . . (the value predicted by duality for the infinite system) for all systems sizes considered here. On the basis of finite-size scaling13), we expect the following asymptotic behaviour for large system size L: CL ~ L2y'r-2 ,

XL ~ L2yh-2 •

(2)

Thus we obtain estimates yT(L) = 1 + ~ ln(CL/CL/2)/In 2

(3)

for the t e m p e r a t u r e exponent from pairs of successive values (CL/2, CL). Similarly, estimates Yh(L) are obtained from pairs (XL/2, XL). These estimates

398

H.W.J. BLOTE et al. TABLE 1 Energylike data for n = 3 models with periodic boundaries and size L x L. The entries at L = 2 and L = 4 are obtained by a numerically exact summation over all spin states. The data for L > 4 are Monte Carlo results obtained with the DISP. The second column gives the number of sweeps (the average number of times that each spin is visited by the Monte Carlo algorithm). The energy /~L in units J = K k T was sampled at intervals of 100 sweeps. The specific heat CL in units k was calculated from the root mean square deviations of the energy. The temperature exponent is estimated according to yT( L ) = 1 + ½ In( CL / CLj2) /ln 2. Statistical uncertainties in the last decimal places are given in parentheses. L

Sweeps

uL

CL

y.r(L)

2 4 8 16 32 64 128

2x 2x 2x 2× 3x

-0.6232 -0.6232 -0.6232(2) -0.6230(2) -0.6231(3) -0.6237(4) -0.6229(4)

0.4478 1.0553 1.9070(13) 3.2520(22) 5.409 (5) 9.016 (22) 15.23 (7)

1.6184 1.4268(5) 1.3850(7) 1.3670(9) 1.3686(19) 1.3782(35)

108 108 108 108 108

Table II Magnetic data for n = 3 models with periodic boundaries and size L × L. The entries at L = 2 and L = 4 are obtained by a numerically exact summation over all spin states. The data for L > 4 are Monte Carlo results obtained with the DISP. The second column gives the number of sweeps of the simulation. The magnetization m L was sampled at intervals of 100 sweeps. The susceptibility XL was calculated from the square magnetization values without correction for the magnetization which is nonzero for these finite systems. The magnetic exponent is estimated according to yh(L) = 1 + ½ In(xL/XL/2)/In 2. Statistical uncertainties in the last decimal places are shown between parentheses. L

Sweeps

mL

A'L

Yh(L)

2 4 8 16 32 64 128

2x 2× 2× 2× 3x

0.6306 0.5126 0.5155(6) 0.4681(10) 0.4316(20) 0.4015(40) 0.3604(63)

2.5789 9.5888 32.113(24) 108.87 (19) 364.5 (13) 1235 (10) 4055 (56)

1.9473 1.8719(5) 1.8807(14) 1.8716(29) 1.8803(62) 1.858 (12)

108 108 108 108 108

a r e a l s o s h o w n in t a b l e s I a n d II. T h e y e x h i b i t s o m e a l t e r n a t i o n , as e x p e c t e d . The absence of apparent convergence with L to a trivial value associated with a T = 0 o r a T = ~ f i x e d p o i n t c o n f i r m s t h a t t h e s y s t e m is i n d e e d c r i t i c a l a t t h e self-dual point. estimated

The

final result f o r t h e

t o b e Yh = 1 . 8 7 ( 1 ) .

magnetic

The temperature

exponent

exponent

is s u b j e c t i v e l y

is d i f f i c u l t t o e s t i -

CRITICAL BEHAVIOUR OF TWO ISING MODELS

399

mate because there is very probably an extremum in YT(L) at L = 32 or 64; the difference between the values found for YT(L) at L = 64 and L = 128 is over two standard deviations. Before proceeding with the analysis of these data, we have to check for possible effects due to the fact that the DISP is unable to simulate systems with sizes that are multiples of 3. To this purpose we have also performed Monte Carlo simulations in software, for L = 6, 9, 12 . . . . . 30, also using periodic boundaries. The length of these simulations was between 5 x 105 and 106 sweeps. The results for the specific heat join smoothly with those obtained by the DISP. There is also agreement between the susceptibilities if the DISP results are corrected for the magnetization which is nonzero for the system sizes given in table II. The data are shown in figs. 1 and 2 on a logarithmic scale. The straight lines in these figures yield the following estimates of the critical exponents: YT = 1.34 -----0.02, and Yh = 1.89 -----0.02, in agreement with the trends of the DISP results for L ~< 32. This confirms our expectation that the question whether L is a multiple of 3 becomes unimportant for large L. Our results for Yh of the n = 3 model are close to ~ , the value of the 4-state Potts model. However, our estimates of YT are much lower than the four-state Potts value 3, though in agreement with earlier investigations. The question arises whether these data are consistent with 4-state Potts universality. In an attempt to answer this question we make use of the renormalization transformation of Nienhuis et al. 9) for the Potts model. Introducing the inverse system size as a relevant scaling field into the recursion relations of Nauenberg and Scalapino14), the finite-size expansion of the specific heat C L of systems in the 4-state Potts universality class can be obtained11):

1.5

I

I

I

I

In C L

1.0

0.5

----~ln L 0 1.5

I

I

I

I

2.0

2.5

3.0

3.5

Fig. 1. Log-log plot of the specific heat of the n = 3 model versus system size. These data were obtained by simulations in software, for system sizes which are multiples of three; periodic boundaries were used.

400

H.W.J. B L O T E et al.

In

I

I

I

i

I

I

I

I

2.0

2.5

3.0

3.5

"

4.0

3.0

2.0

----~ In L 1.5

Fig. 2. L o g - l o g plot of the susceptibility of the n = 3 model versus system size. T h e s e data were obtained by simulations in software, for system sizes which are multiples of three; periodic boundaries were used.

C L = L2yr-2[ao(ln L ) -3/2 + al(ln L ) -5/a + . . . ] .

(4)

Here we have neglected terms with lower powers of L. Furthermore, we have assumed that the marginal scaling field (on which the amplitudes a i depend) is large in comparison with 1 / l n L . Hence eq. (4) does not apply to the Baxter-Wu modelS), for which the marginal scaling field is zero. Analysis of eq. (4) by means of eq. (3) leads to YT(L) = y x - 3 / ( 4 1 n L ) + p 2 / ( l n L ) 2 + . . - ,

(5)

where P2 is a function of a 0 and an: its value is unknown in principle. Eq. (5) converges very slowly with increasing L. It is interesting that it can lead to a behaviour very similar to that of y x ( L ) in table I. For instance, if we choose P2 = 1.5 and neglect the higher terms represented by the dots in eq. (5), we reproduce a very fiat minimum such as in the YT(L) data in table I. We conclude that the data in table I are consistent with eq. (5) which was obtained on the basis of 4-state Potts universality. On the other hand, the existence of an extremum at L = 32 or 64 does not combine very well with the possible absence of logarithmic corrections: this would require very large correction-to-scaling amplitudes. We mention that a finite-size analysis of Monte Carlo results 15) for the specific heats of the simple quadratic Ising model, and of the Baxter-Wu model yielded YT estimates which seemed to converge rapidly and monotonically to the exact values 1 and 3 respectively. Assuming the presence of logarithmic corrections as given by eq. (4), we can

CRITICAL BEHAVIOUR OF TWO ISING MODELS I

I

I

401

[

'l,ll,

/

1.6

J io

tf

J

1.5 I I 128 64

I 32

I 16

I 8

Fig. 3. Estimates y~(L) of the temperature exponent yr of the n = 3 model as a function of L on a (i/In L) 2 scale. As explained in the text, linear behaviour may be expected for large L. The dashed line represents a plausible extrapolation leading to YT= 3. Except for the largest value of L, the numerical uncertainties do not exceed the size of the data points. correct the Y'r estimates for the term p r o p o r t i o n a l to i / I n L , which has a k n o w n amplitude (eq. 5). T h u s we define

y~-(L ) = y,r(L ) + 3 / ( 4 In L ) . Fig. 3 shows this quantity as a function of L on a (In L ) 2 scale, on which it is e x p e c t e d to d e p e n d linearly for large L. T h e data points show that the asymptotic range is not quite reached, since there is still an appreciable curvature. But the value of the expansion p a r a m e t e r I / I n L , which is 0.21 for the largest system size, suggests that the asymptotic range is not far off. A plausible extrapolation (dashed line in fig. 3) to the value YT = 3 is possible. In o u r opinion, one can, with a m o d e s t degree of optimism, estimate y~(~) to be equal to 3 within two times 10 -2 f r o m the data in fig. 3. We conclude that o u r results for the critical e x p o n e n t s of the n = 3 m o d e l agree with those of the 4-state Ports model.

3. Results for the n = 4 model

We have also used the D I S P to p e r f o r m M o n t e Carlo simulations of the n = 4 m o d e l at the self-dual point K 2 = K 4 = K D. We have used again periodic b o u n d a r i e s and system sizes which are powers o f 2. Numerical data and calculated e x p o n e n t s are shown in tables I I I and IV. T h e energies do not seem to c o n v e r g e to the duality value 2 -1/2 In(1 + V ~ ) for systems with a c o n t i n u o u s

402

H.W.J. B L O T E et al. TABLE III Energylike data for the n = 4 model. The entries at L = 2 and L = 4 are numerically exact results. The data for L > 4 are Monte Carlo results obtained by the DISR The energy of these systems was sampled at intervals of 100 sweeps. The temperature exponent is estimated as described in the caption of table I. These estimates do not appear to converge rapidly to the expected discontinuity fixed point value YT=2. Statistical errors are shown between parentheses. L

Sweeps

uL

CL

yx(L)

2 4 8 16 32 64 128 256

108 108 7.5 × 2.5 x 2x 8x

-0.7532 -0.6766 -0.678(2) -0.676(2) -0.669(1) -0.669(2) -0.669(3) -0.669(7)

0.1924 0.4535 1.387(2) 3.232(5) 7.57 (5) 18.7 (4) 54 (5) 155 (25)

1.612 1.806(1) 1.610(2) 1.614(4) 1.65 (2) 1.76 (7) 1.76(13)

107 107 107 106

energy. This already shows that the system undergoes a first order transition at the self-dual point. The magnetization does not differ significantly from zero, in agreement with the symmetry of the Hamiltonian. The magnetic exponent Yh as estimated from consecutive system sizes up to 642 exhibits good apparent convergence to the expected discontinuity fixed point value Yh = 2. However, the same cannot be said of the temperature exponent. It seems that higher values of L are needed to find the true asymptotic behaviour. For this reason, we have executed additional calculations for system sizes up to 5122. A problem with these large systems is that the phase switching times (i.e. the average time needed by the system to evolve from a disordered phase to an ordered one and vice versa) become very long, leading to a decreasing accuracy of the specific heat results. In fact, no phase switches at all took place during the simulations of the 5122 system. The temperature exponents as calculated from the L = 128 and L = 256 systems (see table III) show a tendency to increase towards the expected value YT = 2, but the convergence and accuracy are not yet satisfactory. It appears that surprisingly large values of L are needed before the asymptotic behaviour is reached. This is also evident from the energy distribution functions shown in fig. 4: only at L = 256, two clearly separated peaks appear. Although there is some statistical inaccuracy in this figure, we believe that the small peak at the right-hand side is significant and signals the evolution of the energy distribution function into the sum of two delta functions when the system size grows to infinity. The stability of the ordered and disordered phases of the 5122, n = 4 model at the self-dual temperature allows us to estimate the latent heat za using rather long simulations. The results for the energy averages are: U ( K D + ) / k T D = - 0 . 6 9 6 ( 2 ) and

C R I T I C A L B E H A V I O U R O F T W O ISING M O D E L S

403

TABLE IV Magnetic data for the n = 4 model. T h e entries at L = 2 and L = 4 are numerically exact results. T h e data for L > 4 are M o n t e Carlo results obtained by the D I S E T h e magnetization of these systems was sampled at intervals of 100 sweeps. T h e magnetic exponent is estimated as described in the caption of table II. These estimates seem to converge well to the expected discontinuity fixed point value Yh = 2, but the inaccuracy (shown between parentheses) increases rapidly with the system size. L

Sweeps

mL

XL

Yh(L)

2 4 8 16 32 64

108 108 7.5 × 107 2.5 × 107

0 0 - 0 . 0 0 1 2 (8) -0.0008(20) - 0 . 0 1 4 (11) - 0 . 0 1 5 (16)

1.7071 2.630 10.14(3) 37.6 (3) 146 (6) 594 (44)

1.312 1.974(2) 1.945(6) 1.98 (3) 2.01 (6)

U ( K D - ) / k T D =0.550(2) so that A / k T D =0.146(3). It seems plausible that these energies are close to those of the infinite system. These simulations consisted of 5 × 105 sweeps for the ordered phase and 1 × 105 sweeps for the disordered phase. We have not analyzed the magnetic susceptibilities for the largest system sizes because of the inaccuracy due to the long phase switching times. Instead, we have obtained a second series of estimates of Yh from the magnetization distribution functions. For the larger system sizes, we observe that the magnetization distribution has three rather narrow peaks (see fig. 5): one near zero magnetization due to the disordered and six of the ordered phases, and two symmetrically positioned peaks due to the + + + + and ordered phases. We have determined the positions of the latter peaks as the average of the magnetization values at half maximum (table V). Since the width of the magnetization distribution is expected to scale a s L yh-2, new estimates of Yh follow (see also table V) from consecutive magnetization values. Again, the results agree with the expectation Yh = 2 for the large L limit. The simulation of the ordered L = 512 system allows us to estimate the discontinuity in the spontaneous (sublattice) magnetization of the ordered phases at the transition point: m(Ko+ ) = 0.769(6). The first order character of the transition is also evident from the jump in the energy as a function of temperature at the transition point as illustrated in fig. 6 for a 2562 system. Finally we illustrate the different character of the n = 3 and n = 4 models at the transition temperature by means of spin configurations of 5122 lattices. Starting in a random configuration, the n = 3 model usually evolves rapidly into a state in which one of the ordered states dominates, just as the n = 2 (Ising) model with periodic boundaries. Such a configuration of

H.W.J. B L O T E et al.

404

16

15t

.020

.tO .OtO .05

0.00

• 80

-.60

-.40

0.00~

-.20

-.80

-.60

-.40

-.20

i

64

32 .030 .040 .020 .020 .OtO

O.OOC

)

|

-.80

-;60

-.40

0.000

- .20

/

- .40

- .20

256

.080

.041]

.0s

0.00'

-.60

- •

128 ,12

0.00¢

I

• 00

-.60

-.40

-.20

'

-.BO

.60

z

-.40

-.20

Fig. 4. E n e r g y distributions of the n = 4 model, for system sizes L = 8 up to L = 256. Only for the largest value of L, two clearly separated m a x i m a appear. T h e lengths of these simulations are given in table III.

C R I T I C A L B E H A V I O U R O F T W O ISING M O D E L S

405

16 .075

.20 .050 .10

JL

.025

O.O[ -

.00

- .50

0.00

1

0.000

,50

-1.00

32

- .50

0.00

.50

.12

64

.075

.050 .06

,025

0.000-1.00

1

- .50

0.01

.50

o.oe/',,~ -1,00

JL -.50

0.00

.50

128

.20

256

.20

.10 .10

0.00

-1.00

-.50

0.00

.60

o.oo -1.00

-.50

0.00

,50

Fig. 5. Magnetization distributions of the n = 4 model, for system sizes L = 8 up to L = 256. For the largest system sizes, some peaks at nonzero magnetization are absent because the corresponding ferromagnetic phases did not occur during the simulation. The lengths of these simulations are given in table III.

406

H.W.J. BLOTE et al. TABLE V Analysis of the magnetization distribution. The third column contains the half width m~ of the magnetization distribution function as estimated from the position of the peaks corresponding to the + + + + and . . . . ordered states. The average magnetization value was taken at half maximum of the peaks. For L = 512, the system was in the + + + + phase during the whole simulation. The magnetic exponent is estimated according to yh(L) = 2 +ln(m'Jm'L/2)/ln2. The estimated inaccuracies (shown between parentheses) are smaller than those in table IV.

L

Sweeps

m~

Yh(L)

8 16 32 64 128 256 512

108 108 7.5 × 2.5 x 2× 8× 5x

0.976(5) 0.914(5) 0.870(5) 0.834(5) 0.807(5) 0.783(10) 0.784(15)

1.870(10) 1.929(11) 1.939(12) 1.953(12) 1.96 (2) 2.00 (3)

I

i

107 107 107 106 105

i

i

U(T)/kT o -.50

-.60

-.70

----,. T/I -.80

I

.96

. 8

1.00

I

1.02

I

I

1.04

Fig, 6. Energy of the n = 4 model with linear size L = 256 as a function of the temperature divided by the self-dual temperature T o given by K o = J / k T D = ½In(1 + X/2). These data points consist of two series: one with increasing temperature and T / T D <~ 1, and one with decreasing temperature and T / T D I> 1. Each data point stands for 2 × 104 or 5 × 104 sweeps. This is short in comparison to the time scale on which the different phases evolve into each other at the transition point, but it is long in comparison with the apparent relaxation time within the different phases. The apparent statistical errors do not exceed the size of the data points.

CRITICAL BEHAVIOUR OF TWO ISING MODELS

407

t h e n = 3 m o d e l is s h o w n in fig. 7. T h e b e h a v i o u r o f t h e n = 4 m o d e l is d i f f e r e n t : S t a r t i n g in a r a n d o m c o n f i g u r a t i o n , a d i s o r d e r e d p h a s e d e v e l o p s c o n s i s t i n g o f m a n y p a t c h e s , e a c h c o r r e s p o n d i n g to o n e o f t h e o r d e r e d p h a s e s (fig. 8). Such a d i s o r d e r e d p h a s e was o b s e r v e d to b e s t a b l e o v e r long s i m u l a t i o n s (--3 x 10 5 s w e e p s ) . S i m u l a t i o n s s t a r t i n g in o n e o f t h e g r o u n d states s h o w t h a t also t h e o r d e r e d s t a t e s a r e s t a b l e u n d e r l o n g s i m u l a t i o n s . N o p h a s e c h a n g e s t o o k p l a c e in a r u n c o n s i s t i n g o f 5 x 10 5 s w e e p s . A t y p i c a l c o n f i g u r a t i o n is s h o w n in fig. 9: o n l y i s o l a t e d p a t c h e s o f t h e o t h e r p h a s e s a r e seen. I t can b e e s t i m a t e d f r o m figs. 8

Fig. 7. The n = 3 model with linear size L = 512 at criticality. Plus spins are white, minus spins are black. This picture was taken after a simulation of 2 × 104 sweeps starting from a random configuration. It suggests the presence of excitations on all length scales, as expected for a continuous transition in an infinite system.

408

H.W.J. BLOTE et al.

Fig. 8. The n = 4 model with linear size L = 512 at the transition point. This configuration was obtained after 1.5 x 105 sweeps, starting from a disordered (random) configuration. No transition to an ordered phase was observed: the simulation time to observe such a transition grows very rapidly with system size.

and 9 that the spin-spin correlations have a range of roughly a hundred lattice units. T h e o b s e r v a t i o n t h a t this is m u c h s m a l l e r t h a n t h e l a t t i c e size (512) gives s o m e s u p p o r t to t h e i n t e r p r e t a t i o n o f t h e results f o r t h e d i s c o n t i n u i t i e s in t h e e n e r g y a n d t h e m a g n e t i z a t i o n as t h o s e o f t h e infinite system. O f c o u r s e , in this i n t e r p r e t a t i o n t h e n u m e r i c a l u n c e r t a i n t i e s for t h e i n f i n i t e - s y s t e m d i s c o n t i n u i t i e s e x c e e d t h o s e q u o t e d a b o v e ; w e s u b j e c t i v e l y e s t i m a t e t h a t this is b y a f a c t o r o f ten.

CRITICAL

.

. ~ . i { . . . . . ,

tt5

i

..: .... -

;.

• "

.

,.'.,.;..

.

" .

"'.. ,

-

•

.

"

-

"

. . . . .

'

:i

"

pp

,N ;[ll~l_

.

, --.

.--

H'ilg

"

-

-

"

.0

.

.

.]-;~,f

~

"'~

.

.,

.:" t. Itil~

-

,~ ~""

-

,

"'

..

"..,.

"

: ,-

-

~

;',

...

"l.

409

•""'Y

""

-.-

.,"

-t

~

-

""

"

.

"':

"

,-

•

,

t"

i,"

,,..

. .i ,

"

",&',-~'-" " "

~--

:

~"

•

....

~

-.

'.

.•-

!

,.

,,,.~

.

,,.I..."

"

.

I.

d

"'

•

,ill,.'r..

"

'

""

>..

,

.

a :~."

--.

: .

"

""

" -"

..:

7:,' lQtli -

-'l'~'.'

~ } ~ . " . ' l

..

-ltr

-

' .... •~-:

'

,

.. .....

""

'

'.'-" :LI

..

"4,,;'~-?

.-

oo "sit"

a`'

._''.

,-

:o

~'. ,..'-

.~..illlll~lt,,.."

•

yJI.. . . .

.

- ....

lllli'"

~]11 ~l~l~l~l. • " -

•

¢~"

"i . . . . . .

r,.:--"'"-

-:. '-

K: ~ ~..-.,'g. ""

.

-" .•. "'-

"'

- " iil,n

-

," •

-~•

. ~ ' ~ I B " ~ : • i/illJ[ni

• v=

. . . . . .

":

",~'~',

,~','" .•.-.,-j~,.:, .~,: . ' . - . , . , •: ~ ",, . .....;;, I~l,~ : ~ , , ' . : , . .',~';~---'.:-'~;~ " "<'. ":,:.~." : . ; ~ . , ~ ; ] , ~ .: . . . . "- " ' ' : - ?a '" " ~ " " : l l ~ t i ~ l l ~ ; " "': : :•': :". .... . " ". : ~ - " . :-~: ,-~.•.. n . ,,.. ~. . . . . ..;:.. • Ulpl "" %-,~-

.",;~PlTtl'."

~m

:

•

"" "'

.

.*'qlli~'~'..

,,..,::~, :-

•~

I•. "=,'

"

-

J

. . . . . "

.' . "

-"

,"

"

"..~ .

•

. . . . . . . . . . . . . •, . . • •

.

j .

..:| .

.

"

'

~

-in.,,.:

'~"

~.. ..;"

"

.

,.

Y.

•"

".

;,;.':~.,

- ..........

; .....

;':.

..-'=

"....

,

.

",.

•.

°

;4.

•

.,r ~

;,

"-"

::

.

.

---

.,.

"'"

;..-

:

"'='=',:.'-

-

"

.

"~

..

"

-

-

. . . .

.%

'~" --

.' ~

"~

.'

,.,, "

.,"

"

_ :,

,~L~.

-'.'v,,.,

i. . . .

r ~

~["

"

-'.-

I

-

."

• '

!

-

arrjl-

"

•

.-':," "~'~

fNII

.

:.

- •

"

.

-:

'

.I

~"

,: :

•~- ..';.IY-"7".~".~

"

.

,.~"

"lllam

~il~

"

.".

.

-

. ,.-,~-.m.,

;-

.

,,] ,.

"i',~pe- -

-

-.

.

"'"

---,;-

"~.

~.~ ""."

•

,

- ' " "

!#

•

"^"

..7 -

,"

.

. . , . , - , , . .';".

.'i'~| ,.

-

" . .

.

%

"¢-~.c,•.

,~l:..

.

..

""~

" -.I~!

"

-"

•

,it , ~ .•;'. . . . . . . . . . . "' '%1

"

~".

~ ' ~. . "~ ; . . '.~ ' ~

- _~.:..'..o..r,- ~--...~l,~) ..

-

"~.'..-"~,-a,~"

.

1#11E,

'

.." .'..,: . ".

,v"

" ¢ i ¢"

- ~ . .", - ...... ,,~

." u'*.

-

: -•;;-

~,v, , "-" ,,,,: "

. . . .

....

"-,-"

._",,~

•

" .... .

?t"

-'~.,, .''l~|ll

"'...

"-lit.

. • .-". . .

" "

-

.

•

•

"u, •

,

"

-

--~

.

...ll.m.

. . . . . l; ~i "-

-~r,~

"

=.. : - = " , : ;

- ~ " ".,.~.~ ."

~ . . . .

...

."-

"#g~,

--

•

::.ql~.~,;:'"

|til~d~.t~'u,,f.m~

-'a~.-:'

~Jil"

."

•

""

-,~/i!5~.~. a,~.nT~. . ." . ...~ ..... ..~n~,~. p .t .~i[~j~ ~ •-; .

~

~..~.-.

"

""-

"

.

.=

S'.'Y','l',"r~.

.... ,~-

"'"

.v.~..-..(;.:

.

,

.

.

.. % / , . ~ ~ = ; , , • "~%,~,,-'.'..

" "i~l;l~ ~e~l.

"%;

.-

~.,:'kl~l,--'.':: ..

•,

..

..... - .

-.z .

~ll'

"~

.

"

o.•

.

-""

"-r:k,..:'~,.~ .,~t,,,-~

•

"."

".

."

;

-. "-,-:. "

-~:': ...• ,, • . .,~. . ~"y .... : ' -

.~

. ....,Ct~.4,f!..~lHF,18tili......-.,

."

\;.~.

c - .

"

~ c . ' ~ ."

..

"

-'"

•"

m,~';

"" " "

.

"",,'7":

. . . . .

41Ill''Ill'"

-:.~-"

"

:'~-: ....

~

;

.':-.~nt~,,....

"•,

;"

"~.ti~

- "E-"~

...

•'

- , ~ . , , - : ....o . . . . .

li~"

"

. " - " ~¢..

"'"

•

.,."

:-

".'~/)',,,r..

u,

"

.....

"

.

.

.

.,.,,~•.: :

~:.

•

. . . .

•

:

" ,

"" " " • ..... "' ." "" . |~ ~.•.,r., . . . . . - . . . . . . "~.",:'~'lli'll~. . : . : : . .•.. . " ~- ."~= . . . . . . . . %". . , • ,... :,; :,;:.~e .::[.~:. , ~ ¢ : - . . . . ;..~, . , . . ,....:: .,~@tl~[il~t¢~-~...- . . . . . , , . , , . :. ,..., :~r . . . . ' ~ . : . - . .

"

.af,..... "'" " "'..,,~l" '~ln

"• ~'~"

'.,,-~F~

,"

- ~ , ~ , . ' - " ; , : ', ,~." .~

"•-~ --

"

-.,lt/'_~Ulll~,

~....

,.,ix,

.........

" • ..... ...-" :',',, " :'4,.,..

,.

• "

.......

"

.

",-- ."~,

"."

. . . . ...'~tl,..In/! .

|

:i~a'~

• ""

-

.....

. . . . . . ,,t~,-. .... :: . . . . . .

.~,~,,~#l~'....-. i : . " ' : " ' " . ~ d t l l ~ t P L ....

,;-

-'"

~

".

-~...'i,u"

"~)p' • . -. '; ~'"" --" ~ ,.,., ~,.4t .-, .... -

~,i,,s,I-.,d#dl~l

.. •.,

,"

- ",,:~ ~ ", ~

,

,"

. _ .

""

.

,1" ~"

•

,';""'-""

•

.

. x . . . . ..

. . . . . . . . . . .

'•

:• .

-

-Jill

,

~

4 ) 1 ~

•

•

:

..

"

".

-

-~tl,,,,.

-•o ..-,nlltl2r

:~:;l~'~J~-.."

""

,.

,

,

-. .. . . .. . .., . . ~ .

". ;;-

,-'.. [";.

".

. . . . . . .

.

y,

-

.-

.Y.

-~,...t~,.:.-.~.~r/~h,'l[l~/,V'!!ixt..

"" .

"

•

,

"

-

?,

-" "

. ,:.,...-" • _.-'~ ;"

,.. .

"

....

..,,,'" -." . . . .

" "

.

":

";" ~'._

.......

,,-,,,. -

v.

.

.

"" , _ J l i i ~ l m n l - n

,.. -

"

°

'~'"

.,.~,illlk/illi~

.

.

~-

I..~

~ i

.....

..•

". .....

,.

- , . ' ; I Z~l~... ~,,,-• :.,.

,: ...-7".~.

..l~.

•

" ~

""-.""

".

" *°'"

.....

°

. . . . .

"

•l,.: ".~ -. •..

.~-~-" -

-.

•

.. "

~,

-..,-,-

~

" "

"

I

..

" .~

"'-.'-'n~'~.'-'-."'.-

"

, %"

.

-

. . . . • ..

• ....

....

.....

-

. . ;:

..

.

--i- . . . . . ,t ".'. ~":~il " ~..~.......;.,...=: nlji • . . .',

.

- . . . .

,, ,".

~ l l l , • +

"

. '.. .

~.

."

"-

• -it

--.dr

.

":

,,, .... ~:'~gii/~i"--"'"

., ."

...L.

,~i~Aikl

-

l.ii." :

-,,

-

. ".

",~kil~'.~

. .

.

'

?

•

.

;in

Z,..lttll,'ll'i,'.r"-

• ~.. "~

,

4

~"

". '-

.

-ni4ff~..~..~:."-:'

v,l.,iq~-

°-

~i~.'lii"j.

..

"liI. -% .l~ "~

.¢ fl~l.~,-"-;

.

4," -•{"

.

. . I,,~t..

"

,".

"r.;".

. . . . . . . . " ° ' . • ..

,.-<.. ~ .

"'"

-

, ....

-"

".IL.

OF TWO ISING MODELS

.

.... ~

',;%"" ~.|~,,,~Pl " -r," " .,.''~41~" , in.i .. .

,'

;....

.

- ",,-, ~'aiCsi,. - ~

"~-'"~" :. •

.

BEHAVIOUR

,.4e,,

""

I "

,

Fig. 9. T h e n = 4 m o d e l w i t h l i n e a r size L = 512 at the t r a n s i t i o n point. T h i s c o n f i g u r a t i o n was o b t a i n e d a f t e r 5 x 104 s w e e p s , s t a r t i n g f r o m a n o r d e r e d (all spins plus) p h a s e . A l t h o u g h s m a l l p a t c h e s c o r r e s p o n d i n g w i t h o t h e r o r d e r e d p h a s e s a p p e a r , t h e b a c k g r o u n d ( t h e + + + + p h a s e ) is s t a b l e o v e r l o n g s i m u l t a t i o n s . It is c l e a r t h a t s p i n - s p i n c o r r e l a t i o n s in the h o r i z o n t a l d i r e c t i o n h a v e a l o n g e r r a n g e t h a n t h o s e in t h e v e r t i c a l d i r e c t i o n ; this can b e i n f l u e n c e d by c h a n g i n g t h e r a t i o K 4 / K 2. T h e s a m e effect can b e s e e n in figs• 7 a n d 8.

4. Discussion In section n = 3 model specific

2, we

have

which

is c l o s e

heat data contain

found

a value

for the temperature

exponent

to 3. This result relies on the assumption

logarithmic

corrections

s u c h as p r e s e n t

of the that the

in the 4-state

410

H.W.J. BLOTE et al.

Potts universality class. This assumption was supported a posteriori by the result for YT. The accuracy of the finite-size analysis is comparatively poor: the expansion parameter is In L instead of L as usual. Hence, our analysis is to be compared with an investigation of a model without logarithmic corrections and with finite sizes up to 5! It is also useful to consider other evidence for the universality of the n = 3 model. In the first place, we mention the equivalence of a model resembling the n = 3 model with the B a x t e r - W u model, which was discovered by Horiguchi and Gonqalves16). Indeed, the Hamiltonians of the n = 3 model and the B a x t e r - W u model have a c o m m o n symmetry: spin inversion in two of the three sublattices. Furthermore, there exists a mapping between the n = 3 model and the 4-state Potts model which becomes exact in the anisotropic limit e -2K3 = 8 , K 2 = Ae, e - + 0 , where A -1 is a temperaturelike variable (the limit is different in the case investigated by Igl6i et al. 7) and VanderzandeS): the roles of K 2 and K 3 are interchanged). The proof of this equivalence can be given by showing that the transfer matrices of both models are identical, apart from a trivial prefactor. We consider a transfer matrix which builds up the lattice in the direction of the strong 3-spin interactions, let us say from the left to the right, by adding three consecutive columns of Ising spins. In each row (labelled i) of spins we define a variable o-i which has the value 1, 2, 3 or 4 if the two rightmost spins in the row are ( + + ) , ( + - ) , ( - + ) or ( - - ) , respectively. Using vector notation, the variables %, i = 1, 2, 3 , . . . , are taken together as ~r. This quantity is a suitable index of the transfer matrix T(tr, tr'). The elements of this matrix can be easily calculated in the first two orders of e. The equivalence with the 4-state Potts model follows by a comparison with the Potts transfer matrix with couplings K x = - I n e, Ky = 4Ae. Since it seems rather unlikely that universality depends on anisotropy, this equivalence gives further support to the classification of the n = 3 model in the 4-state Potts universality class. The combined evidence for this classification leaves little room for doubt. In section 3, we have found a clear answer concerning the nature of the phase transition of the n = 4 model: it is confirmed to be first order. However, the finite-size behaviour of the specific heat makes the transition look "almost second order". This is reminiscent to the situation in the two-dimensional Potts model with a number of states that is not much greater than fourl7). On the other hand, the discontinuity in the magnetization is very pronounced, and the different phases become very long lived for the largest system sizes. In conclusion, we have used very long Monte Carlo simulations to investigate the critical behaviour of two Ising models outside the usual Ising universality class. By means of a special-purpose computer, long simulations were feasible, leading to satisfactory answers to some questions of interest.

CRITICAL BEHAVIOUR OF TWO ISING MODELS

411

Acknowledgements One of the authors (H.B.) thanks Professors J.M.J. van Leeuwen and P.W. Kasteleyn for several useful discussions. Part of this research was supported by the Dutch "Stichting voor Fundamenteel Onderzoek der Materie".

References 1) C. Gruber, A. Hintermann and D. Merlini, Group Analysis of Classical Lattice Systems (Springer, Berlin, New York, 1977), p. 25. 2) J.M. Debierre and L. Turban, J. Phys. A16 (1983) 3571. 3) A. Maritan, A.L. Stella and C. Vanderzande, Phys. Rev. B 29 (1984) 519. 4) R.B. Potts, Proc. Cambr. Phil. Soc. 48 (1952) 106. 5) R.J. Baxter and F.Y. Wu, Phys. Rev. Lett. 31 (1973) 1294. 6) K.A. Penson, R. Jullien and P. Pfeuty, Phys. Rev. B 26 (1982) 6334. 7) F. Igl6i, D.V. Kapor, M. Skrinjar and J. S61yom, J. Phys. A16 (1983) 4067. F. Igl6i, D.V. Kapor, M. Skrinjar and J. S61yom, J. Phys. A19 (1986) 1189. 8) C. Vanderzande and F. Igl6i, to be published. 9) B. Nienhuis, A.N. Berker, E.K. Riedel and M. Schick, Phys. Rev. Lett. 43 (1979) 737. 10) H.W.J. B16te, M.P. Nightingale and B. Derrida, J. Phys. A14 (1981) L45. 11) H.W.J. Bl6te and M.P. Nightingale, Physica l12A (1982) 405. 12) A. Hoogland, J. Spaa, B. Selman and A. Compagner, J. Comp. Phys. 51 (1983) 250. A. Compagner and A. Hoogland, J. Comp. Phys, submitted. 13) M.P. Nightingale, J. Appl. Phys. 53 (1982) 7927 and references therein. 14) M. Nauenberg and D.J. Scalapino, Phys. Rev. Lett. 44 (1980) 837. 15) A. Hoogland and H.W.J. Bl6te, unpublished DISP results. 16) T. Horiguchi and L.L. Gonqalves, Physica 133A (1985) 460. 17) C. Rebbi and R.H. Swendsen, Phys. Rev. B 21 (1980) 4094.