We introduce a novel variance-reducing Monte Carlo algorithm for accurate
determination of autocorrelation times. We apply this method to two-dimensional
Ising systems with sizes up to $15 \times 15$, using single-spin flip dynamics,
random site selection and transition probabilities according to the heat-bath
method. From a finite-size scaling analysis of these autocorrelation times, the
dynamical critical exponent $z$ is determined as $z=2.1665$ (12)

A. Linke

B. Bernu

B. Dammann

C. J. Umrigar

C. Münkel

C. P. Yang

D. M. Ceperley

D. Stauffer

F. Wang

G. A. Baker, Jr.

G. F. Mazenko

H. W. J. Blöte

J. R. Heringa

M. P. Nightingale

M. Takahashi

N. Ito

P. Grassberger

R. Matz

S. Kirkpatrick

S. Zhang

W. Feller

W. H. Press

Z. B. Li

English

arXiv

Crossref

Dynamic Exponent of the Two-Dimensional Ising Model and Monte Carlo Computation of the Subdominant Eigenvalue of the Stochastic Matrix

Nightingale, M.P.

Blöte, H.W.J.

NARCIS 

Dynamic exponent of the two-dimensional ising model and Monte Carlo computation of the subdominant eigenvalue of the stochastic matrix

Nightingale, M.P. (author)

Blöte, H.W.J. (author)

TU Delft Repository

We introduce a novel variance-reducing Monte Carlo algorithm for accurate determination of correlation times. We apply this method to two-dimensional Ising systems with sizes up to 15×15, using single-spin flip dynamics, random site selection, and transition probabilities according to the heat-bath method. From a finite-size scaling analysis of these correlation times, the dynamic critical exponent z is determined as z=2.1665(12)

Nightingale, M. P.

Blöte, H. W.J.

DigitalCommons@URI

