On the order of the QCD chiral phase transition for different numbers of quark flavours

The nature of the QCD chiral phase transition in the limit of vanishing quark masses has remained elusive for a long time, since it cannot be simulated directly on the lattice and is strongly cutoff-dependent. We report on a comprehensive ongoing study using unimproved staggered fermions with $N_\text{f}\in[2,8]$ mass-degenerate flavours on $N_\tau\in\{4,6,8\}$ lattices, in which we locate the chiral critical surface separating regions with first-order transitions from crossover regions in the bare parameter space of the lattice theory. Employing the fact that it terminates in a tricritical line, this surface can be extrapolated to the chiral limit using tricritical scaling with known exponents. Knowing the order of the transitions in the lattice parameter space, conclusions for approaching the continuum chiral limit in the proper order can be drawn. While a narrow first-order region cannot be ruled out, we find initial evidence consistent with a second-order chiral transition in all massless theories with $N_\text{f}\leq 6$, and possibly up to the onset of the conformal window at $9\lesssim N_\text{f}^*\lesssim 12$. A reanalysis of already published $\mathcal{O}(a)$-improved $N_\text{f}=3$ Wilson data on $N_\tau\in[4,12]$ is also consistent with tricritical scaling, and the associated change from first to second-order on the way to the continuum chiral limit. We discuss a modified Columbia plot and a phase diagram for many-flavour QCD that reflect these possible features.


Introduction
The order of the QCD chiral phase transition in the limit of massless quarks has been an outstanding question for decades.Many aspects of hadrons and their interactions are remnants of the chiral symmetry of QCD, which nature breaks only weakly by the small u, d-quark masses.One therefore expects the thermal crossover [1], which gradually restores chiral symmetry as a function of increasing temperature, to be similarly related to the phase transition in the chiral limit.Moreover, the phase transition in the chiral limit severely constrains the possibilities for the physical phase diagram at finite baryon densities (see, e.g., ref. 2 and references therein).A definitive, non-perturbative understanding of this transition is therefore mandatory.
Unfortunately, lattice QCD cannot be simulated in the chiral limit because the fermion determinant is singular for vanishing quark masses.Any numerical study necessarily proceeds by extrapolating results from a sequence of finite masses to the chiral limit, which inevitably introduces systematic uncertainties.Moreover, the traditional staggered and Wilson discretisations break chiral symmetry partially or entirely, so that the continuum limit has to be taken before the chiral limit can be approached.On the other hand, lattice actions with chiral symmetry at finite lattice spacing, such as domain wall or overlap fermions, are by an order of magnitude more expensive to simulate.For these reasons, the nature of the chiral phase transition could not be conclusively settled so far.
In this work, we attempt to reduce some of these limitations by performing a novel and comprehensive analysis of the nature of the chiral phase transition as a function of the number of degenerate quark flavours N f P r2, 8s and lattice spacing, N τ P t4, 6, 8u.We use the standard Wilson gauge and staggered fermion actions without improvement terms.On coarse lattices, this discretisation shows clear first-order transitions in the light mass regime, which weaken with increasing bare quark mass to vanish at some critical value.The bare parameter space of our lattice action, consisting of the lattice gauge coupling β, bare quark mass am, the number of flavours N f and the temporal extent of the lattice N τ , then features a three-dimensional region with first-order transitions, bounded by a critical surface separating it from the crossover region, which we systematically map out and trace along decreasing lattice spacing.Our analysis differs from others in exploiting the crucial fact that this surface has to terminate in a tricritical line in the chiral limit of the lattice theory, to which one can extrapolate using tricritical scaling with known exponents.Whenever this applies, it eliminates one important source of systematic uncertainty.
Knowing the order of the phase transition in the bare parameter space of the lattice theory, one can draw conclusions about the continuum phase transition when the continuum and chiral limits are taken in the appropriate order.Surprisingly, our currently available data prefer the continuum chiral phase transition to be second-order for all N f P r2, 6s, and suggest that it might stay second-order up to the onset of the conformal window at 9 À N f À 12.As a check of our analysis in an independent discretisation scheme, we reanalyse already published data for N f " 3 Opaq-improved Wilson fermions.These are fully consistent with tricritical scaling when approaching the lattice chiral limit, which implies the same conclusion of a second-order chiral transition for N f " 3.  2 The chiral phase transition in the continuum and on the lattice

