Order-disorder transition in a two-dimensional associating lattice gas
A. P. Furlan, Tiago J. Oliveira, Jürgen F. Stilck, Ronald Dickman
OOrder-disorder transition in a two-dimensional associating lattice gas
A. P. Furlan, a) Tiago J. Oliveira, J¨urgen F. Stilck, and Ronald Dickman Departamento de F´ısica, ICEx, Universidade Federal de Minas Gerais, C. P. 702, 30123-970 Belo Horizonte,Minas Gerais - Brazil Departamento de F´ısica, Universidade Federal de Vi¸cosa, 36570-900, Vi¸cosa,Minas Gerais - Brazil Instituto de F´ısica and National Institute of Science and Technology for Complex Systems,Universidade Federal Fluminense, Niter´oi, Rio de Janeiro - Brazil Departamento de F´ısica and National Institute of Science and Technology for Complex Systems,ICEx, Universidade Federal de Minas Gerais, C. P. 702, 30123-970 Belo Horizonte,Minas Gerais - Brazil (Dated: 21 March 2019)
We study an associating lattice gas (ALG) using Monte Carlo simulation on the triangular lattice and semi-analytical solutions on Husimi lattices. In this model, the molecules have an orientational degree of free-dom and the interactions depend on the relative orientations of nearest-neighbor molecules, mimicking theformation of hydrogen bonds. We focus on the transition between the high-density liquid (HDL) phaseand the isotropic gas phase in the limit of full occupancy, corresponding to chemical potential µ → ∞ ,which has not yet been studied systematically. Simulation results show a continuous phase transition at τ c = k B T c /γ = 0 . − γ is the bond energy) between the low-temperature HDL phase, witha non-vanishing mean orientation of the molecules, and the high-temperature isotropic phase. Results forcritical exponents and the Binder cumulant indicate that the transition belongs to the three-state Potts modeluniversality class, even though the ALG Hamiltonian does not have the full permutation symmetry of thePotts model. In contrast with simulation, the Husimi lattice results furnish a discontinuous phase transition,characterized by a discontinuity of the nematic order parameter. The transition temperatures ( τ c = 0 . . µ , three possible scenarios arise. The first is that in the limit µ → ∞ the first-orderline ends in a critical point; the second is a change in the nature of the transition at some finite chemicalpotential; the third is that the entire line is one of continuous phase transitions. Results from other ALGmodels and the fact that mean-field approximations show a discontinuous phase transition for the three-statePotts model (known to possess a continuous transition) lends some weight to the third alternative.PACS numbers: 68.35.RhKeywords: Phase transitions, critical phenomena, Monte Carlo simulations, lattice gas, Mean-field approxi-mations I. INTRODUCTION
It is no news that water exhibits quite unusual ther-modynamic behavior, characterized by a set of anoma-lies, among which the anomaly in density is the bestknown . In recent decades, many models were developedwith the purpose of investigating the fundamental mech-anisms which lead to the water anomalies. Among them,lattice models have attracted much attention due to theireasy implementation and low computational cost. Thesemodels usually include soft-core potentials to take careof excluded volume effects, and orientational interactionsto represent hydrogen bonding between molecules. Sofar, these models are only able to reproduce some of theanomalies qualitatively. In spite of this, they exhibit richphase diagrams that provide an ideal environment for thestudy of phase transitions as well as the validation of newcomputational techniques. a) Electronic mail: apfurlan@fisica.ufmg.br
The first orientational lattice model for water was pro-posed by Bell and Lavis (BL) . It is defined on a tri-angular lattice, in which each site can be either vacantor occupied by a molecule. The molecules possess threebonding directions with 120 ◦ between them, resulting intwo orientational states per molecule. The model ex-hibits three phases: gas, low density liquid (LDL) andhigh-density liquid (HDL) at low, intermediate and highchemical potentials, respectively. While the gas-LDLtransition is known to be discontinuous, characterizedby a jump in density, the nature of the liquid-liquidtransition is still controversial. For instance, mean-fieldapproximations (Bethe lattice solutions and cluster-variation methods ) point to a discontinuous phase tran-sition, whereas Monte Carlo simulations show a con-tinuous transition. There is no consensus regarding theuniversality class of this transition. Fiore et al. assertthat it falls the Ising universality class, while ˇSim˙enas etal. report three-state Potts critical exponents.Following the ideas of Bell and Lavis, Henriques andBarbosa introduced a two-dimensional (2D) associating a r X i v : . [ c ond - m a t . s t a t - m ec h ] M a r lattice gas (ALG) . In this model, also defined on a tri-angular lattice, each molecule has four bonding arms andtwo inert ones, the latter taking opposite directions onthe lattice. Two of the bonding arms are proton donorsin hydrogen bonds, while two are receptors, leading toeighteen states per molecule. The model also exhibitsgas, LDL and HDL phases, but in contrast to simulationresults for the BL model, the LDL and HDL phases areseparated by a discontinuous transition line that ends ata bicritical point. There is also a gas-HDL line of contin-uous transitions, which starts at the bicritical point andseems to extend to large values of the chemical potential.The LDL and gas phases are separated by a continuousand a discontinuous transition line, which connect at atricritical point. The former line meets the LDL-HDLand gas-HDL lines at the bicritical point .Variants of the ALG model have been investigated,among them, three-dimensional (3D) and symmetricversions . The 3D model is defined on a body-centered cubic lattice in which each molecule possessesfour bonding, and four inert arms. Despite the differencesin geometry and number of orientational states in rela-tion to the original ALG model , the three-phase (gas,LDL, HDL) behavior is preserved. While a phase dia-gram featuring two tricritical points was suggested byBuzanno et al. based on a cluster-variation approach,simulation results by Szortyka et al. indicate that thereis in fact a tricritical and a bicritical point, similar to the2D case . In the 3D model, however, a gas-LDL coexis-tence line meets the continuous LDL-HDL and gas-HDLtransition lines at the bicritical point .The symmetric ALG model makes no distinction be-tween donor and receptor bonding arms. This leads toa simplification, since the number of states per parti-cle is substantially reduced. Balladares et al. inves-tigated this model on the triangular lattice and foundonly two discontinuous (gas-LDL and a LDL-HDL) tran-sition lines in the phase diagram, each ending at a criti-cal point. As an aside, let us remark that this very samescenario was reported in the former studies of the orig-inal 2D ALG model and of the 3D version . It turnsout that they were incorrect, as demonstrated in morerecent analyses , as discussed above. In fact, thesemi-analytical solution of the symmetric model on aHusimi lattice build with hexagons (which is a mean-fieldapproximation for the triangular lattice) unveils a phasediagram with three coexistence lines (gas-LDL, gas-HDLand LDL-HDL) meeting at a triple point . Moreover,more recent simulations of this model have provided ev-idence that the critical points reported in are actuallytricritical points , so that the thermodynamic behaviorof this model is closer to the other versions. In contrastwith these cases and with the mean-field results, how-ever, a gas-HDL transition has not yet been observed insimulations of the 2D symmetric ALG, so that the ex-istence of a continuous order-disorder transition and itsuniversality class remains unclear for this model. In fact,studies of the critical exponents at the continuous phase transitions of all versions of the ALG model are still miss-ing.Motivated by this issue, we study the phase transi-tion in the full lattice limit of the symmetric 2D ALGmodel , using MC simulations (employing both Wang-Landau and Metropolis algorithms), as well as ob-taining the thermodynamic properties of the model inthe core of Husimi cacti . A fully occupied lattice cor-responds to the limit of infinite chemical potential, forwhich only the HDL and an isotropic, disordered phase(corresponding to what is identified as the gas phase inprevious studies of the AGL) are expected, since the LDLphase becomes metastable already for small values of thechemical potential . Indeed, we find a single tran-sition between the isotropic and the ordered HDL phase,whose loci, nature and universality class will be addressedin detail in the following.The remainder of this paper is organized as follows.In section II we detail the model. In Sec. III we explainthe simulation and analytic methods used. In Secs. IVand V we report our simulation results and Husimi-cactifindings, respectively. Finally in Sec. VI, we discuss ourconclusions and perspectives for future work. II. MODEL
We consider the associating lattice gas (ALG) intro-duced in Ref. in its symmetric version . The model wasproposed in the context of water-like anomalies; despiteits simplicity, it captures some features of liquid water,such as density and diffusion anomalies. The modelis defined on a triangular lattice (coordination number z = 6) in which each site can be either empty or occu-pied by a molecule. Each molecule has six arms, four ofwhich are bonding arms, while the other two are non-bonding (“inert”). In the symmetric version , all thebonding arms interact in the same manner, there beingno distinction between proton donors and receptors. Theinert arms do not interact and assume diametrically op-posed positions, giving rise to three orientational states, η i with i = 1 , and . We denote thegenerators of the triangular lattice by ˆ e , ˆ e and ˆ e , with,ˆ e = i , (1)ˆ e = + 12 i + √ j , (2)and ˆ e = − i + √ j . (3)We use the same set of vectors to label the orientationalstates. Consider, for example, state 1, with bonding armsalong ± ˆ e and ± ˆ e . As illustrated in Fig. 1, we associatethe vector η = ˆ e with state 1, η = ˆ e with state 2 and η = ˆ e with state 3. (Thus η points along one of the nonbonding directions.)In the ALG, interactions are restricted to nearest-neighbor (NN) pairs, so that the separations r i − r j be-tween interacting pairs again fall in the set {± ˆ e j } . Aparticle at site k in state η i has no interaction with itsneighbors at sites r k ± ˆ e i since it has no bonding armspointing toward these sites. On the other hand, for i (cid:54) = j we have | ˆ e i · ˆ e j | = 1 /
2. Thus the interaction between apair of particles at sites i and j , separated by r (a unitvector in the set {± ˆ e i , i = 1 , , } ), can be written so: u ij = − γ (cid:18) (cid:19) (cid:104) − ( η i · r ) (cid:105) (cid:104) − ( η j · r ) (cid:105) , (4)where − γ denotes the energy associated with a NN par-ticle pair with bonding arms pointed toward each other.Using this, the energy of the ALG model may be written: H ( η , r ) = (cid:88) (cid:104) i,j (cid:105) σ i σ j [ ε + u ij ] , (5)where the summation (cid:104) i, j (cid:105) runs over pairs of nearestneighbors, the σ i = 0 , ε > σ i = 1 , ∀ i ), so that the parameter ε is mean-ingless and γ is the only remaining parameter in theHamiltonian. Therefore, from here on we measure en-ergy in units of γ , and define a dimensionless tempera-ture τ = k B T /γ . We remark that when comparing ourresults with those in the literature for the full model ,one should keep in mind that in these papers the particu-lar case γ/ε = 2 was considered, and the parameter ε wasused to define the reduced temperature and chemical po-tential. Thus, in these studies the reduced variables aretwice those used here.At full occupancy, each particle interacts with four ofits six neighbors if all particles have the same orienta-tional state (see Fig. 1). Any configuration such that oneor more pairs of nearest neighbors have distinct statesmust have a higher energy. Thus the ground state isthreefold degenerate, with energy per particle (in unitsof γ ), e = −
2. By contrast, the mean energy per parti-cle in a configuration in which each orientation is chosenat random, independently, is e random = − /
3. Since thegain in entropy per site is ∆ s/k B = ln 3, a crude estimatefor the critical temperature is τ c ≈ / (3 ln 3) (cid:39) . µ →∞ of the full model, and at low temperatures it will bein the HDL phase. Varying the temperature, we expectto observe a transition from the low-temperature orderedphase, with a majority of particles in one of the orien-tational states, to a high-temperature disordered one,with equal populations among the three states. Thisisotropic phase corresponds to the gas phase in the gen-eral model. As noted in the Introduction, continuous gas-HDL phase transitions have already been observed insimulations of the nonsymmetric two-dimensional andthree-dimensional ALG models, while a discontin-uous transition was found in the mean-field approachfor the symmetric model considered here . In the 2Dcase , the transition is observed for reduced chemicalpotentials ¯ µ ≡ µ/γ (cid:38)
1. Szortika et al. proposed anorder parameter θ = 32 (cid:18) max ( n , n , n ) N − (cid:19) (6)where n i is the number of particles in state i and N isthe total number of particles. Evidently, θ = 1 when allparticles are in the same state, and θ = 0 for equallypopulated states.In view of the orientation-dependent interaction andthe orientationally ordered ground state, we can inter-pret the transition as one between nematic and isotropicphases. Analysis of such a transition suggests an alter-native definition of the order parameter, Q = 1 N (cid:113) n + n + n − n n − n n − n n . (7)Of central interest is the nature of the order-disordertransition in the model at full occupancy. (In principle,we would expect this also to describe the transition atlarge, but finite µ .) Since the ground state is threefolddegenerate, it is tempting to suppose that the transitionfalls in the three-state Potts class. While this may in factbe the case, we note that the energy, Eq. (5), does nothave the symmetry of the Potts model, that is, it is notinvariant under permutations among the states. (Theground state is permutation invariant, trivially, but dueto the orientational dependence of the interactions, anarbitrary configuration is not.)In light of the above considerations, we regard it as anopen question whether the model belongs to the 3-statePotts universality class. The simulation results reportedbelow provide some insight on this issue. III. METHODSA. Simulation details
We performed extensive Monte Carlo (MC) simula-tions of the symmetric ALG model using both Wang-Landau (WL) simulations for triangular lattices withlateral sizes 24 ≤ L ≤
72, and Metropolis simulations for larger systems. The simulations are performed atfull occupancy ( N = L ). The WL algorithm is an en-tropic sampling method, designed to estimate Ω ( E ), thenumber of configurations with energy E . Starting fromΩ ( E ) = 1 , ∀ E , the estimates for Ω are gradually refinedin a series of iterations. Each iteration generates a se-quence of configurations as described below. An energy +ˆ e η +ˆ e η +ˆ e η FIG. 1. (Color online) Definition of orientational states η ≡ ˆ e , η ≡ ˆ e and η ≡ ˆ e respectively. The thick black linesrepresent bonding arms, dashed gray lines represent directions on the triangular lattice and the arrows indicate the generatorsof triangular lattice.FIG. 2. A ground-state configuration of the fully occupiedlattice; all molecules are in the same state. histogram H [ E ], records the number of configurationshaving energy E . These simulations use a random initialconfiguration Γ = { η i } with energy E (Γ). To generatea candidate for the next configuration in the sequence(Γ new ), a site i is chosen at random and its state is al-tered ( η i → η (cid:48) i ), generating configuration Γ (cid:48) . This isaccepted as Γ new with probability p [Γ → Γ (cid:48) ] = min (cid:20) , Ω [ E (Γ (cid:48) )]Ω [ E (Γ)] (cid:21) . (8)With complementary probability, the current configu-ration Γ is taken as the new configuration. The his-togram and logarithm of the number of states updated so: H [ E (Γ new )] → H [ E (Γ new )] + 1 and ln Ω [ E (Γ new )] → ln Ω [ E (Γ new )] + ln f . The parameter f , known as themodification factor, is set to e , the base of natural log-arithms, on the first iteration. A sequence of configura-tions is generated by repeating this procedure, until thehistogram is flat , i.e., there is no E such that H ( E ) isless than (or greater than) X % of the mean value of thehistogram over all energies. Following the usual practice,we use X = 20%, that is, the 80% flatness criterion.Once the flatness criterion is satisfied, the current iter-ation ends and a new one is started with the histogram set to zero, but with the Ω [ E ] carried forward from the pre-ceding iteration. The new iteration proceeds as before ex-cept for a smaller modification factor, taken as the squareroot of the previous value (i.e., ln f → ln f / f < − orequivalently after 27 iterations.As is evident from Eq. (8), WL sampling is a kindof Metropolis importance sampling with target probabil-ity distribution P (Γ) ∝ / Ω [ E (Γ)] in place of the usualBoltzmann distribution, P (Γ) ∝ exp[ − βE (Γ)]. Thus ifiterations extended arbitrarily long and there were nosampling noise, the estimates for the Ω ( E ) would con-verge to their true values. For further details about theconvergence (or lack thereof) of WL sampling, see andreferences therein.Recent studies show that it is not necessary to useall 27 iterations of WL sampling. As the modificationfactor is barely greater than unity, estimates for the Ω ( E )are not significantly modified in the final few iterations.We used 27 iterations for system sizes L ≤
32, 22 for32 < L <
72, and 20 for L = 72.Using the estimates for the Ω ( E ), the canonical aver-age of a given thermodynamic observable O ( E ) at tem-perature T can be computed via, (cid:104)O(cid:105) ≡ Z (cid:88) E O ( E )Ω( E ) exp ( − βE ) , (9)where O ( E ) represents the microcanonical average, and Z is the canonical partition function. Microcanonicalaverages are computed as simple averages of property O over all configurations with energy E generated in thefinal iteration of the WL procedure. We compute thecanonical average of the energy E and its variance, andof order parameters θ and Q and their second throughfourth moments.Although WL sampling yields useful results for sys-tems with L ≤
90, it is not effective for larger systems.(The time to achieve a flat histogram becomes excessive.)We therefore use standard Metropolis sampling for sys-tem sizes L = 128 and L = 256. In these studies, we use10 MCS for equilibration followed by 2 × MCS to gen- (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)
FIG. 3. Six-coordinated Husimi tree built with triangles. Twogenerations of triangles are shown. The numbers indicate thethree directions of lattice edges, as defined in Fig. 2. Thehatched triangles should be removed to define a sub-tree. erate data. In both simulation methods, we perform anaverage over 60 independent realizations, starting fromrandomly generated initial configurations.
B. Mean-field approximations
Generally, mean-field approximations of a particularmodel correspond also to the solution of this model ona complete graph, with properly rescaled interactions .Since solutions of a model on the core of a Cayley tree(a graph with no loops, in which each node has the samenumber of neighbors) usually correspond to the Betheapproximation of this model on a regular lattice with thesame coordination number as the tree, Baxter suggestedthat these treelike lattices be called Bethe lattices . Ad-ditional correlations are taken into account considering Husimi trees , built with clusters (polygons or polyhe-drons) rather than single sites, so that short closed pathsare present. Analysis of a model on the core of a Husimitree yields its behavior on a
Husimi lattice .We study the ALG model on such treelike lattices, inthe limit in which each site is occupied by a molecule.Although we could begin with the six-coordinated Bethelattice, one readily verifies that the ground state of themodel on this lattice, without any closed paths, is highlydegenerate. The thermodynamic behavior of the modelon the Bethe lattice is qualitatively distinct from thatfound on the triangular lattice: there is no ordered phase.We therefore study the model on a Husimi tree built withtriangles, in which three triangles meet at each site, asshown in Fig. 3. We also study the model on a moreelaborate Husimi lattice built with hexagons composedof six elementary triangles (see Fig. 4). While the com-plete model, at finite chemical potential, was studied us- FIG. 4. Husimi tree built with hexagons composed by sixtriangles. The central hexagon and the six hexagons of thefirst generation are shown. ing similar techniques on this Husimi lattice , the limit(1 /µ →
0) of interest here was not analyzed in that work.As explained in Ref. , the use of hexagons in contrastto simple triangles as elementary clusters is imperativefor capturing the symmetries of the ground states of thethree phases of the full model (gas, LDL, and HDL) onthe Husimi lattice; in the case of full occupancy, how-ever, a tree of triangles (Fig. 3) already captures thesymmetry.As usual, to solve a model on treelike lattices we de-fine rooted sub-trees and their partial partition functions(ppf’s) for the possible configurations of the site at theroot. For the present model, there are three possible ori-entations of the molecule on the root site of a sub-tree,which may be represented by the variable η assuming thevalues 1 , . The details of the calculations arepresented in the appendix. IV. SIMULATION RESULTS
In this section we report simulation results for themodel at full occupation.The critical temperature τ c is obtained via finite sizescaling (FSS) analysis of the susceptibilities χ θ , χ θ ( τ ) = L τ (cid:104)(cid:10) θ (cid:11) − (cid:104) θ (cid:105) (cid:105) , (10)and χ Q , defined analogously, and of the specific heat c , c ( τ ) = β L (cid:104)(cid:10) E (cid:11) − (cid:104) E (cid:105) (cid:105) . (11)In Eqs. (10) and (11) the terms (cid:104) . . . (cid:105) represent canoni-cal averages, L the linear size of the system and τ thetemperature. (As mentioned in Sec. III, our estimatesfor thermal averages such as θ , χ and c are averages oversixty independent realizations). The expected finite-sizescaling forms for the order parameter, susceptibility, andspecific heat are, θ ≈ L − β θ /ν F (cid:16) tL /ν (cid:17) , (12)and similarly for Q , χ ≈ L γ/ν X (cid:16) tL /ν (cid:17) , (13)and c ≈ c + L α/ν C (cid:16) tL /ν (cid:17) , (14)where β , γ , and α are critical exponents as usually de-fined, and t ≡ ( τ − τ c ) /τ c is the reduced temperature.Additionally we determined Binder’s fourth-order cumu-lant of the order parameters , U ,θ ( T ) = 1 − (cid:10) θ (cid:11) (cid:104) θ (cid:105) . (15) U ,Q is defined analogously. At the critical point, U tends to a universal value, characteristic of the univer-sality class, system shape and boundary conditions em-ployed . A. Determining τ c The results for the order-parameters θ and Q , plottedin Fig. 5, suggest that a continuous phase transition oc-curs at some temperature in the interval 0 . < τ < . ρ i in themajority, minority and intermediate states, for systemsize L = 32, illustrating a continuous variation betweenground state (all particles in the same state) and nearlyequal populations. Relative uncertainties are plotted inthe upper inset. As expected, they are largest in the crit-ical region. In all cases, the relative uncertainty in θ isless than 3%; for system sizes other than L = 64, it is < .
40 0 .
45 0 .
50 0 . τ . . . . . . o r d e r p a r a m e t e r .
40 0 .
45 0 . τ . . . ρ ( η ) ρ ( η ) ρ ( η ) .
47 0 . τ . . σ FIG. 5. Order parameters θ and Q versus temperature forsystem sizes (upper to lower) L = 24 , , , , , , , θ and Q obtained via the WL algorithm. Black crosses repre-sent results for θ obtained via Metropolis sampling ( L = 128and 256); dashed lines represent polynomial fits to the data.Lower inset: particle fractions ρ in the majority (upper), mi-nority (lower) and intermediate states, for L = 32. Upperinset: relative uncertainties for θ in WL simulations; colorsfollow the main plot. The susceptibilities χ θ and χ Q are plotted in Fig. 6.Although these quantities exhibit similar behaviors, thereare slight differences in the critical region; the differencesare more evident for larger systems. The discrepanciesbetween the susceptibilities are related to their struc-tures. While Q (see Eq. 7) takes into account the densi-ties of all three states, θ only involves the majority den-sity (see Eq. 6). This may be why the uncertainties (in-set Fig. 6) in χ Q are approximately one-third those in χ θ .For the sizes shown in Fig. 6, the maximum relative un-certainty χ θ in is about 6%. The specific heat c , shown in .
47 0 .
48 0 .
49 0 . τ χ ( τ ) .
47 0 . τ σ FIG. 6. Susceptibilities versus temperature for system sizes L = 24 , , , , , ,
72 and 128. Solid lines and emptycircles represent, respectively, χ θ and χ Q obtained via the WLalgorithm. Black crosses represent results for χ θ obtained viaMetropolis simulations ( L = 128); dashed lines represent apolynomial fit to the data for L = 128. (For better visibility,the data for L = 256 are not shown.) Inset: uncertainties in χ θ (solid lines) and χ Q (dashed lines). Fig. 7, exhibits behavior consistent with that of the orderparameters and susceptibilities. Relative uncertainties in c are smaller than 4% for 48 ≤ L ≤
74 and smaller than2% for
L < .
47 0 .
48 0 .
49 0 . τ c ( τ ) .
46 0 .
48 0 . τ . . . σ FIG. 7. Specific heat versus temperature for system sizes L = 24 , , , , , , ,
128 and 256. Solid lines repre-sent results from WL simulations and black crosses results ob-tained via Metropolis sampling for L = 128 and 256. Dashedlines represent a polynomial fit to the simulation data. Inset:uncertainties in c . The Binder cumulants of the order parameters areshown in Fig. 8. The crossings of the Binder cumulantsfor θ and Q provide the estimates τ ( Q ) c = 0 . τ ( θ ) c = 0 . c , χ θ , χ Q , and the cumulants, .
45 0 .
47 0 .
49 0 . τ . . . . U ( τ ) .
475 0 . . . .
63 0 .
475 0 . . . .
63 0 .
993 0 . . . . . FIG. 8. Binder cumulants U ,θ (solid lines) and U ,Q (circles) versus temperature for system sizes L =24 , , , , , , ,
128 and 256. Crosses: results for U ,θ using Metropolis sampling for L = 128 and 256; dashed lines:sixth-order polynomial fits to the simulation data. Left in-set: detail of the crossing region of U ,θ ; upper-right inset:detail of the crossing region of U ,Q . Lower-right inset: cu-mulant of the three-state Potts model on the square lattice forsizes L = 32 , , , ,
64 and 128; the dashed line represents L = 128. we estimate the critical temperature, τ c . For c and the χ ’s, we define a size-dependent pseudocritical tempera-ture as the temperature associated with the maximumvalue. Pseudocritical temperatues for the cumulants areidentifying as the crossing temperatures of U between:1) the smallest system size studied, L = 24, and theothers ( L = 32 , . . . , L = 256) and 2) the crossing tem-peratures between a given system size L and the nextsystem size, for example, L = 24 with L = 32, L = 32with L = 40 and so on. The results from the adjacentsizes crossings are represented as U (cid:48) .The pseudocritical temperatures are plotted versus1 /L in Fig. 9. All six sets of pseudocritical temperaturesappear to converge to similar values as L → ∞ . Of noteis the relative insensitivity to system size of the pseud-ocritical temperatures derived from cumulant crossings.The crossings between Binder cumulant of two adjacentsystem sizes for the order parameter θ suffer from largeuncertainties for L ≥
64, affecting the estimate of thepseudocritical temperature. For this reason, we disregardthis property in the calculation of the critical tempera-ture.The resulting estimates for τ c are listed in Table I.From the six estimates for the pseudocritical tempera-tures we derive the global estimate τ c = 0 . τ c was obtained through a weightedaverage of the properties with weights 1 /σ , where σ rep-resents the uncertainty of each quantity. c χ θ χ Q U ,θ U ,Q U (cid:48) ,Q τ c T c obtained fromanalysis of the quantities listed on the first line, while the seventh column contains the global estimate and its uncertainty. .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 . /L . . . ˜ τ ( X ) c X = χ θ X = χ Q X = cX = U ,θ X = U ,Q X = U ,Q FIG. 9. Pseudocritical temperatures associated with χ θ , χ Q , c , U ,θ and U ,Q , with symbols as indicated, versus 1 /L .Broken lines are fits to the data (quadratic, except for the cu-mulant data, for which linear fits are very good). B. Critical exponents
Extrapolating the cumulant crossings to infinite sizeyields two estimates for the critical cumulant: U ∗ ,θ =0 . U ∗ ,Q = 0 . U ∗ = 0 . U ∗ Potts = 0 .
61. (In Ref. , no uncertainty estimateis provided; we may assume it is on the order of 0.01.)In efforts to improve on this estimate we simulated the3-state Potts model on L × L square lattices with peri-odic boundaries. The results for U , Potts , shown in Fig. 8(lower-right insert), yield U ∗ = 0 . β/ν and γ/ν using fits to the data for θ , Q , and the associated suscep-tibilities versus system size (see Fig. 10).From the linear fits shown in Fig. 10 (left) and in-cluding the effect of the uncertainty in τ c we obtain β θ /ν = 0 . β Q /ν = 0 . ln( L ) − . − . − . − . l n [ O ( τ c ) ] O = θ O = Q . . . ln( L ) − r e s i du a l × − ln( L ) O = χ θ O = χ Q . . . ln( L ) − r e s i du a l × − FIG. 10. (Left) Order parameters θ ( Q ) versus system size.The circles(stars) represent simulation data; dashed(dot-dashed) lines are least-squares linear fits of the data. Insetsshow the residuals of the fits. Graphs on the right are analo-gous plots for the associated susceptibilities. the two-dimensional, three-state Potts model is β Potts /ν Potts = 2 /
15 = 0 . . . . . In the worst case( β Q ), the discrepancy is about 4%. Linear fits to ln χ θ and ln χ Q as functions of ln L (see Fig. 10(right)) yield γ θ /ν = 1 . γ Q /ν = 1 . γ θ differs from the value of the three-state Potts model, γ Potts /ν Potts = 26 /
15 = 1 . . . . , by 1 . c as a function of ln L does not yield agood description of the data. We therefore fit the datawith c = aL α/ν + c with a , α/ν and c as adjustableparameters. This form is capable of fitting the simulationdata, as can be seen in Fig. 11. The inset shows that theresiduals are smaller than the uncertainties, and do notexhibit a systematic tendency. A least-squares procedureyields a = 2 . α/ν = 0 . c = − ±
1. Theexponent ratio agrees to within uncertainty with that ofthe three-state Potts model, γ/ν
Potts = 0 . α + 2 β + γ = 2, we have, ν = 2 α (cid:48) + 2 β (cid:48) + γ (cid:48) , (16)where exponents with primes denote the correspondingexponents divided by ν , determined in the finite-size scal-ing analysis discussed above.
20 40 60 80 100 120 L c ( τ c )
50 100 − . . . r e s i du a l FIG. 11. Specific heat versus system size at critical temper-ature. Squares: simulation data with their respective uncer-tainties; dashed line: fit, c = aL α/ν + c . The inset shows theresiduals, with error bars. Using the values for β (cid:48) and γ (cid:48) obtained using the orderparameter θ , equation 16 yields ν θ = 0 . Q furnish ν Q =0 . ν = 0 . ν Potts = 5 / . . . . .As a further test regarding the universality class, weperform a data collapse using the Potts critical expo-nents, as shown in Fig. 12. The quantities θ, Q, χ θ and χ Q exhibit a good collapse using these exponents. Weobserve that the specific heat exhibits the poorest col-lapse. − θ L β / ν − − Q L β / ν − − . . . . χ θ L − γ / ν × − − − . . . . χ Q L − γ / ν × − − − . . . . c L − α / ν [( τ − τ c ) /τ c ] L /ν FIG. 12. Data collapse for θ, Q, χ θ , χ Q and c respec-tively. The symbols (cid:13) , (cid:3) , × , (cid:52) , (cid:53) , (cid:63) and (cid:66) represent L =24 , , , , , ,
72 respectively. Filled black circles ( • )represent L = 128. The results for the exponents α, β , γ and ν , as well asfor the Binder cumulant strongly suggest that the phasetransition studied here belongs to the three-state Pottsmodel universality class. Deviations of the critical ex-ponents from the Potts class values may be due to the flatness criterion, limited sample size, and restricted sys-tem sizes. V. RESULTS ON HUSIMI LATTICES
The thermodynamic behavior of the model on the coreof both treelike lattices will be presented here, after abrief description of the calculations which lead to it.Some more details of these calculations may be foundin the appendix. As mentioned in Sec. III, recursionrelations for the partial partition functions on a rootedsub-tree are obtained, and ratios of these fpp’s usuallyremain finite as the recursion relations are iterated. Thestable (real and positive) fixed points attained iteratingthe recursion relations for the ratios correspond to thethermodynamic limit. The behavior on the complete treemay then be found connecting sub-trees to a central poly-gon. Expected values of densities at this central polygonmay then be calculated, and may be seen as approxima-tions to the results on the triangular lattice.Two stable fixed points are found. At higher temper-atures, only a disordered isotropic fixed point is stable,for which the densities of molecules in the three possi-ble orientations are equal. At low temperatures, threeordered nematic fixed points are stable, where the den-sity of molecules in one of the three orientations is largerthan the ones in the other two orientations. The stabil-ity of these fixed points may be studied using the Ja-cobian of the recursion relations for the ratios. It isfound that there is an interval of temperatures whereboth kinds of fixed point are stable, signalling that thetransition between the phases is discontinuous. Thus, tofind the coexistence temperature the free energy of bothphases should be equated. It should be remembered thatthis free energy corresponds to the model on the core ofthe tree; the free energy of the whole tree is dominatedby the surface . We find a coexistence temperature τ c = 0 . τ c = 0 . Q is shown in Fig. 13 asa function of the temperature for both treelike lattices.We note that the transition temperature is reduced, andthe jump in the order parameter becomes smaller as wemove from the triangle tree to that with hexagons. Also,the coexistence temperature decreases and the temper-ature interval in which both fixed points are stable be-comes narrower. These results are consistent, as the tran-sition temperature on the triangular lattices estimatedfrom the simulations is still lower and the transition iscontinuous.Finally, let us discuss how transition we study here fitsinto the general phase diagram of the model for finitechemical potential. The phase diagram of the completemodel on the hexagon-Husimi tree is shown in Fig. 14, inthe τ × ¯ µ = µ/γ plane. In this diagram we see that thegas-HDL coexistence line, as well as the gas and HDLspinodals, approach vertical asymptotes, and are very0 .
48 0 .
49 0 .
50 0 .
51 0 . τ . . . . . Q FIG. 13. Nematic order parameter Q as a function of τ closeto the coexistence temperature τ c . The upper curve corre-sponds to the triangle-Husimi lattice. The lower curve is forthe hexagon tree, calculating Q at the central site of the cen-tral hexagon, in the limit γ/µ → . . . . . . τ ¯ µ LDLHDL gasCoex. linesSpin. LDLSpin. HDLSpin. gas
FIG. 14. Phase diagram on the hexagon-Husimi tree, for finitechemical potential, in the ¯ µ × τ plane. The three phases (HDL,LDL, and gas) are separated by coexistence lines (full linesin graph), which meet at a triple point. Spinodal lines of thephases are also shown. The red circle represents the criticaltemperature obtained via MC simulations. close to them already for ¯ µ (cid:38)
10, whose τ values arevery close to the ones presented above for the solutionin the full-occupancy case. Hence, our results for theinfinite chemical potential limit correspond to the finalpoint of the coexistence line between the HDL and gasphases. For finite µ the coexisting phases differ in theorder parameter Q , and in the particle density ρ ; thedensity is lower and Q vanishes in the gas phase. In thefull-occupancy limit we consider here, Q still vanishes inthe gas phase and is finite in the HDL phase, but ρ = 1in both phases.Figure 15 shows the discontinuity of the nematic order ¯ µ . . . . ∆ Q . . . / ¯ µ . . . . ∆ Q FIG. 15. Discontinuity of the nematic order parameter at thegas-HDL coexistence line as a function of the reduced chem-ical potential ¯ µ . The insertion shows the same data against1 / ¯ µ . parameter at the gas-HDL coexistence line as a func-tion of the reduced chemical potential ¯ µ . We note that∆ Q becomes smaller as the chemical potential grows, andreaches its minimum value in the limit studied here. Wecan imagine three possible scenarios for this part of thephase diagram on the triangular lattice. One possibilityis that in the limit 1 /µ → µ , and thus the coexistence line would end at a tricriticalpoint. Finally, the whole transition may be continuous.In this case, if the other transition lines (LDL-HDL andgas-LDL), which meet with the gas-HDL transition lineat a triple point on the hexagon tree (see Fig. 14), re-main discontinuous, this point would become a criticalend point. VI. CONCLUSIONS
We investigate the gas-HDL phase transition of theALG model in its symmetric version in the limit1 /µ → τ c = 0 . U ∗ for the three-state Potts model on a square lattice: U ∗ = 0 . U ∗ = 0 . . The criti-cal exponents reported in Sec. IV also support a 3-statePotts model universality class.The mean-field calculations are performed on Husimi1trees constructed with triangles and hexagons. Both ap-proaches revel a discontinuous phase transition with co-existence temperatures somewhat higher than the simu-lation estimate, viz., τ c = 0 . τ c =0 . (cid:46) ,which may explain the difference found here for the ALGmodel.The complete (mean-field) phase diagram shows a gas-HDL coexistence line which extends from ¯ µ ≈ µ → ∞ , along which a reduction of the discontinuityof the nematic order parameter Q is observed as µ in-creases. The existence of the discontinuous line gives ori-gin to three possible scenarios that connect the results ofsimulations and mean-field. The first is that at the limit1 /µ → ∞ the line ends in a critical point, the secondis a change in the order of transition at a finite chemi-cal potential, and the third raises the possibility of thewhole transition be continuous for the triangular lattice.We remark that the ALG model with finite ¯ µ is some-what similar to a diluted three-state Potts system, forwhich a continuous transition is also observed . Thissuggests that the entire gas-HDL line might be contin-uous. Another indication of this is the fact that it iscontinuous in the original 2D ALG model , as well asin the 3D version . In future work, we intend to per-form Wang-Landau sampling, as well as transfer matrixcalculations in order to obtain the complete phase dia-gram and identify which of the three scenarios mentionedabove is correct. VII. ACKNOWLEDGMENTS
The authors acknowledge financial support from theCNPq, CAPES and Fapemig, Brazilian agencies. Wethank to M.C. Barbosa and W. Selke for helpful discus-sions and also Statistical physics Laboratory of Univer-sidade Federal de Minas Gerais for the computationalsupport.
Appendix: Solution of the model on Husimi trees
We show in some detail the calculations which lead tothe thermodynamic properties on a Husimi tree, partic-ularly the one built with triangles, as shown in Fig. 3.As described above, we start by defining partial partitionfunctions (ppf’s) for a sub-tree with a fixed configurationof the molecule on the root site. We thus define nine ppf’s g i,η of rooted sub-trees, where the first index assumes thevalues 1 , η of thenon-bonding arms of the molecule at the root. The recursion relations are obtained by considering theoperation of building a sub-tree with an additional gener-ation ( M + 1) by connecting two pairs of sub-trees (with M generations each) to the vertices which are oppositeto the root vertex of the new root triangle. We call j i theorientation of the molecule placed on the vertex oppositeto the edge of orientation i of the root triangle. The ppf’sof a sub-tree with M + 1 generations are thus obtainedthrough the recursion relations: g ( M +1)1 ,j = (cid:88) j =1 3 (cid:88) j =1 W { j ,j ,j } g ( M )1 ,j g ( M )1 ,j g ( M )2 ,j g ( M )3 ,j , (A.1a) g ( M +1)2 ,j = (cid:88) j =1 3 (cid:88) j =1 W { j ,j ,j } g ( M )2 ,j g ( M )2 ,j g ( M )3 ,j g ( M )1 ,j , (A.1b) g ( M +1)3 ,j = (cid:88) j =1 3 (cid:88) j =1 W { j ,j ,j } g ( M )3 ,j g ( M )3 ,j g ( M )1 ,j g ( M )2 ,j . (A.1c)The function W { j ,j ,j } is the statistical weight as-sociated to the edges of the root triangle, defined as W = exp[2 n b ( j , j , j ) /τ ], where n b is the number ofedges of the triangle with no non-bonding arms on themand τ is the temperature. Thus: n b ( j , j , j ) = n b, ( j , j ) + n b, ( j , j ) (A.2)+ n b, ( j , j ) = δ j , δ j , + δ j , δ j , + δ j , δ j , , where δ i,j is the Kronecker delta. In this expression n b,i stands for the number of edges in direction i with a bondon them, assuming the values 0 or 1, so that, as expected, n b ∈ [0 , M → ∞ ).Usually, the ppf’s diverge in this limit, so it is convenientto define ratios of ppf’s R i,j = g i,j g i, + g i, + g i, . (A.3)These ratios generally converge to finite values in thethermodynamic limit; six of them are independent, since (cid:80) j R i,j = 1. Therefore, the thermodynamic behavior ofthe model will be determined by the stable fixed pointsof the recursion relations for the ratios R i,j , which areobtained from Eqs. (A.1).We find fixed points with two different symmetries.In the isotropic fixed point the values of the 9 ratios,represented as a 3 × R ∗ iso = x (1 − x ) / − x ) / − x ) / x (1 − x ) / − x ) / − x ) / x , (A.4)where the parameter x is a function of the temperature τ ,whose exact expression can be easily found with the helpof an algebra software, but is too long to be presented2here. The three nematic fixed points are: R ∗ nem, = x (1 − x ) / − x ) / x x − x − x x − x − x x , (A.5) R ∗ nem, = x x − x − x (1 − x ) / x (1 − x ) / − x − x x x , (A.6) R ∗ nem, = x − x − x x − x − x x x (1 − x ) / − x ) / x . (A.7)Again the parameters x i , with i = 1 , ,
3, are functions ofthe temperature. We note that the isotropic fixed pointis stable at high enough temperatures, while at lowertemperatures one of the nematic fixed points is reached.We will see below that in the isotropic phase the ori-entation of the non-bonding arms of the molecules hasno preferred direction, while in the nematic phase one ofthe three directions is more probable than the other two.As expected for a discontinuous transition, the isotropicand nematic fixed points are both stable in an interval oftemperatures around the coexistence temperature. Theregion of stability of a given fixed point may be foundfrom the largest eigenvalue of the Jacobian of the recur-sion relations for the ratios: J i,j = (cid:18) ∂R (cid:48) i ∂R j (cid:19) R ∗ , (A.8)where R (cid:48) i and R i denotes the ratios in generations M + 1and M , respectively, and the derivative is calculatedat the fixed point whose stability is being considered.We find that the isotropic fixed point is stable for τ ≥ . τ ≤ . τ ∈ [0 . , . , we assumethat the contribution to the free energy per triangle isdifferent for triangles located at the surface of the tree( φ s ) and the ones in the bulk ( φ b ). In a tree with M generations of triangles, the number of triangles in thebulk is N b = 1 + 6 + 6 · ... + 6 · M − , while the numberof triangles on the surface is N s = 4 · M − . Therefore,if ϕ M is the free energy of a tree with M generations oftriangles, we find that ϕ M +1 − ϕ M = 3 φ b , and since ϕ M = − k B T ln Y M , we have that: φ b = − k B T ln Y M +1 Y M , (A.9)where Y M is the partition function on the whole tree,which can be obtained by connecting six sub-trees to the central triangle. This procedure leads to: Y M = (cid:88) j ,j ,j =1 W { j ,j ,j } g ( M )1 ,j g ( M )1 ,j g ( M )2 ,j g ( M )2 ,j g ( M )3 ,j g ( M )3 ,j . (A.10)In the thermodynamic limit M → ∞ , using the recur-sion relations [Eqs. (A.1)], the ratios [Eqs. (A.3)] at thefixed points, and this expression for the partition function[Eq. (A.10)], the bulk free energy per triangle may beobtained through Eq. (A.9). This yields the coexistencetemperature τ c = 0 . discontinuous transition between the isotropic and ne-matic phases in the solution of the model on this Husimilattice.To further confirm this, we calculate the order param-eter Q as defined in Eq. 7. From the partition function(A.10) it is simple to obtain the mean value of the num-ber of molecules with orientation i of non-bonding armsin the central triangle: (cid:104) n i (cid:105) ( M ) = 13 Y ( M ) 3 (cid:88) j ,j ,j =1 n i W { j ,j ,j } (A.11) × g ( M )1 ,j g ( M )1 ,j g ( M )2 ,j g ( M )2 ,j g ( M )3 ,j g ( M )3 ,j , where this number has been normalized so that 0 ≤(cid:104) n i (cid:105) ≤
1. If we divide the numerator and the denomina-tor by (cid:81) i =1 ( g ( M ) i, g ( M ) i, g ( M ) i, ) we may express this meanvalues in terms of the ratios, which in the thermodynamiclimit assume their fixed point values. The result is: (cid:104) n i (cid:105) = (cid:80) j ,j ,j =1 n i f ( j , j , j , { R ∗ } ) (cid:80) j ,j ,j =1 f ( j , j , j , { R ∗ } ) , (A.12)where: f ( j , j , j , { R ∗ } ) ≡ W { j ,j ,j } (A.13) × R ∗ ,j R ∗ ,j R ∗ ,j R ∗ ,j R ∗ ,j R ∗ ,j . As expected, in the isotropic gas phase, one has (cid:104) n (cid:105) = (cid:104) n (cid:105) = (cid:104) n (cid:105) , so that the order parameter is Q gas = 0.On the other hand, in the nematic HDL phase we find anon-vanishing Q HDL , whose variation with the tempera-ture near the coexistence is displayed in Fig. 13. At thecoexistence point the discontinuity in Q is ∆ Q = 0 . , when all lat-tice sites are occupied by molecules. For brevity, we willnot detail the calculations here, and present the main re-sults only. Since the tree in this case is built with largerclusters, hexagons with one central site and six sites atthe border, it is expected that such a calculation shouldlead to results which are closer to those on the triangularlattice. Again, we find a discontinuous transition, at asomewhat lower temperature, τ c = 0 . Q = 0 . τ ∈ [0 . , . / site or a mean value of the order parameter overall sites of the central hexagon . We have done both cal-culations and found out that these values are quite close.The deviation is maximum at the coexistence tempera-ture and the relative difference is about 0.25% there. M. Chaplin, “Water structure and science,” (2018). G. M. bell and D. A. Lavis, “Two-dimensional bonded latticefluids i. interstitial model,” J. Phys. A: Gen. Phys. , 427 (1970). G. M. Bell and D. A. Lavis, “Two-dimensional bonded latticefluids. ii. orientable molecule model,” J. Phys. A: Gen. Phys. ,568 (1970). D. A. Lavis, “The steam-water-ice system: a two-dimensionalbonded lattice model. the first-order approximation,” J. Phys.C: Solid State Phys. , 1530 (1973). C. E. Fiore, M. M. Szortyka, M. C. Barbosa, and V. B. Hen-riques, “Liquid polymorphism, order-disorder transitions andanomalous behavior: A monte carlo study of the bell-lavis modelfor water,” J. Chem. Phys. , 164506 (2009). M. A. A. Barbosa and V. B. Henriques, “Frustration and anoma-lous behavior in the Bell-lavis model of liquid water,” Phys. Rev.E , 051204 (2008). P. Bruscolini, A. Pelizzola, and L. Casetti, “Comment on “novelphase behavior in a two-dimensional network-forming latticefluid”,” Phys. Rev. Lett. , 089601 (2002). M. ˇSim˙enas, A. Ibenskas, and E. E. Tornau, “Phase transitionproperties of the Bell-lavis model,” Phys. Rev. E , 042124(2014). V. B. Henriques and M. C. Barbosa, “Liquid polymorphism anddensity anomaly in a lattice gas model,” Phys. Rev. E , 031504(2005). M. M. Szortyka, V. B. Henriques, M. Girardi, and M. C. Bar-bosa, “Dynamic transitions in a two dimensional associating lat-tice gas model,” J. Chem. Phys. , 184902 (2009). M. Girardi, A. L. Balladares, V. B. Henriques, and M. C. Bar-bosa, “Liquid polymorphism and density anomaly in a three-dimensional associating lattice gas,” J. Chem. Phys. , 064503(2007). M. M. Szortyka, M. Girardi, V. B. Henriques, and M. C. Bar-bosa, “Dynamic transitions in a three dimensional associatinglattice gas model,” J. Chem. Phys. , 134904 (2010). A. P. Furlan, N. G. Almarza, and M. C. Barbosa, “Lattice modelfor water-solute mixtures,” J. Chem. Phys. , 144501 (2016). A. L. Balladares, V. B. Henriques, and M. C. Barbosa, “Liq-uid polymorphism, density anomaly and h-bond disruption inassociating lattice gases,” J. Phys.: Condens. Matter , 116105(2007). A. P. Furlan, C. E. Fiore, and M. C. Barbosa, “Influence of disordered porous media on the anomalous properties of a simplewater model,” Phys. Rev. E , 032404 (2015). C. Buzzano, E. De Stefanis, and M. Pretti, “Cluster-variationapproximation for a network-forming lattice-fluid model,” J.Chem. Phys. , 024506 (2008). T. J. Oliveira, J. F. Stilck, and M. A. A. Barbosa, “Solutionof an associating lattice-gas model with density anomaly on ahusimi lattice,” Phys. Rev. E , 051131 (2010). F. Wang and D. P. Landau, “Efficient, multiple-range randomwalk algorithm to calculate the density of states,” Phys. Rev.Lett. , 2050–2053 (2001). F. Wang and D. P. Landau, “Determining the density of states forclassical statistical models: A random walk algorithm to producea flat histogram,” Phys. Rev. E , 056101 (2001). K. Husimi, “Note on mayers’ theory of cluster integrals,” J.Chem. Phys. , 682–684 (1950). M. Plischke and B. Bergersen,
Equilibrium Statistical Physics ,3rd ed. (WORLD SCIENTIFIC, 2006). N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H.Teller, and E. Teller, “Equation of state calculations by fastcomputing machines,” J Chem. Phys. , 1087–1092 (1953). R. Belardinelli, V. Pereyra, R. Dickman, and B. Louren¸co, “In-trinsic convergence properties of entropic sampling algorithms,”J. Stat. Mech. , P07007 (2014). A. A. Caparica, “Wang-landau sampling: A criterion for haltingthe simulations,” Phys. Rev. E , 043301 (2014). R. J. Baxter,
Exactly solved models in statistical mechanics (WORLD SCIENTIFIC, London, 1985). P. Gujrati, “Bethe or bethe-like lattice calculations are more reli-able than conventional mean-field calculations,” Phys. Rev. Lett. , 809 (1995). M. E. Fisher and M. N. Barber, “Scaling theory for finite-sizeeffects in the critical region,” Phys. Rev. Lett. , 1516–1519(1972). K. Binder, “Critical properties from monte carlo coarse grainingand renormalization,” Phys. Rev. Lett. , 693–696 (1981). W. Selke, “The critical binder cumulant for isotropic ising modelson square and triangular lattices,” J. Stat. Mech. , P04008(2007). T. Tom and A. Petri, “Cumulants of the three-state potts modeland of nonequilibrium models with c v symmetry,” J. Phys. A:Math. Gen. , 5379 (2002). F. Y. Wu, “The potts model,” Rev. Mod. Phys. , 235–268(1982). Nagai, Okamoto, and Janke, “Crossover scaling in the two-dimensional three-state potts model,” Condensed Matter Physics , 23605 (2013). L. Mittag and M. J. Stephen, “Mean-field theory of the manycomponent Potts model,” J. Phys. A , L109 (1974). B. Nienhuis, A. N. Berker, E. K. Riedel, and M. Schick, “First-and second-order phase transitions in potts models: Renormal-ization group solutions,” Phys. Rev. Lett. , 737 (1979). X. Qian, Y. Deng, and H. W. J. Bl¨ote, “Dilute Potts model intwo dimensions,” Phys. Rev. E72