Wang-Landau simulation for the quasi-one-dimensional Ising model
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] M a r Typeset with jpsj2.cls < ver.1.2 > Full Paper
Wang-Landau Simulation for the Quasi-One-Dimensional Ising Model
Takayuki
Tanabe and Kouichi
Okunishi Graduate School of Science, Niigata University, Igarashi 2, 8050 Niigata 950-2181, Japan. Department of Physics, Faculty of Science, Niigata Univeristy, Igarashi 2, 8050 Niigata 950-2181, Japan.
We revisit the nature of the quasi-one-dimensional Ising model on the basis of Wang-Landausimulation. We introduce the density of states in the two-dimensional energy space correspond-ing to the intra- and interchain directions. We then analyze the interchain coupling dependenceof specific heat of the anistropic two-dimensional Ising model in the context of the density ofstates, and further discuss the size dependence of peak temperature. We also discuss the featureof the phase transition in the three-dimensional case.
KEYWORDS: quasi-one-dimension, Wang-Landau simulation, density of states
1. Introduction
The effect of the interchain coupling in the quasi-one-dimensional(Q1D) spin system has recently attractedmuch attention. In general, the 1D spin system showsno phase transition at a finite-temperature. In the Q1Dsystem, however, the weak interchain coupling inducesa finite-temperature phase transition. In fact, the 3Dlong-range ordering for the Q1D system has been along-standing issue. The recent experimental develop-ments enable the precise investigation of 3D ordering ina wide variety of Q1D systems, such as coupled S = 1 / Moreover, very recently, it has beenshown that the Q1D Ising-like XXZ antiferromagnetBaCu V O exhibits an exotic incommensurate spin or-der in a magnetic field,
3, 4 in which the Ising anisotropyplays an essential role.Motivated by the experimental results above, we reex-amine the phase transition of the Q1D Ising system: H = − J X q S i S j − J ′ X ⊥ S i S j , (1)where S ∈ ± < i, j > q indi-cates spin pairs along the chain direction, and < i, j > ⊥ means the pairs perpendicular to the chain. Thus, J and J ′ respectively represent the coupling constants for theintra- and inter-chain interactions. Of course the univer-sality of the Ising transition due to the Z symmetrybreaking itself is independent of the interchain interac-tion. However, quantitative details of the interchain de-pendence of the phase transition still involve an interest-ing problem, which is essential to resolve experimentalresults. Recently, the universal reduction in the effectivecoordination number has also been reported for the Q1DIsing model. The aim of this paper is to understand therole of interchain coupling in the context of the energydensity of states (DOS) based on Wang-Landau simula-tion.
12, 13
The size dependence of the transition temper-ature of the Q1D Ising system is also discussed.For the quantitative analysis of the phase transitionof the Q1D system, recall that the energy scale of inter-chain coupling is much smaller than that of intrachaincoupling, and thus the transition temperature becomesvery low; the conventional metropolis Monte Carlo simu- lation based on a local spin flip often fails in relaxation tothe appropriate equilibrium state. Recently, an efficientcluster algorithm has been proposed for the Q1D sys-tem. In this paper, however, we employ such a gener-alized ensemble method as multicanonical simulation. In particular, the Wang-Landau simulation
12, 13 enablesus to estimate DOS efficiently through a random walk inthe energy space, and to avoid trapping in a metastablestate. Then we can resolve the contribution of typicalconfigurations at low temperatures, which provides anessential viewpoint of the low-temperature behavior ofthe Q1D system.This paper is organized as follows. In the next section,we briefly explain details of the Wang-Landau simulationfor the Q1D system. In particular, we introduce the DOSof two-dimensional energy space for the intra- and inter-chain directions. In §
3, we discuss the phase transitionin the 2D case in the context of DOS and then analyzethe interchain interaction dependence of the phase tran-sition. For 3D case, we also discuss the nature of thephase transition. In §
4, the summary and discussions aregiven.
2. Simulation Details
The Wang-Landau simulation is based on a randomwalk in the energy space without trapping metastablestate and enables us to estimate DOS. For the Q1Dsystem, however, the energy scale of interchain cou-pling is fairly different from that of intrachain coupling.For the purpose of treating such a highly anisotropicsystem more efficiently, we further introduce the two-dimensional energy space defined by e q ≡ X q S i S j , and e ⊥ ≡ X ⊥ S i S j , (2)where e q and e ⊥ respectively denotes the (dimensionless)unbiased energies for the intra- and inter-chain direc-tions. Then the total energy is given by E = − Je q − J ′ e ⊥ .The Wang-Landau simulation itself is performed for thistwo-dimensional space of the spatially isotropic Isingmodel and then DOS g ( e q , e ⊥ ) is obtained in ( e q , e ⊥ )space. The expectation value for various J ′ can be ob-tained through reweighting. Here, it should be noted that Full Paper
Author NameFig. 1. Two-dimensional DOS of the 2D Ising model of systemsize L=10. the multicanonical simulation for the two-dimensionalparameter space was successfully applied to such com-plex systems as the polymer system and the spin glasssystem. Moreover, the Wang-Landau sampling drasti-cally enhances the accessibility to the multi-dimensionalparameter space in various situations. An interestingpoint in the present case may be that the dimensionalityof the parameter space is directly related to the spatialdimension.The detailed conditions for the Wang-Landau simu-lation are given as the follows. At the start of simula-tion, DOS is unknown, so it is simply set at g ( e q , e ⊥ ) =1 for all possible ( e q , e ⊥ ) . Then we begin a randomwalk in ( e q , e ⊥ ) space with a probability proportional to1 /g ( e q , e ⊥ ) . The transition probability from ( e ′ q , e ′ ⊥ ) to( e ′′ q , e ′′ ⊥ ) isProb(( e ′ q , e ′ ⊥ ) → ( e ′′ q , e ′′ ⊥ )) = min " g ( e ′ q , e ′ ⊥ ) g ( e ′′ q , e ′′ ⊥ ) , , (3)and g ( e q , e ⊥ ) is iteratively updated by the modificationfactor f as ln g ( e q , e ⊥ ) → ln g ( e q , e ⊥ ) + ln f, (4)when the state is visited. At the same time, the his-togram is updated like H ( e q , e ⊥ ) → H ( e q , e ⊥ ) + 1. Whenthe histogram becomes “flat”, the modification factor isreduced, ln f → (ln f ) /
2, and we reset the histogramto zero. Then we perform a random walk again. Tocheck the flatness of the histogram, we use the criterion( H max − H min ) / ( H max + H min ) ≤ . ∼ .
3, where H max and H min are the maximum and minimum histogramcounts, respectively. We end the simulation when themodification factor is smaller than a predefined value (weset ln f final = 10 − ). The initial value of the modificationfactor is ln f = 1. The update of g ( e q , e ⊥ ) and H ( e q , e ⊥ )is performed every N spin flips, where N is the numberof spins in the system.
3. Results
Let us first consider the 2D Ising model on a L × L square lattice, the exact solution of which is well known and very useful for verifying the simulation results. The Fig. 2. Specific heat for J ′ /J = 0 . Hamiltonian of the 2D Ising model is written as H = − J L X i,j S i,j S i +1 ,j − J ′ L X i,j S i,j S i,j +1 , (5)where i and j are the indexes of the intra- and inter-chain directions, respectively. The unbiased energies forthe intra- and inter-chain directions are explicitly givenby e q = P Li,j S i,j S i +1 ,j and e ⊥ = P Li,j S i,j S i,j +1 , re-spectively. Since g ( e q , e ⊥ ) = g ( − e q , e ⊥ ) = g ( e q , − e ⊥ ) = g ( − e q , − e ⊥ ) holds, it is sufficient to perform simulationin the region e q ≧ , e ⊥ ≧
0. The system sizes are L = 10 , , ,
40, and 50. The maximum histogramcount per stage is H max = 3024, and the maximumMonte Carlo steps per stage is 1 . × for the L = 10system. Then the total number of stages is 21, and the to-tal CPU time is 2 minutes with a 2.66GHz Core2Duo pro-cessor. For L = 30, H max = 7981, the maximum MonteCarlo steps per stage is 3 . × . The total number ofstages is also 21 and the total CPU time is 55 hours. InFig. 1, we show the typical result of DOS g ( e q , e ⊥ ) for L = 10.On the basis of DOS g ( e q , e ⊥ ), we calculate specificheat; Figure 2 shows the size dependence of the specificheat for J ′ /J = 0 . J ′ /J = 0 .
025 is given as T c /J = 0 . · · · . The resultclearly shows that the peak of C corresponding to thecritical divergence gradually develops for L = 50 in thevicinity of T c . In addition to the critical point, we can alsosee a small peak in the low-temperature region T /J ∼ .
2. As L increases, the peak temperature of this smallpeak shifts to a higher temperature side and the peakheight itself decreases rapidly.In order to see the origin of the low-temperature peak,we show the DOS g ( e q , e ⊥ ) in the low-energy region inFig. 3. Note that the scale of e q is much smaller than thatof e ⊥ ; Since J ′ ≪ J for the Q1D system, the range of thehorizontal axis in Fig. 3 is adjusted by the ratio: e q /e ⊥ ∼ J ′ /J = 1 /
40. The ground state energy is E g = − ( J + J ′ ) L × L and its configuration is illustrated as Fig. 3(a),which is located at ( e q , e ⊥ ) = (100 , E ( b ) = E g + 4( J + J ′ ). For the Q1D system, however, another . Phys. Soc. Jpn. Full Paper
Author Name 3Fig. 3. Two-dimensional DOS of the 2D Ising model of L = 10near the groundstate(88 ≦ e q ≦ , ≦ e ⊥ ≦ J ′ /J = 0 .
025 with L = 10. The solid squares indicate the distribution function for T/J = 0 .
23, which corresponds to the low temperature peakof the specific heat. The open circles indicate that for the hightemperature peak (
T/J = 0 . important excitation we should discuss is “chain flippedexcitation”, a typical example of which is depicted inFig. 3(c), and its energy is given as E ( c ) = E g + 4 J ′ L .The DOS of the chain flipped configuration is located at e ⊥ = 60 , · · · on the edge of e q = 0. In Fig. 1, we canalso confirm that the DOS of these configurations at theedges deviate from the “bulk” value in the ( e q , e ⊥ ) plane.Moreover, note that the “gap” in DOS at every e ⊥ = 20in Fig. 3 also originates from the chain structure of thelattice. Since L < J/J ′ , we can see that the chain flippedexcitation has a lower energy than the single spin flippedstate, but it becomes the predominant excitation at a lowtemperature. As L increases beyond J/J ′ , the energy ofthe chain flipped excitations shifts to the higher-energyregion, so that the contribution of such configurationsdecreases gradually. Thus the low-temperature peak ofspecific heat in Fig. 2 can be well described by the chainflipped excitations, which is peculiar to the Q1D system.In order to see the weight of each energy state in theequilibrium, we calculate the energy distribution func-tion P ( e q , e ⊥ ) = g ( e q , e ⊥ ) exp [( Je q + J ′ e ⊥ ) /T ] (6) for J ′ /J = 0 . T /J = 0 .
23 is the temperature of the low-temperaturepeak of specific heat and
T /J = 0 .
91 corresponds to thehigh-temperature peak of L = 10. In the figure, the pre-dominant contribution at T /J = 0 .
23 clearly comes fromthe states at the edge of e ⊥ = 100, which implies that thechain flipped state is essential for the small peak; For arelatively small system size ( L < J/J ′ ), the energy of theexcitations actually satisfies 4( J + J ′ ) > J ′ L . On theother hand, P ( e q , e ⊥ ) for T /J = 0 .
91 shows a Gaussian-like shape at approximately ( e q , e ⊥ ) ∼ (90 , T /J = 0 .
23) is the chainflipped states near the ground state. This implies thatthe polarization of the spins in the same chain is basi-cally frozen and the aligned spins in the chain can be-have as a single spin, which forms an effective 1D spinchain through the weak interchain coupling LJ ′ in theinterchain direction. Thus, we can see that the fluctu-ation in the interchain direction is predominant for thelow-temperature peak, while for the high temperaturepeak, the fluctuations in both the intra- and inter-chaindirections give the significant contributions. Of course,the low temperature peak is basically a finite-size effectand vanishes in the bulk limit. However, the region wherethe finite-size effect can clearly appear is up to L ∼ J/J ′ ,which is a certain large number for the Q1D system. Thisimplies that the true critical divergence of specific heatis eventually masked by the analytic contribution of thelow-temperature peak, up to L ∼ J/J ′ . Thus, the finite-size scaling analysis based on the data L < J/J ′ shouldbe performed carefully. In order to extract the proper critical behavior for
L < J/J ′ , we examine the decomposition of specificheat into three parts: the fluctuation along the chain, C q ; the fluctuation in the interchain direction, C ⊥ ; andcross term of the intra- and inter-chain directions, C cross . C = ( h E i − h E i ) /N T = C q + C ⊥ + C cross , (7) C q = J ( h e q i − h e q i ) /N T , (8) C ⊥ = J ′ ( h e ⊥ i − h e ⊥ i ) /N T , (9) C cross = 2 JJ ′ ( h e q e ⊥ i − h e q i h e ⊥ i ) /N T . (10)The intra- and inter-chain fluctuations basically reflectthe 1D nature of the spin fluctuations. Thus, for C , C q ,and C ⊥ , we have to see the critical fluctuation in thebackground of the predominant 1D behavior. On theother hand, the cross term C cross exhibits no such 1Dbehavior and thus we can expect more direct observa-tion of the critical fluctuation.In Fig. 5, we show C ⊥ and C cross for J ′ /J = 0 . C q isnot presented here). In the figure, C ⊥ shows a Schottky-like peak for small system sizes ( L < L increases,the peak position shifts to the high-temperature side andthe peak height itself rapidly decreases. This behavior isconsistent with the fact that the interchain fluctuation isthe predominant contribution for L < J/J ′ . Indeed, we J. Phys. Soc. Jpn.
Full Paper
Author NameFig. 5. Decomposed specific heats for J ′ /J = 0 . C ⊥ , and the crossterm of the chain and interchain directions, C cross .Fig. 6. Size dependences of the peak temperatures for C , C q , C ⊥ ,and C cross : (a) J ′ /J = 0 . J ′ /J = 0 . have verified that the low temperature peak of L = 10can be well fitted by the specific heat of the 1D Isingchain of the effective coupling LJ ′ with L = 10. In ad-dition to the low-temperature peak, a broad peak alsoemerges at approximately T /J ∼ .
6, as L increases;this peak corresponds to the critical divergence in thebulk limit. Thus, the crossover of C ⊥ from the effective1D Ising model behavior to the 2D Ising model clearlyappears at L ∼ J/J ′ . On the other hand, the cross term C cross shows the divergence behavior only near the cor-rect critical temperature T c = 0 . · · · . This suggeststhat C cross captures the critical behavior more effectivelythan the total specific heat C .We further analyze the size dependence of the peaktemperatures of the decomposed specific heats C , C q , C ⊥ , and C cross . Let us write the peak temperature forthe system size L as T c ( L ). Then, in the critical regime,peak temperature is expected to follow the size scaling T c ( L ) − T c = AL − /ν , (11) where A is a nonuniversal constant. First we discuss theresults for J ′ /J = 0 .
1, which are illustrated in Fig. 6(a).In the figure, we can see that C , C q , and C ⊥ graduallyapproach T c , for which no scaling behavior is observed.However, the cross term of the specific heat C cross wellsatisfies eq. (11) within a relatively small system size(1 /L < . C cross is more effective forcapturing the critical behavior than C . Another interest-ing feature of the Q1D system is that T c ( L ) approaches T c from the under side of T c , namely, A >
0, for suffi-ciently large L . This behavior is contrasted to A < C s monotonically approach T c from T > T c .The peak temperatures for J ′ /J = 0 . C q monotoni-cally decreases from the upper side of T c , while C ⊥ clearlyexhibits the crossover behavior. For small L ( < . C ⊥ originates from the chain flip config-uration. However, we can see that it rapidly crossoversto that of the critical behavior at 1 /L ∼ .
04. On theother hand, C cross seems to approach T c smoothly, sug-gesting that C cross is more suitable for finite size analysisof the critical behavior. For J ′ /J = 0 . ν . Let us discuss the 3D Ising model along the same lineof argument as that of the 2D case. The Hamiltonian iswritten as H = − J L X i,j,k S i,j,k S i,j,k +1 − J ′ L X i,j,k [ S i,j,k S i +1 ,j,k + S i,j,k S i,j +1 ,k ] , (12)where k is assumed to run in the chain direction. In ac-tual computations, the maximum histogram count perstage is H max = 6024, and the maximum Monte Carlosteps per stage is 6 . × for L = 6 system. The to-tal number of stages is 21 and the total CPU time is 50minutes. For L = 10, H max = 6748, the maximum MonteCarlo steps per stage is 7 . × , and the total CPUtime is about 6 days. Here, we note that, for the 3D Isingmodel, Wang-Landau simulation in the 2D energy space(2) sometimes does not achieve good convergence nearthe edges of the 2D energy space. In such a case, we haveadditionally performed Wang-Landau simulation for theconventional 1D energy space to obtain specific heats.Figure 7(a) shows the size dependence of the specificheat C for J ′ /J = 0 .
025 up to L = 18. In the figure,we can see the broad maximum of C of L = 6 at ap-proximately T /J ∼ .
0. At the same time, a small peakemerges in the low-temperature region, reflecting the 1Dnature of the system. As L increases, this small peakrapidly merges with a broad peak coming down fromthe higher-temperature side. We can then see that themerged peak rapidly develops into a sharp peak associ-ated with the critical divergence at T /J ∼ . . Phys. Soc. Jpn. Full Paper
Author Name 5Fig. 7. C , C q , C ⊥ , and C cross for the 3D Ising model of J ′ /J =0 . We next resolve the peak structure of the total specificheat using C q , C ⊥ , and C cross . Figure 7(b) illustrates C q , which exhibits a broad peak for a small system size.Since J ≫ J ′ , this broad peak can be attributed to thespin fluctuation in the chain, which is governed by theenergy scale J . As L increases, the peak temperaturedecreases to T ∼ . C ⊥ in Fig. 7(c) shows a clear finite-size peak originatingfrom the interchain fluctuation of the energy scale LJ ′ .As L increases, the peak temperature gradually increasesfrom T /J ∼ . .
8. Then, an interesting point about C ⊥ is that the shape of the peak is almost unchangedduring shifting, in contrast to the 2D case where thepeak considerably reduces its shape. Here, let us recallthat, at sufficiently low temperatures, the effective spinsare frozen in the chain direction form the 2D network.Thus, an important difference between the 2D and 3Dcases is that, for 3D, the effective 2D Ising model in thesmall J ′ limit can involve the quasi-critical divergence,while for 2D, the specific heat of the effective 1D Isingmodel shows no such divergence since there is no phasetransition in the 1D Ising model. In Fig. 7(d), we finallyshow C cross , the peak of which develops near T c and issmoothly connected to the critical divergence.In Fig. 8, we summarize the above size dependences ofthe peak temperatures for the specific heats. In the fig- ure, the horizontal axis indicates the scaled system size L − /ν , where we have used ν = 0 . We can see thatthe round peak of C q comes from the higher temperatureside, but it still does not reach the scaling region. On theother hand, we can see that C ⊥ and C cross are well fit-ted by linear functions, which are respectively shown assolid and broken lines in Fig.8. The straightforward ex-trapolation yields T c ≃ .
85, which is consistent with theprecise estimation T c = 0 .
834 based on the simulation upto the size 10 × ×
80 (the result is not presented here).A similar analysis of susceptibility was also reported inRef. 8, where the peak temperature of the intrachain spinfluctuation behaves similarly to C ⊥ . The present resultis consistent with the previous analysis of susceptibility. As can be seen Fig. 7(a), the divergences of C cross and C ⊥ massively contribute to the critical divergence of thetotal specific heat C within a small system size. Thissuggests that C cross can be expected to be suitable forthe finite-size analysis of the critical behavior as well, al-though C exhibits a rather complicated size dependenceof the peak structure in the 3D case. Fig. 8. Size dependences of the peak temperatures of C , C q , C ⊥ ,and C cross for the 3D Ising model of J ′ /J = 0 . L − /ν with ν = 0 .
4. Summary and Discussion
We have studied the feature of the Q1D Ising modelusing Wang-Landau simulation. In order to treat thedifference between the energy scales for the intra- andinter-chain directions, we have particularly introducedthe two-dimensional energy space ( e q , e ⊥ ) correspondingto the intra- and inter-chain directions. We further de-composed the total specific heat C into contributionsfrom intrachain fluctuation, C q , interchain fluctuation, C ⊥ , and the cross term of the intra- and inter-chain fluc-tuations, C cross . Then the finite-size effect peculiar tothe Q1D system is discussed on the basis of the two-dimensional DOS, and it was demonstrated that thechain flip configuration plays an essential role for the low-temperature peak of specific heat. We have also analyzedthe shift exponent of the peak of the specific heat, andthen found that C cross can capture the critical behaviormore effectively than the total specific heat C , within J. Phys. Soc. Jpn.
Full Paper
Author Name a relatively small system size. We have also discussedthe qualitative difference between 2D and 3D cases inthe low-temperature and small- J ′ limits; in the 2D case,the crossover of C ⊥ from an effective 1D chain to a 2Dmodel occurs rapidly at L ∼ J/J ′ . In the 3D case, theeffective 2D Ising model itself involves the quasi-criticaldivergence, because C ⊥ smoothly crossovers from the ef-fective 2D model into the 3D critical behavior.In this paper, we have analyzed the Q1D Ising model inthe context of the two dimensional DOS. The actual com-putational cost for obtaining the two-dimensional DOSincreases rapidly, with increasing system size, and thenthe cluster algorithm seems to be more efficient for asimulation of a larger system. However, the present de-scription based on the two-dimensional DOS provides theessential insight for the qualitative understanding of thelow-energy excitations in the Q1D system. In addition, tosuppress the finite-size effect peculiar to the Q1D system,the aspect ratio usually follows the ratio of anisotropiccorrelation length in analyzing the critical behavior. We can also see that a possible aspect ratio of the sys-tem is L ′ /L = ( D − J ′ /J , where D is the dimension ofthe system, L is the length of a chain, and L ′ is the sizeof the interchain direction. This is because the scale ofthe single-spin flip excitation and the chain flipped statecan be on the same order 4 J + 4( D − J ′ ≃ D − J ′ L for J ≫ J ′ . Acknowledgments
This work is supported by Grants-in-Aid for Scien-tific Research from the Ministry of Education, Culture,Sports, Science and Technology of Japan (Nos. 18740230and 20340096), and by a Grant-in-Aid for Scientific Re-search on Priority Area “High-Field Spin Science in100T”. One of the authors (K.O.) would like to thankT. Suzuki and M. Kikuchi for valuable discussions.
1) L. J. de Jongh, and A. R. Miedem: Adv. Phys. (1974) 1.2) I. Tsukada, Y. Sasago, K. Uchinokura, A. Zheludev, S. Maslov,G. Shirane, K. Kakurai, and E. Ressouche: Phys. Rev. B (1999) 6601.3) S. Kimura, T. Takeuchi, K. Okunishi, M. Hagiwara, Z. He, K.Kindo, T. Taniyama, M. Itoh: Phys. Rev. Lett. (2008)057202; S. Kimura, M. Matsuda, T. Masuda, S. Hondo, K.Kaneko, N. Metoki, M. Hagiwara, T. Takeuchi, K. Okunishi,Z. He, K. Kindo, T. Taniyama, and M. Itoh: Phys. Rev. Lett. (2008) 207201.4) K. Okunishi and T. Suzuki: Phys. Rev. B (2007) 224411;T. Suzuki, N. Kawashima, and K. Okunishi: J. Phys. Soc. Jpn. , (2007) 123707.5) C. Y. Weng, R. B. Griffiths and M. E. Fisher: Phys. Rev. (1967)475.; M. E. Fisher: Phys. Rev. (1967) 480.6) L. L. Liu and H. E. Stanley: Phys. Rev. Lett. (1972) 927;L. L. Liu and H. E. Stanley: Phys. Rev. B (1973) 2279.7) T. Graim and D. P. Landau: Phys. Rev. B (1981) 5156.8) K. W. Lee: J. Phys. Soc. Jpn. (2002) 25919) S. Todo: Phys. Rev. B. (2006) 104415.10) T. Nakamura: Phys. Rev. Lett. (2008) 210602.11) B. A. Berg and T. Neuhaus: Phys. Rev. Lett. (1992) 9.12) F. Wang and D. P. Landau: Phys. Rev. Lett. (2001) 2050.13) F. Wang and D. P. Landau: Phys. Rev. E. (2001) 056101.14) Y. Iba, G. Chikenji, and M. Kikuchi: J. Phys. Soc. Jpn. (1998) 3327.15) N. Hatano and J. E. Gubernatis: Prog. Theor. Phys. Suppl. (2000) 442.16) H. K. Lee, Y. Okabe, and D. P. Landau: Bull. Comput. Phys.Commun. (2006) 36.17) C. Zhou and R. N. Bhatt: Phys. Rev. E. (2005) 025701(R).18) L. H. Onsager: Phys. Rev. (1944) 117.19) In order to extract the proper critical behavior, the aspectratio of the system may be adjusted to the Q1D system. Inthis paper, however, we discuss the size dependence withinsquare or cubic systems.20) H. W. J. Bl¨ote, E. Luijten, and J. R. Heringa: J. Phys. A (1995) 628921) K. Binder and J.-S. Wang, J. Stat. Phys55