University of Rhode IslandDigitalCommons@URIPhysics Faculty Publications Physics1996Dynamic Exponent of the Two-Dimensional IsingModel and Monte Carlo Computation of theSubdominant Eigenvalue of the Stochastic MatrixM. P. NightingaleUniversity of Rhode Island, nightingale@uri.eduH. W.J. BlöteFollow this and additional works at: https://digitalcommons.uri.edu/phys_facpubsTerms of UseAll rights reserved under copyright.This Article is brought to you for free and open access by the Physics at DigitalCommons@URI. It has been accepted for inclusion in Physics FacultyPublications by an authorized administrator of DigitalCommons@URI. For more information, please contact digitalcommons@etal.uri.edu.Citation/Publisher AttributionNightingale, M. P., & Blöte, H. W.J. (1996). Dynamic Exponent of the Two-Dimensional Ising Model and Monte Carlo Computationof the Subdominant Eigenvalue of the Stochastic Matrix. Physical Review Letters, 76(24), 4548-4551. doi: 10.1103/PhysRevLett.76.4548Available at: http://dx.doi.org/10.1103/PhysRevLett.76.4548VOLUME 76, NUMBER 24 P HY S I CA L REV I EW LE T T ER S 10 JUNE 1996Dynamic Exponent of the Two-Dimensional Ising Model and Monte Carlo Computationof the Subdominant Eigenvalue of the Stochastic MatrixM. P. NightingaleDepartment of Physics, University of Rhode Island, Kingston, Rhode Island 02881H.W. J. BlöteDepartment of Applied Physics, Delft University of Technology, Lorentzweg 1, 2628 CJ Delft, The Netherlands(Received 16 January 1996)We introduce a novel variance-reducing Monte Carlo algorithm for accurate determination ofcorrelation times. We apply this method to two-dimensional Ising systems with sizes up to 15 3 15,using single-spin flip dynamics, random site selection, and transition probabilities according to theheat-bath method. From a finite-size scaling analysis of these correlation times, the dynamic criticalexponent z is determined as z ­ 2.1665s12d. [S0031-9007(96)00379-1]PACS numbers: 64.60.Ht, 02.70.Lq, 05.50.+q, 05.70.JkThe onset of criticality is marked by a divergence ofboth the correlation length j and the correlation timet. While the former divergence yields singularities instatic quantities, the latter manifests itself notably ascritical slowing down. To describe dynamic scalingproperties, only one exponent is required in addition tothe static exponents. This dynamic exponent z links thedivergences of length and time scales: t , jz . In ourcomputation of z we exploit that, for a finite system, jis limited by the system size L, so that t , Lz at theincipient critical point.In this Letter, we focus on the two-dimensional Isingmodel with Glauber-like dynamics. Values quoted in theliterature for z vary vastly, from z ­ 1.7 to z ­ 2.7 [1],but recent computations seem to be converging towardsthe value reported here. Finally, results are beginningto emerge of precision sufficient for sensitive tests offundamental issues such as universality.The numerical method introduced in this Letter is re-lated to Monte Carlo methods used to compute eigenval-ues of Hamiltonians of discrete or continuous quantumsystems [2,3] and transfer matrices of statistical mechani-cal systems [4]. In particular, the current method is suit-able to obtain more than one eigenvalue by adaptation ofthe diffusion Monte Carlo algorithm of Ref. [5].To compute the correlation time of small L 3 L latticeswe exploit the following properties of the single-spin-flip Markov (or stochastic) matrix P [6]. It operatesin the linear space of all spin configurations and itslargest eigenvalue equals unity. The corresponding righteigenvector contains the Boltzmann weights of the spinconfigurations; the left eigenvector is constant, reflectingprobability conservation. The correlation time tL (in unitsof one flip per spin, i.e., L2 single-spin flips) is determinedby the second-largest eigenvalue lL,tL ­ 21L2 lnlL. (1)For a system symmetric under spin inversion, the corre-sponding eigenvector is expected to be antisymmetric.We used two methods to compute lL: exact, numericalcomputation for L # 5 and Monte Carlo for 4 # L # 15.The exact method used the conjugate gradient algorithm[7] and the symmetries of periodic systems. This calcu-lation resembles that in Ref. [8], but currently we realizeGlauber-like dynamics using heat-bath or Yang [9] transi-tion probabilities and random site selection.The Monte Carlo method used a stochastic form of thepower method, as follows [5]. A spin configuration s withenergy Essd has a probabilityexpf2EssdykT gZ;cBssd2Z, (2)where Z is the partition function. The element Pss0jsdof the Markov matrix is the probability of a single-spin-flip transition from s to s0. Since P satisfies detailedbalance,Pˆss0jsd ; 1cBss0dPss0jsdcBssd (3)is symmetric. For an arbitrary trial state j fl an effectiveeigenvalue lstdL is defined bylstdL ­kPˆt11lfkPˆtlf, (4)where k?lf is the expectation value in the state j fl.In the limit t ! ‘, the effective eigenvalue convergesgenerically to the dominant eigenvalue allowed by thesymmetry of j fl. The convergence is exponential in thetime lag t.Given a trial state j fl, standard Monte Carlo methodsuffices to compute the right-hand side of Eq. (4), i.e.,4548 0031-9007y96y76(24)y4548(4)$10.00 © 1996 The American Physical SocietyVOLUME 76, NUMBER 24 P HY S I CA L REV I EW LE T T ER S 10 JUNE 1996the denominator of Eq. (4),N std ; k fjPˆt j fl ­Xs1,...,st11fsst11dPˆsst111jstd · · · Pˆss2js1dfss1d­Xs1,...,st11fss1dfsst11dcBss1dcBsst11dPsst11jstd · · · Pss2js1dcBss1d2 ­ Z¿fss1dfsst11dcBss1dcBsst11dÀP, (5)is an autocorrelation; fssd ; ksj fl and k?lP denotes theaverage with respect to the probabilityPsst11jstd · · · Pss2js1dcBss1d2yZ (6)of finding a configuration s1 in equilibrium and subse-quent transitions to configurations s2 through st11.Similarly, the numerator of Eq. (4),Hstd ; k fjPˆt11j fl­Xs0,...,st11fsst11dPˆsst11jstd · · · Pˆss1js0dfss0d­12Z¿flLss1d 1 lLsst11dgfss1dfsst11dcBss1dcBsst11dÀP(7)is a cross correlation, where the “configurational eigen-value” lLssd of spin configuration s is defined aslLssd ­1fssdXs0fss0dPˆss0jsd . (8)Finally, with Eqs. (5) and (7), one has lstdL ­ HstdyN stdfor the effective eigenvalue.In practice, Hstd and N std are estimated by conventionalMonte Carlo methods. As usual, these estimators involvetime averages of stochastic variables. Thus, on the rightof Eqs. (5) and (7) si is replaced by st01i21 (i ­ 1, . . . , t),and the Monte Carlo average is taken over an appropriatelychosen subset of times t0 after thermal equilibration.In principle, one could choose f ­ mcB, where m isthe magnetization. In that case, the above method reducesto estimating the effective eigenvalue of the Markovmatrix in terms of the magnetization autocorrelationfunction gstd via lstdL ­ gst 1 1dygstd. To estimate gstdone would average over time products of the formmss1dmsst11d. Equation (7) would then yield gst 1 1dby replacing msstd by the conditional expectation value ofthe magnetization at time t 1 1, evaluated explicitly asPst11 msst11dPsst11jstd.The crux is that the estimator of lstdL satisfies a zero-variance principle [5], since Eqs. (5) and (7) contain anoptimizable trial state j fl. In the ideal case, j fl isan exact eigenstate of the symmetrized Markov matrixPˆ, and the “configurational eigenvalue” lLssd equals theeigenvalue independent of s. Then, the estimator ofthe effective eigenvalue lstdL yields the exact eigenvaluewithout statistical and systematic errors at finite t, if careis taken to arrange cancellation of the fluctuating factors inthe estimators of Hstd and N std. It should be noted that thisis true only if the numerator of Eq. (4) is evaluated withEq. (7), in which the change from t to t 1 1 is made byan explicit matrix multiplication, rather than by using theanalog of Eq. (5) with t replaced by t 1 1. In practice, j flis not an exact eigenstate, and this introduces statistical andsystematic errors. However, these errors are kept small bythe zero-variance principle, if the trial states are accurate.Such optimized trial states are constructed prior to themain Monte Carlo run, by minimization of the variancex2 of the configurational eigenvaluex2spd ­ ksPˆ 2 kPˆlfd2lf . (9)As indicated, the variance depends on the parameters pof the trial state. Optimization over p is done followingUmrigar, Wilson, and Wilkins [10]: one samples M con-figurations si , typically a few thousand, with probabilityc2BZ21 and approximates x2spd byx2spd øPMi­1f fssi , pdycBssidg2 flLssi , pd 2 l¯Lspdg2PMi­1f fssi , pdycBssidg2.(10)Here l¯L denotes the weighted average of the configura-tional eigenvalue over the sample, while the modified no-tation explicitly shows dependences on the parameters pof the trial state j fl. Near-optimal values of the parame-ters p can be obtained by minimization of the expres-sion on the right-hand side of Eq. (10) for a fixed sample.Statistical independence in the sample requires that theconfigurations be selected at intervals on the order of thecorrelation time.A guiding principle for the construction of trial statesis that long-wavelength fluctuations of the magnetizationhave the longest decay time. Furthermore, analysis ofthe exact left eigenvectors of the Markov matrix P forsystems with L # 5 shows that the elements dependonly on the magnetization to good approximation. Thissuggests trial functions depending on long-wavelengthcomponents of the Fourier transform of si , the zero-momentum component of which is just the magnetizationm. The formfssd ­ c˜Bssdc s1dssdcs2dssd , (11)where c s6d ! 6c s6d under spin inversion, yields anantisymmetric trial function, as required. The tilde in4549VOLUME 76, NUMBER 24 P HY S I CA L REV I EW LE T T ER S 10 JUNE 1996c˜B indicates that the temperature is used as a variationalparameter, but we found that its optimal value is virtuallyindistinguishable from the true temperature. The c s6dwere chosen asc s1d ­Xkaksm2dms1dk 1 mXkbksm2dms2dk , (12)c s2d ­ mXkcksm2dms1dk 1Xkdksm2dms2dk , (13)where the index k runs through a small set of multiplets offour or fewer long-wavelength wave vectors defining thems6dk , translation and rotation symmetric sums of prod-ucts of Fourier transforms of the local magnetization;the k are selected so that ms2dk is odd and ms1dk is evenunder spin inversion; the coefficients ak, bk, ck, anddk are polynomials of second order or less in m2. Thedegrees of these polynomials were chosen so that no termsoccur of higher degree than four in the spin variables. Weused trial functions dependent on system size only in theoptimal values of the parameters. This yielded a x2 andan error in the variational estimate ls1dL decreasing with L;yet, the relative error in tL increases.Since the probability distribution Eq. (6) is preciselythe one purportedly generated by standard Monte Carlomethod, the sampling procedure is straightforward. TheMonte Carlo algorithm used a random-number generatorof the shift-register type. It was selected with care toavoid the introduction of systematic errors; see discussionand references in Ref. [11]. We used two Kirkpatrick-Stoll [12] generators, the results of which were combinedby a bitwise exclusive or [13]. For test purposes we re-placed one Kirkpatrick-Stoll generator by a linear congru-ential rule, but this did not reveal clear differences [11].For each system size 4 # L # 15, Monte Carlo aver-ages were taken over 8 3 108 spin configurations. ForL ­ 13 15 these were separated by intervals of 16 sweeps(Monte Carlo steps per spin); 8 sweeps for L ­ 11 and 12;2 sweeps for L ­ 5 and 6; and only one sweep for L ­ 4.The simulations of the remaining system sizes consisted ofparts using intervals of 2, 4, or 8 sweeps.The numerical results for the effective second largesteigenvalue lstdL as a function of the projection time t ap-peared to converge rapidly. In agreement with scaled re-sults for L # 5 spectra, we observe that convergence oc-curs within a few intervals as given above. Monte Carloestimates of lL are shown in Table I, as are exact resultsfor small systems. For system sizes L ­ 4 and 5, the twotypes of calculation agree satisfactorily. The small numer-ical errors indicate that the variance-reducing method in-troduced above is quite effective.For finite system size L there are corrections tothe leading scaling behavior tL , Lz . In the two-dimensional Ising model corrections to static equilibriumquantities occur with even powers of 1yL [14]; thus weTABLE I. Second-largest eigenvalue lL of the Markov ma-trix. The first column indicates the method: exact numerical orMonte Carlo.Method L lL ErrorExact 2 0.985 702 260 395 516 0.000 000 000 001Exact 3 0.997 409 385 126 011 0.000 000 000 001Exact 4 0.999 245 567 376 453 0.000 000 000 001Exact 5 0.999 708 953 624 452 0.000 000 000 001MC 4 0.999 245 568 5 0.000 000 009 4MC 5 0.999 708 945 3 0.000 000 006 0MC 6 0.999 865 719 4 0.000 000 004 5MC 7 0.999 929 970 8 0.000 000 003 1MC 8 0.999 960 085 4 0.000 000 002 3MC 9 0.999 975 663 0 0.000 000 001 7MC 10 0.999 984 357 7 0.000 000 001 4MC 11 0.999 989 505 6 0.000 000 001 0MC 12 0.999 992 710 7 0.000 000 000 8MC 13 0.999 994 784 0 0.000 000 000 6MC 14 0.999 996 173 6 0.000 000 000 5MC 15 0.999 997 131 4 0.000 000 000 5expecttL ø LzncXk­0akL22k , (14)where the series was arbitrarily truncated at order nc, butother powers of 1yL might occur as well. Ignoring thelatter, we fitted the correlation times of Table I to thisform. Typical results of such fits are given in Table II.The smallest systems do not fit Eq. (14) well, at least notfor the nc values used. The residuals decrease rapidlyTABLE II. Results of least-squares fits for the dynamicexponent. The first column shows the minimum system sizeincluded, the second the number of correction terms included,and the third column whether ( y) or not (n) numerical exactresults (for L # 5) are included. The last column contains thechi-square confidence index [16].L $ nc Exact z Error Q4 1 n 2.1769 0.0001 0.005 1 n 2.1705 0.0002 0.006 1 n 2.1688 0.0003 0.237 1 n 2.1679 0.0006 0.438 1 n 2.1672 0.0010 0.424 2 n 2.1650 0.0003 0.175 2 n 2.1665 0.0006 0.706 2 n 2.1662 0.0013 0.607 2 n 2.1648 0.0024 0.524 3 n 2.1672 0.0009 0.645 3 n 2.1656 0.0020 0.616 3 n 2.1625 0.0044 0.563 3 y 2.1653 0.0004 0.344 3 y 2.1670 0.0009 0.645 3 y 2.1657 0.0020 0.494550VOLUME 76, NUMBER 24 P HY S I CA L REV I EW LE T T ER S 10 JUNE 1996TABLE III. Comparison of recent results for the dynamicexponent z. Numerical errors are in parentheses.Reference Year ValuePresent work 1996 2.1665 (12)Li et al. [15] 1995 2.1337 (41)Linke et al. [17] 1995 2.160 (5)Grassberger [18] 1995 2.172 (6)Wang et al. [19] 1995 2.16 (4)Baker and Erpenbeck [20] 1994 2.17 (1)Ito [21] 1993 2.165 (10)Dammann and Reger [22] 1993 2.183 (5)Matz et al. [23] 1993 2.35 (5)Münkel et al. [24] 1993 2.21 (3)Stauffer [25] 1993 2.06 (2)when the minimum system size is increased and theconsistency between the results for different nc suggeststhat Eq. (14) captures the essential scaling behavior of tL.From these results we chose the entry for L $ 5 and nc ­2 as our best estimate: z ­ 2.1665 6 0.0012, where weconservatively quote a 2s error. To our knowledge, thisis the most precise estimate of z obtained to date, asevidenced by recent results summarized in Table III. Thetable shows that the mutual consistency of the resultsfor z has tended to improve in recent years. The onlyrecent result that appears inconsistent with ours is dueto Li, Schulke, and Zheng [15]. Its error is copied fromTable I of Li et al. The data in that table display finite-size dependences that exceed the quoted errors, whichmay explain the discrepancy with our result.This research was supported by the (U.S.) NationalScience Foundation through Grant No. DMR-9214669, bythe Office of Naval Research and by the NATO throughGrant No. CRG 910152. This research was conducted inpart using the resources of the Cornell Theory Center,which receives major funding from the National ScienceFoundation (NSF) and New York State, with additionalsupport from the Advanced Research Projects Agency(ARPA), the National Center for Research Resources atthe National Institutes of Health (NIH), IBM Corporation,and other members of the center’s Corporate ResearchInstitute.[1] See, e.g., references in G. F. Mazenko and O. T. Valls,Phys. Rev. B 24, 1419 (1981), and in M.-D. Lacasse,J. Viñals, and M. Grant, Phys. Rev. B 47, 5646 (1993).[2] S. Zhang, N. Kawashima, J. Carlson, and J. E. Gubernatis,Phys. Rev. Lett. 74, 1500 (1995), and references therein.[3] M. Takahashi, Phys. Rev. Lett. 62, 2313 (1989).[4] M. P. Nightingale, E. Granato, and J.M. Kosterlitz, Phys.Rev. B 52, 7402 (1995); M. P. Nightingale and H.W. J.Blöte, Phys. Rev. B (to be published), and referencestherein.[5] D.M. Ceperley and B. Bernu, J. Chem. Phys. 89, 6316(1988); B. Bernu, D.M. Ceperley, and W.A. Lester,J. Chem. Phys. 93, 552 (1990).[6] See, e.g., W. Feller, An Introduction to Probability Theoryand its Applications (John Wiley and Sons, New York,1968), Vol. 1; H. Haken, Synergetics: an Introduction:Nonequilibrium Phase Transitions and Self-Organizationin Physics, Chemistry, and Biology (Springer-Verlag, Ber-lin, 1978).[7] M. P. Nightingale, V. S. Viswanath, and G. Müller, Phys.Rev. B 48, 7696 (1993).[8] M. P. Nightingale and H.W. J. Blöte, Physica (Amster-dam) 104A, 352 (1980).[9] C. P. Yang, Proc. Symp. Appl. Math. 15, 351 (1963).[10] C. J. Umrigar, K. G. Wilson, and J.W. Wilkins, Phys. Rev.Lett. 60, 1719 (1988); in Computer Simulation Studies inCondensed Matter Physics: Recent Developments, editedby D. P. Landau, K. K. Mon, and H. B. Schüttler, SpringerProceedings of Physics Vol. 33 (Springer, Berlin, 1988);C. J. Umrigar, Int. J. Quant. Chem. Symp. 23, 217 (1989).[11] H.W. J. Blöte, E. Luijten, and J. R. Heringa, J. Phys. A28, 6289 (1995).[12] S. Kirkpatrick and E. P. Stoll, J. Comput. Phys. 40, 517(1981).[13] J. R. Heringa, H.W. J. Blöte, and A. Compagner, Int.J. Mod. Phys. C 3, 561 (1992).[14] See, e.g., H.W. J. Blöte and M. P.M. den Nijs, Phys. Rev.B 37, 1766 (1988), and references therein.[15] Z. B. Li, L. Schulke, and B. Zheng, Phys. Rev. Lett. 74,3396 (1995).[16] W.H. Press, B. P. Flannery, S.A. Teukolsky, and W.T.Vetterling, Numerical Recipes (Cambridge UniversityPress, Cambridge, 1986 and 1992), Sec. 6.2.[17] A. Linke, D.W. Heermann, P. Altevogt, and M. Siegert,Physica (Amsterdam) 222A, 205 (1995).[18] P. Grassberger, Physica (Amsterdam) 214A, 547 (1995).[19] F. Wang, N. Hatano, and M. Suzuki, J. Phys. A 28, 4543(1995).[20] G. A. Baker, Jr. and J. J. Erpenbeck, in Computer Simu-lation Studies in Condensed Matter Physics, edited byD. P. Landau, K. K. Mon, and H. B. Schüttler, SpringerProceedings in Physics Vol. 78 (Springer-Verlag, Berlin,1994), Vol. VII.[21] N. Ito, Physica (Amsterdam) 196A, 591 (1993).[22] B. Dammann and J. D. Reger, Europhys. Lett. 21, 157(1993).[23] R. Matz, D. L. Hunter, and N. Jan, J. Stat. Phys. 74, 903(1993).[24] C. Münkel, D.W. Heermann, J. Adler, M. Gofman, andD. Stauffer, Physica (Amsterdam) 193A, 540 (1993).[25] D. Stauffer, J. Phys. A 26, L599 (1993).4551

https://digitalcommons.uri.edu/cgi/viewcontent.cgi?article=1238&amp;context=phys_facpubs