Deconfinement critical point of lattice QCD with N f =2 Wilson fermions
Francesca Cuteri, Owe Philipsen, Alena Schön, Alessandro Sciarra
TThe deconfinement critical point of lattice QCD with N f “ Wilson fermions
Francesca Cuteri, ˚ Owe Philipsen,
1, 2, : Alena Schön, ; and Alessandro Sciarra § Institut für Theoretische Physik, Goethe-Universität FrankfurtMax-von-Laue-Str. 1, 60438 Frankfurt am Main, Germany John von Neumann Institute for Computing (NIC) GSI, Planckstr. 1, 64291 Darmstadt, Germany
The SU p q pure gauge theory exhibits a first-order thermal deconfinement transition due to spon-taneous breaking of its global Z center symmetry. When heavy dynamical quarks are added, thissymmetry is broken explicitly and the transition weakens with decreasing quark mass until it dis-appears at a critical point. We compute the critical hopping parameter and the associated pionmass for lattice QCD with N f “ degenerate standard Wilson fermions on N τ P { , , } lattices,corresponding to lattice spacings a “ .
12 fm , a “ .
09 fm , a “ .
07 fm , respectively. Significantcut-off effects are observed, with the first-order region growing as the lattice gets finer. While cur-rent lattices are still too coarse for a continuum extrapolation, we estimate m cπ « with aremaining systematic error of „ . Our results allow to assess the accuracy of the LO and NLOhopping expanded fermion determinant used in the literature for various purposes. We also providea detailed investigation of the statistics required for this type of calculation, which is useful forsimilar investigations of the chiral transition. PACS numbers: 12.38.Gc, 05.70.Fh, 11.15.HaKeywords: QCD phase diagram
I. INTRODUCTION
For physical quark mass values, the thermal QCD tran-sition is known to be an analytic crossover [1]. Its con-tinuation to small baryo-chemical potentials has beenstudied in detail by means of Taylor expansion or an-alytic continuation from imaginary chemical potential,with so far no hints of a non-analytic phase transition [2].Since a severe sign problem precludes lattice QCD simu-lations at finite baryon density, one way to constrain thephase diagram is to study the thermal transition withthe QCD parameters varied away from their physical val-ues. Non-analytic phase transitions related to the spon-taneous breaking of the global center and chiral sym-metries, respectively, are explicitly seen in simulationsemploying unimproved fermion discretizations on coarselattices, in the heavy and light quark mass regime, asindicated in Figure 1. One can then study how thesecritical structures evolve when a chemical potential isswitched on, for an overview see Ref. 3.In this work we focus on the heavy mass corner, whosethermodynamics can be addressed also at finite baryonchemical potentials, either by means of effective latticetheories obtained from hopping expansions [4–6], or by ef-fective Polyakov loop theories in the continuum [7, 8]. Inorder to assess the accuracy of these approaches and theirpossible extensions to light quarks, reliable benchmarksare warranted. Moreover, the phase transitions in thisparameter regime are interesting in their own right. Inthe infinite quark mass limit, the theory reduces to SU p q ˚ [email protected] : [email protected] ; [email protected] § [email protected] pure gauge theory with its first-order transition [9], whichis caused by the spontaneous breaking of the Z centersymmetry. Recently, the latent heat associated with thistransition has also been determined [10]. When dynami-cal quarks are added to the theory, the center symmetryis broken explicitly and the phase transition weakens,i.e, the latent heat decreases until it vanishes at a criticalquark mass. The value of the critical quark mass andthe latent heat are thus intimately related by the dy-namics of the deconfinement transition, so that valuablenon-perturbative insights are possible in this parameterregion.An early lattice investigation for N f “ Wilsonfermions was restricted to N τ “ and, for large vol-umes, required mapping to an auxiliary model [11]. Thefirst systematic studies in the standard Wilson discretiza-tion [12, 13] on N τ P { , , } lattices are based on an ef-fective potential for the order parameter determined bya histogram method [5], which is reweighted with a hop-ping expanded quark determinant from quenched config-urations. This allows for a flexible investigation and in-terpolation between the N f P { , , } cases. In this workwe focus on N f “ only, but simulate the full fermiondeterminant with no approximation beyond the latticediscretization. Preliminary results for N τ P { , } havebeen reported in Refs. 14 and 15. Here we significantlyincrease statistics and add a third lattice spacing witha set of simulations on N τ “ lattices. In agreementwith earlier work, a shift of the deconfinement criticalpoint towards smaller bare quark masses is observed asthe lattice is made finer, Figure 1 (right). On the otherhand, on the finer lattices we observe a quantitative de-viation regarding the location of the critical point withrespect to the hopping expanded results [13]. It is inter-esting to observe that the direction of the cut-off effectin the bare parameter space is the same for the critical a r X i v : . [ h e p - l a t ] S e p a Ñ N τ Ñ 8
Crossover st order 1 st order Z Z m u,d “ m phys. u,d m u,d “ 8 Figure 1. Left: Columbia plot for Wilson fermions on coarse lattices. The red circle in the heavy mass region at N f “ indicates the critical point which we determine on N τ P { , , } lattices. Right: Qualitative behavior of the second-orderboundary point with increasing N τ , i.e., decreasing lattice spacing. deconfinement boundary in the heavy quark mass regimeas for the critical chiral boundary in the light quark massregime with both Wilson [16, 17] and staggered [18–22]discretizations. The shift of the chiral Z boundary im-plies an increase in the simulation cost while the contin-uum is approached, which is particularly drastic in thechiral transition region. This further motivates to at-tempt a continuum limit in the heavy quark mass regimefirst, where it should be more feasible.After devoting Section II to the description of our lat-tice setup and of the relevant symmetry at work in theinfinite mass limit, we discuss our finite size scaling anal-ysis in Section III. Details on the simulations and theanalysis strategy are provided in Section IV, followed bya critical appraisal of the growing statistics requirementfor decreasing lattice spacings in Section V. Finally, ourresults for the deconfinement critical point are reportedin Section VI and conclusions are drawn in Section VII. II. LATTICE ACTION AND CENTERSYMMETRY
We work with the standard Wilson gauge action S g “ β (cid:88) n (cid:88) µ ď ν Re (cid:110) Tr (cid:2) ´ P µν p n q (cid:3)(cid:111) , (1)with the lattice coupling β “ N c { g , the plaquette P µν p n q and n labeling the lattice sites. The standardWilson fermion action for N f mass-degenerate quarks isdefined as S f “ a N f (cid:88) f “ (cid:88) n ,n ¯ ψ f p n q D p n | n q ψ f p n q (2) with the fermion matrix D p n | n q “ δ n ,n ´ κ ˘ (cid:88) µ “˘ (cid:2) p ´ γ µ q U µ p n q δ j ` ˆ µ,n (cid:3) , (3)where γ ´ µ ” ´ γ µ . The bare fermion mass m is adjustedvia the hopping parameter κ “ p am ` q . (4)The lattice coupling β controls the lattice spacing a p β q ,and temperature is defined as T “ a p β q N τ . (5)We do not use any improvements on the Wilson dis-cretization to make sure the phase structure does not getmodified by any unknown effects of improvement terms.The Polyakov loop L p n q is defined by the product ofall temporal links at space point n , closing through theperiodic boundary, L p n q “
13 Tr C (cid:34) N τ ´ (cid:89) n “ U p n , n q (cid:35) . (6)Physically, it describes the (euclidean) time evolution ofa static quark.At finite temperature, the periodic boundary inthe temporal direction permits topologically non-trivialgauge transformations, which are twisted by a global cen-ter element of the group when crossing the boundary, g p τ ` N τ , n q “ z g p τ, n q , g P SU p q , z “ e i πk , (7)with k P { , , } . The pure gauge action is invariantunder such transformations, S G r U g s “ S G r U s , while thePolyakov loop picks up the twist factor, L g p n q “ z ´ L p n q . (8)In the pure gauge theory and in the thermodynamic limit,the expectation value (cid:104) L (cid:105) is therefore an order parame-ter for center symmetry breaking, with non-zero valuessignalling deconfinement [23]. In the presence of dynam-ical quarks, S F r U g s ‰ S F r U s and the center symmetry isexplicitly broken, i.e., (cid:104) L (cid:105) ‰ always, with { m playingthe role of the symmetry-breaking field. Nevertheless, arapid change accompanied by large fluctuations of (cid:104) L (cid:105) still signals the deconfinement transition, which weak-ens with decreasing quark mass until a critical point isreached.Being the endpoint of a first-order transition, this crit-ical point will be in the same universality class as 3dliquid-gas transitions, i.e., the 3d-Ising model. This isfully confirmed by our data. In the effective 3d-IsingHamiltonian governing the vicinity of the critical point,the Polyakov loop will then be the dominant contribu-tion to the magnetization-like variable, coupling to thesymmetry-breaking magnetic field { m . It is interestingto contrast this situation with the chiral transition in thelight quark regime, where the Polyakov loop contributesto the energy-like variable in the corresponding effectiveHamiltonian, since it is invariant under chiral transfor-mations [24]. III. FINITE SIZE SCALING ANALYSIS
To determine the order of the phase transition as afunction of the quark mass, we compute standardizedmoments of X ” | L | , the absolute value of the Polyakovloop, B n p X ; α , α , . . . q “ (cid:10) p X ´ (cid:104) X (cid:105) q n (cid:11)(cid:10) p X ´ (cid:104) X (cid:105) q (cid:11) n , (9)and study their behavior as a function of the physicalvolume. Here n P { , } and { α n } is a set of relevantphysical parameters, in our case β, κ, N τ , N s . (In the fol-lowing, the dependency of B n on X will be understood).The phase boundary at a pseudo-critical coupling β c is defined by the vanishing of the skewness B p β c q “ .We use finite size scaling of the fourth moment, the kur-tosis B , trivially linked to the Binder cumulant [25], tolocate the critical point κ Z . For V Ñ 8 , B evaluatedon the phase boundary takes the form of a step function,assuming values that are characteristic for the respectivenature of the phase transition (c.f. Table I). On finite vol-umes, B is an analytic curve that gradually approachesthe step function with increasing volume. In the vicinityof the critical point κ Z , universality implies it to be afunction of the scaling variable x ” p κ ´ κ Z q N { ν s only, Crossover st order nd order Z B . ν - { . p q γ - 1 . p q α - - . p q Table I. Values of the kurtosis at the transition and of therelevant critical exponents [26]. which can be expanded around κ “ κ Z in a Taylor series.For sufficiently large volumes, the series can be truncatedafter the linear term, B p β c , κ, N s q “ B p β c , κ Z ,
8q ` a x ` O p x q , (10)and fitted to extract the critical parameters.This simple picture holds for asymptotically large vol-umes, and previous studies show that these can be pro-hibitively expensive to attain, cf. [27, 28]. Moreover, wefind the required aspect ratios in the heavy mass regionto be even significantly larger than those in our lightquark studies. The reason cannot be related to the light-est state in the spectrum, the glueball in this case, whosecorrelation length is shorter than that associated with thepion in the light quark mass regime. Instead, a possibleexplanation is that regular (non-divergent) contributionsto the fluctuations grow towards the heavy mass region,so that it takes larger volumes for the diverging terms todominate.It is therefore expedient to also consider the leading fi-nite volume corrections to Eq. (10). Since for finite quarkmasses the center symmetry is explicitly broken, thePolyakov loop is no longer a true order parameter, buta mixture of the magnetization- and energy-like observ-ables entering the effective 3d-Ising Hamiltonian whichgoverns the vicinity of the critical point. Magnetization-and energy-like observables scale with different expo-nents, and only for asymptotically large volumes thelarger of them dominates to produce the simple expres-sion Eq. (10). With the leading mixing correction takeninto account, the kurtosis takes instead the form B p β c , κ, N s q “ (cid:2) B p β c , κ Z ,
8q ` a x ` O p x q (cid:3) ˆ (cid:2) ` B N p α ´ γ q{ ν s (cid:3) , (11)where B is another non-universal parameter to be ex-tracted from a fit. For a detailed derivation and expla-nation in the context of the light quark mass regime,see [17]. As a consequence of this correction, the B curves are no longer fitted by straight lines, and differentvolumes intersect at different points, with the pairwiseintersection approaching the universal value with grow-ing volumes. This is shown schematically in Figure 2. IV. SIMULATION AND ANALYSIS DETAILS
We study cut-off effects by simulating three differenttemporal lattice extents N τ P { , , } , corresponding ´ . ´ . ´ . . . . B p β c , κ Z , x B L L ą L L ą L L ą L L ą L Increasinglatticesize ´ . ´ . ´ . . . . B p β c , κ Z , x B L L ą L L ą L L ą L L ą L Increasinglatticesize
Figure 2. Qualitative behavior of the kurtosis B for spatial volumes too small to be described by Eq. (10). The effect of theleading finite size correction, Eq. (11), is to shift B at κ Z (i.e. for x “ ) to larger values and approach the universal valuewith growing volumes. The enlargement shows that different volumes do not cross in one point. Pairwise intersections convergeto B p β c , κ Z , in the thermodynamic limit. − − . . n σ = 1 . n σ = 1 . n σ = 0 . B .
84 5 .
842 5 .
844 5 .
846 5 .
848 5 . . . . n σ = 1 . n σ = 3 . n σ = 1 . β B .
84 5 .
842 5 .
844 5 .
846 5 .
848 5 . β Raw data from multiple chains Raw (merged) and reweighted data
ReweightedRaw merged B ( β c ) ≡ B ( β c ) β c Figure 3. Example of the data analysis for B , B at κ “ . on a ˆ lattice. In the left plots different Markov chainsare displayed, with data points slightly shifted horizontally for better visibility. The label n σ denotes by how many standarddeviations the maximally incompatible pair is apart. On the right, the merged raw and reweighted data for B (top) and B (bottom) are shown. The determination of β c and B p β c q is depicted in red. to different lattice spacings at the respective critical cou-plings. For every value of N τ , the critical κ Z in theheavy quark mass region was determined (i.e. the posi-tion of red circle in Figure 1) for to values of κ . Ateach value of κ up to four spatial volumes have been sim-ulated, corresponding to aspect ratios of N s { N τ P r , s .For N τ “ , also simulations at N s “ have been done to give better insights into the position of κ Z . The de-confinement transition has been located simulating at to different values of β around the pseudo critical cou-pling. Configurations have been produced using a stan-dard HMC algorithm [29] with unit trajectories and theacceptance rate tuned to stay between 75% and 85%.At least five thousand thermalization trajectories havealways been discarded. Observables (the plaquette andthe Polyakov loop) have then been measured on the flyfor all Monte Carlo steps. At each value of β , between56k and 800k trajectories have been accumulated, mak-ing sure to always have over « independent eventsper β . In total , . and . millions trajectorieshave been produced for N τ “ , N τ “ and N τ “ ,respectively. Details about simulation parameters andstatistics are listed in Tables IV to VI.All the finite temperature simulations, except fromthose on the ˆ lattices, have been performed using v1.0 of our publicly available [30] OpenCL based CL QCDcode [31], which is optimized to run on GPUs. TheL-CSC supercomputer [32] at GSI in Darmstadt has beenused for this set of runs. The few simulations on aspectratio lattices have been performed with openQCD-FASTSUM [33], a highly MPI-optimized software thathas been run on the Goethe-HLR supercomputer. In allcases, monitoring and handling of thousands of jobs hasbeen efficiently automatized using the software package BaHaMAS [34, 35], whose new version v0.3.1 can also beused in conjunction with the openQCD-FASTSUM code.In order to faster accumulate statistics and to gainbetter control over statistical errors, for each parameterset we simulated four different Markov chains, except atthe outermost β values for aspect ratio where onlythree Markov chains were run. In this way a criterion todecide when the gathered statistics is large enough canbe established: All replicas included in the final analysiswere run until B was compatible within at most stan-dard deviations between all of them. An example of onedataset can be found in Figure 3 (left). At the beginningof the simulations, for very low statistics, this criterioncan be trivially fulfilled due to very large statistical er-rors. Therefore, to ensure a proper bracketing of β c andto avoid to stop the simulations too early, it has been al-most always demanded that B is incompatible with zeroat one standard deviation at the smallest and largest β on every chain. Once both these guidelines are satisfied,the chains were merged to the raw data points shown inFigure 3 (right).In order to have precision on β c , the multi-histogrammethod [36], also known as reweigthing, was used to in-terpolate between our measurements. This interpolationis repeated for the B data, thus pinning down B p β c q .Carrying out the same procedure for each simulated valueof κ , the resulting B p β c , κ, N s q data points can be plot-ted as a function of κ and fitted according to Eqs. (10)and (11), as it will be discussed in Section VI.Finally, to set the physical scale for our simulationpoints we use the Wilson flow parameter ω { a , whichwe determined and converted using the publicly avail-able software described in Ref. 37. To this end, we pro-duced independent zero-temperature configurationson ˆ lattices for each κ at the value of the crit-ical coupling β c . As a physical quantity to parametrizethe determined critical point κ Z , we then computed thepseudo-scalar meson mass m π corresponding to those .
11 0 .
115 0 .
12 0 .
125 0 .
13 0 . κ N s “ N s “ N s “ N s “ N s “ κ Z τ i n t N τ “ (a) Fixed N τ . ´ . ´ . ´ .
005 0 0 .
005 0 .
01 0 .
015 0 .
02 0 . N s N τ “ κ ´ κ Z τ i n t N τ “ N τ “ N τ “ (b) Fixed aspect ratio. Figure 4. The integrated autocorrelation time τ int of the skew-ness of the order parameter is shown for the simulated β clos-est to β c for different values of κ . In (a) N τ is kept fixed andthe spatial volume is varied ( κ Z is marked by the dashedline). In (b) τ int is plotted against κ ´ κ Z for different N τ atfixed aspect ratio. bare parameter values. Note that for all lattice spacingsconsidered here, am π ą at and around κ Z . Therefore,our pion mass estimates are afflicted by large cut-off ef-fects. This problem naturally reduces as N τ is increased.All critical couplings β c , the corresponding lattice spac-ings a , m π and critical temperatures T c are summarizedin Table III for each value of κ and N τ . V. STATISTICS REQUIREMENTS TOWARDSTHE CONTINUUM
A crucial parameter to judge the statistical quality ofthe analysed ensembles is the integrated autocorrelationtime τ int of the observables. Qualitatively speaking, τ int parametrizes the memory that the simulated system hasof its dynamics, which in our case strongly depends onthe order of the phase transition. Consider a first-ordertransition to start with and let us briefly summarize whathappens to the behaviour of the order parameter. In afinite volume, far away from the phase transition, thesystem stays in a given phase. The order parameter willfluctuate around its mean value and its probability dis-tribution will be approximately Gaussian, with some as-sociated τ int . As soon as β gets sufficiently close to β c ,the system will explore both phases and the order pa-rameter will jump from time to time from fluctuatingaround its mean value in one phase to fluctuating aroundits mean value in the other one. These fluctuations be-tween the phases are slower than the fluctuations withinone phase, since tunneling between different phases isexponentially suppressed by a potential barrier growingwith volume [38]. As a result, the system has a muchlonger-term memory since its dynamics now occurs ona larger timescale. τ int is related to the average tun-neling rate between the two phases and thus increasessignificantly. This in turn implies that sufficiently manytunneling events in a simulation are needed to reliablyestimate τ int . At a crossover transition, instead, the dis-tribution of the order parameter does not show a two-peak structure for any β values, even around β c , whereonly its variance slightly increases. Therefore, τ int is ex-pected to be relatively small in this case. Finally, at asecond-order phase transition the critical slowing downphenomenon will play an important role and the diverg-ing correlation length of the system will be reflected in anincrease of τ int , too. In a finite volume this effect will beonly partially felt, though, because the correlation lengthcannot literally diverge.We estimate τ int for both B and B using a Pythonimplementation of the Γ -method [39], and distribute ourdata in bins of size τ int to remove autocorrelations. Ta-bles V and VI give an overview of τ int for the skewnessand the kurtosis (averaged over the simulated β values)and the number of statistically independent events ob-tained from the binning procedure. The qualitative ex-pectation outlined above is completely met in Figure 4,where for each p κ, N s q pair, τ int is displayed at the simu-lated β closest to β c . In Figure 4(a) the autocorrelationtime is shown as a function of κ and spatial volume for N τ “ . For each volume, a maximum is seen around κ Z , showing critical slowing down near a second-ordertransition (that maximum would be sharper if τ int hadbeen evaluated exactly at β c ). The more drastic effect,however, is the increase of τ int with increasing volume.In Figure 4(b) we compare different N τ values at fixedaspect ratio. As expected when approaching the contin-uum limit, we observe another increase of τ int around thecritical as well as the first order region. The statistics for N τ “ is effectively smaller compared to N τ P { , } ,which explains the larger error bars and possibly an over-estimate in the first order region.Note also, that the observed rise in the autocorrelation time feeds back into the practical organization of the sim-ulation. The choice of β -values to be simulated is an op-timization of having a sufficiently narrow spacing neededfor reweighting, and covering a large enough interval tobracket β c , with as few values as possible. This requiresmonitoring and frequently analysing running simulationsin order to optimally adjust the parameters. Since a givenstatistics which turned out to be sufficient on a coarserlattice, will have to be increased on a finer one, even thisprocess of parameter tuning requires more and more tra-jectories as the lattices get finer. Altogether, our studyof the autocorrelation time illustrates the crucial impor-tance and formidable difficulties of obtaining sufficientstatistics for a reliable determination of the order of thephase transition as the continuum is approached. VI. RESULTS AND DISCUSSION
We are now ready to extract our desired results fromfits of our data to Eqs. (10) and (11), measuring thefit quality by the reduced chi-square χ d.o.f and the Q -parameter. The latter amounts to the probability of get-ting a value of chi-square larger than χ d.o.f , assumingthat the data have Gaussian noise, and gives a measureof the quality of the fit (its optimal value is ). In theideal situation of negligible finite size effects, we expectto obtain a good fit to Eq. (10) and, simultaneously, agood fit to Eq. (11) with a consistent κ Z and B compat-ible with zero. The general strategy then is to comparethe fits performed with and without the correction termdescribed in Eq. (11), with the goal to isolate the lead-ing terms. For this purpose, it is in principle enough tosuccessively exclude points from the smallest volumes inthe fit to Eq. (11), until the coefficient B is compatiblewith zero and check for consistency of the fit parameters.However, we generally observe that the inclusion of allvolumes corresponding to the smallest κ -values, i.e., theones deepest in the first order region, significantly deteri-orates the fit quality, irrespective of the Ansatz. Becausetunneling between the phases is exponentially suppressedby volume in the first-order region, these points are in-creasingly difficult to determine accurately, as discussedin the previous section. Moreover, the entire B curveis asymmetric about its inflection point, because its twoasymptotes are not symmetrically displaced with respectto B p β c , κ Z , . As the entries in Table III illustrate,the masses associated with the smaller κ -values rise veryquickly and are thus quite far from the critical quarkmass in which we are interested. This means that fit-ting points corresponding to the smallest κ -values wouldprobably require higher order terms in the Taylor ex-pansion of the kurtosis, which we neglect in Eqs. (10)and (11). On the other hand linear fits over ranges in κ that only cover the crossover region can result in biasedestimates for κ Z (cf. the shift in κ Z comparing the two e. . and i. . fits at N τ “ represented in Table II).Consequently, we made sure, that for each N τ we have .
075 0 .
08 0 .
085 0 .
09 0 .
095 0 . .
105 0 . . . . . . . B p β c , κ Z , κ B N s “ N s “ N s “ κ Z (a) Fit e. . at N τ “ . .
11 0 .
115 0 .
12 0 .
125 0 .
13 0 . . . . . . . B p β c , κ Z , κ B κ Z N s “ N s “ N s “ N s “ N s “ (b) Fit e. . at N τ “ . .
115 0 .
12 0 .
125 0 .
13 0 .
135 0 . . . . . . . B p β c , κ Z , κ B N s “ N s “ N s “ κ Z (c) Fit i. . at N τ “ . Figure 5. Final B fits for N τ P { , , } . These refer to the coloured fields in Table II. Shaded points have not been includedin the fits. included at least one κ in the fits, which belongs to thefirst-order region.Distinctive features for κ values in the first-order re-gion are the reversed N s ordering of the correspondingkurtosis (central) values with respect to what happens inthe crossover region (the value of the kurtosis decreaseswith N s in the first order region), as well as a more pro-nounced two-peak structure in the distribution of the or-der parameter. For N τ “ and κ “ . , a parameterset in doubt, we simulated aspect ratio specifically tocheck these features. As can be seen in Figure 5(b), nohints for a first order phase transition were found whichis consistent with the fit indicating κ Z ă . . In ourgeneral strategy it is important that the fit gives a κ Z larger than the smallest simulated κ value, and that the κ Z extracted from the fit is cross-checked against evi-dence for a first order phase transition for all simulated κ ă κ Z . Some representative fits are shown in Table II, wherethe best ones chosen as our final result are highlighted.Note that for N τ “ , we are able to find good fitswith B “ , as well as several consistent ones includingadditional data points and B ‰ . This indicates that allsimulation data are indeed described by Eq. (11) and thecoefficient B is fully controlled. Our data for B togetherwith the best fits highlighted in the table are also showngraphically in Figure 5.Our final results for the critical hopping parameter κ Z as a function of lattice spacing are collected in Figure 6(top). Strong cut-off effects are apparent. To comparewith other approaches and get a feeling for the physi-cal scales involved, Figure 6 (bottom) shows the criticalcouplings converted to pseudoscalar meson masses. Asalready indicated in Section IV, these numbers have to betaken with care because am π ą in all cases. This prob-lem could in principle be circumvented by heavy quark Fit label Included/Excluded κ per N s κ Z a B χ
NDF Q r % s e. . . p q . p q ´ .
87 9 55 i. . . p q . p q . p q .
98 8 45 .
075 0 .
085 0 .
09 0 . . e. . . p q . p q ´ .
49 7 85 i. . . p q . p q . p q .
54 6 78 .
075 0 .
085 0 .
09 0 . . e. . . p q . p q ´ .
05 13 39 i. . . p q . p q . p q .
14 12 32 .
11 0 .
115 0 .
12 0 .
125 0 .
13 0 . e. . . p q . p q ´ .
88 8 53 i. . . p q . p q ´ . p . q .
99 7 43 .
11 0 .
115 0 .
12 0 .
125 0 .
13 0 . e. . . p q . p q ´ .
36 9 20 i. . . p q . p q . p . q .
25 8 26 .
115 0 .
12 0 .
125 0 .
13 0 .
135 0 . e. . . p q . p q ´ .
59 7 77 i. . . p q . p q . p . q .
41 6 87 .
115 0 .
12 0 .
125 0 .
13 0 .
135 0 . Table II. An overview of the outcome of the final fit analysis is presented here. For each value of N τ the effect of excludingsome data points is shown. The fits are labeled by x.y.z , where x refer to the fit ansatz, y is the value of N τ and z simply acounter. For x “ e the fit has been done according to Eq. (10) ( excluding the correction term, B “ ” ´ ” ), while for x “ i Eq. (11) has been used and the correction term has been included . z “ always represents the fit with the least excluded datapoints, for z ą , more and more data points get excluded. The rows with blue background contain the best fit as a compromisebetween all parameters. Sub-tables in the second column show which κ (columns) have been included ( (cid:51) ) or excluded ( (cid:55) ) forwhich N s (rows). Gray cells denote simulations that were not done. κ Z N τ = 6 N τ = 8 N τ = 10 .
07 0 .
08 0 .
09 0 . .
11 0 . a { fm } m π { M e V } Figure 6. Above: The κ Z for three different lattice spacings.Below: Pion masses corresponding to the critical parameters. effective theory methods [40]. However, it is apparentthat at least two finer lattice spacings are needed beforea continuum extrapolation can be attempted. For thesethe mass in lattice units might well be small enough, sowe postpone this issue until such data are available. Nev-ertheless, the shift in the critical pion mass is significantlyreduced between N τ “ , compared to N τ “ , .A linear extrapolation using the last two points (unim-proved Wilson fermions have O p a q effects) would thenpredict m π « GeV, where twice the shift to the extrap-olated value should amount to a conservative estimate ofthe remaining systematic error of „ .It is now interesting to compare our results with thoseof [13], where the critical points for N f P { , , } Wil-son fermions were determined on N τ P { , } by meansof reweighting with a next-to-leading order hopping ex-panded fermion determinant and a histogram techniqueinstead of the cumulant analysis. For N f “ on N τ “ , N s “ , Ref. 13 reports κ Z “ . p q at NLO hop-ping expansion, which translates to m cπ { T c « . . Com-pared with our κ Z “ . p q or m cπ { T c « . , oneobserves a large discrepancy of «
50% in the criticalpion mass. Note that, since the same lattice action isemployed, this discrepancy is not related to the cut-off error on the determination of m cπ . As discussed in [13],the histogram method is applied at a fixed spatial volumeand does not include an extrapolation to the thermody-namic limit. Indeed, the authors report a reduction of κ Z by « when increasing volume to N s “ and by « on N s “ when going from LO to NLO in thehopping expansion, i.e., both truncations have a system-atic error in the same direction. The comparison to ourresult highlights the necessity for either refined approxi-mations or a full calculation in order to achieve a reliablecontinuum limit. VII. CONCLUSIONS
While qualitatively well understood, the heavy mass(top right) corner of the Columbia plot in Figure 1 (left)still lacks a quantitative determination, in the continuumlimit, of the location of the deconfinement critical bound-ary. In this work we focused our attention on the N f “ deconfinement critical point and studied its location onprogressively finer lattices simulating at N τ P { , , } .Our results show that the continuum limit is not yetwithin reach given the extent of the observed cut-off ef-fects. It is, indeed, apparent that for a continuum ex-trapolation to become feasible at least two finer latticespacings will be needed. Nevertheless, this work docu-ments the progress that has been made both in refiningthe fitting strategy for the finite size scaling analysis andin appraising the growing statistics requirements towardsthe continuum limit. Concerning the former, we showedhow a correction term can be used as a probe for finitesize effects, which allowed us to isolate the leading termsin the linear kurtosis expansion around κ Z and identifythe required aspect ratios for the linear regime. Withregards to the latter, our detailed analysis of the inte-grated autocorrelation times of the relevant observableillustrates the alarming prospects in terms of the statis-tics needed to reliably establish, as the continuum limitis approached, the order of the phase transition and thelocation of κ Z . Furthermore, our results allow for quan-titative comparisons with those obtained using hoppingexpansions and thus to assess their systematic error. ACKNOWLEDGMENTS
We thank C. Czaban for collaboration during the earlystages of this project [14, 15] and R. Kaiser for providinga very handy GUI for the fitting analysis. The authorsacknowledge support by the Deutsche Forschungsgemein-schaft (DFG) through the grant CRC-TR 211 “Strong-interaction matter under extreme conditions”. We alsothank the computing staff of the Goethe-HLR and L-CSCclusters for their support.0 [1] Y. Aoki, G. Endrödi, Z. Fodor, S. D. Katz, and K. K.Szabo, Nature , 675 (2006), arXiv:hep-lat/0611014[hep-lat].[2] C. Ratti, PoS
LATTICE2018 , 004 (2019).[3] O. Philipsen, in (2019) arXiv:1912.04827 [hep-lat].[4] M. Fromm, J. Langelage, S. Lottini, and O. Philipsen,JHEP , 042 (2012), arXiv:1111.4953 [hep-lat].[5] H. Saito, S. Ejiri, S. Aoki, K. Kanaya, Y. Nakagawa,H. Ohno, K. Okuno, and T. Umeda, Phys. Rev. D89 ,034507 (2014), arXiv:1309.2445 [hep-lat].[6] G. Aarts, F. Attanasio, B. Jäger, and D. Sexty, JHEP , 087 (2016), arXiv:1606.05561 [hep-lat].[7] C. S. Fischer, J. Luecker, and J. M. Pawlowski, Phys.Rev. D , 014024 (2015), arXiv:1409.8462 [hep-ph].[8] P. M. Lo, B. Friman, and K. Redlich, Phys. Rev. D ,074035 (2014), arXiv:1406.4050 [hep-ph].[9] G. Boyd, J. Engels, F. Karsch, E. Laermann, C. Lege-land, M. Lutgemeier, and B. Petersson, Nucl. Phys. B , 419 (1996), arXiv:hep-lat/9602007.[10] M. Shirogane, S. Ejiri, R. Iwami, K. Kanaya, andM. Kitazawa, Phys. Rev. D , 014506 (2016),arXiv:1605.02997 [hep-lat].[11] C. Alexandrou, A. Borici, A. Feo, P. de Forcrand,A. Galli, F. Jegerlehner, and T. Takaishi, Phys. Rev.D , 034504 (1999), arXiv:hep-lat/9811028.[12] H. Saito, S. Ejiri, S. Aoki, T. Hatsuda, K. Kanaya,Y. Maezawa, H. Ohno, and T. Umeda (WHOT-QCD), Phys. Rev. D84 , 054502 (2011), [Erratum: Phys.Rev.D85,079902(2012)], arXiv:1106.0974 [hep-lat].[13] S. Ejiri, S. Itagaki, R. Iwami, K. Kanaya, M. Ki-tazawa, A. Kiyohara, M. Shirogane, and T. Umeda(WHOT-QCD), Phys. Rev. D , 054505 (2020),arXiv:1912.10500 [hep-lat].[14] C. Czaban and O. Philipsen, PoS
LATTICE2016 , 056(2016), arXiv:1609.05745 [hep-lat].[15] F. Cuteri, C. Czaban, O. Philipsen, and A. Sciarra, EPJWeb Conf. , 07032 (2018), arXiv:1710.09304 [hep-lat].[16] O. Philipsen and C. Pinke, Phys. Rev.
D93 , 114507(2016), arXiv:1602.06129 [hep-lat].[17] X.-Y. Jin, Y. Kuramashi, Y. Nakamura, S. Takeda,and A. Ukawa, Phys. Rev. D , 034523 (2017),arXiv:1706.01178 [hep-lat].[18] F. Karsch, E. Laermann, and C. Schmidt, Phys. Lett. B520 , 41 (2001), arXiv:hep-lat/0107020 [hep-lat].[19] P. de Forcrand and O. Philipsen, Nucl. Phys.
B673 , 170(2003), arXiv:hep-lat/0307020 [hep-lat].[20] C. Bonati, P. de Forcrand, M. D’Elia, O. Philipsen,and F. Sanfilippo, Phys. Rev.
D90 , 074030 (2014),arXiv:1408.5086 [hep-lat].[21] F. Cuteri, O. Philipsen, and A. Sciarra, Phys. Rev. D , 114511 (2018), arXiv:1711.05658 [hep-lat].[22] F. Cuteri, O. Philipsen, and A. Sciarra, PoS LAT-TICE2018 , 170 (2018), arXiv:1811.03840 [hep-lat]. [23] L. D. McLerran and B. Svetitsky, Phys. Rev. D , 450(1981).[24] D. A. Clarke, O. Kaczmarek, F. Karsch, A. Lahiri, andM. Sarkar, (2020), arXiv:2008.11678 [hep-lat].[25] K. Binder, Zeitschrift für Physik B Condensed Matter , 119 (1981).[26] A. Pelissetto and E. Vicari, Phys. Rept. , 549 (2002),arXiv:cond-mat/0012164.[27] O. Philipsen and C. Pinke, Phys. Rev. D89 , 094504(2014), arXiv:1402.0838 [hep-lat].[28] C. Czaban, F. Cuteri, O. Philipsen, C. Pinke,and A. Sciarra, Phys. Rev.
D93 , 054507 (2016),arXiv:1512.07180 [hep-lat].[29] S. Duane, A. D. Kennedy, B. J. Pendleton, andD. Roweth, Phys. Lett.
B195 , 216 (1987).[30] M. Bach, C. Pinke, A. Sciarra, et al. , “CL QCD v1.0 ,” https://github.com/AG-Philipsen/cl2qcd (2018).[31] M. Bach, V. Lindenstruth, O. Philipsen, andC. Pinke, Comput. Phys. Commun. , 2042 (2013),arXiv:1209.5942 [hep-lat].[32] D. Rohr, M. Bach, G. Nešković, V. Lindenstruth,C. Pinke, and O. Philipsen, “Lattice-CSC: Optimiz-ing and Building an Efficient Supercomputer for Lattice-QCD and to Achieve First Place in Green500,” in HighPerformance Computing: 30th International Conference,ISC High Performance 2015, Frankfurt, Germany, July12-16, 2015, Proceedings , edited by J. M. Kunkel andT. Ludwig (Springer International Publishing, Cham,2015) pp. 179–196.[33] J. R. Glesaaen and B. Jäger, “openQCD-FASTSUM,” http://fastsum.gitlab.io/ (2018).[34] A. Sciarra,
Proceedings, 35th International Symposiumon Lattice Field Theory (Lattice 2017): Granada, Spain,June 18-24, 2017 , EPJ Web Conf. , 09003 (2018),arXiv:1710.08831 [hep-lat].[35] A. Sciarra et al. , “BaHaMAS: A Bash Han-dler to Monitor and Administrate Simulations,” https://github.com/AG-Philipsen/BaHaMAS (2020).[36] A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. , 1195 (1989).[37] S. Borsanyi et al. , JHEP , 010 (2012), arXiv:1203.4469[hep-lat].[38] J. Langer, Annals Phys. , 258 (1969).[39] U. Wolff (ALPHA), Comput. Phys. Commun. , 143 (2004), [Erratum: Comput. Phys. Com-mun.176,383(2007)], arXiv:hep-lat/0306017 [hep-lat].[40] R. Sommer, in Les Houches Summer School: Session 93:Modern perspectives in lattice QCD: Quantum field the-ory and high performance computing (2010) pp. 517–590,arXiv:1008.0710 [hep-lat]. Appendix: Tables with scale setting, statistics and integrated autocorrelation time N τ κ β c a m π a { fm } m π { GeV } T c { MeV } .
075 5 . . p q . p q . p q p q .
085 5 . . p q . p q . p q p q . . . p q . p q . p q p q .
09 5 . . p q . p q . p q p q . . . p q . p q . p q p q .
11 5 . . p q . p q . p q p q .
11 6 . . p q . p q . p q p q . . . p q . p q . p q p q .
115 6 . . p q . p q . p q p q .
12 6 . . p q . p q . p q p q .
125 5 . . p q . p q . p q p q .
13 5 . . p q . p q . p q p q .
135 5 . . p q . p q . p q p q .
115 6 . . p q . p q . p q p q .
12 6 . . p q . p q . p q p q . . . p q . p q . p q p q .
125 6 . . p q . p q . p q p q .
13 6 . . p q . p q . p q p q .
135 6 . . p q . p q . p q p q .
14 5 . . p q . p q . p q p q Table III. Outcome of the pion mass and the scale setting simulations, performed on ˆ lattices, accumulating 800independent configurations. The values of β c and the critical temperature T c of the deconfinment phase transition are alsoincluded to provide a more complete overview. a and m π have also been measured at the κ Z -values obtained from the fits andare reported on gray background. The value of β c at κ Z where to set the scale and measure the pion mass has been obtainedvia interpolation of β c -values at simulated κ nearby. N τ κ β c | Total statistics per N s | Number of simulated β valuesAspect ratio Aspect ratio Aspect ratio Aspect ratio . – 5.88884 | 1.6M | 2 5.88895 | 1.6M | 2 5.88933 | 1.6M | 2 . – 5.88407 | 1.6M | 2 5.88448 | 1.6M | 2 5.88452 | 1.6M | 2 . – 5.88097 | 2.4M | 3 5.88104 | 2.4M | 3 5.87985 | 2.4M | 3 . – 5.86865 | 1.6M | 2 5.86762 | 1.6M | 2 5.86758 | 1.6M | 2 . – 5.84677 | 1.6M | 2 5.84624 | 2.4M | 3 5.84623 | 2.4M | 38 . . . . . . . . . – 6.13558 | 1.8M | 3 6.13558 | 1.2M | 2 – . . – 6.05851 | 1.8M | 3 6.05758 | 1.6M | 4 – . N τ P { , , } . For N τ “ , N s “ we simulated at 3 different β values, β c “ . , andhave an overall statistics of 2.0M. N τ κ Average τ int p B q ¨ ´ | Average number of independent eventsAspect ratio Aspect ratio Aspect ratio Aspect ratio . – . p q | 437 . p . q | 222 p q | 95 . – . p . q | 397 p q | 187 p q | 104 . – . p . q | 381 p q | 243 . p . q | 351 . – . p q | 437 . p . q | 382 . p . q | 432 . – . p q | 790 . p q | 825 . p q | 8648 .
11 10 . p q | 419 . p . q | 239 p q | 139 p q | 45 .
115 12 . p . q | 343 . p . q | 236 p q | 146 – .
12 9 . p q | 429 . p . q | 271 . p . q | 158 – .
125 7 . p q | 550 . p q | 468 . p . q | 202 – .
13 5 . p q | 731 . p q | 683 . p q | 677 – .
135 3 . p q | 1293 . p q | 1273 – –10 .
115 12 . p . q | 258 – – – .
12 13 . p . q | 233 p q | 109 p q | 65 – . – . p . q | 147 p q | 87 – .
13 13 . p . q | 230 . p . q | 147 . p . q | 138 – . – . p q | 379 . p q | 213 – .
14 4 . p q | 765 . p q | 820 – –Table V. Average of the integrated autocorrelation time τ int and of the number of independent events for the skewness of theorder parameter. The average has been done among the merged chains at the simulated β . N τ κ Average τ int p B q ¨ ´ | Average number of independent eventsAspect ratio Aspect ratio Aspect ratio Aspect ratio . – . p q | 1313 . p . q | 292 p q | 131 . – . p q | 1293 . p . q | 323 p q | 139 . – . p . q | 792 p q | 308 . p q | 499 . – . p q | 857 . p q | 540 . p q | 524 . – . p q | 1399 . p q | 1216 . p q | 12678 .
11 5 . p q | 1049 . p q | 507 . p . q | 294 p q | 53 .
115 6 . p q | 845 . p . q | 362 . p . q | 331 – .
12 3 . p q | 1105 . p q | 592 . p q | 281 – .
125 3 . p q | 1116 . p q | 739 . p . q | 288 – .
13 3 . p q | 1508 . p q | 1128 . p q | 1102 – .
135 1 . p q | 2381 . p q | 1945 – –10 .
115 3 . p q | 892 – – – .
12 5 . p q | 576 p q | 215 p q | 103 – . – . p . q | 307 p q | 165 – .
13 4 . p q | 722 . p . q | 261 . p . q | 244 – . – . p q | 596 . p q | 336 – .
14 2 . p q | 1393 . p q | 1429 – –Table VI. Average of the integrated autocorrelation time τ int and of the number of independent events for the kurtosis of theorder parameter. The average has been done among the merged chains at the simulated ββ