PHYSICAL REVIEW B VOLUME 43, NUMBER 1 1 JANUARY 1991 Numerical determination of the order of phase transitions in Ising systems with multispin interactions Nestor Caticha* Instituto de Fisica e Quimica de Sao Carlos da Universidade de Sao Paulo, 13560 Sao Carlos, Sao Paulo, Brazil Jorge Chahine Insti tuto de Biociencias, I.etras e Ciencias Exatas da Uni Uersidade Estadual Paulista, 15100Sao Jose do Rio Preto, Sao Paulo, Brazil J. R. Drugowich de Felicio Instituto de Fisica e Quimica de Sao Carlos da Universidade de Sao Paulo, 13560Sao Carlos, Sao Paulo, Brazil (Received 14 June 1990) The method of the fourth-order cumulant of Challa, Landau, and Binder is used together with the Monte Carlo histogram technique of Ferrenberg and Swendsen to study the order of the phase tran- sitions of two-dimensional Ising systems with multispin interactions in the horizontal direction and two-body interactions in the vertical direction. Drastic modifications in the behavior of statistical- mechanical models can occur when multispin interac- tions are present. The best known cases where this hap- pens are the Baxter and Ashkin-Teller models. There are, however, other instances of these types of models that have deserved attention in the past decade. ' These include three-, four-, and more-body terms in d 2. The usual Landau expansion lore holds that three-body terms in dimension larger than two will lead to first-order transitions. The case, however, is not total- ly clear in some two-dimensional (2D) problems. Methods of study have included mean-field theories, ' Monte Carlo simulations, ' and continuous-time formu- lations. As an example, Turban and Debierre' have looked at an anisotropic Ising-like model in two dimen- sions (5 =+1 ), with a Hamiltonian given by m —1 As shown in Ref. 1, the above Hamiltonian is self-dual for any m, which proves that as long as it has a single transition the critical point will be located at the self-dual point, which is given by the known relation sinh2K sinh2K„= 1, independent of the value of m. From the fact that the ground state is 2 ' times degen- erate, as can be seen from the spin reversal invariance of the Hamiltonian, as well as from an analysis of the energy of domain walls, one can expect that these models will be on the same universality class as the q-state Potts model whenever q =2 ', so that the m = 3 and the q =4 Potts should be on the same universality class, as well as the Baxter-Wu model. To check these conjectures, at least at the numerical level, a precise determination of the order of the transitions has to be achieved, a question that we will address in what follows for the cases m =3 and 4. The first attempts to study these problems used mean- field theory' and finite-size scaling. ' Further improve- ments of the mean-field methods and Monte Carlo simu- lations helped in pointing to a first-order transition in the m =4 case, by detecting quite clearly the energy di6'erences between the two di6'erent phases. This issue was, however, not totally settled in a definite way for the m = 3 case, presumed to be of second order. In this paper we intend to study these problems by us- ing the method proposed by Ferrenberg and Swendsen, which combines their histogram method with the cumu- lant technique of Challa, Landau, and Binder. The cumulant technique to determine the order of phase transitions takes advantage in a smart way of the fact that the energy probability distribution of the finite- size system (L XL), PL(E), has two peaks for the first- order transition, while it only has one in the second-order case. As the size of the system increases, these probabili- ty distributions tend to two thin Gaussians at E+ and E in the first case, while they tend to a single one, cen- tered at, say, Eo in the second case. It follows then, that for a second-order transition, and in a region around it, the energy expected values satisfy (E )=(E ) =Eo in the thermodynamic limit. Although this relation is also satisfied away from a first-order transition point in the thermodynamic limit, because basically only one of the peaks of the distribution contributes, the two-valued na- ture of PL is felt exactly at the transition point. Thus, ( E4) (E+4+E —4) &1 for I (E ) (E 2++E —2)2 We consider the fourth-order cumulant, (E'), VL =1 3(E2)2 43 1173 1991 The American Physical Society 1174 BRIEF REPORTS 43 (3) where JV, the ratio of the partition functions at the two temperatures, ensures the proper normalization and can be calculated from that condition. In general, a multidi- mensional histogram has to be stored if any complex ex- pectations are to be evolved in parameter space. But, for the problem in hand, there are only two couplings (K,K, ), and we stored a double histogram to obtain the distribution PI (E~,EH) of the vertical (two-body) and horizontal (m-body) parts of the energy, so that K, and E could be varied independently. The simulations were performed at the self-dual symmetrical point. Figures 1(a) and 1(b) show typical histograms for the m =3 and 4 which tends to —, in the thermodynamic limit for any tem- perature at and around the transition for the second- order case and to —', also in the case of the first order, everywhere except at the transition point, where a small- er nontrivial limit is expected. The numerical determination of such a cumulant is particularly well suited for the Monte Carlo histogram method, since only low-dimensional arrays have to be stored. As a reminder of the histogram method we men- tion very brieAy that during a usual Monte Carlo run (e.g. , heat bath) a histogram of the values of the energy is stored. Once normalized, this gives an approximation to the energy probability distribution, Pr (E, T), at the tem- perature T where the simulation was performed. It can be thought of arising from two different contributions. The first is a hard to determine geometrical (entropic) part, and the second is a temperature-dependent part that is essentially a Boltzmann-type exponential. It can be evolved in temperature by using the well-known formula PI (E, T')=APL(E, T) exp[E(T ' —T' ')], cases, respectively. We simulated these models in L XL lattices with L from 12 up to 36, with periodic boundary conditions. A heat bath, multispin coding algorithm was used, and runs were typically (5—4) X 10 Monte Carlo Steps (MCS). Thermalization was assumed after looking at the correlation times, and that lead to discarding the initial (1—3) X 10 MCS. Once in possession of the histo- grams the cumulants were then calculated quite easily over a range of temperatures. Of course simulations are done at finite size and then the cumulants tend to —', quite far from the transition point. The cumulants graphs for several system sizes L appear in Figs. 2(a) and 2(b) for m =3 and 4, respectively. The evolution of the histo- gram data is, of course, not very trustworthy on a very large interval of temperatures. We believe that this might be the cause of the additional structure, which ap- pears in the form of shoulders to the left of the minima, and we attribute no physical cause to them. It is hard to tell just from these pictures if there is any difference at all between the phase transitions of the two models. Howev- er, this difference can be established by looking at the se- quence of minima of Vl as a function of L. These values of the minima were fitted by using a power law as sug- gested by finite-size scaling of the form 0.67 VL 0.65 0.63 P(E) I 906 608 ~y ~ ~ ~ ~ ~ ~,'\ ~ ~ 0 ~ ~ ~~ ~ ~ ~~ ~ 0 ~ ~ 0 ~ ~ ~ ~ ~y ~ ~ ~ ~ ~ ~ ~ ~~ ~ ~0 ~ ~g ~ 0 I I 96 VL 0.6 I 0.59 I I I 0.20 0.25 0.30 0.35 0.67 "x 0.40 I I I 0.45 0.50 0.55 P IE) I 537 ~ ~~ ~ ~ ~ ~ ~ ~ ~ ~ ~0 ~ ~ 4 ~~ ~ ~ ~ ~ 0 ~ ~ %~ ~ 0 ~ ~ 0 ~ ~~ ~ ~ ~ ~ ~ ~ ~ ~y ~~ ~ ~ ~ ~ ~ +ye 0.65 0.63 L= 28 L= 24 L= 20 0.6 I L= I6 ~oP I ~+ 'f 532 4 ~ I 248 0.59 0.20 I I I I I 0.25 0.30 0,35 0.40 0.45 Kx L = l2 I 0.50 0.55 FIG. 1. Energy probability distribution: (a) m =3 and {b) m =4. FIG. 2. Fourth-order cumulants: {a) m = 3 and (b) m =4. 43 BRIEF REPORTS 11'75 VL=V (1 a—L ) . (4) In both cases T =0.4409 is the best value obtained by a least-squares fitting in almost perfect agreement with the exact value T, =0.4406867. . . . This fitting also fur- The best value obtained for V characterizes the order of the transition to be second for m =3, V =0.667= 3, quite consistently within the estimated errors which affect only the last digit. For m =4 we can surely assert that the transition is of first order, since nontrivially V„=0.649%—', , with the same errors. Also one can ex- trapolate the value of K at which the minima occur, by supposing as well a power-law convergence to the infinite volume transition temperature of the form TI =T (1—a'L ) . (5) nishes the value of b'=1.66 (1.94) for m =3 (m =4) to be compared with 1/v= l. 5 in the first case and d =2 in the second one. Whereas these are not terrific numbers, they are quite good when considered to be spinoffs of the method. One can finally extract a latent heat I =E+—F. , for the m =4 case. A small error is made if one assumes that the two peaks of Fig. 1(b) contribute equally to the aver- age energy, E =(E++E )/2, but then we can estimate E+ and E to determine l =0.18+0.02 in good agree- ment with the previously reported value of 0.20+0.05. In conclusion we have shown the feasibility of using the cumulant method in association to the Monte Carlo histogram technique to determine the order of the phase transitions occurring in a class of two-dimensional spin systems with multispin interactions. Present address: Instituto de Fisica, Universidade de Sao Pau- lo, Brazil. L. Turban and J. M. Debierre, J. Phys. A 16, 3571 (1983). R. Bidaux and N. Boccara, Phys. Rev. B 34, 4881 (1986}. ~F. Alcaraz and M. N. Barber, J. Phys. (N.Y.) A 20, 179 (1987). J. R. Heringa, H. %. J. Blote, and A. Hoogland, Phys. Rev. Lett. 63, 1546 (1989). 5K. A. Penson, R. Julien, and P. Pfeuty, Phys. Rev. B 26, 6334 (1982); F. Igloi, D. V. Kapor, M. Skrinjar, and J. Solyom, J. Phys. (N.Y.) A 16, 7067 (1983). ~A. Maritan, A. Stella, and C. Vanderzande, Phys. Rev. B 29, 519 (1984). 7F. C. Alcaraz, Phys. Rev. 8 34, 4885 (1986). SA. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 63, 1195 (1988). 9M. S. Challa, D. P. Landau, and K. Binder, Phys. Rev. B 34, 1841 (1986).