Compelling evidence for the theory of dynamic scaling in first-order phase transitions
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] A p r Compelling evidence for the theory of dynamic scaling in first-order phase transitions
Fan Zhong
State Key Laboratory of Optoelectronic Materials and Technologies and School of Physics,Sun Yat-sen University, Guangzhou 510275, People’s Republic of China (Dated: July 9, 2018)Matter exhibits phases and their transitions. These transitions are classified as first-order phasetransitions (FOPTs) and continuous ones. While the latter has a well-established theory of therenormalization group, the former is only qualitatively accounted for by classical theories of nucle-ation, since their predictions often disagree with experiments by orders of magnitude. A theory tointegrate FOPTs into the framework of the renormalization-group theory has been proposed butseems to contradict with extant wisdom. Here we show first that classical nucleation and growththeories alone cannot explain the FOPTs of the paradigmatic two-dimensional Ising model drivenby linearly varying an externally applied field. Then we offer compelling evidence that the transi-tions agree well with the renormalization-group theory when logarithmic corrections are properlyconsidered. This unifies the theories for both classes of transitions and FOPTs can be studied usinguniversality and scaling similar to their continuous counterpart.
Matter as a many-body system exists in various phasesand/or their coexistence and its diversity comes fromphase changes. It thus exhibits just phases and theirtransitions. These transitions are classified as first-orderphase transitions (FOPTs) and continuous ones. Al-though the phases can be studied by a well-establishedframework and the continuous phase transitions havea well-established theory of the renormalization group(RG) that has predicted precise results in good agree-ment with experiments, the FOPTs gain a different sta-tus in statistical physics. They proceed through eithernucleation and growth or spinodal decomposition [1–3]. Although classical theories of nucleation [4–14] andgrowth [15–17] correctly account for the qualitative fea-tures of a transition, even an agreement in the nucle-ation rate of just several orders of magnitude betweentheoretical predictions and experimental and numericalresults is considered as a feat [11–14, 18–20]. A lot ofimprovements have thus been proposed and tested in thetwo-dimensional (2D) Ising model whose exact solutionis available. One theory of nucleation, called FT here-after, considers field theoretic corrections to the classicaltheory [21, 22]. Its field dependence was quantitativelyverified for a constant applied magnetic field H that di-rects oppositely to the equilibrium magnetization M eq at a temperature T below the critical temperature T c by Monte Carlo simulations [23]. By employing the re-sults of such relaxation processes, FT was also shown toaccurately produce numerical results of hysteresis loopareas in a single droplet (SD) regime in which only asingle droplet nucleates and grows quickly throughoutthe system [24]. So was in a multidroplet (MD) regimewhere many droplets nucleate and grow even in the caseof a sinusoidally varying H by using Avrami’s growthlaw [15, 25] and an adiabatic approximation [26]. In thisregime, an adjustable parameter was needed to match thearea of just one frequency but then yielded good resultsfor others [26]. Another theory, referred to as BD below, adds appropriate corrections to the droplet free energy ofBecker and D¨oring’s nucleation theory [6]. Such a theorywas found to accurately predict nucleation rates for the2D Ising model without adjustable parameters [27, 28].However, it is well-known that classical nucleation the-ories are not applicable in spinodal decompositions inwhich the critical droplet for nucleation is of the lat-tice size and thus no nucleation is needed [1]. Althoughsharply defined spinodals that divide the two regimes ofthe apparently different dynamic mechanisms do not ex-ist for systems with short-range interactions contrary tothe mean-field case which has long-range interactions [1–3], it is generally believed that there exists a crossoverregion between them at least at the early stage of anFOPT for systems with short-range interactions [1–3].One may then characterize this crossover by fluctuation-shifted mean-field spinodals and expand near such in-stability points below T c of a usual φ theory that de-scribes the critical behavior of the Ising model. Thisresults in a φ theory for the FOPT due to the lack ofthe up–down symmetry in the expansion [29, 30]. AnRG theory for the FOPT can then be set up in parallelto that for the critical phenomena, giving rise to univer-sality and dynamic scaling characterized by “instability”exponents corresponding to the critical ones. The pri-mary qualitative difference is that the nontrivial fixedpoints of such a theory are imaginary in values and arethus usually considered to be unphysical, though the in-stability exponents are real. Yet, it is later shown thatcounter-intuitively imaginariness is physical in order forthe φ theory to be mathematically convergent, since atthe instability points, the unstable degrees of freedom ofthe system flows to the fixed points upon coarse grain-ing [31]. Moreover, the other degrees of freedom thatneed finite free energy costs for nucleation are coarse-grained away with the costs and are thus irrelevant tothe transition [31]. This indicates that nucleation is ir-relevant to the scaling. Although no clear evidence ofan overall power-law relationship was found for the mag-netic hysteresis in a sinusoidally oscillating field in two di-mensions [24, 26, 32], recently, with properly logarithmiccorrections a dynamic scaling near a temperature otherthan the equilibrium transition point T was found for thecooling FOPTs in the 2D Potts model [33]. This resultshows that spinodal-like dynamic scaling does exist forFOPTs in systems with short-range interactions if loga-rithmic corrections are properly considered. However, inthat case only one hysteresis exponent found numericallyis consistent with a similar theory [34].Here we first compare results arsing from both FT [24,26] and BD [27, 28] and numerical simulations of the2D Ising model. We see that both the theories agreequite well generally with the numerical results. However,the slight but systematic deviations for different sweep-ing rates of the external driving indicate that the theoriesalone cannot explain such a driven transition. Then wefind good agreement with the RG theory of FOPTs in-cluding instability exponents and even scaling forms aswell as existence of finite instability points for two dif-ferent T below T c after account of additional logarithmiccorrections. This offers compelling evidence for the the-ory and thus one can study the universality and scalingof FOPTs similar to their continuous counterpart. Finite-time scaling
Crucial in our analysis is the the-ory of finite-time scaling (FTS) [35, 36]. We drive theFOPT by linearly rather than sinusoidally varying H .This linear driving is a direct implementation of theFTS [35, 36], whose essence is a constant finite timescale associating with the given sweeping rate R of thefield. This single externally imposed time scale can thusprobe effectively the transition when it is of the orderof the nucleation time. In contrast, the sinusoidal driv-ing has two controlling parameters, the field amplitude H and the frequency ω , and thus complicates and con-ceals the essence of the process [37]. In particular, at afixed H , for ω →
0, the hysteresis loop area is governedby H ω , which is equivalent to R , and increases with ω ;while for ω → ∞ , the area is determined by H /ω inmean field and vanishes [38]. At least these two mecha-nisms compete and produce an area maximum at some ω [24, 26, 38, 39]. In addition, for high ω , the hysteresisloops are rounded and even not close and thus their ar-eas are not well defined [24]. This shortcoming does notcontaminate the linear driving [40, 41]. Deficiency of nucleation theories for driving
InFT [24, 26], if a positive constant H is applied against − M eq , the field-theoretically corrected nucleation rate I ( T, H ) per unit time and volume is given by [21, 22] I = B ( T ) H K e − F c /k B T = B ( T ) H K e − Ξ /H (1)with Ξ = Ω σ / M eq k B T (see Supplemental material fordetails), where F c is the free-energy cost for the criticalnucleus, B ( T ) is a parameter, K = 3 for the 2D kinetic Ising model [21–23, 42], Ω d ( T ) is a shape factor in a d -dimensional space, σ is the surface tension along a prim-itive lattice vector, and k B is Boltzmann’s constant.In the MD regime, Avrami’s growth law [15] gives themagnetization M at time t as [15, 25, 26] M ( t ) = 1 − ( − Ω d Z t I (cid:20)Z tt n v ( t ′ ) dt ′ (cid:21) d dt n ) , (2)where v ( t ) is the interface velocity of a growing droplet. v ≈ gH θ with θ = 1 and a constant proportionality g inthe Lifshitz-Allen-Cahn approximation [1, 43, 44].For a time-dependent field H ( t ) = Rt , by assumingan adiabatic approximation in which the constant field issimply replaced with its time dependent one [26], Eqs. (1)and (2) then result in Γ( − , x ) /x − Γ( − , x ) /x +Γ( − , x ) = 4 R ln 2 / [Ω g B ( T )Ξ ] with x ≡ Ξ /H c intwo dimensions, where the coercivity H c is the field at M = 0 and Γ is the incomplete gamma function. Anidentical equation has been derived for the sinusoidaldriving in the low frequency approximation [26] in which R = H ω ≡ πH / [ τ ( H , T ) R ] with τ ( H , T ) being theaverage lifetime of the metastable state at H and T [26].In the SD regime [23], by neglecting the growth time fora supercritical nucleus to occupy half the system volume L d compared with the nucleation time, the probabilityfor the system to make the transition by time t is [24] P ( t ) = 1 − exp (cid:20) − L d Z t I ( T, H ) dt (cid:21) . (3)Accordingly, H c is approximately given by the time t c atwhich P ( t c ) = 1 /
2. Using again the adiabatic approxi-mation for I , one obtains in this regime in two dimen-sions [24] Γ( − , x ) /x = R ln 2 / [ B ( T ) L Ξ ] ≡ CR .In BD, on the basis of the Becker-D¨oring theory ofnucleation [6, 27, 28], the nucleation rate can also be castin the form of Eq. (1) but with a complicated B ( T, H )that is H dependent (see Supplemental material). H c inthe MD and SD regimes can then be found similar to FT.An asymptotic form H c ∼ [ − ln( CR )] − can be foundby expanding Γ( a, x ) in large x in the SD regime [24, 26].This was argued to be the leading behavior for small R [32]. However, it has been shown that such a behaviorif exists could only be detected for extremely low R [24,26], as seen by the curves marked asymptotic logarithmin Fig. 1(a). We shall thus not pursue it.Figure 1 shows the simulation results (see Supplemen-tal material for detailed method) along with theoreti-cal ones from solving numerically the relevant equationsand their BD counterparts. Using the values of H c at R = 200 in the linear driving, we find B ( T ) = 0 . . -4 -2 -0.10.0 (b) sin FT BD BD-FT 256 H c (a) linear sin FT BD H c DS MFS S D A sy m p t o t i c l oga r i t h m (c) -3 -2 -1 Linear FT BD 256 (d) H c FIG. 1. (Color online) (a) H c versus scaled sweep rate R .Linear and sin indicate the data obtained numerically fromthe 2D Ising model using a linearly and a sinusoidally vary-ing external field, respectively. Note that the “error bars” givethe standard deviations of the distributions of the transitioninvolved [26]. The three curves around SD are theoretical re-sults for the single-droplet regime [one BD and two FT curveswith B ( T ) = 0 . B ( T ) = 69 .
73 for thelower] and the two lower curves are results of the asymptoticlogarithmic approximation [the results of the larger B ( T ) arefar smaller and absent]. The horizontal lines with arrows indi-cate the dynamic spinodal (DS) and the mean-field spinodal(MFS) [23, 47]. (b) Differences in H c . BD-FT denotes thedifferences of the two theories, while the others are the dif-ferences to the linear driving. 256 symbols the results aboutthe 256 lattices. (c) and (d) Finite time effects of κ and H c ,respectively. Each curve is obtained by successively omittingthe datum with the smallest R and plotting the results atthe remaining smallest R . Different curves start with dif-ferent largest R . The widths of the distributions have notbeen included into the fits, since their inclusion only slightlychange the results of large R for large ranges. For clarity, weplot only every other curve for the theories. Lines connectingthe symbols are only a guide to the eye. lowest rate, we find B ( T ) = 69 .
73, larger by more thantwo thousand times. On the other hand, BD yields goodresults even remarkably in the SD regime without anyadjustable parameters, though they are slightly smalleras seen in Fig. 1(b) and the H range is far larger than0 .
01 to 0 .
13 studied in Refs. [27, 28] for a constant field.Even though Fig. 1(a) appears to demonstrate both FTand BD are quite good generally, comparing with othercurves in Fig. 1(b), one sees that both theories exhibitsystematic deviations from the numerical results. Thiscan be clearly seen from Figs. 1(c) and (d), where weshow the results of a systematic fits to the simple powerlaw [40, 41], H c = H c + aR − κ , with constants H c , a ,and κ . For the theories, both κ and H c change contin-uously with the range of R that is used to find them,even if we change θ and K to give better agreement withthe numerical results, conforming to the expectation thatthe results described by such theories exhibit no scal-ing [24, 26]. However, the simulation results are quali- tatively distinct. If we include the theoretical data fromthe SD regime into the fits, we see a similar upturn near R = 10 and a descent at larger R . This would indicatethat the feature of the simulation results were related toa crossover from the MD to the SD regimes. However,deviations from the theoretical upturn are large (see Sup-plemental material for details). If we neglect in Figs. 1(c)and (d) the two rightmost data, we see monotonic vari-ations roughly up to the 12th curve (light cyan). Thisimplies that the theories might be valid within the rangefrom R = 0 . H c between the theories and the numerical results even inthe reduced range, though κ may agree. Note that thislarge gap cannot be removed by adjusting parameters like B ( T ), because bigger H c leads to bigger κ and thus thegap transfers to κ . Moreover, such possible adjustmentshave only a negligible effect since the differences in H c between the theories and the numerical results are small. Evidence for the RG theory
We next show that the φ theory can explain the results. Within the theory,scaling exists similar to the critical phenomena. For ex-ample, the scaling form for M is [29, 30], M ln m t = M s + R β/rν ( − ln R ) m f [( H ln n t − H s ) R − βδ/rν ( − ln R ) n ],where β , δ , ν , and r = z + βδ/ν are instability exponentsfor M , H , the correlation, and R , respectively, with z being the dynamic exponent, each corresponding to itscritical counterpart [29, 30], and f is a scaling function.When n = m = 0, H s and M s compose simply the insta-bility point around which the theory is expanded and arethus finite, in sharp contrast with the critical phenomena.In the presence of the special logarithmic corrections in t ,the point appears effectively at H s ln − n t and M s ln − m t ,which are scale dependence in consistent with previousstudies [48–50]. The ln n t term with n = d/ ( d −
1) wasargued to arise from the interplay between the exponen-tial time in tunneling between the two phases and dropletformations in the low- T phase in the Potts model [33].In that case, the field is replaced by T − T . The curvesof normalized energies versus ( T − T ) ln t for variouscooling rates cross at a finite value, which was suggestedto show a dynamic transition with spinodal-like singular-ity [33]. Figure 2(a) shows that this crossing does appearfor the Ising model studied here at T = 0 . T c . However,it is absent at T ≈ . T c . This indicates that the mech-anism can not be dominated generally, as varying T andvarying H cannot change the mechanism. We thus regard n as an adjustable parameter and introduce generally theother exponents for the logarithmic corrections.Our task is to show that the scaling form can indeedaccount for the data. This demands that there exist a -4 -3 -4 -3 -2 -0.90-0.85-0.12-0.0810 -4 -3 -2 -4 -4 -4 -3 -0.38-0.370.20.40.60.810 -4 -3 (h) H s (e) / r (i) M s R (f) / r R (g) H s R / r (d) (k) / r (n) H s R (l) / r (o) M s R (j) / r (m) H s R -1 0 1 0.0000.0010.002 c (c) ( < M >< M > ) R . ( < M >< M > ) R . (Hln -1/3 t 0.0915)R -0.62 ln -1.5 R -20 -10 0 100.00000.00020.00040.00060 10 20-101 ~0.6T c (a) M Hln t c -1 0 1 02040 ~0.6T c c (b) ( M l n - / t + . ) R . l n R (Hln -1/3 t 0.0915)R -0.62 ln -1.5 R -20 -10 0 100204060 (Hln t 6.93)R -0.60 ln -2 R ( M + . ) R . l n R (Hln t 6.93)R -0.60 ln -2 R FIG. 2. (Color online) (a) M versus H ln t for nine R about from 0 . . T = 0 . T c and from 0 . . T ≈ . T c . (b) Rescaled of those curves in (a). (c) Rescaled of h M i − h M i . Note that only the rising parts of the curves are expected to collapse after being rescaled in line with (b). In (b)and (c), the arrows indicate the bottom-left and top-right axes used for the two T . (d) to (o) Finite time effects of βδ/rν and H s and β/rν and M s fitted out of Eqs. (4) and (5), respectively, for the four n and m given in (b). Different from Fig. 1, foreach curve, starting from the rightmost data point that represents the fit to the R it stands and five others which are largerthan it, each connected successive point denotes the fit of its R and all the foregoing ones. The panels on the middle columnzoom in on the corresponding panels on the left. (d) to (i) are results of T = 0 . T c and (j) to (o) are the corresponding ones of T ≈ . T c . (d) and (g) [(j) and (m)] along with their respective enlarged ones (e) and (h) [(k) and (n)] are results of the fits ofthe H at M s = − .
883 [ − .
37] for various R at T = 0 . T c [ T ≈ . T c ] and (f) and (i) [(l) and (o)] are those of the fits of the M at H s = 6 .
93 [0 . H s is roughly 6 .
93 [0 . M s is about − .
883 [ − .
37] for thecorresponding data. They are thus self-consistent in that at M s the curves produce H s and just at this H s they give back to M s correctly. Errorbars are not shown in (f), (i), (l), and (o) since they are relatively large possibly due to the negative powerof the logarithms, though the fits are good for a not-large R range. Also the fits in these four panels appear not so approachingone another or level off as the others show, possibly because sub-leading contributions and corrections to scaling are strongerfor M . Lines connecting symbols are only a guide to the eye. In (a) to (c), the data points are dense and only lines connectingthem are displayed. single point, ( H s , M s ), such that at the particular M s H ln n t = H s + a R βδ/ ( rν ) ( − ln R ) − n , (4)while at the corresponding H s , M ln m t = M s + f (0) R β/ ( rν ) ( − ln R ) m , (5) self-consistently, where a is a constant satisfying f ( a ) =0. In order to reduce the parameters to be fitted and liftprecision, we choose the values of the four n and m asinput. We find this condition is highly restrictive fortheir choices. For example, if all n and m are set tozero, the condition cannot be satisfied. Neither can theseemingly plateau in Fig. 2(o). In addition, since we havenot considered sub-leading contributions and correctionsto scaling, Eqs. (4) and (5) are not expected to hold fora large range of R . Nevertheless, we require that theexponents obtained should somehow not depend on R ina certain range.Figure 2(d) to (o) show the results. Except for (f) and(i), all other figures show that the fitted results exhibitjumps from large to small R values. It is remarkable thatwhen the self-consistent H s and M s are reached, the fit-ted results minimize their variations with R and approacheach other for some R ranges. For example, at other M s ,the three lowest curves in Figs. 2(k) and (n) tilt and sep-arate from each other. For T ≈ . T c , n = m = − / − . − .
45, with βδ/rν and β/rν varying from 0 .
589 to 0 . − .
077 to − . M and its fluctuation h M i − h M i . The latter is rescaled just by R ( βδ − β ) / ( rν ) rather than follows the susceptibility ∂M/∂H , thoughthe exponents for the two functions are identical. Thisarises from the violation of fluctuation-dissipation theo-rem in the nonequilibrium driving [37]. The collapses asdisplayed in Figs. 2(b) and (c) are reasonably quite good,noting that only the leading behavior is considered, thusconfirming the results. Note however that data collapsesare sometimes deceptive. We show in Supplemental ma-terial an example in which the collapse appears perfectbut unreasonable.Besides the existence of the single finite H s and M s , themost striking result is that the estimated exponents andtheir deviations from results of both T , βδ/rν ≈ . β/rν ≈ − . .
575 and − . β in two dimensions [30].Moreover, although why the two T data take on n and m values of opposite signs and their consequences haveyet to be explored, a possible reason being the proximityof the high T to T c , the scaling functions appear to beuniversal up to a proper overall displacement and scalingas seen in Figs. 2(b) and (c). These therefore provide acompelling evidence for the RG theory.I thank Professor Per Arne Ridvold for his informationand Shuai Yin, Baoquan Feng, Yantao Li, Guangyao Li,and Ning Liang for their useful discussions. This workwas supported by National Natural Science Foundationof China (Grant Nos. 10625420 and 11575297). [1] J. D. Gunton, M. San Miguel, and P. S. Sahni, in PhaseTransitions and Critical Phenomena , eds. C. Domb andJ. L. Lebowitz Vol. 8, 267 (Academic, London, 1983).[2] A. J. Bray, Adv. Phys. , 357 (1994).[3] K. Binder and P. Fratzl, in Phase Transformations in Materials , ed. G. Kostorz, 409 (Wiley, Weinheim, 2001).[4] M. Volmer and A. Weber, Z. Phys. Chem. (Leipzig) ,277 (1926).[5] L. Farkas, ibid. , 236 (1927).[6] R. Becker and W. D¨oring, Ann. Phys. (Leipzig) , 719(1935).[7] Ya. B. Zeldovich, Acta Physicochim. USSR , 1 (1943).251 (1990).[8] F. F. Abraham, Homogeneous Nucleation Theory (Aca-demic, New York, 1974).[9] P. Debenedetti,
Metastable Liquids (Princeton UniversityPress, Princeton, NJ 1996).[10] D. Kashchiev,
Nucleation: Basic Theory with Applica-tions (Butterworth-Heinemann, Oxford, 2000).[11] D. W. Oxtoby, Acc. Chem. Res. , 91 (1998).[12] J. D. Gunton, J. Stat. Phys. , 903 (1999).[13] S. Auer and D. Frenkel, Annu. Rev. Phys. Chem. , 333(2004).[14] R. P. Sear, J. Phys.: Condens. Matter , 033101 (2007).[15] A. N. Kolmogorov, Bull. Acad. Sci. USSR, Class Sci.,Math. Nat. , 355 (1937)[16] W. A. Johnson and P. A. Mehl, Trans. Am. Inst. Min.Metall. Eng. , 416 (1939).[17] M. Avrami, J. Chem. Phys. , 1103 (1939).[18] S. Auer and D. Frenkel, Nature (London) , 1020(2001).[19] T. Kawasaki and H. Tanaka, Proc. Nat. Acad. Sci. ,14036 (2010).[20] L. Filion, R. Ni, D. Frenkel, and M. Dijkstra, J. Chem.Phys. , 134901 (2011).[21] J. S. Langer, Ann. Phys. (N. Y.) , 108 (1967).[22] N. J. G¨unther, D. J. Wallace, and D. A. Nicole, J. Phys.A , 1755 (1980).[23] P. A. Rikvold, H. Tomita, S. Miyashita, and S. W. Sides,Phys. Rev. E , 5080 (1994).[24] S. W. Sides, P. A. Rikvold, and M. A. Novotny, Phys.Rev. E , 6512 (1998).[25] R. A. Ramos, P. A. Rikvold, and M. A. Novotny, Phys.Rev. B , 9053 (1999).[26] S. W. Sides, P. A. Rikvold, and M. A. Novotny, Phys.Rev. E , 2710 (1999).[27] S. Ryu and W. Cai, Phys. Rev. E , 030601(R) (2010).[28] S. Ryu and W. Cai, Phys. Rev. E , 011603 (2010).[29] F. Zhong and Q. Z. Chen, Phys. Rev. Lett. , 175701(2005).[30] F. Zhong, Front. Phys. , 126402 (2017)[arXiv1205.1400 (2012)].[31] F. Zhong, Phys. Rev. E , 022104 (2012).[32] P. B. Thomas and D. Dhar, J. Phys. A , 3973 (1993).[33] A. Pelissetto and E. Vicari, Phys. Rev. Lett. , 030602(2017).[34] N. Liang and F. Zhong, Phys. Rev. E , 032124 (2017).[35] S. Gong, F. Zhong, X. Huang, and S. Fan, New J.Phys. , 043036 (2010).[36] F. Zhong, in Applications of Monte CarloMethod in Science and Engineering , 144103(2016).[38] M. Rao, H. R. Krishnamurthy, and R. Pandit, Phys. Rev. B , 856 (1990).[39] For a review, see K. Chakrabarti and M. Acharyya, Rev.Mod. Phys. , 847 (1998).[40] F. Zhong, J. X. Zhang, and X. Liu, Phys. Rev. E ,1399- (1995).[41] F. Zhong, Phys. Rev. B , 060401(R) (2002).[42] C. K. Harris, J. Phys. A , L143 (1984).[43] I. Lifshitz, Sov. Phys. JETP , 939-942 (1962) [Zh.Eksp. Teor. Fiz. , 1354 (1962)].[44] S. Allen and J. Cahn, Acta Metall. , 1084 (1979).[45] V. A. Shneidman, K. A. Jackson, and K. M. Beatty, J.Chem. Phys. , 6932 (1999). [46] X. Huang, S. Gong, F. Zhong, and S. Fan, Phys. Rev. E , 041139 (2010).[47] H. Tomita and S. Miyashita, Phys. Rev. B , 8886(1992).[48] C. Billoted and K. Binder, Z. Phys. B , 195 (1979).[49] K. Kawasaki, T. Imaeda, and J. D. Gunton, in Perspec-tives in Statistical Physics , ed. H. J. Ravech´e, 201 (NorthHolland, Amsterdam, 1981).[50] K. Kaski, K. Binder, and J. D. Gunton, Phys. Rev. B29