The Columbia plot
To motivate our analysis by the general picture, we briefly summarise current knowledge about the chiral phase transition.The thermal QCD transition with physical quark masses has been known for some time to be an analytic crossover [1].Away from the physical point, the nature of the N f " 2`1 QCD thermal transition as a function of quark masses is usually summarised in a so-called Columbia plot [3], figure 1. 1 In the quenched limit QCD reduces to SU p3q Yang-Mills theory in the presence of static quarks, and shows a first-order phase transition in the continuum limit [4], associated with the spontaneous breaking of the Z 3 centre symmetry.Finite quark masses break this symmetry explicitly and weaken the firstorder phase transition, until it disappears along a critical line in the 3D Z 2 universality class.This region can be simulated directly [5,6], and the currently finest lattices predict the critical pseudoscalar mass for N f " 2 to be about m P S " 4 GeV [6].
In the chiral limit, the situation is more complicated, because it cannot be simulated.For a long time expectations have been based on an analysis using the epsilon expansion about " 1 applied to a linear sigma model in three dimensions [7].It predicts the chiral transition to be first-order for N f ě 3, whereas the case of N f " 2 is found to crucially depend on the fate of the anomalous U p1q A symmetry: If the latter remains broken at T c , the chiral transition should be second order in the Op4q-universality class, whereas its effective restoration would enlarge the chiral symmetry and push the transition to firstorder.For non-zero quark masses, chiral symmetry is explicitly broken.A first-order chiral  phase transition then weakens to disappear at a Z 2 second-order critical boundary, whereas a second-order transition disappears immediately.This results in the two scenarios for the Columbia plot depicted in figure 1.These early findings were confirmed by numerical simulations of the 3D sigma model [8] and a high-order perturbative analysis [9].Yet another analysis of the renormalisation group flow in a 3D effective theory [10] finds a possible symmetry breaking pattern to be U p2q L bU p2q R Ñ U p2q V in the case of a restored U p1q A , which would amount to a second-order transition in a different universality class.Finally, a recent functional renormalisation group analysis applied to U p1q A restoration in QCD, which is in good agreement with lattice susceptibilities at finite quark masses [11], favours the Op4q scenario in the chiral limit [12].
A lot of effort has been devoted to disentangle which of these situations is realised in lattice QCD.There, the bare parameter space of the Columbia plot is enlarged by a third axis for the lattice spacing a.The chiral critical line then traces out a chiral critical surface whose shape is discretisation-dependent, while all valid discretisations must of course merge to the same critical line in the continuum limit.Early simulations on coarse N τ " 4 lattices seemed to confirm the expected first-order transition for N f " 3, both with unimproved staggered [13,14] as well as standard [15] and Opaq-improved [16] Wilson fermions.A narrower first-order region is also seen for N f " 2 with unimproved staggered [17,18] and unimproved Wilson [19] fermions on N τ " 4.However, the location of the boundary line varies widely between these, indicating large cutoff effects.On the other hand, simulations with an improved staggered action (HISQ) do not see a first-order region on N τ " 6 lattices even for N f " 3 down to m P S Á 50 MeV [20].
The emerging picture then is as sketched in figure 2 for any fixed strange quark mass, with a shrinking first-order region as the lattice spacing a is decreased.How much of a first-order region remains in the continuum limit, i.e. what is the value of the critical quark mass m c u,d pm s q of the second-order boundary?It is obvious that any answer to this question strongly depends on the details of extrapolations.For the direct approach of following the chiral critical line along finer lattices, this has been demonstrated in the case of N f " 3 with Opaq-improved Wilson fermions [16,21,22].While on the finest N τ " 12 lattice the pion mass on the chiral critical point is still m c π « 110 MeV, continuum extrapolations give either finite or zero critical masses, depending on their functional form.Even for N f " 4 the situation is similar, both for unimproved staggered [23] and Opaq-improved Wilson quarks [24].For improved staggered fermions, no first-order region is seen at all at the N f values simulated so far.Recently, the critical temperature of the chiral phase transition was determined using HISQ fermions tuned to the physical strange quark mass, simulating a sequence of decreasing light quark masses down to " 58 MeV in the crossover region [25].These were extrapolated to the chiral limit employing the well-known scaling relation Note that this formula applies in the continuum for infinite volume.The quoted errors cover the systematic uncertainties obtained by interchanging the order of continuum and chiral extrapolations, as well as using either Op2q or Z 2 values for the critical exponents β, δ [25].This statement implies conversely, that an unambiguous distinction between the two scenarios based on this scaling equation is impossible at foreseeable accuracies.
A similar investigation with twisted mass Wilson fermions based on slightly larger pion masses confirms this behaviour [26].

