Optimal modularity for nucleation in network-organized Ising model Hanshuang Chen and Zhonghuai Hou∗ Hefei National Laboratory for Physical Sciences at Microscales & Department of Chemical Physics, University of Science and Technology of China, Hefei, 230026, China (Dated: January 19, 2011) We studynucleation dynamicsof Ising model in a topology that consists of two coupled random networks,therebymimickingthemodularstructureobservedinreal-worldnetworks. Byintroducing 1 a variant of a recently developed forward flux sampling method, we efficiently calculate the rate 1 and elucidate the pathway for nucleation process. It is found that as the network modularity 0 becomesworsethenucleationundergoesatransitionfromtwo-steptoone-stepprocess. Interestingly, 2 the nucleation rate shows a nonmonotonic dependency on the modularity, in which a maximal nucleation rate occurs at a moderate level of modularity. A simple mean field analysis is proposed n toqualitatively illustrate thesimulation results. a J PACSnumbers: 89.75.Hc,64.60.Q-,05.50.+q 8 1 I. INTRODUCTION in three-dimensional lattice has also been studied using ] h transitionpathsamplingapproach[25]. Inaddition,Ising c model has been frequently used to test the validity of e Inthe lastdecade,criticalphenomenain complexnet- classical nucleation theory (CNT) [26–30]. However, all m works have received an enormous amount of attention these studies are limited to regular lattices in Euclidean in the field of statistical physics and many other disci- - space. Since many real systems can be properly mod- t plines (see, for example, a recent review [1]). Extensive a eledbynetwork-organizedstructure,itisthusnaturalto t research interests have focused on the onset of critical s ask how the topology of a networked system affects the behaviors in diverse network topology, which included a . nucleation process of Ising model. t wide range of issues: percolationphenomenon [2–5], epi- a m demic thresholds [6, 7], order-disordertransitions [8–12], In a recent work [31], we have studied nucleation dy- synchronization[13,14],self-organizedcriticality[15,16], namicsonscale-freenetworks,inwhichwefoundthatnu- - d nonequilibrium pattern formation [17], etc. However, cleation starts from, on average, nodes with more lower n thereismuchlessattentionpaidtothedynamics/kinetics degrees, and the rate for nucleation decreases exponen- o ofphasetransitionitselfincomplexnetwork,suchasnu- tially with network size and the size of critical nucleus c cleation in a first-order phase transition. increase linearly with network size, implying that nucle- [ Nucleation is a fluctuation-driven process that initi- ation is relevant only for a finite-size network. Herein, 1 we want to study nucleation dynamics of Ising model ates the decay of a metastable state into a more stable v in a modular network composed of two coupled random one [18]. A first-order phase transition usually involves 0 networks. It is known that many real-world networks, 3 the nucleation and growthof a new phase. Many impor- asdiverseasfromsocialnetworkstobiologicalnetworks, 4 tant phenomena in nature, including crystallization[19], havefoundtoexhibitmodularitystructures[32,33],that 3 glass formation [20], and protein folding [21], etc., are is,linkswithinmodulesaremuchmoredenserthanthose . associated with nucleation. Despite its apparent impor- 1 between modules. Many previous studies have revealed 0 tance, many aspects of nucleation process are still un- that such modular structures have a significant impact 1 clear and deserve more investigations. The Ising model, on the dynamics taking place on the networks, such as 1 which is a paradigm for many phenomena in statistical : physics, has been widely used to study the nucleation synchronization [34, 35], neural excitability [36], spread- v ing dynamics [37, 38], opinion formation [39, 40], and i process. Despite itssimplicity,Isingmodelhasmadeim- X Ising phase transition [41–43]. In particular, for major- portantcontributionstothe understandingofnucleation ity model [39] and Ising model [41–43] in modular net- r phenomena in equilibrium systems and is likely to yield a works, it has been shown that there exists a region in a important insights also for nonequilibrium systems. In discontinuous transition where modular order phase and two-dimensionallattices, forinstance,shearcanenhance globalorderphasecoexist. However,thesestudiesmainly the nucleation rate and the rate peaks at an intermedi- focused on phase diagrams in parameters space, and did ate shear rate [22], a single impurity may considerably notmakedetailedinvestigationforthetransitionprocess enhance the nucleation rate [23], and the existence of from a phase to another that may undergo a nucleation a pore may lead to two-stage nucleation and the overall process. nucleationratecanreachamaximumlevelataninterme- diate pore size [24]. Nucleation pathway of Ising model Sincenucleationisanactivatedprocessthatoccursex- tremely slow, brute-force simulation is prohibitively ex- pensive. Toovercomethisdifficulty,wewilluseavariant of a recently developed simulation method, forward flux ∗Electronicaddress: [email protected] sampling (FFS) [44]. This method allowsus to calculate 2 nucleationrateanddeterminethepropertiesofensemble spins pointing up. We are interested in the pathwayand towardnucleationpathways. Wefoundthatasthedegree rate for this nucleation process. ofnetworkmodularitydecreasesnucleationgoesthrough FFS method has been usedto calculaterate constants a transition from two-step to one-step process, and the and transition paths for rare events in equilibrium and rate exhibits a maximum at an intermediate degree of nonequilibriumsystems[22–24,44,47,48]. This method modularity. Free energy profiles for different modularity uses a seriesofinterfacesin phase spacebetween the ini- obtained by umbrella sampling (US) [45] and a simple tial and final states to force the system from the initial mean-field (MF) analysis help us understand the FFS stateAtothe finalstateB inaratchet-likemanner. Be- results. fore the simulation begin, an order parameter λ is first This paper is organized as follows. In Sec.II, we de- defined, such that the system is in state A if λ<λ and 0 scribethedetailsofourmodelandthesimulationmethod it is in state B if λ > λ . A series of nonintersecting M applied to this system. In Sec.III, we present the results interfaces λ (0 < i < M) lie between states A and B, i forthenucleationrateandpathwayinmodularnetworks. such that any path from A to B must cross each inter- InSec.IV,asimplemeanfieldanalysisisusedtoqualita- facewithoutreachingλ beforeλ . Thealgorithmfirst i+1 i tivelyillustratethesimulationresults. Atlast,discussion runs a long-time simulation which gives an estimate of and main conclusions are addressed in Sec.V. the fluxΦ¯ escapingfromthe basinofAandgenerates A,0 a collection of configurations corresponding to crossings of interface λ . The next step is to choose a configura- 0 II. MODEL AND SIMULATION DETAILS tion from this collection at random and use it to initiate a trial run which is continued until it either reaches λ 1 Consider a network consisting of N nodes arranged orreturns to λ0. If λ1 is reached,storethe configuration into two modules with N and N nodes. For simplicity, of the end point of the trial run. Repeat this step, each 1 2 we only consider the case ofN =N =N/2throughout time choosing a random starting configuration from the 1 2 this paper. The connection probability between a pair collectionatλ0. Thefractionofsuccessfultrialrunsgives of nodes belonging to the same module is ρi, while that an estimate of of the probability of reaching λ1 without for nodes belonging to different modules is ρo. The pa- going back into A, P(λ1|λ0). This process is repeated, rameter σ = ρo ∈ [0,1], defined as the ratio of inter- to stepbystep,untilλM isreached,givingtheprobabilities intra-modularρcionnectivity,measures the degreeof mod- P(λi+1|λi)(i=1,··· ,M−1). Finally,wegetthetransi- ularity. The higher degree ofmodularity ofa networkis, tionrateRfromAtoB,whichisthe productoftheflux the smaller value of σ is. As σ → 0, the network be- Φ¯A,0 andtheprobabilityP (λM|λ0)=QMi=−01P(λi+1|λi) comes two isolated clusters,while as σ →1, the network of reaching λM from λ0 without going into A. The de- approachesaErdo¨s–R´enyi(ER)randomnetwork. When tailed descriptions of FFS method see Ref.[49]. σ is varied the total number of links of the network is However, conventional FFS method will become very kept unchanged, Nhki, where hki is the average degree. inefficient if one intermediate metastable state exists be- 2 This restriction leads to ρ = 2hki and ρ = 2hkiσ . tween initial state and final state, as a two-step nucle- i N(1+σ) o N(1+σ) ation process demonstrated in Fig.1. This is because Eachnode is endowedwithanIsing spinvariables that i that sampling paths will be trapped in these long-lived can be either +1 (up) or −1 (down). The Hamiltonian metastable states so that they hardly return to initial of the system is given by state. To solve this problem, we will perform instead two-stepsamplings from initial down-spin state to inter- H =−JXi<jaijsisj −hXsi, (1) mediatemetastablestate,andthentofinalup-spinstate, i givingthe two-steprates, R andR , respectively. Since 1 2 the total mean time for nucleation is simply the sum of whereJ(>0)isthecouplingconstantandhistheexter- the meantime ofthe two-stepprocess,the totalratecan nalmagneticfield. Theelementsoftheadjacencymatrix ofthenetworktakeaij =1ifnodesiandj areconnected be expressed as R = (cid:0)R1−1+R2−1(cid:1)−1. To determine the and a =0 otherwise. location of the intermediate state, during FFS we moni- ij Our simulation is performed by Metropolis spin-flip tor the sampling time for the probability P(λi+1|λi) be- dynamics[46],inwhichweattempttoflipeachspinonce, tween two neighboring interfaces. If the sampling time on average, during each Monte Carlo (MC) cycle. In spent on between interfaces i and i + 1 is much more eachattempt, arandomlychosenspinisflippedwiththe than its previous step and the probability P(λi+1|λi) is probability min(1,e−β∆E), where β = 1/(k T) and k nearly one, we consider the ith interface as the location B B istheBoltzmannconstantandT isthetemperature,and ofthe intermediate state. If suchconditionsdo notmeet ∆E is the energy difference due to the flipping process. during the whole sampling, the intermediate metastable Here, we set J = 1, h > 0, T < T (T is the critical state does not exist, meaning that nucleation is a one- c c temperature),andstartwithametastablestateinwhich stepprocess. Note thatthe methodis straightforwardto s = −1 for most of the spins. The system will stay in generalize to a multi-step nucleation process. i that state fora significantly longtime before undergoing Here,wedefinetheorderparameterλasthetotalnum- a nucleation transition to a more stable state with most berofupspinsinthenetworks. The spacingbetweenin- 3 III. RESULTS d 400 To begin with, in Fig.1 we exhibit typical time evolu- e 300 tions of the number of up spins λ corresponding to two different values of network modularity σ = 0.011 and n σ = 0.051 via brute-force simulations, with relevant pa- 200 rameters N = 400, hki = 6, T = 2.0, and h = 1.2. It is c clearlyobservedthatthesystemundergoesatwo-stepnu- 100 a b cleation process for σ =0.011 and a one-step nucleation process for σ = 0.051. We also plot several represen- 0 tative configurations in Fig.2, corresponding to different phases of the system, respectively. Before the nucleation 6 6 6 6 6 6 0 1x10 2x10 3x10 4x10 5x10 6x10 happening, the system lies in a metastable state, where MC step most of the nodes are in down-spin state (indicated by blue circles), as shown in Fig.2(a) and Fig.2(d). When FIG. 1: (Color online) Typical time evolutions of the num- ber of up spins λ. It is shown that the system undergoes a the network modularity is very good, the system enters two-step nucleation process for σ =0.011 and a one-step nu- into an intermediate metastable state via the first-step cleation process for σ = 0.051. The representative network nucleation, where nodes in one of modules are in up- configurations at different moments indicated by arrows are spin state (indicated by red triangles), while nodes in shown in Fig.2. Other parameters are N = 400, hki = 6, theothermodulearestillindown-spinstate,asshownin T =2.0, and h=1.2. Fig.2(b). When the network modularity becomes worse, suchanintermediatemetastablestatedisappearssothat the nucleation becomes a one-step process. Finally, the system will enter into the most stable state, where al- most all spins are in up-spin state, as shown in Fig.2(c) andFig.2(e). Moreover,wenotethatthatthenucleation processtypicallytakestheorderof106 ormoreMCsteps thatiscomputationallycostly. Therefore,inwhatfollows we will give the results obtained by FFS method. The nucleation rate R as a function of σ is plotted in Fig.3(a), with relevant parameters being the same as those in Fig.1 except for h=1.0. One can see that as σ increasesRreachesamaximumR atσ ≃0.031andthen c decreases. Obviously, there exists a maximal nucleation rate that occurs at a moderate degree of network modu- larity. In Fig.3(b), we plot the results of the nucleation rates, R and R , for two-steps process as functions of 1 2 σ. As σ increases, R seems to exponentially decreases 1 withσ, whileR increasesmonotonicallyuntil σ =0.051 2 is reached. For σ > 0.051, nucleation becomes one-step processsothatR cannotbewelldefinedandtheoverall 2 FIG. 2: (Color online) Five representative network configu- nucleationrate is only determined by R . FromFig.3(b) 1 rations at different moments indicated in Fig.1, where down- one can find that R is much lower than R when the 2 1 spinnodesandup-spinnodesaredenotedbybluecirclesand value of σ is relatively small, so that R is dominantly red triangles, respectively. (a)-(c) correspond to the case of determined by the second step nucleation. While for σ=0.011 and (d)-(e)correspond to the case of σ=0.051. σ >0.031, R is almost determined only by the first step nucleation. Thus,thereexistsaregion0.001<σ <0.031 whereRisdeterminedbybothR andR . Notethatwe 1 2 have also made extensive simulations for other param- eters such as h = 0.7,1.2 and T = 1.5,1.8, and found terfaces is fixed at 3 up spins, but the computed results that the qualitative features of the above results do not do not depend on this spacing. We perform 1000 trials change (results now shown). per interface for each FFS sampling, from which at least To further understand the above results, we calculate 100configurationsaresavedateachinterfacein orderto freeenergyofthesystembyUSmethod,inwhichwehave investigatethestatisticalpropertiesofanensembleofre- used a bias potential 0.1k T(λ−λ¯)2, with λ¯ being the B active pathways to nucleation. The results are obtained centerofeachwindow. Thefreeenergy∆F asafunction by averagingover10 independent FFS samplings and50 ofλforthreedifferentvaluesofσaredepictedinFig.4(a). different network realizations. For σ = 0.001 and σ = 0.031, there are two free-energy 4 -27 (a) 0 -28 -50 (a) F -29 -100 nR-30 -150 l -200 -31 -250 -32 -300 -33 0 50 100 150 200 250 300 350 0.00 0.02 0.04 0.06 0.08 0.10 95 90 (b) 260 -12 (b) -16 85 250 lnR -20 lnR1 80 240 -24 lnR 2 -28 0.00 0.02 0.04 0.06 0.08 0.10 -32 24 20 0.00 0.02 0.04 0.06 0.08 0.10 * 212 (c) 15 * 2 F F 20 10 FIG. 3: (Color online) (a) The logarithm of nucleation rare lnR as a function of the degree of modularity σ. (b) lnR1 18 5 (squares) and lnR2 (circles) as a function of σ, and dotted line indicates theoverall rate lnR. The parameters are same 16 0 as in Fig.1 except for h=1.0. 0.00 0.02 0.04 0.06 0.08 0.10 maximums, occurring at the locations of critical nucleus FIG. 4: (Color online) (a) Free energy ∆F as a function of of λ = λ∗ and λ = λ∗, respectively. This picture is λ for three different σ = 0.001,0.031,0.101. For smaller σ 1 2 two free-energy barriers are clearly observed, while for larger consistentwiththetwo-stepnucleationprocessdescribed before. For a larger σ = 0.101, just the first free-energy σ thesecond onevanishes. (b)Thesizeof critical nucleusλ∗1 and λ∗2, and (c) the free-energy barriers ∆F1∗ and ∆F2∗, as barrier is present, implying that the nucleation becomes functionsof σ. Theotherparameters arethesame asFig.2. one-stepprocess. Withtheincrementofσ,λ∗ movestoa 1 largervaluewhilethevalueofλ∗getssmaller,asshownin 2 Fig.4(b). Fig.4(c)showsthatthefirstfree-energybarrier ∆F∗, defined as the difference between the free energy and all the nodes in the other module (module II) are 1 at λ∗ and the first minimum in free energy (λ = 23), in down-spins. The energy change due to the spin-flip 1 nearly increases linearly with σ, while the second free- of these λ nodes can be expressed as the sum of two energybarrier∆F2∗ (definition issimilarto thatof∆F1∗, parts ∆U1 = −2hλ+ 2JN1in, where the first part de- and the second minimum in free energy is an increasing notes the energy loss due to the creation of λ up-spins, function of σ, within the range λ ∈ [78,93]) decreases which favors the growth of the nucleus, while the sec- monotonically with σ until ∆F∗ vanishes at σ > 0.051, ond part denotes the energy gain due to the forma- which is in agreement with the2result of Fig.3(b). tion of N1in new interfacial links between up and down spins, which disfavors the growth of the nucleus. Ac- cording to MF approximation, Nin can be written as 1 Nin =ρ λ N −λ +ρ λN, wherethe first partandthe IV. MEAN FIELD ANALYSIS 1 i (cid:0)2 (cid:1) o 2 second part arise from interfacial links inside module I and between modules, respectively. For the second-step Inordertounveilpossiblemechanismbehindtheabove nucleation,we assume that allthe nodes inmodule I are phenomenon, we will present an analytical understand- inup-spins, andλ nodes are inup-spins andthe remain- ing by CNT and simple MF approximation. Firstly, let ing nodes are in down-spins in module II. This process us assume, for the first-step nucleation, that λ nodes creates new interfaciallinks inside module II, and at the are in up-spins and the remaining nodes are in down- sametimeremovesoldinterfaciallinksbetweenmoduleI spins in one of modules (say module I for convenience), andmoduleII.Thus,theenergychangeforthisprocessis 5 ofFig.4. FromFig.5,onecanseethatwiththeincrement 100 ofσ the size ofthe firstcriticalnucleus andthe heightof (a) 270 thefirstfree-energybarrierincreasealmostlinearly,while 95 the size of the second critical nucleus and the height of 90 the second free-energy barrier decrease until σ ≃ 0.13 is 260 reached. This implies that the analysis also predicts the 85 extinction of the second nucleation stage, but this pre- 250 dictionobviouslyoverestimatesthe transitionvalueofσ. 80 75 240 V. CONCLUSIONS 0.00 0.05 0.10 0.15 0.20 80 In conclusion, we have studied nucleation dynamics of Ising model in modular networks consisted of two ran- 110 dom networks. Using FFS method, we found that as (b) 60 the networkmodularity graduallybecomesworseatran- 100 F F sition occurs from one-step to two-step nucleation pro- 40 cess. Interestingly, the nucleation rate is a nonmono- 90 tonic function of the degree of modularity and a max- 20 imal rate exists for an intermediate level of modularity. 80 UsingUSmethod,weobtainedfreeenergyprofilesatdif- 0 ferent network modularity, from which one can see that 70 two free-energy barriers exist for very good modularity 0.00 0.05 0.10 0.15 0.20 and the second one vanishes when the network modu- larity becomes worse. This picture further confirms the FIG. 5: The results of mean field analysis. (a) The size of FFS results. Finally, a mean field analysis is employed ∆criFt1i∗caalnndu∆cleFu2∗s,λas∗1 faunndctλio∗2n,saonfdσ(.bT)htheeotfhreeer-epnaerragmyebtearrsriaerres to understand the nature of nucleation in modular net- worksandthesimulationresults. Sincestochasticfluctu- thesame as Fig.2. ationandthecoexistentofmulti-statesareubiquitousin social and biological systems, our study may shed valu- ∆U =−2hλ+2JNinwhereNin =ρ λ N −λ −ρ λN able insights into fluctuation-driven transition phenom- is t2he net number2of the int2erfaciail l(cid:0)in2ks. (cid:1)Theoen2- enatakingplaceinnetwork-organizedsystemswithmod- tropy changes for the two nucleation processes are both ular structures. ∆S = −kBN 2λln 2λ + 1− 2λ ln 1− 2λ . Then, 2 (cid:2)N (cid:0)N(cid:1) (cid:0) N(cid:1) (cid:0) N(cid:1)(cid:3) the changes of free energy for the two-step processes are ∆F = ∆U −T∆S (i = 1,2). In Fig.5 we give the ana- Acknowledgments i i lyticalresults ofthe critical nucleus and free-energybar- riersasfunctionsofthenetworkmodularity. Clearly,the This work was supported by NSFC (Nos. 20933006, analysis qualitatively agrees with the simulation results 20873130). [1] S.N.Dorogovtsev,A.V.Goltseve,andJ.F.F.Mendes, [10] S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Rev.Mod. Phys. 80, 1275 (2008). Phys. Rev.E 66, 016104 (2002). [2] R. Cohen, K. Erez, D. ben Avraham, and S. Havlin, [11] M. Leone, A.Va´zquez, A. Vespignani, and R. Zecchina, Phys.Rev.Lett. 85, 4626 (2000). Eur. Phys. J. B 28, 191 (2002). [3] D. S. Callaway, M. E. J. Newman, S. H. Strogatz, and [12] C. P. Herrero, Phys. Rev.E 69, 067109 (2004). D.J. Watts, Phys.Rev.Lett. 85, 5468 (2000). [13] T. Nishikawa, A. E. Motter, Y.-C. Lai, and F. C. Hop- [4] M. E. J. Newman, Phys.Rev.Lett. 89, 208701 (2002). pensteadt, Phys. Rev.Lett. 91, 014101 (2003). [5] R. Cohen, D. ben Avraham, and S. Havlin, Phys. Rev. [14] J. G´omez-Garden˜es, Y. Moreno, and A. Arenas, Phys. E 66, 036113 (2002). Rev. Lett.98, 034101 (2007). [6] R. Pastor-Satorras and A. Vespignani, Phys. Rev. Lett. [15] K.-I.Goh,D.-S.Lee,B.Kahng,andD.Kim,Phys.Rev. 86, 3200 (2001). Lett. 91, 148701 (2003). [7] C. Castellano and R. Pastor-Satorras, Phys. Rev. Lett. [16] A. E. Motter and Y. C. Lai, Phys. Rev. E 66, 065102 105, 218701 (2010). (2002). [8] A.Aleksiejuk, J.A.Holysta,and D.Stauffer,PhysicaA [17] H.NakaoandA.S.Mikhailov,Nat.Phys.6,544(2010). 310, 260 (2002). [18] D. Kashchiev, Nucleation: basic theory with applications [9] G. Bianconi, Phys. Lett.A 303, 166 (2002). (Butterworths-Heinemann, Oxford, 2000). 6 [19] L. Gr´ana´sy and F. Iglo´i, J. Chem. Phys. 107, 3634 Buldu´, S. Havlin, and S. Boccaletti, Phys. Rev. Lett. (1997). 101, 168701 (2008). [20] G.Johnson,A.I.Mel´cuk,H.Gould,W.Klein,andR.D. [36] C. Zhou, L. Zemanova´, G. Zamora, C. C. Hilgetag, and Mountain, Phy.Rev.E 57, 5707 (1998). J. Kurths,Phys. Rev.Lett. 97, 238103 (2006). [21] A. R. Fersht, Proc. Natl. Acad. Sci. USA 92, 10869 [37] Z. Liu and B. Hu,Europhys.Lett. 72, 315 (2005). (1995). [38] L. Huang, K. Park, and Y.-C. Lai, Phys. Rev. E 73, [22] R. J. Allen, C. Valeriani, S. Tanase-Nicola, P. R. ten 035103 (2006). Wolde, and D. Frenke, J. Chem. Phys. 129, 134704 [39] R. Lambiotte and M. Ausloos, J. Stat. Mech. p. P08026 (2008). (2007). [23] R.P. Sear, J. Phys. Chem. B 110, 4985 (2006). [40] R. Lambiotte, M. Ausloos, and J. A. Hol yst, Phys. Rev. [24] A. J. Page and R. P. Sear, Phys. Rev. Lett. 97, 065701 E 75, 030101 (2007). (2006). [41] R. K. Pan and S. Sinha, Europhys. Lett. 85, 68006 [25] A.C.PanandD.Chandler,J.Phys.Chem.B108,19681 (2009). (2004). [42] S. Dasgupta, R.K.Pan, and S.Sinha,Phys.Rev.E80, [26] M. Acharyya and D. Stauffer, Eur. Phys. J. B 5, 571 025101(R) (2009). (1998). [43] K. Suchecki and J. A. Hol yst, Phys. Rev. E 80, 031110 [27] V. A. Shneidman, K. A. Jackson, and K. M. Beatty, J. (2009). Chem. Phys.111, 6932 (1999). [44] R. J. Allen, P. B. Warren, and P. R. ten Wolde, Phy. [28] S. Wonczak, R. Strey, and D. Stauffer, J. Chem. Phys. Rev. Lett.94, 018104 (2005). 113, 1976 (2000). [45] J.S.vanDuijneveldtandD.Frenkel,J.Chem.Phys.96, [29] K. Brendel, G. T. Barkema, and H. van Beijeren, Phys. 15 (1992). Rev.E 71, 031601 (2005). [46] D. P. Landau and K. Binder, A Guide to Monte Carlo [30] S.RyuandW.Cai, Phys.Rev.E81, 030601(R) (2010). Simulations in Statistcal Physics (Cambridge University [31] H. Chen, C. Shen, Z. Hou, and H. Xin, arXiv:1008.0704 Press, Cambridge, 2000). (2010). [47] C. Valeriani, R. J. Allen, M. J. Morelli, D. Frenkel, and [32] M.E.J.Newman,Proc.Natl.Acad.Sci.USA103,8577 P. R. ten Wolde, J. Chem. Phys.127, 114109 (2007). (2006). [48] R. J. Allen, D. Frenkel, and P. R. ten Wolde, J. Chem. [33] S.Fortunato, Phys. Rep.486, 75 (2010). Phys. 124, 024102 (2006). [34] A. Arenas, A. D´ıaz-Guilera, and C. J. P´erez-Vicente, [49] R. J. Allen, C. Valeriani, and P. R. ten Wolde, J. Phys.: Phys.Rev.Lett. 96, 114102 (2006). Condens. Matter 21, 463102 (2009). [35] D.Li,I.Leyva,J.A.Almendral,I.Sendin˜aNadal,J.M.