Locations of Roberge-Weiss transition endpoints in lattice QCD with N_f=2 improved KS quarks
aa r X i v : . [ h e p - l a t ] D ec Locations of Roberge-Weiss transition endpoints in lattice QCD with N f = 2 improvedKS quarks Liang-Kai Wu ∗ Faculty of Science, Jiangsu University, Zhenjiang 212013, Peoples Republic of China
Xiang-Fei Meng
National Supercomputer Center in Tianjin, Tianjin, 300457, Peoples Republic of China (Dated: October 12, 2018)Result on the locations of the tricritical points of N f = 2 lattice QCD with imaginary chemicalpotential is presented. Simulations are carried out with Symanzik improved gauge action and Asqtadfermion action. With imaginary chemical potential iµ I = iπT , previous studies show that theRoberge-Weiss (RW) transition endpoints are triple points at both large and small quark masses,and second order transition points at intermediate quark masses. The triple and second orderendpoints are separated by two tricritical ones. Our simulations are carried out at 7 values of quarkmass am ranging from 0.024 to 0.070 on lattice volume 12 × , × , ×
4. The susceptibilityand Binder cumulant of the imaginary part of Polyakov loop are employed to determine the natureof RW transition endpoints. The simulations suggest that the two tricritical points are within therange 0 . − .
026 and 0 . − . PACS numbers: 12.38.Gc, 11.10.Wx, 11.15.Ha, 12.38.Mh
I. INTRODUCTION
The phase diagram of QCD has significantly phe-nomenological implications. It is relevant to the earlyUniverse, compact stars and heavy ion collision experi-ments. Reviews on the study of phase diagram can befound in Refs. [1, 2] and references therein. While sub-stantial lattice simulation has focused on the phase ofQCD at finite density, a great amount of study centresaround QCD with imaginary chemical potential. QCDwith imaginary chemical potential has a rich phase struc-ture, and it not only deserves detailed investigations in itsown right theoretically, but also has significant relevanceto physics at zero or small real chemical potential [3–11].The Z(3) symmetry which is present in the pure gaugetheory is explicitly broken at the presence of dynamicalquarks. However, Ref. [12] shows that the Z(3) sym-metry is restored when imaginary chemical potential isturned on and Z(3) transformation can be compensatedby a shift in µ I /T by 2 π/
3, so the partition function ofQCD with imaginary chemical potential has periodicityin µ I /T with period 2 π/N c as well as reflection symmetryin µ = iµ I .Different Z(3) sectors are distinguished by the phaseof Polyakov loop. At high temperature, the spontaneousbreaking of Z(3) symmetry implies transition betweenadjacent Z(3) sectors in µ I and this transition is of firstorder, while at low temperature, unbroken Z(3) symme-try guarantees the transition is analytic. The first ordertransition takes place at those critical values of imagi-nary chemical potential µ I /T = (2 n + 1) π/ ∗ Corresponding author. Email address: [email protected] a transition line which necessarily ends at an endpoint T RW when the temperature is decreased sufficiently low.Recent numerical studies [3–5, 15, 16] show that theRW transition endpoints are triple points for small andheavy quark masses, and second order points for interme-diate quark masses. So there exist two tricritical pointsseparating the first order transition points from the sec-ond ones. Moreover, it is pointed out [3, 10, 11] that thescaling behaviour at the tricritical points may shape thecritical line which separate different transition region forreal chemical potential, and thus, the critical line for realchemical potential is expected to be qualitatively consis-tent with the scenario suggested in Refs. [17, 18] whichshows that the first order transition region shrinks withincreasing real chemical potential. In addition, Ref. [19]employs the scaling behaviour at the tricritical point todetermine the nature of 2 flavour QCD transition in thechiral limit.So far, the investigation for the Roberge-Weiss transi-tion endpoints are implemented through standard gaugeand fermion actions. In this paper, we aim to inves-tigate the endpoints of N f = 2 QCD with one-loopSymanzik-improved gauge action [20–23] and Asqtad KSaction [24, 25]. These actions have discretization errorof O ( α s a , a ) and O ( α s a , a ), respectively. These im-provements are significant on N t = 4 lattice where thelattice spacing is quite large. Standard KS fermions sufferfrom taste symmetry breaking at nonzero lattice spacing a [26]. This taste symmetry breaking can be illustratedby the smallest pion mass taste splitting which is compa-rable to the pion mass even at lattice spacing a ∼ . f m [27]. Asqtad KS action has good taste symmetry and freedispersion relation by introducing fattened links and theso-called ”Naik terms” [28, 29].The paper is organized as follows. In Sec. II, we definethe lattice action with imaginary chemical potential andthe physical observables we calculate. Our simulationresults are presented in Sec. III followed by discussionsin Sec. IV. II. LATTICE FORMULATION WITHIMAGINARY CHEMICAL POTENTIAL
After introducing pseudofermion field Φ, the partitionfunction of the system can be represented as: Z = Z [ dU ][ d Φ ∗ ][ d Φ] e − S g − S f , where S g is the Symanzik-improved gauge action, and S f is the Asqtad quark action with the quark chemicalpotential µ . Here µ = iµ I . For S g , we use S G = β C P X x ; µ<ν (1 − P µν ) + C R X x ; µ = ν (1 − R µν )+ C T X x ; µ<ν<σ (1 − T µνσ ) ! , with P µν , R µν and T µνσ standing for 1 / ×
1, 1 × × × P µν = 13 ReTr ✲ ✻✛❄ ,R µν = 13 ReTr ✲ ✲ ✻✛✛❄ ♣♣♣♣♣♣♣♣♣ ,T µνσ = 13 ReTr ✲ ✟✟✯✟✟✻✛✟✟✙✟✟❄ ♣ ♣ ♣ ♣ ♣ ♣♣♣♣♣♣♣♣♣♣ ♣ ♣ ♣ ♣ ♣ ♣ ♣ ♣ ♣ . The coefficents C P , C R , C T are tadpole improved [27], C P = 1 . ,C R = − u (1 − (0 . − . n f ) ln( u )) ,C T = 1 u (0 . − . n f ) ln( u ) . The Asqtad action with pseudofermion field Φ is S f = D Φ (cid:12)(cid:12)(cid:12)(cid:2) M † [ U ] M [ U ] (cid:3) − n f / (cid:12)(cid:12)(cid:12) Φ E , where the form of M x,y [ U ] = 2 m x,y + D x,y ( U ) reading2 mδ x,y + P ρ =1 η x,ρ (cid:16) U Fx,ρ δ x,y − ˆ ρ − U F † x − ˆ ρ,ρ δ x,y +ˆ ρ (cid:17) + η x, (cid:16) e iaµ I U Fx, δ x,y − ˆ4 − e − iaµ I U F † x − ˆ4 ,µ δ x,y +ˆ4 (cid:17) + P ρ =1 η x,ρ (cid:16) U Lx,ρ δ x,y − ρ − U L † x − ρ,ρ δ x,y +3ˆ ρ (cid:17) + η x, (cid:16) e i aµ I U Lx, δ x,y − − e − i aµ I U L † x − ˆ4 ,µ δ x,y +3ˆ4 (cid:17) , where U Fx,ρ stands for the fattened link which is producedby Fat7 smearing and U Lx,ρ stands for the naik term.ˆ ρ, ˆ4 are the unit vector along ρ − direction,4 − direction,respectively. η x,µ is the staggered fermion phase.We carry out simulations at θ = µ I /T = π . As it ispointed out that the system is invariant under the chargeconjugation at θ = 0 , π , when θ is fixed [9]. But the θ -odd quantity O ( θ ) is not invariant at θ = π under chargeconjugation. When T < T RW , O ( θ ) is a smooth functionof θ , so it is zero at θ = π . Whereas when T > T RW , thetwo charge violating solutions cross each other at θ = π . Thus the charge symmetry is spontaneously brokenthere and the θ -odd quantity O ( θ ) can be taken as orderparameter . In this paper, we take the imaginary part ofPolyakov loop as the order parameter.The expression of Polyakov loop L is defined as thefollowing: h L i = * L s L t X x Tr " L t Y t =1 U ( x , t ) ,L s , L t are the spatial, time extent of lattice, respectively.To simplify notation, we use X to represent the imagi-nary part of Polyakov loop Im( L ). The susceptibility ofimaginary part of Polyakov loop Im( L ) is defined as χ = L s (cid:10) ( X − h X i ) (cid:11) , which is expected to scale as: [4, 5] χ = L γ/νs φ ( τ L /νs ) , (1)where τ is the reduced temperature τ = ( T − T RW ) /T RW .This means that the curves χ/L γ/νs at different latticevolume should collapse with the same curve when plottedagainst τ L /νs . In the following, we employ β − β RW inplace of τ = ( T − T RW ) /T RW . The critical exponentsrelevant to our study are collected in Table. I [5, 30].We also consider the Binder cumulant of imaginarypart of Polyakov loop which is defined as the following: B = (cid:10) ( X − h X i ) (cid:11) / (cid:10) ( X − h X i ) (cid:11) (2)with h X i = 0. In the vicinity of the RW transition lineendpoints, B with the finite size correction is a functionof x = ( β − β RW ) L /νs and can be expanded as a series [3,10, 11], B = B ( β c , ∞ ) + a x + a x + · · · . (3) -60 -40 -20 0 20 40 600.0000.0020.0040.0060.0080.0100.0120.0140.016 am=0.024 L S =12 L S =16 L S =20 / L S / ( - C ) L S1/ -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.00.100.150.200.250.300.35 am=0.024 L S =12 L S =16 L S =20 / L S / ( - C ) L S1/
FIG. 1: Scaling behavior of the susceptibility of imaginary part of Polyakov loop according to first order critical index (leftpanel), and to 3D Ising critical index (right panel) at am = 0 . -60 -40 -20 0 20 40 600.000.010.020.030.040.05 am=0.038 L S =12 L S =16 L S =20 / L S / ( - C ) L S1/ -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.00.200.250.300.350.400.450.500.550.600.65 am=0.038 L S =12 L S =16 L S =20 / L S / ( - C ) L S1/
FIG. 2: Scaling behavior of the susceptibility of imaginary part of Polyakov loop according to first order critical index (leftpanel), and to 3D Ising critical index (right panel) at am = 0 . B ( β c , ∞ ) ν γ γ/ν
3D ising 1.604 0.6301(4) 1.2372(5) 1.963tricritical 2 1/2 1 2first order 1.5 1/3 1 3crossover 3 - - -TABLE I: Critical exponents relevant to our study.
In the thermodynamic limit, the critical index ν and B ( β c , ∞ ) takes on the corresponding value summarizedin Table. I. However, on finite spatial volumes, the stepsof B ( β c , ∞ ) are smeared out to continuous functions. III. MC SIMULATION RESULTS
Before presenting the simulation results, we describethe simulation details. Simulations are carried out atquark mass am = 0 . , . , . , . , . , . , . M + ( U ) M ( U )] − n f / for the pseudofermion field. Forthe heat bath updating and for computing the actionat the beginning and end of the molecular dynamicstrajectory 10’th rational function is used to approxi-mate [ M + ( U ) M ( U )] n f / and [ M + ( U ) M ( U )] − n f / , re-spectively. The step is chosen to ensure the acceptancerate is around 80% − -60 -40 -20 0 20 40 600.000.010.020.030.040.050.060.070.08 am=0.040 L S =12 L S =16 L S =20 / L S / ( - C ) L S1/ -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.00.00.10.20.30.40.50.60.70.80.91.01.11.21.31.41.51.61.71.8 am=0.040 L S =12 L S =16 L S =20 / L S / ( - C ) L S1/
FIG. 3: Scaling behavior of the susceptibility of imaginary part of Polyakov loop according to first order critical index (leftpanel), and to 3D Ising critical index (right panel) at am = 0 . -60 -40 -20 0 20 40 600.010.020.030.04 am=0.070 L S =12 L S =16 L S =20 / L S / ( - C ) L S1/ -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.00.10.20.30.40.50.60.70.8 am=0.070 L S =12 L S =16 L S =20 / L S / ( - C ) L S1/
FIG. 4: Scaling behavior of the susceptibility of imaginary part of Polyakov loop according to first order critical index (leftpanel), and to 3D Ising critical index (right panel) at κ = 0 . ration are taken as warmup form a cold start. In orderto fill in observables at additional β values, we employthe Ferrenberg-Swendsen reweighting method [36].The critical coupling β RW ’s on various spatial volumeat different quark mass am are summarized in Table. II.These β RW ’s are determined from the locations of peaksusceptibility of imaginary part of Polyakov loop.We present the rescaling susceptibility of imaginarypart of Polyakov loop χ/L sγ/ν as a function of ( β − β RW ) L /νs at am = 0 .
024 in Fig. 1. From Fig. 1, wecan find that χ/L sγ/ν according to the first order transi-tion index collapses with the same curve, while χ/L sγ/ν according to 3D index does not.The rescaling susceptibility of imaginary part ofPolyakov loop χ/L sγ/ν as a function of ( β − β RW ) L /νs at am = 0 .
038 is depicted in Fig. 2. From Fig. 2, we
TABLE II: Results of critical couplings β RW on different spa-tial volume at different κ . am
12 16 200.024 6 . . . . . . . . . . . . . . . . . . can find that χ/L sγ/ν according to the first order tran-sition index or 3D index does not collapse with the same -60 -40 -20 0 20 40 600.0100.0120.0140.0160.0180.0200.0220.0240.0260.0280.0300.0320.0340.0360.0380.040 am=0.050 L S =12 L S =16 L S =20 / L S / ( - C ) L S1/ -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.00.20.30.40.50.60.70.80.9 am=0.050 L S =12 L S =16 L S =20 / L S / ( - C ) L S1/
FIG. 5: Scaling behavior of the susceptibility of imaginary part of Polyakov loop according to first order critical index (leftpanel), and to 3D Ising critical index (right panel) at κ = 0 . am=0.024 L S =12 L S =16 L S =20 B -60 -40 -20 0 20 40 602.002.052.102.152.202.252.302.352.402.452.502.552.602.652.702.752.80 am=0.024 L S =12 L S =16 L S =20 B ( - C ) L S1/
FIG. 6: Binder cumulants as a function of β on various spatial volume intersect at one point (left panel), and as a function of( β − β c ) L /νs with values of β c , ν from Table. III collapse (right panel) at am = 0 . curve. We cannot determine the nature of Roberge-Weisstransition endpoint at am = 0 .
038 from χ/L sγ/ν .The behaviour of rescaling susceptibility of imaginarypart of Polyakov loop χ/L sγ/ν at am = 0 .
040 and am =0 .
070 are presented in Fig. 3, and Fig. 4 respectively.Form Fig. 3 and Fig. 4, we can find that The rescalingsusceptibility of imaginary part of Polyakov loop χ/L sγ/ν at am = 0 .
040 and am = 0 .
070 have similar behaviourto the that at am = 0 . χ/L sγ/ν as a function of ( β − β RW ) L /νs at am = 0 .
050 is depicted in Fig. 5. From Fig. 5, wecan find that χ/L sγ/ν as a function of ( β − β RW ) L /νs atlattice 12 × × χ/L sγ/ν and ( β − β RW ) L /νs in Fig. 5,the first order transition index may be the better choice. χ/L sγ/ν as a function of ( β − β RW ) L /νs at am = 0 . am = 0 .
050 which tendsto be in favour of first order transition index.In order to discern the scaling behaviour, we turn toinvestigate Binder cumulant B as defined in Eq. (2)whose scaling behaviour is described in Eq. (3). B de-creases with the increase of β , and at one fixed quarkmass am , B as a function of β on various spatial vol-ume is expected to intersect at one point. The intersec-tion gives an estimate of accurate location of β RW . Byfitting to Eq. (3), we can extract critical index ν, β RW and B ( β c , ∞ ). The results are collected in Table. III.We present B as a function of β at am = 0 .
024 inthe left panel of Fig. 6, and B as a function of ( β − am=0.026 L S =12 L S =16 L S =20 B -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.201.551.601.651.701.751.801.851.901.952.00 am=0.026 L S =12 L S =16 L S =20 B ( - C ) L S1/
FIG. 7: Binder cumulants as a function of β on various spatial volume intersect at one point (left panel), and as a function of( β − β c ) L /νs with values of β c , ν from Table. III collapse (right panel) at am = 0 . am=0.040 L S =12 L S =20 B -0.3 -0.2 -0.1 0.0 0.1 0.2 0.31.71.81.92.02.12.22.32.42.52.62.7 am=0.040 L S =12 L S =20 B ( - C ) L S1/
FIG. 8: Binder cumulants as a function of β on various spatial volume intersect at one point (left panel), and as a function of( β − β c ) L /νs with values of β c , ν from Table. III collapse (right panel) at am = 0 . β RW ) L /νs in the right panel of Fig. 6 with ν taken tobe the extracted value through fitting procedure. FromTable. III, we find that the critical index ν = 0 . am = 0 .
024 can explain the behaviour of B as a functionof ( β − β RW ) L /νs , especially, on lattice L s = 16 ,
20. Thisbehaviour implies that the transition endpoint at am =0 .
024 belongs to first order transition.We also present B as a function of β at am = 0 . B as a function of ( β − β RW ) L /νs in the right panel of Fig. 7 with ν taken tobe the extracted value through fitting procedure. Wefind that the critical index ν = 0 . am = 0 . B as a function of ( β − β RW ) L /νs . ν = 0 . am = 0 .
026 is of 3D transition nature. At am = 0 . B as a functionof β on lattice L s = 12 ,
20 intersects at one point. B as a function of β and as a function of ( β − β RW ) L /νs at am = 0 .
040 are depicted in the left, right panel ofFig. 7, respectively. The extracted value ν = 0 . am = 0 .
040 is of 3D transition nature. At other val-ues of am , B as a function of β and as a function of( β − β RW ) L /νs have similar behaviour. For clarity, theyare not presented.From the behaviour of χ/L sγ/ν and B , we concludethat the nature of endpoint transition at am = 0 . , . , . , .
070 is of first order, while at am = 0 . , . , . , the endpoint transition nature is of 3D Isingclass. This conclusion suggests that the two tricritical TABLE III: Results of critical couplings β RW and the critical index ν by fitting Eq. (3) to data on different spatial volume. Iferrors are very small, we take them to be zero. am L s β RW ν B ( β c , ∞ ) a a r-square0.024 12 , ,
20 6 . . . − . . , ,
20 6 . . . − . . ,
20 6 . . . − . . ,
20 6 . . . − . . , ,
20 6 . . . − . . ,
20 6 . . . − . − ,
20 6 . . . − . − points are between 0 . < am tricl < .
026 and 0 . We have studied the nature of critical endpoints ofRoberge-Weiss transition of two flavor lattice QCD withimproved KS fermions. When iµ I = iπT , the imaginarypart of Polayakov loop is the order parameter for study-ing the transition from low temperature phase to hightemperature one.Our simulations are carried out at 7 values of quarkmass am on L t = 4 lattice on different 3 spatial volumes.Our central result is that the two tricritical points arebetween 0 . < am tricl < . 026 and 0 . < am tricl < . ν is ex-pected to change smoothly, while our simulation showsthat index ν changes rapidly within a narrow quark massinterval.In Ref. [5], the two locations of tricritical point for N f = 2 QCD are am = 0 . , . N f = 3 QCD, Ref. [3] concludes that the two tri-critical points are between 0 . < am tricl < . . < am tricl < . 5. Comparing with those results, thesecond transition region from our simulation is narrow.Apart from monitoring the behaviour of susceptilityof imaginary part of Polayakov loop Im( L ), we also lookinto the change of Binder cumulant of Im( L ). In or-der to fill in observables at additional β values, theFerrenberg-Swendsen reweighting method [36] is em-ployed. It is noted that when applying Ferrenberg-Swendsen reweighting method, the number of β pointstaken to calculate susceptility is not completely the sameas the number taken to calculate Binder cumulant.In our simulations, the behaviour of susceptility ofimaginary part of Polayakov loop Im( L ) at am = 0 . 024 can give us clear signal to determine the nature of transi-tion, while at other quark mass, it is difficult to determinethe nature of transition.The values of B ( β c , ∞ ) extracted through fitting pro-cedure are not in consistent with what are expected. Thisis because logarithmic scaling corrections will be presentnear the tricritical point [3, 37], and our simulations arecarried out on finite size volume on which large finitesize corrections are observed in simpler spin model [38].However, the critical exponent ν is not sensitive to finitesize corrections [3]. So index ν extracted through fittingprocedure can provide us information to determine thetransition nature.In our simulation, we can find that the values of B on lattice with spatial volumes 12 , , inter-sect approximately at one point at quark masses am =0 . , . , . B ’s from three spatialvolumes. It is expedient to determine the intersectionpoint from two spatial volumes as indicated in Table. III.Taking what mentioned above into account, furtherwork along this direction which can provide crosscheck isexpected, especially simulations with larger time extentwhich is being under our consideration. Acknowledgments We thank Philippe de Forcrand for valuable helps. Wemodify the MILC collaboration’s public code [39] to sim-ulate the theory at imaginary chemical potential. We usethe fortran-90 based multi-precision software [40]. Thiswork is supported by the National Science Foundationof China (NSFC) under Grant Nos. (11347029). Thework was carried out at National Supercomputer Cen-ter in Wuxi, We appreciate the help of Qiong Wang andZhao liu when carrying out the computation. [1] K. Fukushima and T. Hatsuda, Rept. Prog. Phys. [3] P. de Forcrand and O. Philipsen, Phys. Rev. Lett. ,152001 (2010) [arXiv:1004.3144 [hep-lat]].[4] M. D’Elia and F. Sanfilippo, Phys. Rev. D , 111501(2009) [arXiv:0909.0254 [hep-lat]].[5] C. Bonati, G. Cossu, M. D’Elia and F. Sanfilippo, Phys.Rev. D , 054505 (2011) [arXiv:1011.4515 [hep-lat]].[6] M. D’Elia, F. Di Renzo and M. P. Lombardo, Phys. Rev.D , 114509 (2007) [arXiv:0705.3814 [hep-lat]].[7] Y. Sakai, H. Kouno and M. Yahiro, J. Phys. G , 105007(2010) [arXiv:0908.3088 [hep-ph]].[8] G. Aarts, S. P. Kumar and J. Rafferty, JHEP , 056(2010) [arXiv:1005.2947 [hep-th]].[9] H. Kouno, Y. Sakai, K. Kashiwa and M. Yahiro, J. Phys.G , 115010 (2009) [arXiv:0904.0925 [hep-ph]].[10] O. Philipsen and P. de Forcrand, PoS LATTICE ,211 (2010) [arXiv:1011.0291 [hep-lat]].[11] C. Bonati, P. de Forcrand, M. D’Elia, O. Philipsenand F. Sanfilippo, PoS LATTICE , 189 (2011)[arXiv:1201.2769 [hep-lat]].[12] A. Roberge and N. Weiss, Nucl. Phys. B , 734 (1986).[13] P. de Forcrand and O. Philipsen, Nucl. Phys. B , 290(2002).[14] M. D’Elia and M. -P. Lombardo, Phys. Rev. D , 014505(2003) [hep-lat/0209146].[15] O. Philipsen and C. Pinke, Phys. Rev. D ,no. 9, 094504 (2014) doi:10.1103/PhysRevD.89.094504[arXiv:1402.0838 [hep-lat]].[16] L. K. Wu and X. F. Meng, Phys. Rev. D ,no. 9, 094506 (2014) doi:10.1103/PhysRevD.90.094506[arXiv:1405.2425 [hep-lat]].[17] P. de Forcrand and O. Philipsen, JHEP , 077 (2007)[hep-lat/0607017].[18] P. de Forcrand and O. Philipsen, JHEP , 012 (2008)[arXiv:0808.1096 [hep-lat]].[19] C. Bonati, P. de Forcrand, M. D’Elia, O. Philipsen andF. Sanfilippo, Phys. Rev. D , no. 7, 074030 (2014)doi:10.1103/PhysRevD.90.074030 [arXiv:1408.5086 [hep-lat]].[20] K. Symanzik, Nucl. Phys. B , 187 (1983).doi:10.1016/0550-3213(83)90468-6[21] M. Luscher and P. Weisz, Phys. Lett. , 250 (1985).doi:10.1016/0370-2693(85)90966-9 [22] G. P. Lepage and P. B. Mackenzie, Phys. Rev.D , 2250 (1993) doi:10.1103/PhysRevD.48.2250[hep-lat/9209022].[23] M. G. Alford, W. Dimm, G. P. Lepage, G. Hockneyand P. B. Mackenzie, Phys. Lett. B , 87 (1995)doi:10.1016/0370-2693(95)01131-9 [hep-lat/9507010].[24] T. Blum et al. , Phys. Rev. D , 1133 (1997)doi:10.1103/PhysRevD.55.1133 [hep-lat/9609036].[25] K. Orginos et al. [MILC Collaboration], Phys. Rev.D , 054503 (1999) doi:10.1103/PhysRevD.60.054503[hep-lat/9903032].[26] A. Bazavov et al. , Phys. Rev. D , 054503 (2012)doi:10.1103/PhysRevD.85.054503 [arXiv:1111.1710 [hep-lat]].[27] A. Bazavov et al. [MILC Collaboration], Rev. Mod. Phys. 238 (1989).[29] C. W. Bernard et al. [MILC Collaboration], Phys. Rev.D , 549 (2002)[cond-mat/0012164].[31] M. A. Clark and A. D. Kennedy, Nucl. Phys. Proc. Suppl. 850 (2004).[32] M. A. Clark and A. D. Kennedy, Phys. Rev. D 272 (2003).[36] A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. , 1195 (1989).[37] I.D. Lawrie and S. Sarbach, in Phase transitions and crit-ical phenomena , eds. C. Domb and J.L.Lebowitz, vol.9,1 (1984).[38] A. Billoire, T. Neuhaus and B. Berg, Nucl. Phys. B ,779 (1993) [hep-lat/9211014].[39] http://physics.utah.edu/~detar/milc/ [40][40]