Finite size and cut-off effects on the Roberge-Weiss transition in N f =2 QCD with Staggered fermions
FFinite size and cut-off effects on the Roberge-Weiss transition in N f “ QCD withStaggered fermions
Owe Philipsen ˚ Institut für Theoretische Physik - Goethe-Universität, GermanyMax-von-Laue-Str. 1, 60438 Frankfurt am Main andJohn von Neumann Institute for Computing (NIC) GSI, Planckstr. 1, 64291 Darmstadt, Germany
Alessandro Sciarra : Institut für Theoretische Physik - Goethe-Universität, GermanyMax-von-Laue-Str. 1, 60438 Frankfurt am Main (Dated: September 27, 2019)In the absence of a genuine solution to the sign problem, lattice studies at imaginary quarkchemical potential are an important tool to constrain the QCD phase diagram. We calculate thevalues of the tricritical quark masses in the Roberge-Weiss plane, µ “ ıπT { , which separate massregions with chiral and deconfinement phase transitions from the intermediate region, for QCD withN f “ unimproved staggered quarks on N τ “ lattices. A quantitative measure for the quality offinite size scaling plots is developed, which significantly reduces the subjective judgement requiredfor fitting. We observe that larger aspect ratios are necessary to unambiguously determine theorder of the transition than at µ “ . Comparing with previous results from N τ “ we find a „ % reduction in the light tricritical pion mass. The heavy tricritical pion mass stays roughlythe same, but is too heavy to be resolved on N τ “ lattices and thus equally afflicted with cut-offeffects. Further comparison with other discretizations suggests that current cut-off effects on thelight critical masses are likely to be larger than „ %, implying a drastic shrinking of the chiralfirst-order region to possibly zero. PACS numbers: 12.38.Gc, 05.70.Fh, 11.15.HaKeywords: QCD phase diagram
I. INTRODUCTION
The theoretical prediction of the QCD phase diagramas a function of temperature T and baryon chemical po-tential µ B has proved to be a difficult challenge for severaldecades. Because of the non-perturbative nature of thestrong interactions on hadronic scales, a first principlesapproach such as lattice QCD is required. On the otherhand, because of the severe sign problem of lattice QCDat finite µ B , standard Monte Carlo simulations are lim-ited to addressing small densities, µ B ă T , only [1, 2].Even at zero baryon density, there remain open ques-tions. While the thermal transition from a hadron gas toa quark gluon plasma is well established to be an analyticcrossover for physical quark masses [3], the universalityclass of the transition in the chiral limit of the u, d -quarksis still not settled since it cannot be simulated directly.For these reasons, it is useful to study the dependenceof the thermal transition on QCD parameters like quarkmasses, numbers of flavors, imaginary chemical potential,for which there is no sign problem, as well as on the latticespacing. The current knowledge of the nature of the QCDthermal transition as a function of the three light quarkmasses and imaginary chemical potential, as obtained oncoarse lattices with unimproved actions, is sketched in ˚ [email protected] : [email protected] Figure 1. For large and small quark masses, there areregions with first-order deconfinement and chiral phasetransitions, which in the infinite and zero mass limits areassociated with the breaking and restoration of the centerand chiral symmetries, respectively. These are separatedby surfaces of second order transitions from a large regionwhere the transition is merely an analytic crossover, towhich also QCD with physical parameters belongs [3, 4].Note that this qualitative picture is the same for unim-proved staggered [5, 6] and unimproved [7] as well asimproved [8] Wilson discretizations, whereas the preciselocation of the boundary at µ “ differs significantlybetween them, indicating large cut-off effects. These arealso observed for N f “ staggered fermions without root-ing [9]. By contrast, simulations with improved staggeredactions do not see any region of first-order chiral transi-tions within the available mass range, neither at zero [10]nor imaginary chemical potential [4, 11] thus providingupper bounds on the critical mass values.In the present work we continue earlier studies us-ing the unimproved staggered discretization at imaginarychemical potential on finer lattices. In particular, refer-ring to Figure 1, we investigate how the (red) tricriticalpoints on the N f “ line in the Roberge-Weiss-plane(bottom plane at p µ { T q “ ´p π { q ) move as the lat-tice spacing is reduced to „ { of its previous values.Together with similar investigations at µ “ , this estab-lishes the behavior of the critical surfaces when approach-ing the continuum. Such studies are complementary toones with improved actions, where no non-analytic chi- a r X i v : . [ h e p - l a t ] J a n T r i c r i t i c a l T r i c r i t i c a l Figure 1. Qualitative sketch of the three-dimensionalColumbia plot realized on coarse lattices. Whether the chiral,first-order triple region in the Roberge-Weiss plane shrinkson finer lattice enough to make an O p q -region appear in the m u,d “ plane remains unclear. ral transition is seen, and necessary, if all discretizationsare to be understood in the same manner, with expectedagreement in an eventual continuum limit. As a by-product of our study, we develop a new analysis of thefinite size scaling of cumulants, which significantly re-duces the amount of subjective judgement required forfitting.In order to render the paper self-contained, we brieflysummarize the main features of QCD at imaginary chem-ical potential in Section II. We then proceed to describeour numerical methodology in Section III and our novelanalysis method in Section IV. Our numerical results aregiven in Section V before we conclude in Section VI. II. QCD AT IMAGINARY CHEMICALPOTENTIAL
Because of charge conjugation symmetry and its ex-plicit breaking by a non-vanishing baryon density, theQCD partition function is an even function of quarkchemical potential, Z p µ q “ Z p´ µ q . For purely imagi-nary chemical potential, µ “ ıµ i , µ i P R , it is furthermoreperiodic [12], Z p µ i { T q “ Z p µ i { T ` πk { N c q , k “ , . . . N c ´ , (1)and we use N c “ colors for the QCD gauge group.These symmetries imply the phase structure shown inFigure 2, with three different Z p q center sectors, which End or meeting points
First order lines
Figure 2. QCD phase diagram in the T ´ ˆ µ i plane. The dashedline depicts the chiral/deconfinement transition whose naturedepends on the quark masses. The orange lines representthe Roberge-Weiss (RW) transitions. The black dots, wherethe first-order lines terminate, can be first-order triple points,tricritical points or second-order endpoints. are periodically repeated for higher µ i . Physical observ-ables, and in particular the thermodynamic functions,are invariant under a change of sectors, which are char-acterized by different phases of the Polyakov loop L p x q “
13 Tr N τ ´ (cid:89) τ “ U p τ, x q ” | L p x q| e ´ iϕ , (2)with (cid:104) ϕ (cid:105) “ kπ { , k P { , , } . At high temperatures,there are first-order phase transitions between the cen-ter sectors, whereas at low temperatures they are ana-lytically connected. The dotted line represents the an-alytic continuation of the thermal transition, whose or-der depends on the quark masses. For large and smallquark masses, these lines represent first-order deconfine-ment and chiral transitions, respectively, whereas for in-termediate quark masses they correspond to an analyti-cal crossover. Consequently, there are three possibilitiesfor the end-point of the Roberge-Weiss transition: forlarge and small quark mass it is a first-order triple point,where the thermal first-order transition lines meet thatof the center transition. For intermediate quark masses,the thermal transition is only a crossover and the cen-ter transition ends in a critical end-point in the D Isinguniversality class. At the boundaries between these situ-ations, corresponding to specific quark mass values, theend-point is tricritical and corresponds to the red bound-ary points in the Roberge-Weiss plane of Figure 1. Thepurpose of the present work is to locate these tricriticalmasses on N τ “ lattices with N f “ and compare theirvalues with previous determinations on a coarser N τ “ lattice [13], as well as with those of other discretizationschemes. Crossover st triple Tricritical D Ising B . . ν ´ { { . p q γ ´ . p q Table I. Critical values of ν , γ and B ” B p X, . . . q for someuniversality classes [14]. III. NUMERICAL SETUP
We consider the QCD partition function of N f “ mass-degenerate quarks with a purely imaginary chemi-cal potential. After integration over the fermionic fieldsit can be written as Z p T, ˆ µ i q “ (cid:90) D U (cid:0) det D r U, ˆ µ i s (cid:1) { e ´ S g r U s , (3)where S g is the gauge part of the action and D is thefermion matrix. For our investigation we used the stan-dard Wilson gauge action and the standard staggereddiscretization of dynamical fermions. Denoting the lat-tice gauge coupling by β “ { g , with the continuumgauge coupling g , and an elementary plaquette by P , wehave S g r U s “ β (cid:88) P (cid:110) ´ (cid:60) (cid:2) Tr C P (cid:3)(cid:111) . (4)The fermion matrix reads D i,j “ ˆ m u,d δ i,j `` (cid:88) ν “ η i,ν (cid:16) ˜ U i,ν δ i,j ´ ˆ ν ´ ˜ U : i ´ ˆ ν,ν δ i,j ` ˆ ν (cid:17) , (5)where ˆ m u,d “ am u,d is the quark bare mass in latticeunits, a is the lattice spacing, i, j refer to lattice sites, η i,ν are the staggered phases, ˆ ν is a unit vector on the latticeand ˜ U i,ν are the gauge links, which include the purelyimaginary chemical potential ˆ µ i “ aµ i in the temporaldirection, ˜ U i,ν “ (cid:40) U i,ν ν P { , , } e i ˆ µ i U i,ν ν “ . (6)The temperature is specified by the inverse euclideantime extent of the lattice, T “ a p β q N τ . (7)In order to locate a phase transition and study its na-ture, we calculate standardized cumulants B n p X, β, ˆ m u,d , ˆ µ i 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 , (8) constructed from an (exact or approximate) order para-meter X . In particular, a non-trivial zero of the skewnessof the X -distribution, B p β c q “ , (9)determines at which value of β “ β c a thermal transitiontakes place, while the value of the kurtosis B in thethermodynamic limit, evaluated at the critical coupling,will determine the order of the phase transition (referto Table I for common kurtosis values). Note that theso-called Binder cumulant [15], U p X q ” ´ (cid:10) p X ´ (cid:104) X (cid:105) q (cid:11) (cid:10) p X ´ (cid:104) X (cid:105) q (cid:11) “ ´ B , (10)is trivially related to the kurtosis (of the same observable)and contains the same information.We fix µ i { T “ π since in this case the imaginary partof the Polyakov loop is an exact order parameter, X “ L Im ” N s (cid:88) x Im p L p x qq . (11)Referring to Eq. (2), it distinguishes between the low T disordered phase and the high T ordered phase with two-state coexistence, (cid:40) (cid:104) ϕ (cid:105) “ ñ (cid:104) L Im (cid:105) “ low T (cid:104) ϕ (cid:105) ‰ ñ (cid:104) L Im (cid:105) ‰ high T , (12)with the advantage of knowing its mean value exactly.On a N τ “ lattice, we thus set ˆ µ i “ π { . Since thisis the boundary between two Roberge-Weiss sectors forall temperatures, B p L Im q “ for any value of β andwe cannot use Eq. (9) to locate the Roberge-Weiss end-point. However, the kurtosis B is expected to vary fromvalues close to (crossover) at low T to values close to (first order) at high T . Although it becomes a non-analytic step function in the V Ñ 8 limit, it is a smoothfunction on finite volumes, with the curves for differentvolumes crossing at a universal value for B at the criticalpoint β “ β c , provided that the spatial lattice extent islarge enough. This crossing provides the location of theRoberge-Weiss end-point. In the neighborhood of thecritical point β c , the kurtosis shows a well-defined finitesize scaling behavior as a function of the scaling variable x ” p β ´ β c q N { ν s . (13)Its Taylor expansion around the critical point x “ is B p β, x q “ B p β c ,
8q ` a x ` a x ` O p x q . (14)Sufficiently close to the thermodynamic limit, the coef-ficient B p β c , and the critical exponent ν take theiruniversal values depending on the type of transition.In order to locate the two tricritical points in theRoberge-Weiss plane, we performed simulations at dif-ferent values of ˆ m u,d and different values of β around .
435 5 .
445 5 .
455 5 .
465 5 . . . . . . Correct crossing point β B ( L I m ) N s = 12 N s = 18 N s = 24 (a) m “ . .
770 5 .
775 5 .
780 5 .
785 5 .
790 5 .
795 5 . . . . . . Correct crossing point β B ( L I m ) N s = 12 N s = 18 N s = 24 N s = 30 N s = 36 (b) m “ . Figure 3. Kurtosis of the imaginary part of the Polyakov loop as function of β at two different values of the quark mass. Theplot at m “ . is a typical example of what can happen when finite size effects are too large. Clearly, the data at N s “ have not been included in the finite size scaling analysis. the critical temperature. Evaluating the kurtosis in thecritical region and fitting it to Eq. (14), considering thelinear term only, gives B p β c , , a , β c and ν for everyvalue of ˆ m u,d . The change of ν as a function of ˆ m u,d thenpermits to locate the light and heavy ˆ m tric u,d values.Although our main quantitative analysis is based onthe kurtosis of the order parameter, we also calculatedthe susceptibility of | L | , χ p| L |q ” N s (cid:10) p| L | ´ (cid:104) | L | (cid:105) q (cid:11) . (15)which is expected to scale around β c according to χ “ N γ { ν s f p t N { ν s q . (16)Here t ” p T ´ T c q{ T c is the reduced temperature and f is a universal scaling function. Comparing the collapseplots obtained by fixing the critical exponents γ and ν to the first-order or second-order values, and by plotting χ { N γ { ν s evaluated on different lattice sizes against t N { ν s also provides information about the nature of the ther-mal transition and serves as a cross-check of the kurtosisanalysis. A similar cross-check using the chiral conden-sate (cid:104) ¯ ψψ (cid:105) was occasionally performed in the small-massregion, leading to consistent conclusions.We investigated values of ˆ m u,d in the intervals r . , . s and r . , . s . For each value of ˆ m u,d ,three to five spacial lattice sizes have been used, keep-ing N τ “ and ˆ µ i “ π { fixed. This corresponds toaspect ratios N s { N τ P r , s . Larger spatial volumes thaninitially chosen were added whenever the kurtosis of theorder parameter on different volumes was not crossingat the same point (an example is reported in Figure 3).For every lattice size, between three and seven values of β around the critical temperature have been simulated.In between those, the observables have been evaluated at additional β -values using the Ferrenberg-Swendsen mul-tiple histogram method (also known as reweighting tech-nique) [16] to increase resolution (see the discussion atthe end of this and in the next section).Configurations were generated by a standard RHMCalgorithm [17], producing four different Monte Carlochains per β with unit-length trajectories. Where ad-vantageous ( ˆ m u,d ă . ), the multiple-pseudofermionstechnique [18] has been used. The algorithm acceptancehas been tuned to be not lower than . At least k trajectories were always discarded as thermalization andafterwards observables of interests (i.e. the plaquette,the Polyakov loop and the chiral condensate for smallmasses) have been computed for every trajectory. Weincreased statistics until the standard deviation of thekurtosis B p X q decreased below „ . , and B p X q wasthe same on all the four chains at the same β within two(three) standard deviations in the large (small) mass re-gion. For this reason the collected statistics per β is notuniform, detailed information is given in the Appendix.In order to satisfy these strict requirements, millions oftrajectories per volume and almost half a billion in totalhave been produced. This large statistics is necessary dueto the large autocorrelation times shown by L Im , espe-cially when entering the first-order regions (cf. Table IIIin the Appendix). We always insisted on having at least independent events per β in the analysis.In order to determine the lattice spacing and the pionmass, also zero-temperature simulations have been per-formed. We produced independent configurations on ˆ lattices for each value of ˆ m u,d . Using the pub-licly available code described in Ref. 19, the scale wasset by the Wilson flow parameter w . Pion masses, in-stead, were measured with standard spectroscopy tech-niques [20, 21].All our numerical simulations (except those for scalesetting purposes) have been run using the publicly avail-able [22] OpenCL based code CL QCD [23], which is opti-mized for GPUs. The L-CSC [24] supercomputer at GSIin Darmstadt has been used, and the thousands of jobsneeded in the study have been efficiently handled usingthe simulation monitoring package
BaHaMAS [25].Our quite intricate fitting procedure used to extractthe critical exponent ν is completely analogous to theone previously described in Appendix B of Ref. 26. Foreach value of the quark mass, nearly all possible fits ofthe data to the linear part of Eq. (14) are performed, anda filtering procedure is applied afterwards in order to pickthe best fits. This is needed because the range in which B p L Im q can be considered linear for each N s value is notknown a priori. However, here we differ in a technical de-tail from the previous study in Ref. 26. As already men-tioned, reweighting was used not only to smoothen thesignal, but also to supply additional β -points for the fit.This approach allows to reduce the number of requiredsimulations, provided there are clear criteria accordingto which such points are added. Reweighted points in-troduce a correlation with the others, and too many ofthem would render the fits unreliable. Hence, we addedmore reweighted points between simulated points onlyif with lower resolution it was not possible to obtain agood fit. Moreover, another important aspect should beconsidered in choosing the reweighting resolution. The β -region where the kurtosis is linear shrinks on largervolumes. Thus, choosing the same resolution in β on dif-ferent N s would imply to include fewer points from largervolumes in the fit and, consequently, to enhance finite sizeeffects. Therefore, we increased the reweighting resolu-tion in β on larger N s , making the information comingfrom the smallest volume systematically less important,see Table V for detailed overview. This aspect is evenmore explicit looking at how many points per N s havebeen included in the fit. IV. AN ALTERNATIVE TO FITTING:QUANTITATIVE COLLAPSE PLOTS
Any fitting procedure, however careful, relies on a fewsubjective decisions like the number of reweighted pointsand the filtering parameters to judge a good fit. We nowpropose an alternative procedure to independently de-termine the critical exponent from scaling/collapse plots,which can then be compared with the results of the fittingprocedure.A collapse plot is obtained if an observable, which dis-plays universal finite size scaling, is plotted as a functionof its scaling variables, such that the curves for differentvolumes fall on top of each other provided the volumes aresufficiently large to represent the thermodynamic limit.There are several common observables used for this pur-pose, here we focus on B p L Im q as function of the scalingvariable x defined in Eq. (13). An example is shown inFigure 5. − − − −
10 0 10 20 301 . . . . . ν = 13 β − β c β c N /ν s B ( L I m ) N s = 12 N s = 18 N s = 24 (a) First-order critical exponents − . − . − . . . . . . . . ν = 0 . β − β c β c N /ν s B ( L I m ) N s = 12 N s = 18 N s = 24 (b) Second-order critical exponents Figure 4. Collapse plot of the kurtosis of the imaginary part ofthe Polyakov loop for ˆ m u,d “ . using known values of thecritical exponent ν . Since this mass is far from the tricriticalpoints, a by-eye judgement is enough to rule out the first orderas possible type of phase transition. Whenever the lattice volume is not large enough, scal-ing is violated and no good collapse is obtained, evenwhen the known critical values (listed in Table I) areused. Finite size corrections are responsible for thatand, in principle, a better collapse can be obtained usingdifferent (non-universal) values of the exponents. Thequality of the collapse is usually judged by eye, which ismostly sufficient to distinguish between a first and a sec-ond order phase transition with known exponents (likein Figure 4). However, a more rigorous method is clearlyneeded in a situation where the scaling exponents changefrom first to second-order. On any finite volume, thiswill lead to intermediate values of the exponents, whichshould be determined unambiguously together with anassociated error. For this purpose we now construct aquantitative measure of the collapse of our data. β c ∆ xf p x q “ a ` a xa “ ˘ a “ -0.000096 ˘ ν ∆ xf p x q “ a ` a xa “ ˘ a “ ˘ Figure 5. Example of linear extrapolation ∆ x Ñ in the quantitative collapse plot analysis for ˆ m u,d “ . for β c and ν . Considering how we judge a collapse plot by eye, themeasure of the quality has to be related to the distancebetween different points at the same value of the uni-versal scaling variable. Inspired by the method used byBarkema and Newman [27] for the thermal random-fieldIsing model (see also a similar analysis in Appendix Aof Ref. 28), we associate a quantitative quality to thecollapse of B p L Im q by estimating the average varianceof the data as Q p ¯ β c , ¯ ν q ” x (cid:90) x max x min (cid:40) N V N V (cid:88) i “ (cid:104) B (cid:0) x p ¯ β c , ¯ ν, V i q (cid:1)(cid:105) (17) ´ (cid:34) N V (cid:88) i “ B (cid:0) x p ¯ β c , ¯ ν, V i q (cid:1)(cid:35) (cid:41) d x . Here, ∆ x ” x max ´ x min is the considered range in thescaling variable, V ” N s , N V is the number of consideredlattice volumes, while ¯ β c and ¯ ν are fixed values for thecritical temperature and for the critical exponent ν , re-spectively. A normalisation factor N ´ V was neglected infront of the expression, since it is irrelevant for the esti-mate of the critical exponent. It is now possible to obtainan estimate for β c and ν by minimizing Q as function ofthese two variables. Nevertheless, this is a non-trivialtask and there are some problems to be addressed.The integration in Eq. (17) must be done numerically,since the exact functional form of B p x q is unknown.This, in principle, would not be a problem, if only we hadthe kurtosis at the same x on different volumes. How-ever, in lattice simulations the kurtosis B is measuredat fixed values of β , and the mapping (13) between β and x depends on the unknown parameters β c and ν . There-fore, it is not possible to have simulated data uniformlyspaced in x for any β c and ν . On the other hand, themeasured data can be interpolated in β using the multi-ple histogram method. Hence, it is possible to reweightthe kurtosis in β in such a way that its values at the same x are available for all the volumes. After this step thecalculation of Q is trivial. In practice, this implies aninterpolation for each pair p ¯ β c , ¯ ν q at which Q has to beevaluated, which is too costly if a precise determinationof the final value of the critical exponent is desired. A cheaper alternative is to use the reweighting techniqueto obtain the kurtosis as an approximately continuousfunction of β , i.e., to add a large number of points be-tween two simulated temperatures. The numerical inte-gration to obtain Q can then be performed with negligi-ble additional error. However, due to the particular formof the map x p β q , sometimes, especially for small valuesof ν , the number of interpolated points needed to havea sufficiently precise numerical integration can becomevery large and, therefore, the reweighting very costly. Asmarter approach is then required.As can be seen in Figure 3, the kurtosis of the imagi-nary part of the Polyakov loop is a quite regular func-tion of β , in the sense that no sudden variations arepresent. This means that a numerical interpolation ofthe data which does not take into account the physics –as the multiple histogram method does – will probablystill find the correct value of the kurtosis. Clearly, thisis true under the assumption that the resolution of thedata to be interpolated is high enough. For example,the simulated data are usually too distant in tempera-ture to be correctly interpolated without the reweightingtechnique. But after having applied the multiple his-togram method to the data, a second interpolation canbe done very cheaply. In practice, we used the software Mathematica to obtain an interpolated function out of aset of points and perform numeric operations on it. Theadvantage of having a kurtosis as a function makes thecalculation of Q straightforward. Furthermore, it is thenpossible to automatically minimize Q p β c , ν q as functionof two variables.Next, we need to estimate the statistical error on β c and on ν , which has to contain the error on thereweighted points and the statistical error of our sim-ulations. An error on reweighted points is often obtainedusing the bootstrap method. This means that, in thereweighting procedure, N boot sets of reweighted kurtosisvalues are calculated, and the bootstrap errors are ex-tracted from them. Now, instead of using these sets tocompute errors on the kurtosis, they can be used to min-imize Q , obtaining N boot different estimates of β c and of ν , which will give the desired final error. Since, typically,the number of bootstrap samples is of the order of somehundreds, it is clear that the minimization of Q shouldnot take too much time .Finally, let us discuss how x min and x max should bechosen. Clearly, no extrapolation outside the simulatedinterval in β should be done. Thus, the largest ∆ x is theinterval in x around where data from all volumes areavailable. Since x c “ , we have x min ă and x max ą and, in order to have a symmetric Taylor expansioninterval [26], we chose | x min | “ | x max | . (18)Using too large an interval of integration is, in general,not correct, since it assumes data collapse possibly out-side the critical region. On the other hand, the width ofthe scaling region is not known a priori. A clever solu-tion to this problem, successfully applied in [27], consistsof repeated analyses for successively decreasing ∆ x fol-lowed by an extrapolation of the resulting parameters to ∆ x Ñ . An example is reported in Figure 5. V. NUMERICAL RESULTS AND DISCUSSION
As mentioned in Section III, our strategy to locate thetricritcal points is to measure the critical exponent ν fordifferent quark masses and see where it changes from itsfirst-order value { for small and large masses to the D Ising value p q for intermediate masses. Thechanges approach a step function in the thermodynamiclimit but remain smooth as far as finite lattice volumesare used to extract ν . The critical exponent ν is prefer-able over B p β c , , since it is known to suffer less fromfinite volume corrections [26, 29, 30]. The main result ofour investigation is reported in Figure 6 (more detailedinformation about the displayed data is available in Ta-bles IV and V in the Appendix). For each value of thequark mass, the critical exponent ν is extracted, bothwith the fit analysis used already in Ref. 26 and withthe new collapse strategy introduced in Section IV. Theagreement between the two methods to extract the crit-ical exponent ν is evident in Figure 6. The quantitativecollapse analysis has systematically smaller errors on ν ,though. This, together with the fact that no arbitrarydecision in the analysis may affect the outcome, shouldmake this method preferable.In the large mass region, the signal is quite smoothand ν changes monotonically from the second-order tothe first-order value. However, approaching and enteringthe first-order region, the minimal aspect ratio N s { N τ needed to extract ν increases significantly from to compared to studies at µ “ . This has been remarkedalready in previous studies [13, 26]. Presumably this is In Mathematica , for example, it is possible to use the
NargMin function, but a user implemented minimization based on a scanin β c and in ν is more efficient, though less precise. due to the fact that in the Roberge-Weiss plane we aredealing with the more complex three-state coexistenceand its coalescence in a tricritical point.The same behavior is expected also in the small-massfirst-order region, where simulations with an aspect ratiolarger than are too costly. Therefore, in Figure 6, themass region ˆ m u,d ď . has been marked with a graybackground to stress that larger volumes are required topolish the result. However, we have reasons to believethat the tricritical point is already located. It is alwayspossible to compare the critical exponent ν extracted us-ing only part of the available volumes, while leaving thesmallest or largest out of the analysis. In this way finitesize effects are made visible by checking whether ν driftstowards first-order or second-order values upon inclusionof larger volumes. This is shown for ˆ m u,d “ . in Fig-ure 6 where a clear decrease in ν is visible adding N s “ and removing N s “ in the analysis. Another aspectthat made us confident to be entering the first-order re-gion for ˆ m u,d À . is the typical “Binder bump” be-havior discussed in detail in Ref. 26. At ˆ m u,d “ . the kurtosis of the order parameter starts to overshootthe value for β À β c , which is due to the coexistenceof three states and thus clearly signals entering the first-order region.A few points in Figure 6 were accepted to be obtainedfrom two spatial volumes only. For the two smallestquark masses which were simulated, it was clear fromthe crossing point of the kurtosis on N s P { , , } that N s “ was too far away from the thermodynamiclimit. On the other hand, to add a larger spatial extentwould have been very costly without the guarantee to besufficient for a conclusive statement. About ˆ m u,d “ . ,instead, we considered the outcome of the analysis withN s P { , } satisfactory, since the kurtosis of the orderparameter reaches values larger than . for β À β c andthe bump shrinks and get larger increasing N s , behaviortypical of the first-order region.After these considerations, our estimates of the tricrit-ical bare quark masses are ˆ m triclight “ . ` . ´ . ˆ m tricheavy “ . p q , (19)where the conservative choice of having an asymmetricerror in the chiral region is to stress that further investi-gation would be needed in the chiral limit to polish themeasurement.In order to asses how much the results are affected bycut-off effects, we measured both the lattice spacing a andthe pion mass m π for all simulated bare quark masses,by running T “ µ i “ simulations at the β c found inthe Roberge-Weiss plane. The outcome is reported inTable II. Having fixed the scale, it is possible to expressEq. (19) in terms of pion masses in physical units, m tric π, light “ ` ´ MeV m tric π, heavy “ ` ´ MeV . (20) .
003 0 .
005 0 .
007 0 .
009 0 . . . . . .
34 34 345 345 234 234 234 234234 ν From fit analysis .
20 0 .
40 0 .
60 0 .
234 234 234 345 3456 456 456 456 456 67 567 N s Nτ From quantitative collapse ˆ m u,d Tricritical nd order st triple Figure 6. Critical exponent ν as function of the bare quark mass ˆ m u,d obtained with two different methods. The quantitativecollapse plot points have been slightly horizontally shifted to avoid superposition (refer to Tables IV and V for the the data).At ˆ m u,d “ . the outcome of the data analysis without the largest spatial value (shaded points) has been included to guidethe discussion in the text. For each mass, a badge containing the aspect ratios N s { N τ used in the analysis has been drawn.The horizontal colored lines are the critical values of ν for some universality classes. The mass axis has been broken and twodifferent scales have been used in order to improve readability. Points in the grey-background region have to be taken with apinch of salt, since it is reasonable to believe that finite size effects are for them dominant. Our results in physical units are given in Figure 7, wherethe critical exponent ν obtained with our new analysisstrategy is plotted as a function of m π . It is importantto stress that in the large-mass region the lattices usedare still too coarse to correctly resolve the pion and wehave am π ą , implying sizeable cut-off effects on thisvalue.We now compare with the previous results obtained onN τ “ lattices [13]. However, there, only the tricriticalbare quark masses and a rough estimate of m tric π, light werereported. We therefore improve the determination of thelatter, by performing additional scale setting simulations,whose outcome can be found in Table II. In particular,we measured m π and a for three values of the quark baremass, corresponding to the light tricritical point quotedin Ref. 13 (the central value and at one standard devi-ation apart from it). The value of β has been chosenusing a polynomial interpolation of the β c obtained bythe authors at the simulated masses. Taking as error onthe tricritical pion mass the difference between its valueand m π resulting from the neighboring bare masses, weobtain m tric , N τ “ π, light “ ` ´ MeV . (21) We thus conclude that, in the light region, a shift ofaround % is found when moving from a N τ “ toa finer N τ “ lattice. In the heavy mass region, only arough comparison is possible, since no pion mass is re-ported in Ref. 13 and in any case on N τ “ the pion isresolved even less. However, it is possible to compare thedimensionless ratio ˆ m u,d { T at the tricritical point, m tric , N τ “ π, heavy T “ . p q m tric , N τ “ π, heavy T “ . p q , (22)which turn out to be compatible. VI. DISCUSSION AND CONCLUSIONS
Since we are still far from being able to perform a con-tinuum extrapolation, it is instructive to compare withother discretizations. The results obtained in Ref. 26with Wilson fermions on N τ “ lattices, i.e. with similar N τ ˆ m u,d β w { a a m π a { fm } m π { MeV } T { MeV }4 0.038 5.356 .
755 05 p q . p q . p q .
759 17 p q . p q . p q .
763 56 p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q . p q T “ simulations have been performed on ˆ lattices always collecting independent configurations. w { a has been determined and converted to physical scales using the publicly available codedescribed in [19]. For the pion mass determination, point sources per configuration have been used. The table also containsthe lattice spacing, the pion mass and the temperature of the corresponding finite temperature ensemble in physical units. m π { MeV } nd order Tricritical st triple m tric1 m tric2
240 280 320 360 4000 . . . . . ν c o ll a p s e Figure 7. Critical exponent ν as function of the pion mass m π . This plot is similar to that in Figure 6, but the values of thecritical exponent are only those obtained with the strategy described in Section IV. The regions with light-red backgrounddenote the position of the tricritical masses within one standard deviation. N f = 2 + 1 N f = 2 N f = 2 m π { MeV } N τ [30][26][13] and this studyThis study[4][11] Ref.Unimproved WilsonUnimprovedStoutHISQ S t agg e r e d Figure 8. Overview of chiral tricritical values of the pion mass in the Roberge-Weiss plane. lattice spacing, appear to have considerably larger cut-offeffects. For example, comparing am tric π, heavy “ . p q from Ref. 26 with our am tric π, heavy “ . p q , the pion-resolution problem is milder in the present study. It isalso interesting to compare the position of the tricriticalpoints in physical units, m tric, Wilson π, light “ ` ´ MeV m tric, Staggered π, light “ ` ´ MeV (23)and m tric,Wilson π, heavy “ ` ´ MeV m tric, Staggered π, heavy “ ` ´ MeV . (24)The large differences between discretizations again implybeing far from the continuum limit, where results fromall discretizations have to merge. The observed trendis consistent with the findings of simulations with im-proved staggered actions, where the tricritical points canonly be bounded to be at much smaller masses, as indi-cated in Figure 8, as well as with the analogous findingsat zero chemical potential (see discussion in the introduc-tion). In particular the comparison across discretizationsimplies enormous cut-off effects in the critical masses,which could end up being over „ τ “ to N τ “ lattices. The aspect ratios and statistics required to ex-tract the correct order of the phase transition are foundto be larger in the Roberge-Weiss plane than at µ “ .We find the cut-off effect on the tricritical masses to besmaller but qualitatively the same as that observed withWilson fermions, and consistent with results for both dis-cretizations at zero chemical potential. This implies inparticular, that the entire chiral critical surface depictedin Figure 1 is shifted significantly towards smaller (andpossibly zero) light quark masses, as the lattice spac-ing decreases, which is also consistent with results fromimproved staggered actions. Unfortunately, our studyalso implies that much finer lattices at inevitably smallerquark masses are necessary, before one can hope the re-sults of the light tricritical mass to stabilize in a contin-uum limit. ACKNOWLEDGMENTS
We thank Francesca Cuteri for useful discussions andinput for Figure 8. The authors acknowledge supportby the Deutsche Forschungsgemeinschaft (DFG) throughthe grant CRC-TR 211 “Strong-interaction matter un-der extreme conditions” and by the Helmholtz Interna-tional Center for FAIR within the LOEWE program ofthe State of Hesse. We also thank the computing staff ofL-CSC for their support.1 [1] O. Philipsen, (2010), arXiv:1009.4089 [hep-lat].[2] C. Ratti,
Proceedings, 36th International Symposium onLattice Field Theory (Lattice 2018): East Lansing, MI,United States, July 22-28, 2018 , PoS
LATTICE2018 ,004 (2019).[3] Y. Aoki, G. Endrodi, Z. Fodor, S. Katz, and K. Szabo,Nature , 675 (2006), arXiv:hep-lat/0611014 [hep-lat].[4] C. Bonati, E. Calore, M. D’Elia, M. Mesiti, F. Negro,F. Sanfilippo, S. F. Schifano, G. Silvi, and R. Tripiccione,Phys. Rev.
D99 , 014502 (2019), arXiv:1807.02106 [hep-lat].[5] C. Bonati, P. de Forcrand, M. D’Elia, O. Philipsen,and F. Sanfilippo, PoS
LATTICE2011 , 189 (2011),arXiv:1201.2769 [hep-lat].[6] F. Cuteri, O. Philipsen, and A. Sciarra, Phys. Rev.
D97 ,114511 (2018), arXiv:1711.05658 [hep-lat].[7] O. Philipsen and C. Pinke, Phys. Rev.
D93 , 114507(2016), arXiv:1602.06129 [hep-lat].[8] X.-Y. Jin, Y. Kuramashi, Y. Nakamura, S. Takeda,and A. Ukawa, Phys. Rev.
D96 , 034523 (2017),arXiv:1706.01178 [hep-lat].[9] P. de Forcrand and M. D’Elia,
Proceedings, 34th In-ternational Symposium on Lattice Field Theory (Lattice2016): Southampton, UK, July 24-30, 2016 , PoS
LAT-TICE2016 , 081 (2017), arXiv:1702.00330 [hep-lat].[10] A. Bazavov, H. T. Ding, P. Hegde, F. Karsch, E. Laer-mann, S. Mukherjee, P. Petreczky, and C. Schmidt,Phys. Rev.
D95 , 074505 (2017), arXiv:1701.03548 [hep-lat].[11] J. Goswami, F. Karsch, A. Lahiri, and C. Schmidt,
Pro-ceedings, 36th International Symposium on Lattice FieldTheory (Lattice 2018): East Lansing, MI, United States,July 22-28, 2018 , PoS
LATTICE2018 , 162 (2018),arXiv:1811.02494 [hep-lat].[12] A. Roberge and N. Weiss, Nucl.Phys.
B275 , 734 (1986).[13] C. Bonati, G. Cossu, M. D’Elia, and F. Sanfilippo,Phys.Rev.
D83 , 054505 (2011), arXiv:1011.4515 [hep-lat].[14] A. Pelissetto and E. Vicari, Phys.Rept. , 549 (2002),arXiv:cond-mat/0012164 [cond-mat]. [15] K. Binder, Z.Phys.
B43 , 119 (1981).[16] A. M. Ferrenberg and R. H. Swendsen, Phys.Rev.Lett. , 1195 (1989).[17] A. D. Kennedy, I. Horvath, and S. Sint, Lattice FieldTheory. Proceedings: 16th International Symposium,Lattice ’98, Boulder, USA, Jul 13-18, 1998 , Nucl. Phys.Proc. Suppl. , 834 (1999), [,834(1998)], arXiv:hep-lat/9809092 [hep-lat].[18] M. A. Clark and A. D. Kennedy, Phys. Rev. Lett. ,051601 (2007), arXiv:hep-lat/0608015 [hep-lat].[19] S. Borsanyi et al. , JHEP , 010 (2012), arXiv:1203.4469[hep-lat].[20] M. F. L. Golterman, Nucl. Phys. B273 , 663 (1986).[21] N. Ishizuka, M. Fukugita, H. Mino, M. Okawa, andA. Ukawa, Nucl. Phys.
B411 , 875 (1994).[22] M. Bach, C. Pinke, A. Sciarra, et al. , “CL QCD v1.0 ,” https://github.com/CL2QCD/cl2qcd .[23] O. Philipsen, C. Pinke, A. Sciarra, and M. Bach, PoS LAT2014 , 038 (2014), arXiv:1411.5219 [hep-lat].[24] D. Rohr, M. Bach, G. Neskovic, V. Lindenstruth,C. Pinke, and O. Philipsen, in
High Performance Com-puting (LNCS) , Vol. 9137 (2015).[25] 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].[26] C. Czaban, F. Cuteri, O. Philipsen, C. Pinke,and A. Sciarra, Phys. Rev.
D93 , 054507 (2016),arXiv:1512.07180 [hep-lat].[27] M. E. J. Newman and G. T. Barkema, Phys. Rev. E ,393 (1996).[28] J. Houdayer and A. K. Hartmann, Phys. Rev. B ,014418 (2004).[29] P. de Forcrand and O. Philipsen, Phys.Rev.Lett. ,152001 (2010), arXiv:1004.3144 [hep-lat].[30] O. Philipsen and C. Pinke, Phys. Rev. D89 , 094504(2014), arXiv:1402.0838 [hep-lat].[31] U. Wolff (ALPHA), Comput. Phys. Commun. , 143 (2004), [Erratum: Comput. Phys. Com-mun.176,383(2007)], arXiv:hep-lat/0306017 [hep-lat].
Appendix: Simulation details
It is known that different powers of the same observable have different integrated autocorrelation times τ int , whichcan be estimated using the Wolff algorithm [31]. This is important to be taken into account when it comes to measurestandardized cumulants, like the kurtosis of a given observable. Binning, i.e. substituting a block of data with itsaverage, allows to obtain uncorrelated data from the correlated ones. This is true if the size of a block is at leasttwice τ int , though. It is then possible to understand how many independent measurement of the quantity of interestare available in a Monte Carlo simulation, just by dividing the number of the trajectories produced by τ int . Clearly,the larger this number is the more accurate the result will be. However, simulations in full QCD are costly and acompromise is needed. We always had at least independent events for B p L Im q in the merged chain obtained byputting together the four independent Markov chain that we produced for each β value. A detailed overview of thecollected statistics is presented in Table III.Tables IV and V contain, instead, the detailed outcome of our analysis, whose data were plotted in Figure 6.2 ˆ m u,d β range Total statistics per spatial lattice size N s (cid:0) of simulated β | ¯ τ int | n eventsmin per chain (cid:1)
12 36 18 42 24 300.004 5.425 - 5.437 1.56M (5 | 88 | 301) 1.56M (4 | 247 | 118) 1.32M (4 | 339 | 64) -0.005 5.427 - 5.442 1.00M (5 | 86 | 225) 3.04M (5 | 199 | 153) 2.12M (5 | 310 | 63) -0.006 5.430 - 5.445 2.28M (6 | 91 | 363) 1.92M (6 | 212 | 84) 1.08M (4 | 410 | 61) 1.00M (4 | 527 | 23)0.007 5.420 - 5.460 0.64M (5 | 69 | 174) 1.72M (5 | 174 | 135) 1.10M (5 | 256 | 63) 1.56M (4 | 420 | 39)0.008 5.430 - 5.470 0.92M (5 | 66 | 217) 0.86M (5 | 175 | 68) 1.46M (4 | 331 | 57) -0.009 5.430 - 5.470 1.26M (5 | 72 | 266) 1.28M (5 | 182 | 91) 1.60M (4 | 325 | 66) -0.010 5.430 - 5.480 1.20M (6 | 65 | 297) 0.60M (4 | 143 | 71) 1.84M (5 | 263 | 48) -0.011 5.430 - 5.490 0.68M (4 | 57 | 229) 1.24M (5 | 180 | 82) 1.92M (4 | 336 | 72) -0.150 5.590 - 5.720 1.08M (5 | 62 | 275) 5.80M (7 | 242 | 203) 5.28M (7 | 362 | 153) -0.250 5.600 - 5.760 4.20M (7 | 81 | 555) 2.00M (4 | 190 | 167) 4.60M (6 | 409 | 85) -0.350 5.720 - 5.780 6.44M (7 | 130 | 510) 3.40M (5 | 279 | 131) 5.88M (6 | 442 | 130) -0.400 5.750 - 5.790 - 9.00M (5 | 305 | 214) 10.60M (6 | 574 | 116) 7.80M (5 | 917 | 95)0.450 5.760 - 5.810 18.40M (4 | 1343 | 287) 4.72M (5 | 330 | 189) 20.00M (5 | 636 | 511) 20.00M (5 | 1010 | 279)0.500 5.780 - 5.820 10.00M (4 | 1507 | 107) 9.40M (5 | 325 | 255) 10.00M (5 | 598 | 112) 11.20M (5 | 889 | 236)0.550 5.760 - 5.840 6.20M (5 | 1355 | 38) 10.80M (6 | 300 | 173) 6.56M (5 | 606 | 77) 7.80M (5 | 917 | 57)0.600 5.812 - 5.827 12.00M (4 | 1720 | 152) - 19.80M (6 | 766 | 404) 14.00M (6 | 1207 | 94)0.650 5.817 - 5.837 15.40M (4 | 1759 | 74) 16.40M (5 | 501 | 673) 20.40M (5 | 780 | 295) 15.40M (5 | 1301 | 192)0.750 5.830 - 5.846 15.20M (4 | 1976 | 79) 5.60M (3 | 2507 | 48) - 16.80M (5 | 1240 | 93)0.850 5.840 - 5.856 19.60M (4 | 1792 | 111) 18.00M (5 | 2066 | 71) - 16.40M (4 | 1288 | 136)
Table III. Overview of the statistics accumulated in all the simulations (red entries are preliminary). Since the resolution in β is not the same at different ˆ m u,d , the number of simulated β has been reported per each range. The accumulated statistics per β varies because of the criterion adopted to stop to increase the statistics on the chains. Therefore we reported here the totalnumber of trajectories produced per given N s . For each N s , the number of simulated β , the average integrated autocorrelationtime and the smallest number of independent events per chain of B p L Im q can be found in the brackets next to the totalstatistics. Observe that ¯ τ int and n eventsmin are not connected. The former is an average among all the different chains run atone fixed spatial lattice extent, while the latter is the effective length of the shorter chain for that given N s . The number ofindependent events is obtained as ratio between the number of produced trajectories and the bin size, which is roughly τ int . ˆ m u,d N s β extr.c ν extr. .
432 61 p q . p q .
436 48 p q . p q .
439 17 p q . p q .
442 60 p q . p q .
445 29 p q . p q .
448 34 p q . p q .
451 68 p q . p q .
453 888 p q . p q ˆ m u,d N s β extr.c ν extr. .
647 61 p q . p q .
712 18 p q . p q .
755 59 p q . p q .
773 73 p q . p q .
787 83 p q . p q .
800 70 p q . p q .
810 96 p q . p q .
819 99 p q . p q .
827 77 p q . p q .
841 14 p q . p q .
851 35 p q . p q Table IV. Result of the quantitative collapse analysis. The critical temperature β c and the critical exponent ν have been foundminimizing Q p ¯ β c , ¯ ν q as defined in Eq. (17) for several decreasing values of ∆ x . β extr.c and ν extr. are the outcome of a linearextrapolation for ∆ x Ñ . Note that the reweighting resolution in β used to add new points between simulated ones variedbetween . and . and it has been chosen in order to have around values of the kurtosis to be later interpolated. ˆ m u,d N s points δβ ¨ β c ν B p β c , a χ d.o.f. Q(%) Ω min Ξ min .
432 41 p q . p q . p q ´ . p q .
436 49 p q . p q . p q ´ . p q .
439 18 p q . p q . p q ´ . p q .
442 47 p q . p q . p q ´ . p q .
445 13 p q . p q . p q ´ . p q .
448 28 p q . p q . p q ´ . p q .
451 54 p q . p q . p q ´ . p q .
453 50 p q . p q . p q ´ . p q .
647 87 p q . p q . p q ´ . p q .
711 76 p q . p q . p q ´ . p q .
755 51 p q . p q . p q ´ . p q .
773 64 p q . p q . p q ´ . p q .
787 83 p q . p q . p q ´ . p q .
800 37 p q . p q . p q ´ . p q .
810 91 p q . p q . p q ´ . p q .
820 12 p q . p q . p q ´ . p q .
827 85 p q . p q . p q ´ . p q .
841 12 p q . p q . p q ´ . p q .
851 24 p q . p q . p q ´ . p q ν (results on grey background are preliminary). The fits have been performed according to Eq. (14),considering the linear term only. The N s column contains the spatial lattice extents that have been included in the fits. Ω min and Ξ min are respectively the minimumoverlap percentage and the minimum symmetry percentage as defined in Eqs. (B3) and (B4) of Ref. 26. In the third and the fourth column, the number of fitted pointsper N s and the reweighting resolution in ββ