Three-state coexistence and tricritical scaling
Some of these obstacles can be circumvented by an alternative analysis, which applies whenever the transition in the massless limit changes between first-order and second-order as a function of some parameter of the theory.As an example, consider the scenario in figure 1(b).A change in the order of the chiral phase transition between N f P t2, 3u implies the existence of a tricritical point on the strange quark mass axis.This is because the first-order transitions for m s ă m tric s in the limit m u,d " 0 represent a triple line, with a three-state coexistence (at the critical temperature) characterised by positive, negative and vanishing chiral condensate, cf.figure 3.As the strange quark mass is increased, the latent heat of the corresponding first-order transition diminishes, until it vanishes in a tricritical point.This point also marks the confluence of two ordinary critical lines separating crossover and first-order regions for the transitions between ˘x ψψy ‰ 0, respectively, and x ψψy " 0. These critical lines enter the tricritical point as a function of the symmetrybreaking scaling field m u,d with a known mean field exponent, since the upper critical dimension for a tricritical point is three.For a general overview of tricritical scaling in and beyond a mean field treatment, see ref. 27.The demonstrated existence of a tricritical point would therefore imply the orders of both N f P t2, 3u chiral phase transitions (but not their universality class).
While the chiral symmetry of staggered fermions is reduced to Op2q compared to Op4q in the continuum, it preserves invariance under rotations x ψψy Ñ ´x ψψy when am u,d " 0, and hence reproduces a possible triple line and tricritical point also at finite lattice spacing.In a first attempt with N f " 2 `1 on coarse N τ " 4 lattices, the chiral critical line was found to be consistent with tricritical scaling [28].Unfortunately, this is inconclusive for the same reasons as described in the last section: A finite portion of the critical line can always be fitted in terms of different polynomial forms, so that a presently impossible accuracy would be required close to the chiral limit in order to get a compelling distinction between the different versions of the Columbia plot reported in figures 1(a) and 1(b).

The chiral phase transition for N f mass-degenerate flavours
The way out is to exploit tricritical scaling in a setup, where a tricritical point is guaranteed to exist.In such a case the scaling form and its exponents are fixed, and one is only concerned about the location of the tricritical point.Such a situation emerges from a slight change of perspective and variables, as we suggested previously [18].We now consider degenerate quark masses only, with continuum partition function Instead of tuning the strange quark mass, an alternative interpolation between N f P t2, 3u, which generalises to larger N f , is achieved by an analytic continuation of N f to continuous, non-integer values.In the lattice formulation with rooted staggered fermions, whose determinant is raised to the power N f {4 in order to describe N f mass-degenerate quarks, this is implemented straightforwardly.The Columbia plot scenario figure 1(b) then translates to the version shown in figure 4, where the tricritical strange quark mass is replaced by a tricritical number of flavours, 2 ă N tric f ă 3, and the N f -axis to the right of it corresponds to the new triple line.The crucial advantage in this modified parameter space is that, since there is no chiral transition for N f " 1, a tricritical point N tric f ą 1 is guaranteed to exist as soon as there is a first-order region for any N f ą 1.In particular, the first-order scenario from figure 1(a) now also features a tricritical point, 1 ă N tric f ă 2. When a third axis for finite lattice spacing a is added to this picture, there must be a tricritical line N tric f paq in the plane m " 0, which represents the chiral limit of the Z 2 -critical surface separating lattice parameter regions with first-order transitions from crossover.
The principle of the analysis is now clear: Starting with the already known first-order transitions for N f P t3, 4u on N τ " 4 lattices, map out the Z 2 boundary lines until the tricritical scaling region is reached, and extrapolate to the chiral limit, In this way, N tric f « 1.7 was obtained on N τ " 4 lattices [18], implying the first-order scenario for N f " 2. As a powerful check of the continuation of N f as well as tricritical scaling, the critical quark mass for N f " 2, N τ " 4 obtained when keeping N f fixed and  varying (imaginary) chemical potential [17] is consistent, based on tricritical scaling, with results at N f P r2.1, 2.2s [18].In that case the quark mass is again the symmetry breaking scaling field, but with N f Ñ pµ{T q 2 in eq.(2.4).In the present work, we systematically extend our study from ref. 18 to larger numbers of flavours and finer lattices.
The chiral phase transition of QCD with many flavours is also of interest in the context of physics beyond the Standard Model.There the main focus is on the "conformal window" f , which denotes a range of theories with an infrared fixed point at a finite value of the gauge coupling.For N f ą N a.f.f asymptotic freedom is lost, while N f is the critical number of flavours marking the endpoint of the existence of a chiral condensate, and hence of chiral phase transitions, as indicated schematically in the possible phase diagram figure 5.A precise non-perturbative determination of this point is difficult for various technical reasons, but is expected in the range 9 À N f À 12.For an overview of the related strong interaction dynamics and its lattice investigations, see refs.29-31.For our present purpose, we stay at N f ă N f .Note that the decrease of the critical temperature in figure 5 is predicted by functional renormalisation group methods [32,33] and lattice simulations [34], which in the chiral limit implies N f to represent a quantum phase transition at T " 0. An interesting feature of the critical temperature T c pN f q is its approximately linear behaviour observed in refs.32, 33.It can be explained by expressing dimensionful quantities in terms of Λ QCD pN f q, whose perturbative expression can be expanded in N f , The perturbative second factor is indeed approximately linear due to the smallness of ε « 0.11, and appears to dominate the remaining N f -dependence of the dimensionless first factor until the conformal window is approached.There, non-perturbative scaling associated with the quantum critical point takes over and leads to a bending of the critical line [32,33].By contrast, the order of the chiral transition as a function of N f is not yet known.A change from second-order to first-order with increasing N f implies the already mentioned tricritical point N tric f .One would then also expect N f to be tricritical, as it represents yet another endpoint of a triple line.

Lattice simulations and analysis
For our numerical investigation, we work with the standard unimproved Wilson gauge and staggered fermion actions.All numerical simulations have been performed using the publicly available OpenCL-based code CL 2 QCD, which is optimised to run efficiently on AMD GPUs and contains an implementation of the RHMC algorithm for unimproved rooted staggered fermions.In particular, version v1.0 [35] has been employed for simulations on smaller N τ on the L-CSC supercomputer, while version v1.1 [36] has been run on the newer Goethe HLR supercomputer to run the most costly simulations.To effectively handle the thousands of necessary simulations, the BaHaMAS software [37] has been extensively used.
Regarding variable numbers of flavours, there is one aspect that is worth mentioning.It is well known that in the RHMC algorithm, when the condition number of the Dirac matrix is dominated by a few isolated tiny eigenvalues (e.g. in the chiral limit), the multiple pseudofermions technique [38] is advantageous, since it reduces the extra noise introduced into the system by the stochastic estimate of the fermion determinant.Using multiple pseudofermions, the maximum fermion force, which might trigger the molecular dynamics integrator instabilities (requiring in turn a smaller integration step size and hence making algorithm more costly), is reduced.However, when increasing the fermion mass at fixed N f , the condition number of the Dirac matrix decreases and using more than one pseudofermion turns out to be disadvantageous.On the other hand, in the present study, we increase the simulated quark mass range and the number of flavours simultaneously.Since a larger N fvalue also implies larger forces in the RHMC algorithm, it is worth wondering whether using multiple pseudofermions can still help for simulations at larger mass values.Therefore, the condition number of the Dirac matrix has been monitored in all simulations and it actually turned out that the increase of the fermion force produced by larger N f values prevails by far over its decrease due to larger fermion mass values.Indeed, it sometimes pays off to even increase the number of multiple pseudofermions when simulating at larger N f , despite am being larger.
Our goal is to determine the location and order of chiral phase transitions in the four-dimensional space spanned by the dimensionless parameters of our lattice action: The lattice gauge coupling β, the bare quark mass in lattice units am, the number of degenerate quark flavours N f , and the number of time-slices N τ .For any fixed value of N τ and N f , we achieve this by making use of two particular standardised moments, are due to the low number of simulated mass values and/or to the still low accumulated statistics.For N f P t6, 7u at N τ " 8, the available data are not yet sufficient to obtain an acceptable fit, and the quoted critical mass represents the middle point of simulated masses at which the order of the transition was clear from the ordering of the kurtosis on increasing spatial lattice size, i.e. towards the thermodynamic limit.
where the chiral condensate has been chosen as observable, O " ψψ, as it becomes the order parameter of the thermal phase transition in the chiral limit.In particular, to extract the order of the transition as a function of the quark mass, we evaluate the kurtosis B 4 pβ c , am, N σ q [39] of the sampled x ψψy distribution, where β c denotes the (pseudo-) critical coupling of the phase boundary, for which the zero-skewness condition B 3 pβ " β c , am, N σ q " 0 holds.In the thermodynamic limit N σ Ñ 8, the kurtosis B 4 pβ c , am, N σ q takes the values of 1 for a first order transition and 3 for an analytic crossover, respectively, with a discontinuity when passing from a first order region to a crossover region via a second order point; for the 3D Ising universality class of interest here, it takes the value 1.604 [40].On finite, increasing volumes this discontinuity is smoothed out and approached gradually with a rate characteristic of the universality class in question, B 4 pβ c , am, N σ q « 1.604 `c pam ´am c q N 1{0.6301σ with c P R . (3.2) Data have been produced and analysed in a completely analogous way to that explained in refs.18, 41 and, in particular, the critical mass am c has been extracted at fixed N τ and N f by fitting the kurtosis data according to this finite size scaling formula.However, some steps in the analysis procedure have been improved and will be discussed next.This is also the reason why older results are reported here slightly changed.First of all the evaluation of centred moments of observables that are stochastically estimated might be biased, if carried out simply using all stochastic sources for all powers of the observable.A bias reduction technique consisting of not mixing estimates coming from the same stochastic source when calculating powers of ψψ can be applied in order to figure out how much the naive estimate is biased.This has been implemented for an arbitrary number of stochastic estimates and has been used in the analysis of all data. 2nterestingly enough, significant corrections to B 3 and B 4 have been found on single chains, but these are smoothed out when it comes to analyse the merged chains.
Since the subsequent analysis presented in section 4 heavily relies on the outcome of the B 4 -fits, we decided to improve the error estimate on am c using a more accurate procedure.Values of B 4 pβ c , am, N σ q are obtained using the multiple-histogram method [42], and their error is in turn calculated out of hundreds of bootstrap estimators.Each of these can be used as new central value of B 4 , so that a new B 4 -data set can be built and fitted to extract a new (bootstrap estimator of) am c . 3 Combining all these estimators of am c , a bootstrap error can be calculated and compared to that coming from the χ 2 -minimisation.This method has been run on all data sets and it often leads to smaller errors on the critical Z 2 mass.The outcome of all fits can be found in table 1, where also the simulated mass range has been included.The uncertainties on am c are those obtained with our new bootstrap-fitting technique.In appendix A a detailed overview of the simulations can be found.To give an idea of the numerical effort: Over 600 values of β have been simulated in total, producing about 120 millions of trajectories.
Finally, we explain how the value of β c at am c reported in the last column of table 1 has been extracted.Rather than running further costly simulations at the extracted critical masses, the values β c found on the largest simulated N σ for each mass value have been linearly interpolated in order to extract β c at am c .In order to attribute an uncertainty to this interpolated value, the data have been also interpolated at am c ˘σamc , namely at the critical mass plus or minus one standard deviation.

The chiral critical surface in the parameter space of staggered fermions
For the interpretation of the lattice data, it is useful to recall the structure of the fourdimensional parameter space, tβ, am, N f , N τ u in the limit of infinite spatial volume.The discussion is facilitated by assuming some analytic continuation of both N f and N τ to real values, so that all parameters can be regarded as continuous and treated on equal footing.The phase boundary associated with chiral symmetry restoration, which we define by vanishing skewness, B 3 " 0, then represents a three-dimensional parameter region.That subspace consists of a region with first-order transitions and a crossover region, separated by a chiral critical surface, on which the kurtosis assumes its 3D Ising value, B 4 " 1.604.The identification of this chiral critical surface in the bare parameter space of the standard staggered lattice action is the central result of our simulations.In the following we analyse this surface projected onto different variable pairings.

Plane of gauge coupling and quark mass
The chiral critical surface projected onto the pβ, amq-plane is shown in figure 6.As in the case of the Columbia plot, every point of figure 6 represents a point of the phase boundary, with the fourth parameter implicitly determined by the other three.In figure 6(a), we show slices of the phase boundary corresponding to different values of N τ .The curve for every given N τ represents the critical bare parameter combinations with second-order phase transitions, and is parametrised by different values of N f , increasing from left to right.The region below the curve corresponds to first-order transitions, whereas the region above the curves correspond to analytic crossover behaviour.Obviously, this qualitative picture is unchanged when it is restricted to a grid with the physically meaningful integer values for N f , as in figure 6(b).However, in that case the underlying functional behaviour is less constrained and difficult to extract.By filling in non-integer values of N f , a convincing picture of tricritical scaling is obtained: For each given N τ , the lower β-axis represents a triple-line of first-order transitions in the chiral limit of the lattice action, which weakens as β is increased.This first-order chiral transition disappears in a tricritical point, from which the line of ordinary Zp2q critical points, β c pamq, emerges.The critical curve is excellently described over a wide mass range by a next-to-leading order fit in the scaling variable, cf.table 2, Figure 6(b) shows the same chiral critical surface, this time sliced at fixed N f , so that the curves are implicitly parametrised by N τ , increasing from right to left.Again, for each N f the parameter region below the curves describes first-order transitions and the region above it analytic crossovers.However, with presently only three N τ P t4, 6, 8u, a functional behaviour is difficult to assess.This further illustrates the advantage that the additional non-integer N f values offer in figure 6(a).Note that a finer resolution in N τ could in principle be achieved by working with anisotropic lattices, but this would come at the price of introducing a second lattice gauge coupling, and thus another dimension to the parameter space of the lattice theory.
Clearly, figures 6(a) and 6(b) show the same chiral critical surface, presented by slicing the three-dimensional subspace of the phase boundary in different directions.In particular, the bare quark mass am is the symmetry breaking field and, sufficiently close to the chiral limit, behaves as a tricritical scaling variable in any plane that contains a tricritical point and is not orthogonal to the am-axis.This is guaranteed to be the case for every curve in figure 6(a), where all N f are included in every fixed N τ curve, but not in figure 6(b), since for a given fixed N f there need not be a tricritical point.The chiral extrapolations shown in figure 6(a) then lead to a one-dimensional tricritical line, With a non-perturbative scale setting for the lattice spacing apβq, β tric could be converted into a tricritical temperature T tric pN τ q " T tric paq, corresponding to the tricritical point in figure 5, displaced due to the finite lattice spacing.This is not easily feasible in practice, because a scale would have to be set for arbitrary values of N f in the chiral limit.Fortunately, we do not need this temperature for our present purpose.Here, we conclude that there is unambiguous evidence for the chiral critical surface terminating in a tricritical line in the chiral limit of the lattice theory, which it approaches as a function of the scaling field pamq 2{5 .With this knowledge, we can now analyse the chiral critical surface in different pairings of variables.Part of the data for N τ " 4 [18] and N τ " 6 [41] have been published in similar form previously, here we have added several N f on N τ P t4, 6u and all of the N τ " 8 data. 4he parameter region below the curves corresponds to first-order transitions and above to analytic crossovers.The N f -axis represents a triple line for larger values, whose latent heat decreases with diminishing N f until it vanishes in a tricritical point N tric f .The critical boundary line must therefore approach this point with scaling behaviour

Plane of quark mass and number of flavours
Since the statistical error of our calculations is on the critical quark mass while N f is exact, we consider the inverted series instead, where trivially am tric pN τ q " 0 in all cases.It was already demonstrated in ref. 18 that for N τ " 4 the two lowest mass points in the present work and the N f " 2 result from ref. 17 are consistent with a leading-order scaling curve.However, as figure 7(a) shows, the scaling window is rather narrow in this variable pairing, and even a next-to-leading order fit breaks off the data rather quickly, both for N τ P t4, 6u.The reason is that scaling gets quickly superseded by a nearly linear behaviour over a large range of N f , analogous to T c pN f q observed in the chiral limit of the continuum theory, cf.figure 5.The emerging qualitative picture is thus: With increasing N τ , i.e. decreasing lattice spacing, the chiral phase transition observed on the lattice at fixed N f weakens, so that the first-order region gets narrower.This is in line with all earlier observations discussed in the introduction.Because of the smallness of the scaling region, it is difficult to determine the value of N tric f quantitatively.On N τ " 4, the N f " 2 theory has a first-order chiral transition, in agreement with a result based on an extrapolation from imaginary chemical potentials [17].On N τ " 6 the N tric f value shifts towards larger values and for N τ " 8 there are not enough data yet to constrain the scaling.Nevertheless, the size of the observed cutoff effects strongly suggests further movement of the critical line on its way to the continuum.One question is whether the critical quark masses scaled by the temperature will settle at finite values, which we checked is not the case yet.However, figure 7(a) raises the additional question whether the intercept N tric f slides further to the right and possibly beyond N f " 3.More data and larger N τ are required for a definite answer, but we can get a hint by analysing yet another variable pairing.

Plane of quark mass and temperature
We now turn to the pam, aT " N ´1 τ q-plane.Figure 7(b) shows a sequence of curves for different fixed N f .The region below the lines corresponds to first-order transitions, that above to analytic crossovers.For this variable pairing, the near-linear N f -dependence through Λ QCD does not dominate the behaviour of the individual curves, but rather shows up in the spacing between them.However, a tricritical point is not guaranteed to exist for any given value of N f , so we have to test for the functional behaviour.
Conventional wisdom predicts a definite first-order transition in the continuum for any N f ě 3 [7,9].This implies a finite Z 2 -critical bare quark mass m c in physical units, with the usual Opa, a 2 , . ..q lattice corrections, and am c pN τ q Ñ 0 in the limit N τ Ñ 8.The critical lines in lattice units then approach the continuum as am c pN τ , N f q " F1 pN f q aT `F 2 pN f q paT q 2 `F 3 pN f q paT q 3 `O`p aT q 4 ˘.(4.5) In figure 8(a) and table 4 we fit our N f " 5 data to two variations of this functional form.
Clearly our data are incompatible with eq.(4.5), and the same holds for N f P t6, 7u.The other possibility is for the critical lines to end on the aT -axis, which implies a tricritical point.This is described by a polynomial in the scaling variable, which we invert as in the previous section, Unfortunately, with only three data points per line at most, meaningful next-to-leading order fits with non-vanishing intercept are impossible.However, data in figure 7(b) with three N τ -values show only a slight deviation from leading-order scaling, with an unambiguous curvature that can be covered by a next-to-leading order term as before.In the following, we take this as an indication for tricritical scaling, for which a much stronger justification will be provided in section 5.2, and attempt to constrain the location of the tricritical points.
A next-to-leading order scaling interpolation, without error or χ 2 -estimate, can be obtained from eq. (4.7) and is shown as the central line in figure 8(b).We estimate the associated uncertainty to its x-axis crossing by alternatively interpolating our data with the inversion of eq.(4.6) with a linear term only, or fitting them to the inversion of eq.(4.6) with a quadratic term only, `aT ´aT tric pN f q ˘. (4.9) The two x-axis crossings of eqs.(4.8) and (4.9) bracket the extrapolation using eq.( 4.7) at NLO. 5 Eq.(4.8) is used to interpolate the data on the two coarsest lattices only, which corresponds to the most conservative choice.While this procedure leaves a sizeable error band to be improved upon in the future, it still gives an interesting estimate of where the chiral critical surface ends in the massless limit.In particular, the lower bounds on the tricritical values aT tric pN f q, or equivalently upper bounds on N tric τ pN f q, can only be avoided by a change of sign in the curvatures of pam c pN τ qq 2{5 for N τ ě 10.   4 pertaining to a first-order scenario, eq. ( 4.5), do not describe the data.The data are instead consistent with a second-order scenario.
In (b) the green solid line interpolates the data according to eq. ( 4.7) at truncated at next-toleading order.Dashed and dotted lines represent interpolations according to eqs.(4.8) and (4.9), respectively.
In the absence of such a drastic change, the near-tricritical scaling of our data suggests that for any fixed N f À 6 there is a maximal value N tric τ pN f q, for which the chiral transition on the lattice can be of first order.Following the qualitative behaviour in figure 7(b), this value keeps increasing with N f .There will then be a special (not necessarily integer) value N 0 f , for which the tricritical point hits the origin of the figure, aT tric pN 0 f q " 0, and this is the situation when a tricritical point exists in the continuum.Only for N f ą N 0 f can a first-order phase transition persist all the way to the continuum.This part of our analysis is consistent with N 0 f ą 6.

Plane of number of flavours and temperature
Finally one can also slice the chiral critical surface at fixed values of the quark mass am and obtain the critical line separating first-order transitions and crossovers in the pN f , N ´1 τ q-plane, which reduces to a grid for integer variable values.Such a diagram relies on the interpolation/extrapolation functions discussed in the previous sections.The most interesting plane, and goal of our investigation, is the chiral limit, am " 0, shown in figure 9.In this case the chiral phase transition must be non-analytic, and the tricritical "line" shown in the figure separates a region of first-order phase transitions above it from a region of second-order phase transitions below it.The two rightmost points correspond to proper extrapolations from table 3, where the existence of the tricritical points is guaranteed and only the accuracy may vary.By contrast, the error band to the left of it reflects the bounds estimated in the last subsection, and hinges upon better data description by tricritical scaling only.We will come back to this issue in section 5.2.
Taking these chiral extrapolations at face value, it is intriguing to speculate how this tricritical line might continue beyond the N f -range covered by our simulations.In fact, the entire first-order region corresponds to an area of triple points, which must be bounded by either a tricritical line or the boundaries of the am " 0 parameter space.Since our error band is again consistent with an approximately linear N f -dependence of T tric pN f q on the lattice in the range N f P r2, 7s, and in the absence of any other special point, we conjecture the tricritical line to hit the N f -axis in the conformal point on the N f -axis, i.e., N 0 f " N f .

Implications for the chiral phase transition in the continuum
Knowing the phase diagram in the bare parameter space of the lattice theory, we can draw conclusions for the continuum limit.Here some care is in order.A necessary condition to avoid problems with the rooting of the staggered fermion determinant requires the continuum limit to be taken before the chiral limit (for reviews of the staggered fermion discretisation and its potential problems, see refs.44, 45 and the references therein).A procedure to approach the continuum chiral transition for given N f would then be: (i) Find the phase boundary β c pN τ , amq for a chosen bare quark mass and N τ ; (ii) Determine the order of the transition; (iii) Set the scale for the lattice spacing, compute the pseudoscalar mass m P S pβ c , amq and the critical temperature T c pam, N τ q; (iv) Repeat steps (i) to (iii) for a sequence of increasing N τ along a line of constant physics, i.e., tuning the bare quark mass such that m P S pβ c , amq in physical units stays fixed; (v) Once N τ is large enough that the order of the transition does not change anymore, continuum extrapolation (β Ñ 8, N ´1 τ Ñ 0q gives T c pm P S q; (vi) Repeat steps (i) to (v) for a sequence of decreasing pseudoscalar masses.The implications of our results are easiest to see in the pam, aT q plots and illustrated schematically in figure 10, where the continuum limit is located in the lower left corner, am " aT " 0. Continuum approaches for different physical masses m P S correspond to different lines of constant physics ending in the origin.For every N f with a finite N tric τ pN f q, simulations with N τ ą N tric τ will show exclusively crossover behaviour for any tuning of the bare masses ampβq ‰ 0, figure 10(b).Consequently, for these N f the chiral limit in the continuum corresponds to an isolated point with a second-order phase transition.By contrast, a first-order chiral phase transition in the continuum extends over a finite mass range to a Z 2 -point at some critical quark mass m c , such as in the scenario figure 1(a).This requires the chiral critical line to terminate in the origin as N τ Ñ 8, cf.figure 10(a).In such a case, there is no tricritical scaling and one expects a different functional form of the critical line.Our data for N f P r2, 6s clearly prefer the first option.
Of course, any conclusion has to be consistent between all parameter pairings.For example, in figure 6(b) the continuum limit is on the upper left corner, pβ " 8, am " 0q.To obtain a continuum first-order transition requires the critical line for every N f in question to bend upwards, so as to include the entire β-axis in the first-order region.Finally, in figure 9, a continuum limit with first-order transitions for N f ě 3 requires the tricritical line to hit the N f -axis below N f " 3.
Based on the available data only, we thus have to conclude that for 2 ď N f ď 6, and possibly 2 ď N f ď N f , the regions of first-order transitions displayed in the bare parameter space of the standard staggered lattice theory are not continuously connected to the chiral limit in the continuum, and hence represent lattice artefacts.Note that this conclusion is independent of whether or not additional unphysical phases, like a spontaneously broken staggered shift symmetry [46,47], exist in the current discretisation.Such phases cannot be detected by the chiral condensate alone and represent lattice artefacts as well.Like the first-order regions described here, the system enters such phases when the chiral limit is approached before the continuum limit, while they are avoided by the correct ordering of limits.
Lattice artefacts are discretisation dependent, i.e., the location of the chiral critical surface investigated here will be different for different lattice actions.This would explain the fact that no first-order transition region is seen in simulations with highly improved staggered actions for N f " 3 [20].It would be interesting to confirm this by investigating how the chiral critical lines shift under tuneable improvement, such as an increasing amount of smearing steps in the action.One would expect an effect similar to that of decreasing the lattice spacing, as was indeed observed in an earlier study with a stout-smeared staggered action for N f " 3 [48].In particular, also the tricritical boundary line in figure 9 should then shift upwards, while its left end should remain anchored at the conformal point.This is because at a second-order point in the continuum the correlation length diverges, so that the UV-differences between different discretisations are suppressed and all actions with a correct continuum limit must reproduce that point.

Wilson fermions
In view of our surprising conclusions for N f ě 3, and in order to rule out potential problems with the rooting procedure of staggered fermions, it is desirable to cross check our analysis in Wilson's discretisation scheme.In that case, matters are complicated by the fact that the Wilson term breaks chiral symmetry completely at any finite lattice spacing.As a consequence there is an additive renormalisation of the quark mass, and the naive chiral condensate is always non-zero.This raises the question whether an analysis based on threestate coexistence and tricritical scaling is applicable at all.In this section we suggest how this could be realised in the bare parameter space of Wilson fermions, and find compelling evidence for a tricritical point in already published N f " 3 Wilson data [22].

Parameter space of Wilson fermions for fixed N f
In the standard formulation of Wilson's fermion action with N f degenerate flavours, the bare quark mass am is traded for the hopping parameter, κ " 1 2pam `4q . (5.1) For the tree-level action, κ c " 1{8 represents a "chiral limit" defined by a vanishing bare quark mass.In the interacting theory, this critical parameter value gets shifted as a function of the gauge coupling β, leading to an additive renormalisation of the quark mass, A subtracted and renormalised chiral condensate can be defined based on an axial Ward-Takahashi identity [49], x ψψy R " 2pam q qZ ÿ x xπpxqπp0qy , (5.In (b) N τ is implicitly tuned so that each point pκ, βpN τ qq is on the phase boundary.The line κ c pβ, N τ pβqq denotes the chiral limit at finite temperature and is a triple line, which ends in a tricritical point.
with πpxq an interpolating operator for the pseudoscalar meson and a renormalisation factor Z. Approaching the chiral limit, the pseudoscalar meson mass and the quark mass are related as in the continuum, am 2 P S 9 am q . (5.4) It is therefore customary to define κ c pβq by the vanishing of the pseudoscalar meson mass in the vacuum, i.e., am P S pκ c pβq, βq " 0 at N τ " 8.This is shown schematically as a dashed line in figure 11(a).Towards the strong coupling region, this line meets the parityflavour violating Aoki phase [50,51], which ends in a cusp [52,53] whose location depends on the lattice action and the value of N τ .Around κ c pβq, Wilson chiral perturbation for the theory predicts a metastability region corresponding to a first-order bulk transition between positive and negative quark mass, while the meson mass stays finite everywhere, both for untwisted and twisted mass [54,55].A metastability region has been identified numerically at zero temperature [56] as well as at finite temperature [57,58], but its location and extent depend strongly on the chosen action and N τ [59].The series of N f " 3 data [16,21,22], which we re-analyse below, is based on the RG-improved Iwasaki gauge action [60] and a non-perturbatively Opaq-improved Wilson clover fermion action [61].We are not aware of a dedicated study of the bare phase diagram pertaining to the precise action and parameter tunings used in those simulations, besides determining the line κ c pβ, N τ " 8q.However, a previous study using the same action with a mean field tuning of the clover coefficient [62] reports a phase diagram as sketched by the dashed lines in figure 11, with no additional structures besides an Aoki phase in the strong coupling region, so we will base our discussion on this situation.
First, it has to be emphasised that for studies of the thermal phase transition we need the lines κ c pβ, N τ q for the finite N τ under consideration, and not κ c pβ, N τ " 8q, which is only needed to set the scale.The former marks the vanishing of the pseudoscalar screening mass in the low temperature phase, and is related to the latter by an expansion in powers of N ´1 τ " aT , In the literature the difference between the two is often dismissed, being of Opaq, whereas in fact it is qualitatively crucial.The partition function at finite N τ has no singularities on the line κ c pβ, 8q (except at its crossings with the thermal transition).Furthermore, the subtracted chiral condensate has finite values with different signs across κ c pβ, N τ q, which should therefore mark a first-order transition. 6Following this line with increasing β at fixed N τ , the thermal chiral phase transition is reached at some critical coupling.From this point the thermal transition lines κ t pβ, N τ q branch off into the positve and negative quark mass directions, respectively, along which the chiral transition weakens to end in a critical point.At the branching point the line κ c pβ, N τ q should terminate, since on the large-β-side of the thermal transition the Matsubara modes " 2πT produce an always non-zero screening mass and the subtracted chiral condensate can pass through zero smoothly.The branching point κ c pβ, N τ q " κ t pβ, N τ q then represents the Wilson version of the continuum triple point discussed in figure 3.For a sequence of increasing N τ 1 ă N τ 2 ă . .., the lines κ c pβ, N τ q gradually approach κ c pβ, 8q, the triple points trace out a triple line and the critical endpoints form a chiral critical line separating first-order transitions for light quarks from the crossover region.This is shown in the Columbia-like figure 11(b), which gives the analogue of the critical quark mass lines at fixed N f in the staggered case, figure 6(b).There are then two possibilities for the situation at larger β-values: Either the triple line extends all the way to the continuum limit, or it terminates in a tricritical point, where it changes to a second-order transition, as sketched in figure 11(b 12 to eq. (5.6).

Numerical analysis for N f " 3
We now apply our analysis to the sequence of N f " 3 clover-improved Wilson fermion simulations with N τ P t4, 6, 8, 10, 12u, whose data are published in refs.16, 21, 22.In these works, the pseudoscalar meson mass was evaluated along the chiral critical endpoints in the lattice parameter space, which we reproduce in figure 12.The region below the critical line corresponds to first-order transitions, the region above to crossover.To elucidate tricritical scaling with the quark mass, we employ eq. ( 5.4) and convert the scaling relation for the quark mass as a function of aT " N ´1 τ , eq. (4.7), to a scaling relation for the pseudoscalar mass, Figure 12 shows three different fits of the data to this scaling form, as detailed in table 5.
The first fit demonstrates that the data for the smallest critical masses, corresponding to the finest N τ P t8, 10, 12u lattices, are excellently described by leading-order tricritical scaling.By contrast, they are incompatible with am c P S " N ´1 τ pχ 2 d.o.f.ą 100q and in tension with am c P S " N ´2 τ pχ 2 d.o.f.ą 3q, one of which would be the expected behaviour in the absence of a tricritical point.Moreover, a next-to-leading-order scaling fit including N τ " 6 is too good to be properly constrained, while fitting the entire range N τ P r4, 12s produces tension but only mild shifts the extracted parameters.Finally, interpolating the next-to-leading order scaling form to N τ P t4, 6, 8u data only, as we did in our staggered analysis in section 4.3, aT tric is only " 25% off its value from the best fit.
Altogether the available Wilson data thus show: (i) a clear preference for tricritical scaling and (ii) the stability of the analysis when restricted to the coarser lattices N τ P t4, 6, 8u, which we employed in section 4.3.
Avoiding the existence of a finite N tric τ would require a change of curvature in figure 12 for N τ ą 12.In the absence of such a drastic change, we conclude that any sequence of simulations with N τ ą N tric τ with N tric τ " paT tric q ´1 « 24 will inevitably approach the continuum chiral limit from the crossover region, implying a second-order transition, in full accord with our staggered results.Conversely this would imply, however, that  simulations of the phase transition with N τ ă 24 cannot reproduce the correct continuum physics despite the improvement of the action.This has repercussions for, e.g., studies of the U p1q A -anomaly with the Wilson action [63,64].By contrast, unimproved staggered fermions should reproduce the second order of the N f " 3 transition for N τ Á 10 already.

Conclusions
In summary, we have conducted a comprehensive analysis of the finite temperature chiral transitions observed in the bare parameter space tβ, am, N f , N τ u of lattice QCD with unimproved staggered fermions.In particular, we have mapped out the chiral critical surface, which separates the first-order transition region from the crossover region.In the plane of vanishing bare quark mass, this surface terminates in a tricritical line N tric τ pN f q, which for our presently available data is consistent with N tric τ ă 8 for all N f P r2, 6s.
The necessity to take the continuum limit before the chiral limit then enforces any appropriate series of simulations to approach the combined limits from the crossover region of the bare phase diagram.This implies the chiral phase transition in the massless limit to be of second order for all N f ď 6, and possibly up to the conformal window N f ď N f .As a crosscheck for our findings, we have reanalysed already published data from simulations with N f " 3 Opaq-improved Wilson fermions, which are equally consistent with tricritical scaling and a tricritical point at a finite N tric τ .Hence, this entirely different discretisation is consistent with the continuum chiral phase transition to be of second order as well.
Taking our results seriously, the continuum Columbia plot for N f " 2 `1 would have to be modified to look as in figure 13, with a second-order line all along the m s -axis.Our analysis has nothing to say about the universality class of this second-order line.However, chiral symmetry being different in the limiting N f P t2, 3u cases, one would expect the set of critical exponents associated with these transitions to smoothly cross from one universality class to the other as the strange quark mass varies.Our findings for larger N f suggest the phase diagram of many flavour QCD shown in figure 14, with a second-order chiral phase transition in the massless limit, which ends in a quantum tricritical point at the onset of the conformal window. 7an our conclusions for N f ě 3 be reconciled with the opposite predictions from 3D sigma-models [7][8][9]?These works investigate 3D effective Landau-Ginzburg-Wilson φ 4 theories with the desired global symmetry and symmetry-breaking patterns, as well as a 'tHooft term added for the anomaly.However, in mean field theory it is well known that φ 4 theories cannot accommodate tricritical points, which requires a φ 6 -term.A φ 6 -term is consistent with renormalisability in 3D, and indeed 3D φ 6 theories display non-trivial infrared stable fixed points also in functional renormalisation group [66] and lattice [67] analyses.Further studies of their connection to QCD would now be interesting.
While our data should be further supplemented for definitive conclusions concerning the continuum, an analysis employing tricritical scaling offers a new extrapolation tool for all discretisations observing first-order chiral transitions.

A Detailed overview of simulations
In tables 6 to 8, a detailed overview of the accumulated statistics is offered.Further information about previous simulations can be found in ref. 18.
The total statistics refers to the whole spatial volume, and the statistics per chain can be inferred using the number of simulated β values, considering that four Monte Carlo chains have been always run at each β value. 8Simulating multiple chains per fixed set of parameters helps in keeping under control the amount of accumulated statistics, and in particular in taking decisions about when to stop simulations.
The value of β c is reported without uncertainty and refers to the reweighted point which gives the closest skewness of the order parameter to zero.
The quantity n c pB 3 q « 0 denotes the number of chains at the outermost β values (i.e. the furthest β from β c ), which are compatible with zero.Ideally, this number should be 0, which would ensure the boundaries of the simulated β-range to belong to the two different phases of the system, while not being too far away from β c .This in turn would guarantee to be correctly cornering the phase transition within the given statistics.Values different from zero require attention and can be accepted if small (e.g. 1) and rare.Larger values signal a to-be-continued simulation, as it is still the case for several N τ " 8 sets.
Finally, n max σ pB 3 q denotes by how many standard deviations the skewness deviates between the two maximally separated chains.This quantity should be as small as possible, and we usually aim at having it not much larger than 3.

Figure 1 .
Figure 1.Possible scenarios for the order of the thermal QCD transition as a function of the quark masses.Every point of the plot represents a phase boundary, with an implicitly associated (pseudo-)critical temperature T c pm u,d , m s q.

Figure 2 .
Figure 2.A Columbia-like plot for the (a,m u,d )-plane.The chiral phase transition for any fixed value of m s weakens as the lattice spacing is decreased.

Figure 3 .
Figure 3. Sketch of the T ´mu,d phase diagram.The temperature axis with T ă T c corresponds to a first-order coexistence line of ˘x ψψy ‰ 0, and T c pm u,d " 0q represents a triple point.

Figure 4 .
Figure 4. Columbia plot for mass-degenerate quarks.Every point represents a phase boundary and has an implicitly associated T c pm, N f q.

1 s t o r d e r t r i p l eFigure 5 .
Figure 5. Schematic representation of a possible scenario for the T ´m ´Nf phase diagram with different numbers of light mass-degenerate flavours.

Figure 7 .
Figure 7.The chiral critical surface projected onto different planes.Every point represents a phase boundary with an implicitly tuned β c pam, N f , N τ q.Lines in (a) represent next-to-leading order scaling fits to eq. (4.4) as in table 3. Solid lines in (b) are the analogue of the solid line in figure 8(b); refer to section 4.3 for more details.

Figure 7 (
Figure 7(a) shows the chiral critical surface projected into the plane of bare quark mass and number of flavours for different values of N τ .This is the lattice version of figure 4.Part of the data for N τ " 4[18] and N τ " 6[41] have been published in similar form previously, here we have added several N f on N τ P t4, 6u and all of the N τ " 8 data.4The parameter region below the curves corresponds to first-order transitions and above to analytic crossovers.The N f -axis represents a triple line for larger values, whose latent heat decreases with diminishing N f until it vanishes in a tricritical point N tric f .The critical boundary line must therefore approach this point with scaling behaviour Fit in a second-order scenario.

Figure 8 .
Figure 8. Chiral critical line for N f " 5. Fits from table4pertaining to a first-order scenario, eq.(4.5), do not describe the data.The data are instead consistent with a second-order scenario.In (b) the green solid line interpolates the data according to eq. (4.7) at truncated at next-toleading order.Dashed and dotted lines represent interpolations according to eqs.(4.8) and (4.9), respectively.

Figure 9 .
Figure 9. Upper and lower bounds on the location of the tricritical line N tric f pN τ q for am " 0. The lower left region represents second-order chiral transitions, the region above the band represents first-order transitions.These have no continuous connection to the continuum theory at aT " 0 and are discretisation artefacts.

Figure 10 .
Figure10.Different scenarios of the approach to the continuum chiral limit along several lines of constant physics.

1
st o r d e r t r i p l e (b) Columbia-like plot for Wilson fermions.

Figure 11 .
Figure 11.Suggested phase structure of finite temperature lattice QCD with Wilson fermions.In (b) N τ is implicitly tuned so that each point pκ, βpN τ qq is on the phase boundary.The line κ c pβ, N τ pβqq denotes the chiral limit at finite temperature and is a triple line, which ends in a tricritical point.

Figure 13 .
Figure13.The Columbia plot in the continuum, as predicted by our analysis.

Figure 14 .
Figure 14.Suggested phase diagram for the chiral phase transition as a function of N f .

3 q
N f am β c | Total statistics | Number of β values | n c pB 3 q « 0 | n max σ pB Aspect

Table 1 .
Overview of the finite size scaling analysis and of the β c interpolation.Large values of χ 2 d.o.f.
The chiral critical surface, separating 3D parameter regions with a first-order transition from crossover regions, projected onto the pβ, amq-plane.Every point represents a phase boundary with an implicitly tuned β c pam, N f , N τ q.In (a) the lines represent next-to-leading order fits in the scaling variable, cf.table 2.

Table 4 .
Fits of the critical lines in figure8(a) to eq. (4.5).
Tricritical scaling of the critical pseudoscalar mass for N f " 3 clover-improved Wilson fermions.The data are taken from Table IV in ref. 22.

Table 5 .
). Fits of the data in figure