Lattice Constraints on the QCD Chiral Phase Transition at Finite Temperature and Baryon Density

The thermal restoration of chiral symmetry in QCD is known to proceed by an analytic crossover, which is widely expected to turn into a phase transition with a critical endpoint as the baryon density is increased. In the absence of a genuine solution to the sign problem of lattice QCD, simulations at zero and imaginary baryon chemical potential in a parameter space enlarged by a variable number of quark flavours and quark masses constitute a viable way to constrain the location of a possible non-analytic phase transition and its critical endpoint. In this article I review recent progress towards an understanding of the nature of the transition in the massless limit, and its critical temperature at zero density. Combined with increasingly detailed studies of the physical crossover region, current data bound a possible critical point to $\mu_B>3T$.


Introduction
Many salient features of the spectrum of hadrons and their interactions observed in nature are determined by the near-chiral symmetry of the QCD Lagrangian and its spontaneous breaking by the QCD vacuum.The fact that nature breaks chiral symmetry also explicitly by non-vanishing quark masses can be taken into account in terms of chiral perturbation theory [1,2], due to the smallness of the uand d-quark masses.
In a hot and/or dense medium, chiral symmetry is expected to be gradually restored once the temperature exceeds T ą " 160 MeV or the baryon chemical potential µ B ą " 1 GeV.Associated with this change of the realised symmetry, one also expects a change of dynamics and its underlying degrees of freedom which, asymptotically, should become quarks and gluons.Similar to the vacuum properties of the theory, one also expects the properties of this thermal transition to be closely related to those of the theory in the massless limit.
The low energy scales of the transition region demand a non-perturbative first principles approach like lattice QCD.Unfortunately, a severe sign problem prohibits simulations by importance sampling for non-vanishing chemical potential, for introductions see [3][4][5].Despite tremendous efforts over several decades, no genuine solution to this problem is available to date, and knowledge of the QCD phase diagram by direct calculation remains scarce.Nevertheless, during the last years there has been considerable progress in studying the properties of the thermal QCD transition in an enlarged parameter space, where the number of flavours and quark masses is varied away from their physical values.Besides giving theoretical insight into the interplay of symmetries and dynamics in controlled situations, such studies also provide constraints on the physical phase diagram, which are beginning to become phenomenologically relevant.In this article, I collect some recent results allowing to better understand and remove lattice artefacts, resulting in a consistent picture that also Symmetry 2021, 1, 0 2 of 28 accommodates apparently conflicting statements based on earlier simulations.In particular, I will focus on developments suggesting a modified version of the Columbia plot and constraints on the location of a possible critical point.

Some Lattice Essentials
For readers less familiar with lattice QCD, it may be useful to recall a number of basic definitions and issues that are important for the interpretation of the following lattice results.For details, see the corresponding textbooks [3,6,7].Conventionally one considers Euclidean space-time discretised by hypercubic lattices with N 3 s ˆNτ points separated by a lattice spacing a.This represents a box with spatial length L " N s a and Euclidean time extent aN τ " 1{T, which corresponds to inverse temperature as in continuum thermal field theory, once (anti-)periodic boundary conditions for (fermionic) bosonic fields have been applied.When a quantum field theory is formulated on the lattice, its action (and all observables) differ from those in the continuum.As an example, consider Wilson's SUpN c q pure gauge action which is constructed from so-called plaquette variables Uppq.These are the fundamental squares of the lattice consisting of four gauge links U µ pxq " expp´igaA µ pxqq, which are group elements of SUpN c q.For sufficiently small lattice spacing the plaquette can be expanded, and the classical action ( 1) is seen to reproduce the continuum Yang-Mills action in the limit a Ñ 0. However: for any finite lattice spacing a the action differs by terms of Opa 2 q from the continuum action.Such terms disappearing in the continuum are called "irrelevant", at finite lattice spacing they constitute lattice artefacts.An improvement of lattice actions in the sense of reducing such artefacts [8] can be achieved by adding irrelevant terms to the lattice action, which subtract those already present.Quite generally, lattice actions are then not unique but may differ by arbitrary irrelevant terms.All such differences between valid actions and observables must disappear in the continuum limit.
Approaching the continuum limit is a highly non-trivial affair, as one has to take a Ñ 0 while keeping other dimensionful quantities, like masses or temperature, fixed along a socalled line of constant physics in the lattice parameter space.The lattice spacing is tuned indirectly through the running gauge coupling evaluated at the cutoff scale, g " gpΛ " a ´1q, so that by asymptotic freedom gpa Ñ 0q Ñ 0. This implies that the lattice gauge coupling βpa Ñ 0q Ñ 8.For thermodynamics, keeping temperature constant during the continuum approach implies N τ Ñ 8 in the continuum limit.Thus, while the precise value of the lattice spacing is determined by some form of scale setting, finer thermodynamics lattices imply larger N τ , and lattice artefacts can be parametrised in units of temperature as powers of the dimensionless aT " N ´1 τ .The most problematic lattice artefacts for QCD by far are those afflicting the fermion formulation.In particular, it is impossible to put chiral fermions on a lattice without artificial doubler degrees of freedom while maintaining locality [9].It is clear that this implies severe caveats for any investigation of chiral symmetry breaking.In a nut shell, the current choices are:

•
The Wilson formulation: By adding irrelevant mass terms " a ´1, the doubler degrees of freedom become heavy and decouple in the continuum limit.However, for any finite lattice spacing chiral symmetry is broken completely by these terms.

•
The staggered formulation: Spin and flavour degrees of freedom are distributed to different lattice sites, which can be done so as to reduce the number of doublers.The remaining ones are removed by taking an appropriate root of the fermion determinant, which can only be fully valid in the continuum.In this formulation the original chiral symmetry is reduced from SUp2q L ˆSUp2q R Ñ Op2q.

•
Formulations with a lattice version of full chiral symmetry exist in the form of domain wall or overlap fermions, and are expected to eventually supersede the previous formulations.However, they require complicated non-local constructions and currently are more expensive to simulate by over an order of magnitude.
Due to the enormous numerical cost of thermodynamical investigations for light fermions, most lattice results covered here are based on different versions of the staggered or Wilson discretisations.In order to avoid the damage these do to chiral symmetry and its spontaneous breaking, it is mandatory to first take the continuum limit to remove the lattice artefacts, and only then approach the chiral limit.
A general difficulty in investigating phase diagrams is the fact that non-analytic phase transitions only exist in infinite volume [10,11], while numerical simulations are obviously limited to finite boxes.A change of dynamics associated with a phase transition is signalled by a rapid change of the expectation values of suitable observables Opxq, accompanied by maximal fluctuations, i.e., peaks in their susceptibilities, This allows to determine the location of a transition as a function of the bare parameters of the theory.In a finite box this is always an analytic crossover, and the behaviour on different volumes has to be compared to see whether a non-analytic phase transition builds up in the thermodynamic limit.The peak height will diverge at a second-order transition with " V σ , where σ is a combination of critical exponents appropriate for the observable and the universality class of the transition.For a first order transition the divergence is " V, whereas for a crossover the peak height will saturate at a finite value.More intricate (and reliable) finite size scaling studies are based on the first three generalised cumulants (n " 2, 3, 4 ) of the fluctuation of an observable, δO " O ´xOy, which besides the peak height also quantify the skewness and the tails of the statistical distribution of the observable, thus giving access to various critical exponents.Finally, the infamous sign problem at finite baryon density appears because of a property of the fermion determinant, defined with complex chemical potential for baryon number (in continuum notation), detp D `m ´µγ 0 q " det ˚p D `m `µ˚γ 0 q . (5) Thus, the fermion determinant is complex for real chemical potential, but stays real for purely imaginary chemical potential.For real chemical potential, all imaginary parts cancel out of the partition function exactly, but the real part is oscillatory and averages to the so-called average sign, evaluated by the partition function with µ " 0, Symmetry 2021, 1, 0 4 of 28 This is also the simplest reweighting factor multiplying any observable, when a µ " 0 ensemble is reweighted to µ ‰ 0. Its size is governed by the free energy difference of the systems with and without chemical potential, and is exponentially suppressed for large volumes.Hence all expectation values of observables are lost in the statistical errors sooner or later.It must be stressed that this is only a problem for numerical evaluations by importance sampling, other evaluations and the physics itself have no issue with chemical potential.For more details, see introductions like [3][4][5].So far there is no genuine algorithmic cure for the sign problem, for an overview of earlier attempts see [12] and the proceedings of the annual lattice conferences for updates.

The Columbia Plot at Zero Baryon Density
As a point of departure, consider the thermal QCD transition at zero baryon density, µ B " 0. The nature of this transition for physical quark masses has been known for some time to be an analytic crossover [13].Away from the physical point, the order of the QCD thermal transition with N f " 2 `1 quarks as a function of quark masses is usually specified in a so-called Columbia plot [14], for which two possible realisations, to be discussed below, are shown in Figure 1.

The Deconfinement Transition
Even though it is far from the physical point or the chiral limit, the heavy mass corner is interesting in its own right and promises insights, which will also be useful for light quarks.In the quenched limit QCD reduces to SUp3q Yang-Mills theory in the presence of static quarks, whose propagators are Polyakov loops which represent an order parameter for the global Zp3q center symmetry related to confinement [15].Its spontaneous breaking at some critical temperature proceeds by a first-order deconfinement phase transition [16], whose latent heat has recently also been computed [17].In the presence of dynamical quarks, the center symmetry is explicitly broken by 1{m q and the first-order phase transition weakens with diminishing quark masses, until it disappears along a second-order critical line in the 3d Zp2q universality class.Contrary to the chiral limit, this entire region can be simulated directly and a complete non-perturbative understanding should be attainable in the near future.The location of the critical line is under investigation [18][19][20][21].While continuum extrapolations are not yet available, the presently finest lattices predict the pseudo-scalar mass evaluated on the critical point for N f " 2 to be about m c PS " 4 GeV [20].Together with the latent heat, these quantities should permit a detailed understanding of the confinement-deconfinement transition and its interplay with screening by dynamical quarks, which will be valuable for the physics interpretation of the crossover region at the physical point.Finally, these quantities can serve as benchmarks to test and tune truncations involved in functional renormalisation group or Dyson-Schwinger approaches [22], or effective lattice theories for the heavy quark region, which are employed to describe finite density physics [23,24].

Figure 1.
Possible scenarios for the nature of the thermal QCD transition as a function of the quark masses, with a first-order (left) or a second-order (right) transition in the chiral limit of the light quarks.Every point of the plot represents a phase boundary, with an implicitly associated (pseudo-)critical temperature T c pm u,d , m s q. (Here and in the following, "phase boundary" refers to a (pseudo-)critical combination of parameters irrespective of the nature of the transition, which can be first order, second order or crossover.)

The Chiral Transition at Zero Baryon Density
The situation is far more difficult in the opposite limit of massless quarks, because it cannot be simulated directly.For several decades expectations have been based on an analysis of three-dimensional sigma models, augmented by a 't Hooft term for the Up1q A anomaly, which represent Landau-Ginzburg-Wilson effective theories for the chiral condensate as the order parameter of the transition.The renormalisation group flow based on the epsilon expansion [25] predicts the chiral phase 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 Up1q 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 first-order.For non-zero quark masses, chiral symmetry is explicitly broken.A first-order chiral phase transition then weakens to disappear at a Zp2q second-order critical boundary, just like the deconfinement transition in the heavy mass regime, whereas a secondorder transition disappears immediately.This results in the two scenarios for the Columbia plot depicted in Figure 1.The results of the epsilon expansion were essentially confirmed by numerical simulations of three-dimensional sigma models [26] and a perturbative expansion in fixed dimensions [27].A later high-order perturbative analysis for the case with an effectively restored Up1q A at the transition temperature finds a possibility for second-order transitions also in this case, but with a symmetry breaking pattern Up2q L b Up2q R Ñ Up2q V signalling a different universality class [28].Finally, a recent functional renormalisation group analysis applied to Up1q A restoration in QCD, which is in good agreement with lattice susceptibilities at finite quark masses [29], favours the Op4q scenario in the chiral limit [30].
A large number of lattice simulations has been devoted to disentangle which of these situations is realised in QCD.When the Columbia plot is considered on the lattice with different lattice spacings, its parameter space is enlarged by an additional axis.The chiral critical line then traces out a chiral critical surface whose shape is discretisation-dependent, while for all valid discretisations it must of course emerge from the same critical line in the continuum limit.The following subsections collect the current evidence on the location of the chiral boundary line in the Columbia plot.
Early simulations (without sophisticated finite size scaling analyses) on coarse N τ " 4 lattices were fully consistent with the expected scenarios shown in Figure 1: A first-order region could be clearly seen for N f " 3, whereas the smallest available masses were consistent with a continuous transition or crossover for N f " 2, both for unimproved staggered [14] and unimproved Wilson [33] fermions.More recently and using finite size scaling of cumulants, a narrow first-order region was also identified for N f " 2, again for unimproved staggered [34,35] and unimproved Wilson [36] fermions on N τ " 4.However, the location of the Zp2qboundary varies widely between these, indicating large cutoff effects.
Recent investigations using the Highly Improved Staggered Quark (HISQ) action start at the N f " 2 `1 physical point and then gradually reduce the light quark mass value, until they correspond to m PS « 45 MeV, on lattices N τ " 6, 8, 10 [31,32].No sign of a first-order transition is detected.Thus either a Zp2q-critical point bounding a very narrow first-order region is approached, or a second-order transition in the chiral limit.The analysis employs a renormalisation group invariant combination of chiral condensates as order parameter representing the magnetisation-like variable, and the light quark masses m l " m u,d in units of the strange quark mass as symmetry breaking field, Near a critical point the magnetic observables are dominated by universal scaling functions with a scaling variable z " t{h 1{βδ expressed in terms of the reduced temperature and external field, t " t 0 pT ´T0 c q{T 0 c , h " H{h 0 , which contain the unknown critical temperature in the chiral limit, T 0 c , and two non-universal parameters t 0 , h 0 .: Chiral susceptibility at the physical strange quark mass for a range of decreasing light quark masses on N τ " 8 lattices with HISQ action, from [31].(Right): Magnetic equation of state at the pseudo-critical temperature approaching the chiral limit.Lines represent fits to OpNq-(N " 2, 4 second-order scenario) or Zp2q-scaling with a finite critical quark mass (first-order scenario), from [32].
Figure 2 (left) shows the chiral susceptibility as a function of temperature for a series of decreasing light quark masses on N τ " 8.A considerable reduction of the pseudo-critical temperature, defined by the peak location of the susceptibility, is apparent as the light quark mass is reduced.Moreover, the peak height is growing as expected when approaching a true phase transition.To test for the universality class of the approached transition, the ratio M{χ M is shown in Figure 2 (right), which near a critical point is dominated by universal scaling functions.While it is difficult to distinguish between Op4q and Zp2q, it is apparent that any finite critical quark mass bounding a first-order region has to be excessively small to be consistent with the observed behaviour.
There is a similar scaling expression for the approach of the pseudo-critical crossover temperature to the critical temperature in the chiral limit, T c pm l " 0q " 132 `3 ´6 MeV .(10) In [31] the variation between the possible sets of critical exponents is observed to be very small, so that an extrapolation makes sense even without knowledge of the true universality class.Moreover, the extrapolations were checked to be stable under an exchange of the order of the continuum and chiral extrapolations, leading to the critical temperature as a first result for the chiral phase transition in the massless limit.Simulations with Wilson fermions do not yet reach such small quark masses, but are fully consistent with this picture down to the physical pion mass.In particular, extrapolations of the pseudo-critical temperature to the chiral limit have also been attempted by an investigation using N f " 2 `1 `1 Opaq-improved twisted mass Wilson fermions, with strange and charm quark masses close to their physical values [37].In these simulations lattice spacings are held fixed in a range a{fm P r0.062, 0.082s, and temperature is varied by N τ P r12, 20s.The pseudo-critical temperature is defined in three different ways: By the maximum of the chiral susceptibility χ as above, and as the inflection point of polynomial fits to the chiral condensate (∆) and a subtracted version without additive renormalisation p∆ 3 q.The resulting temperatures as a function of pseudo-scalar (pion) mass are shown in Figure 3 (left), together with extrapolations to the chiral limit, employing either Op4q exponents or Zp2q-exponents and a critical pseudo-scalar mass up to m π " 100 MeV.Again, it is not possible to distinguish between these scenarios.As in the previous case, the extrapolated critical temperature in the chiral limit is therefore robust under changes of the critical exponents and quoted as in remarkable agreement with the staggered result.Figure 3 (right) shows an investigation of sections of the chiral critical line using Opaq clover-improved Wilson fermions [38].Starting point are the data for N f " 3 to be discussed separately below, and on N τ " 6 further points at larger strange quark masses have been added.The critical line is then fitted assuming a tricritical strange quark mass as explained in Section 3.5 plus polynomial corrections.Note that this discretisation features a much wider first-order region, which even contains the physical point on the coarser lattices.This must be a lattice artefact, and the first-order region rapidly shrinks as N τ is increased.
Pseudo-critical temperature of the crossover defined by the chiral susceptibility χ, the inflection point of the chiral condensate ∆ or an additively renormalised chiral condensate ∆ 3 , for N f " 2 `1 `1 twisted mass Wilson fermions close to the continuum.Lines represent chiral extrapolations according to the Op4q second-order or finite critical Zp2q-mass scenario.From [37].
(Right): Columbia plot expressed in η, π-masses in units of the Wilson flow parameter t 0 .Critical points have been determined using an Opaq-improved Wilson action.The first-order region includes the physical point on coarse lattices, but shrinks drastically as N τ is increased.From [38].
Several conclusions can be drawn from these results.Firstly, the width of a potential first-order region as in Figure 1 (left) is bounded to a small fraction of the physical light quark (or pion) masses.Second, the numerical proximity of the critical exponent combinations 1{pβδq for the 3D Op2q, Op4q and Zp2q universality classes appears to allow for a robust extrapolation of the chiral transition temperature to the massless limit with remarkably small uncertainties.Conversely this statement implies, however, that it is impossible to firmly identify the universality class in this way, which would require exponentially accurate data.This problem might be avoided by looking at the scaling of energy-like variables, which are governed by the critical exponent α that changes sign between the Op2q, Op4q and the Zp2q universality classes.It was shown that the Polyakov loop behaves as an energy-like observable, but unfortunately a firm distinction between universality classes would require a further substantial reduction of the light quark mass [39].Finally, note that the value of T c pm l " 0q is " 25 MeV lower than the pseudo-critical temperature at the physical point.This constrains the possible location of a chiral critical point at finite baryon density, as will be discussed in Section 5.1.

N f " 3 and N f " 4
The N f " 3 theory, represented by the diagonal of the Columbia plot, Figure 1, is particularly interesting because a sufficiently large first-order transition and its associated Zp2qcritical endpoint are accessible by direct simulations, rather than by extrapolation.However, the situation has become more complicated since the early results, as shown by the compilation of critical pseudo-scalar masses obtained by different actions and lattice spacings in Table 1.The largest value in the table differs by almost an order of magnitude from the lowest bound!Unless any of the employed lattice actions is fundamentally flawed, all of them must eventually converge to the same continuum limit, which is not in sight yet.However, the general trend is for the critical mass to shrink when either a lattice with fixed action is made finer (larger N τ ), or when improved actions are employed.This points to enormous cutoff effects, which quite generally appear to increase the first-order region, i.e., make the transition stronger.Another observation is for the first-order region to be considerably larger in the case of Wilson fermions, which is consistent with stronger cutoff effects due to the complete violation of chiral symmetry.As an example, consider Figure 4 ( left), where the critical pseudo-scalar mass is shown from a comprehensive long term study with Opaq-improved Wilson fermions for different lattice spacings.The critical meson mass shrinks by nearly a factor of two when the lattice is changed between N τ " 6 and N τ " 12.Its continuum value depends strongly on the extrapolation function and could even be zero.The same effect is observed for N f " 4 unimproved staggered fermions [48], as shown in Figure 4 (right).One expects the transition in this case to be more strongly first-order than for N f " 3, since the spontaneously breaking symmetry is larger.This is borne out both for Wilson [49] and staggered fermions.Note that for N f " 4 staggered fermions no rooting is applied, which therefore cannot cause any unphysical behaviour.Moreover, the effect of reducing the lattice spacing is similar for N f " 2, 3, 4 staggered fermions, i.e., with and without rooting.
We can then conclude that Wilson and staggered discretisations show the same qualitative behaviour, with a monotonic decrease of the critical pseudo-scalar mass m c PS pN τ q as the continuum is approached by increasing N τ .Quantitative contradictions between different actions are only avoided, if the final answer m c PS p8q is below the lowest available value m PS pN τ q.In summary, since finer lattices and improved actions have become available, the situation for N f " 3, 4 is quite similarly ambiguous to the previously discussed N f " 2 and N f " 2 `1.

Three-State Coexistence and Tricritical Scaling
The uncertainty about the size of a potential first-order region can be removed by a different analysis, which leads to a consistent picture covering the results of all discretisations reported here.Starting point is the observation that a change of the chiral phase transition in the massless limit from first order to second order, brought about by variation of a continuous parameter, necessarily proceeds by a tricritical point [50].To see this, assume the scenario in Figure 1 (right), and consider the corresponding phase diagram for a fixed m s ă m tric s , Figure 5 (left).For non-vanishing quark masses, the first-order transition weakens to disappear in a critical endpoint, and the phase diagram is symmetric under a chiral rotation m u,d Ñ ´mu,d .Hence, in the massless limit the critical temperature marks a triple point of three-state coexistence: Two states with broken chiral symmetry, ˘x ψψy ‰ 0, and one with chiral symmetry restored, x ψψy " 0. When m s is increased, these points trace out a triple line (corresponding to the green m u,d " 0 line in Figure 1) along which the latent heat of the first-order transition decreases, until it vanishes in a tricritical point.At the same time, the two wings of finite mass first-order transitions shrink, so that the tricritical point marks the confluence of two critical endpoints.
This opens an additional path of investigation.Quite generally, tricritical points come with their own set of critical exponents, which take mean field values since their upper critical dimension is three.For a general introduction to the theory of tricritical scaling, see [51].In particular, the functional form of the critical wing lines entering the tricritical point as a function of the symmetry breaking field, which is the bare light quark mass in our case, is governed by a critical exponent, This allows for an extraction of the location of a tricritical point with fixed, known exponents, which is much easier than having to determine the exponents themselves.Knowing the location of the tricritical point then fixes the order of the transition on both sides of it, albeit without information about the universality class.Attempts in this direction were made in [38,52], but they are not yet conclusive because the lattices are coarse, whereas very high precision data very close to the chiral limit would be required on sufficiently fine lattices in order to unambiguously distinguish tricritical scaling from an ordinary polynomial.However, this problem can be circumvented by a change of variables in the bare parameter space of the theory.3.6.The Chiral Phase Transition as a Function of N f and N τ In [35] it was suggested to consider a different version of the Columbia plot with degenerate quarks only, in which the variable strange quark mass m s is traded for a variable N f analytically continued to non-integer values.Starting with integer N f , after Grassmann integration the QCD partition function reads which can now be formally viewed as a statistical system of gauge field variables depending on a continuous parameter N f .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 Figure 1 (right) then translates to the version shown in Figure 5 (right), 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.At finite lattice spacing there is an additional parameter, and the Zp2q-critical boundary traces out a chiral critical surface which terminates in a tricritical line in the chiral limit of the lattice theory.Since there is ample evidence for firstorder transitions at coarse lattice spacings, the task is now reduced to track their boundaries in the lattice bare parameter space and determine the location of the tricritical line in the lattice chiral limit.
Numerical results obtained with unimproved staggered fermions for a wide range of flavours and N τ " 4, 6, 8 [35,53] are shown in Figure 6.The left plot represents the lattice version of the Columbia plot Figure 5 (right).One observes a summary of the previous discussion, namely a first-order transition getting stronger with N f and with increasing lattice spacing.The cutoff effects are growing with N f , for N f " 7 the critical bare quark mass is reduced by a stunning factor of five when the lattice spacing is reduced (N τ is increased) by a factor of two!No continuum critical mass m c is discernible yet for any N f .On the other hand, the intercepts of the critical lines at m " 0 constitute the sought-after tricritical line N tric f pN τ q.
Tricritical scaling can be appreciated in the rescaled Figure 6 (right), where the data approach a leading order scaling relation [35].This allows for an extrapolation The chiral critical surface for unimproved staggered fermions, projected onto the pm{T, N f q plane [53] (left) and onto the pam, N f q plane for N τ " 4 [35] (right).Every point of the plane represents a phase boundary with an implicitly tuned βpam, N f , N τ q.Regions above the Zp2q-lines represent crossover, those below first-order transitions.Fits on the left show next-to-leading order tricritical scaling, on the right leading-order tricritical scaling.
The scaling relation has been inverted here, because in simulations N f is fixed whereas the critical quark mass is computed with statistical errors.Note also that the N f " 2 data point in Figure 6 (right) has been obtained by a tricritical extrapolation in imaginary chemical potential at fixed N f " 2 [34].This is an independent confirmation of the bare quark mass as a tricritical scaling field near its chiral limit.Scaling fits in the left figure take into account the next-toleading order terms to predict N tric f pN τ " 4q « 1.71p3q and N tric f pN τ " 6q « 2.20p8q.While no extrapolation for N τ " 8 is available yet, the critical line further flattens and appears to shift to the right, implying a line N tric f pN τ q.One concludes that, for unimproved staggered fermions, the N f " 2 massless theory shows a first-order transition on N τ " 4, but a second-order transition on all finer lattices and in the continuum.The question now is what happens to the N f ě 3 theories.
This can be seen in a different variable pairing in Figure 7 (left), where the same critical bare quark masses are rescaled by the tricritical exponent, and plotted as a function of N τ for different fixed N f -values.Only a slight curvature is exhibited by those N f -lines with three data points, which thus are compatible with next-to-leading order scaling and a tricritical point at some finite N tric τ pN f q, obtained by the corresponding extrapolations.

Tricritical Scaling for Wilson Fermions
What about Wilson fermions, which exhibit the widest first-order region among the discretisation schemes discussed here?In this case the conceptual situation is more difficult, because chiral symmetry is broken completely so that the bare quark mass and chiral condensate receive additive renormalisation as well as the ordinary multiplicative one.Moreover, the bare parameter phase diagram is expected to show lattice artefacts such as parity broken phases [54,55] or first-order bulk transitions [56][57][58], whose location depends on the details of the lattice action and N τ .This may considerably complicate the bare parameter phase diagram while, of course, all actions should eventually agree in the continuum limit.
Here we consider an RG-improved Wilson gauge and non-perturbatively O(a)-improved Wilson fermion action, for which the Aoki phase is limited to the strong coupling regime and no other unphysical phases have been observed [59], and which was used for the sequence of N f " 3 simulations [45][46][47] displayed in Figure 4 (left).These data have recently been reanalysed in order to test for tricritical scaling [53], in analogy to the staggered case above.
To this end, one can exploit that in chiral perturbation theory which also holds for Wilson fermions provided that an additively renormalised quark mass is used.The data from Figure 4 are replotted in Figure 7 (right), with the vertical axis rescaled to be proportional to the renormalised quark mass as scaling variable.Perfect triritical scaling is observed, with extrapolations that are robust under variation between leading and next-toleading order, as well as when using only the coarsest lattices N τ " 4, 6, 8. Nf " 3 LO Nτ P r8, 12s NLO Nτ P r6, 12s NLO Nτ P r4, 8s Figure 7.The chiral critical surface projected onto the pam, aT " N ´1 τ q plane for unimproved staggered fermions [53] (left) and for Opaq-improved Wilson fermions [47] (right), with tricritical scaling fits to both.From [53].

Conclusions for the Continuum Limit
Knowing the structure of the bare parameter lattice phase diagram, its implications for continuum physics can be stated.The continuum limit is represented by the origin of the pam, N ´1 τ q plane, which can be approached along different lines of constant physics, representing theories with different values for, e.g., the vacuum pseudo-scalar mass.For any N f with a finite N tric τ pN f q, the parameter region of first-order transitions is not continuously connected to the continuum limit, and thus represents a mere cutoff effect.In other words, for all N τ ą N tric τ , simulations will only show crossover behaviour for any vacuum pion mass, such that the transition in the continuum chiral limit is approached from the crossover, and corresponds to an isolated second-order point.By contrast, a first-order transition in the chiral continuum limit implies a finite m c in physical units, which would require the Zp2q-critical line am c pN τ q to terminate in the origin of the plot.Moreover, there would be no tricritical point and no associated tricritical scaling, so that am c pN τ q should have an ordinary Taylor expansion about zero with the usual cutoff corrections " a, a 2 etc.In [53] it was found that polynomial next-to-leading order fits are unable to accommodate the N τ " 4, 6, 8 data for N f " 5, 6, 7 unimproved staggered fermions.For the N f " 3 Opaq-improved Wilson fermions such fits are possible but have worse χ 2 {dof than the tricritical fits.We must then conclude that the staggered action predicts all massless theories with N f ď 7 to have second-order transitions in the continuum limit, and that the Wilson clover-action agrees with this for N f " 3. Note that this conclusion is perfectly consistent with the absence of any first-order region in all simulations with improved staggered actions so far.
In summary: • Qualitatively, the discretisation effects on the chiral phase transition are the same for unimproved staggered fermions and either unimproved or Opaq-improved Wilson fermions, making the transition stronger compared to the continuum.

•
In both discretisations, the boundaries of the first-order transitions at finite lattice spacings exhibit tricritical scaling and extrapolate to a finite N tric τ pN f " 3q.This implies that the first-order region is not connected to the continuum limit.Thus, when the continuum limit is taken before the chiral limit, as is necessary to avoid lattice artefacts, both predict a second-order transition.• Quantitatively, the cutoff effects are larger for Wilson fermions, resulting in a larger N tric τ pN f q than in the case of staggered fermions.This might be explained by the fact that the respective discretisations break chiral symmetry fully or only partially.
If these conclusions withstand the test of further investigations, they require a modification of the Columbia plot, as shown in Figure 8, with a second-order transition all along the m u,d " 0 limit, independent of the strange quark mass.Note that the universality class of these second-order transitions remains open.Since the chiral symmetry is different for the N f " 2, 3 boundary cases, one might expect critical exponents to smoothly cross from one set of values to the other along the line.Numerical simulations are no mathematical proof, and the reported results are in contradiction with the predictions from [25,27].Given the presently available data, what would it take to avoid this conclusion and render the chiral phase transition for N f ě 3 first order?This would require future data points for am c pN τ q (or am PS pN τ q) on larger N τ to extrapolate to zero as a polynomial in N ´1 τ , i.e. without tricritical scaling, independent of the discretisation scheme.The tricritical scaling observed for the presently available N τ would then have to be "accidental".On the other hand, the question arises whether the three-dimensional sigma models investigated in [25,27] represent the most general case compatible with QCD.It is well-known that φ 4 -theories do not permit tricritical points, which requires φ 6 -terms.In three dimensions, φ 6 theories are renorrmalisable and moreover show stable infrared fixed points in non-perturbative studies [60,61], which could reconcile effective theories with the current lattice results.It should be possible to resolve these questions in the near future and finally settle the fate of the thermal QCD phase transition with massless quarks.

The Columbia Plot with Chemical Potential
Once a chemical potential µ for quark (or µ B " 3µ for baryon) number is switched on, the Columbia plot in the continuum (or at fixed lattice spacing viz.N τ ) gets extended into a third dimension.In this case the Zp2q-critical boundary lines sweep out critical surfaces, as sketched schematically in Figure 9 (left).Labelling the third axis by pµ{Tq 2 , both real and imaginary chemical potential can be discussed.Imaginary chemical potential, µ " iµ i with µ i P R, is unphysical, but it does not induce a sign problem and thus can be simulated without difficulty [62].This can be exploited to extract various properties of the phase diagram at sufficiently small real µ by analytic continuation [63,64].Two exact symmetries are important here.Because of CP-invariance, the QCD partition function is an even function of the chemical potential Zpµq " Zp´µq.Furthermore, for arbitrary fermion masses it is periodic in imaginary chemical potential because of the global Roberge-Weiss (or center) symmetry [65], Thus, as imaginary chemical potential is increased, the partition function periodically cycles through the N c center sectors, which are distinguishable by the phase of the Polyakov loop, whereas the thermodynamic functions are invariant under a shift between sectors.The boundaries of these sectors, located at pµ{Tq c " ˘ip2n `1qπ{3 with n " 0, 1, 2, . .., are marked by first-order transitions for high temperatures and crossover for low temperatures, see Figure 9 (right).Varying T at fixed imaginary chemical potential, the analytic continuation of the QCD thermal transition is crossed, whose order depends on N f and the quark masses, as discussed in the previous sections.For a first-order chiral or deconfinement transition, the transition lines meet up in a triple point, which marks a three-phase coexistence between two phases of the Polyakov loop and their average.For a thermal crossover the Roberge-Weiss transition ends in a critical endpoint with 3D Ising universality.The boundary between these cases, corresponding to specific quark mass values, is marked by a tricritical point.In the Columbia plot for the first Roberge-Weiss plane at µ{T " iπ{3, these trace out tricritical lines, in which the deconfinement critical surface ends in the Roberge-Weiss plane [66].These structures have been established explicitly for unimproved staggered [66,67] as well as unimproved Wilson [68] fermions and generalise to any discretisation respecting the center symmetry.
The thermal phase transition at imaginary µ < l a t e x i t s h a 1 _ b a s e 6 4 = " O 5 q q L o 1 E t j G q V E E X O K R l 0 E S / y Y w = " > A A A B 6 n i c d V B N S w M x E M 3 W r 1 q / q h 6 9 B I v g a c l u q 2 1 v R S 8 e K 9 p a a J e S T b N t a J J d k q x Q S n + C F w + K e P U X e f P f m G 0 r q O i D g c d 7 M 8 z M C

The Deconfinement Transition
Once again it is instructive to also consider the more easily accessible heavy mass corner, and to study how the deconfinement transition changes in the presence of a baryon chemical potential.Moreover, for heavy quarks one can expand the fermion determinant in inverse quark mass and, after an additional resummed strong coupling expansion, an effective lattice theory can be constructed that reproduces the Zp2q-boundary on N τ " 4 lattices at µ " 0., but in addition can be solved at finite chemical potential, both imaginary and real [23].The tricritical line in the Roberge-Weiss plane dictates the deconfinement critical surface to emanate from it by tricritical scaling, where now the deviation of the chemical potential from the Roberge-Weiss plane is the center symmetry breaking scaling field.This scaling behaviour is confirmed numerically and, even at leading order, found to extend far into the range of real chemical potentials.The region of first-order deconfinement transitions is thus found to continuously shrink with real chemical potential [23,66].This observation implies that the chiral and deconfinement transition regions remain separated for all chemical potentials, i.e., the respective critical surfaces do not join.The chiral critical surface must then form a closed volume with the boundaries of the quark mass and chemical potential parameter space, since one cannot get from a first-order transition to a crossover without passing through a critical point.

The Chiral Transition
For unimproved Wilson and staggered discretisations on N τ " 4, 6, the 3D Columbia plot looks as in Figure 9 (left), with the region of chiral phase transitions getting wider in the imaginary µ direction.Expanding the Zp2q-critical light mass for small chemical potential for any fixed m s , one may conclude that c 1 ă 1 and by analytic continuation the first order region shrinks in the real-µ direction.As an independent check on these calculations at imaginary chemical potential, the curvature of the chiral critical surface can be computed directly at µ " 0. For staggered fermions one finds c 1 ă 0 both at m s " m l and at m s " m phys s , and the next coefficient is negative as well [40,52,69].Thus the chiral transition strengthens with imaginary and weakens with real chemical potential.This is opposite to a scenario with a chiral critical point close to the temperature axis, which would require the chiral transition at the physical point to strengthen with real µ.The only investigation finding a positive curvature so far is for Opaq-improved Wilson fermions [70].However, we have seen in the previous sections the enormous cutoff effects on the location of the chiral critical surface at µ " 0, and the question is how these vary with µ for a given discretisation.The fact that for µ " 0 the critical surface moves to the chiral limit makes it practically impossible to compute its curvature on fine lattices.
Investigations in the Roberge-Weiss plane on finer lattices reveal the same trend as seen at µ " 0, namely the chiral tricritical line moving towards smaller quark masses, both for unimproved staggered [71] and Wilson [72] quarks, Figure 10 (left).For stout-smeared staggered [73] and HISQ [74] actions on N τ " 4, even the larger first-order region in the Roberge-Weiss plane cannot be detected when starting from the physical point and reducing the pion masses down to m π « 50 MeV, as Figure 10 (right) shows by demonstrating secondorder scaling for the Roberge-Weiss endpoint.[24] [11] and this study This study [4] [28] Ref. lattice spacing, appear to have considerably larger cut-off effects.For example, comparing am tric ⇡, heavy " 2.2302p2q from Ref. 24 with our am tric ⇡, heavy " 1.7260p3q, the pionresolution problem is milder in the present study.It is also interesting to compare the position of the tricritical points in physical units,

Unimproved Wilson
and The large differences between discretizations again imply being far from the continuum limit, where results from all discretizations have to merge.The observed trend is consistent with the findings of simulations with improved staggered actions, where the tricritical points can only be bounded to be at much smaller masses, as indicated in Figure 8, as well as with the analogous findings at zero chemical potential (see discussion in the introduction).In particular the comparison across discretizations implies enormous cut-off effects in the critical masses, which could end up being over "100% of an eventual continuum limit.We remark that cut-off effects in the critical temperatures are much milder.At present, there is no theoretical explanation as to why the discretization effects on critical quark masses in the Columbia plot are so strong.
In conclusion, we have determined the shift of the tricritical points in the Roberge-Weiss plane of unimproved staggered fermions by changing from N⌧ " 4 to N⌧ " 6 lattices.The aspect ratios and statistics required to extract the correct order of the phase transition are found to be larger in the Roberge-Weiss plane than at µ " 0. We find the cut-off effect on the tricritical masses to be smaller but qualitatively the same as that observed with Wilson fermions, and consistent with results for both discretizations at zero chemical potential.This implies in particular, that the entire chiral critical surface depicted in Figure 1 is shifted significantly towards smaller (and possibly zero) light quark masses, as the lattice spacing decreases, which is also consistent with results from improved staggered actions.Unfortunately, our study also implies that much finer lattices at inevitably smaller quark masses are necessary, before one can hope the results of the light tricritical mass to stabilize in a continuum limit.[71].(Right): Finite size scaling with 3d Zp2q exponents for stout-smeared staggered fermions with m PS « 50 MeV and pm l {m s q phys on N τ " 4 in the RW-plane [73].
Note that the order parameter for the center sectors is the imaginary part of the Polyakov loop, even in the light quark regime.It is then interesting to study the interplay between the center and chiral dynamics in the Roberge Weiss plane, since they are independent symmetries.In particular, in the chiral limit there must be a true phase transition for any chemical potential, since the latter does not break chiral symmetry.Figure 11 shows the renormalisation group invariant chiral condensate from Equation ( 7) and the disconnected part of the corresponding susceptibility as functions of temperature, as observed with HISQ fermions on N τ " 4 [74].Indeed, as the quark mass decreases a chiral transition is building up, just as at µ " 0, and it is particularly interesting where the transition happens.The location of the emerging peak in the susceptibility is, at the current accuracy, consistent with the yellow bar in the figure marking the critical temperature of the Roberge-Weiss endpoint, so that the transition lines in the chiral limit would indeed be connected as in Figure 9 (left).Note that there is no a priori reason for this, since we have independent center and chiral symmetries, with a Zp2q universality of the Roberge-Weiss endpoint and a larger, OpNq-like universality along the chiral transition line.If these lines really do meet in the chiral limit, their junction should be multicritical.In summary, from the µ " 0 and imaginary chemical potential results together we conclude that the entire chiral critical surface is shifting drastically towards the chiral limit as the lattice spacing is decreased, and it is an open question whether any first-order transition remains in the continuum limit at any µ.At the same time, this implies that the crossover at the physical point for µ{T " µ B {p3Tq ă " 1, which is where an analytic continuation of Equation ( 18) holds (see Section 5.1), is very soft and does not display any singular behaviour due to the chiral transition.

QCD at the Physical Point
Having discussed the order of the thermal QCD transition in a generalised parameter space, the aim is to see how this constrains the physical quark mass configuration.The qualitative relation between the phase diagram in the chiral limit of the light quarks and their physical values has been worked out some time ago by mean field methods [75,76], as shown in Figure 12.In the chiral limit, there must be a non-analytic phase transition for any value of the baryon chemical potential, since the latter respects chiral symmetry.We have seen that lattice results are indeed converging towards a second-order transition at µ B " 0. By contrast, the first-order nature of the chiral transition at T " 0 is based on expectations from NJL, Gross-Neveu and sigma models in mostly mean field investigations, with no firm QCD predictions available yet.If this assumption holds, there must be a tricritical point marking the change between first-order and second-order behaviour, which might even have an influence on physical QCD if the light quark masses happen to be in the scaling window of the tricritical point [77].At finite quark masses chiral symmetry is broken explicitly, the second-order line turns into crossover and the first-order transition weakens to disappear in a Zp2q-critical line, which emanates from the tricritical point by tricritical scaling [76], This implies an ordering of the critical temperatures to be exploited below, For completeness, we need to also discuss an alternative scenario, where the chiral phase transition in the massless limit is second order all the way to T " 0. At least from a lattice perspective, this is not excluded so far, but crucially depends on whether there is any nontrivial m c u,d pµq-dependence in the continuum limit.Moreover, a recent investigation of the chiral nucleon-meson and chiral quark-meson models finds the phase transition for m " 0 at T " 0 to turn second order, once fluctuations are included [78].In such a scenario there is no tricritical point and no first-order transition anywhere.Instead, non-vanishing quark masses remove the entire second-order line and the chiral transition would be analytic crossover exclusively for physical quark masses.[75,76].(Right): If the entire chiral transition line in the massless limit is of second order, the transition at the physical point is crossover everywhere.

The Crossover at Small Baryon Densities
There are several methods that have been used so far to extract information about the phase structure at the physical point for small baryon density.All of them introduce some approximation which can be controlled as long as µ{T ă " 1: (i) Reweighting [79], (ii) Taylor expansion in µ{T [80] and (iii) analytic continuation from imaginary chemical potential [63,64].When the QCD pressure is expressed as a series in baryon chemical potential, the Taylor coefficients are the baryon number fluctuations evaluated at zero density, which can also be computed by fitting to untruncated results at imaginary µ B .This permits full control of the systematics between (ii) and (iii).These coefficients are presently known up to 2n " 8 on N τ " 16 lattices, Figure 13 (left), and in principle also observable experimentally.
For a review of the equation of state relating to heavy ion phenomenology, see [81,82].Note also, that this low density regime appears to be accessible by complex Langevin simulations without recourse to series expansions, albeit not yet for physical quark masses [83].This offers an additional cross check between different methods.8 , calculated within CEM-CD (red stars).Lattice QCD data of Wuppertal-Budapest [20] and HotQCD [18,19] collaborations are shown by the blue and green nds/symbols, respectively.

Baryon number susceptibilities
The baryon number susceptibilities Leading order baryon number susceptibilities at µ B = 0 have recently been computed in lattice QCD [16, , 18, 19, 20].A comparison with these lattice data can test the predictive power of the CEM.
Figure 2 depicts the temperature dependence of B 2 , B 4 / B 2 , B 6 / B 2 , and B 8 , calculated in CEM and mpared to the lattice data of Wuppertal-Budapest [20] and HotQCD collaborations [18,19].The CEM lculations use the Wuppertal-Budapest data [11] for b 1 (T ) and b 2 (T ) as an input and are therefore labeled M-LQCD in Fig. 2. CEM results are in quantitative agreement with the lattice data for B 2 and B 4 / B 2 .e CEM is also consistent with the lattice data for B 6 / B 2 and B 8 , although these data are still preliminary d have large error bars.One interesting qualitative feature is the dip in the temperature dependence of / B 2 , where this quantity is negative.It was interpreted as a possible signature of chiral criticality [21].iven that this behavior is also present in CEM (see red stars in Fig. 2c), i.e. in a model which has no critical int, we conclude that the negative dip in B 6 / B 2 cannot be considered as an unambiguous signal of chiral iticality.

Reconstructing the Fourier coe cients b 1 and b 2 from susceptibilities
All baryon number susceptibilities at a given temperature are determined in the CEM by two parameters the leading two Fourier coe cients b 1 and b 2 .One can now consider a reverse prescription -assuming e validity of the CEM ansatz one can extract the values of b 1 and b 2 at a given temperature from two dependent combinations of baryon number susceptibilities by reversing Eq. ( 6).We demonstrate this considering the lattice QCD data of the HotQCD collaboration for B 2 and B 4 / B 2 .The temperature pendence of the b 1 and b 2 coe cients, reconstructed from the HotQCD collaboration's lattice data on e basis of CEM [Eq.( 6)], is shown in Fig. 3 by the green symbols.The extracted values agree rather ell with the imaginary µ B data of the Wuppertal-Budapest collaboration, shown in Fig. 3 26) by the WB collaboration, and reverse engineered using CEM from HotQCD baryon number fluctuations.From [88].

by the blue
An important quantity is the pseudo-critical temperature marking the "phase boundary" between the chirally broken and restored regimes.Since the chiral transition at the physical point corresponds to an analytic crossover with a non-zero order parameter everywhere, there are no truly distinct "phases" and no unambiguous definition of a transition temperature exists.In general, definitions based on different observables will give different pseudo-critical temperatures, even in the thermodynamic limit, contrary to the unique locations of singularities for true phase transitions.While this is an issue when comparing with an experimental situation, for theoretical investigations it is convenient to stick to the observables representing the true order parameter in the appropriate limit, i.e., the susceptibility of an appropriately normalised chiral condensate in this case.Following as an implicitly defined function from the partition function, the pseudo-critical temperature can be similarly expressed as a power series in chemical potential, with T pc p0q " 156.5p1.5qMeV [87].Continuum extrapolated results for the leading coefficient are collected in Table 2, the sub-leading coefficient κ 4 is compatible with zero at the current accuracy.This is a remarkable result telling us that up to µ B ă " 3T the dependence of thermodynamic quantities on chemical potential is rather weak and can be accurately described by a truncated leading-order Taylor series in chemical potential.(20) imag.µ, stout-smeared staggered [85] 0.0145 (25) Taylor, stout-smeared staggered [85,86] 0.016 (5) Taylor, HISQ [87] We now have the necessary information to obtain a conservative bound on the location of a possible critical point, which according to Figure 12

The Search for a Critical Point
For any power series of a function with a given domain of analyticity in its complex argument, the radius of convergence gives the distance between the expansion point and the nearest singularity.This implies that the location pT c , µ c B q of a non-analytic QCD phase transition constitutes an upper bound on the radius of convergence of the series Equation ( 21), or that of any other thermodynamic function.Turning this around one may search for a critical point: If a finite radius of convergence can be extracted from the pressure series for real parameter values, it should signal a phase transition.The simplest estimator is the ratio test of consecutive coefficients, whose extrapolation yields the radius of convergence, In practice, however, only the first few coefficients are available and an extrapolation is not feasible.For a compilation of available results, see [89].Moreover, the ratio estimator is inappropriate for series with irregular signs and complex singularities, where it fails to converge.This can be illustrated by modelling the lattice data in a spirit similar to the hadron resonance gas descriptions, such that higher coefficients become available and different scenarios can be tested for compatibility with the data.As an example, consider the fugacity expansion of baryon number density.At imaginary chemical potential, this is a Fourier series whose coefficients can be computed on the lattice without sign problem, In [90] a cluster expansion model (CEM) was proposed, which takes the first two coefficients as input from a lattice calculation [91], and expresses all higher coefficients in terms of these, The α SB k are T-independent and fixed to reproduce the Stefan-Boltzmann limit.This recursion implies that only two-body interactions are included, and corresponds to a truncated virial expansion which is valid for sufficiently dilute systems.The model now predicts all coefficients b kě3 and allows for an all-order closed expression, All existing lattice data in the crossover regime are reproduced with excellent accuracy, as the examples in Figure 13 (left) show.Moreover, Figure 13 (right) compares the coefficient b 1 computed directly from its defining equation, Equation (26), with a calculation from a combination of χ B i , by inverting the CEM.Note that all orders of the fugacity expansion enter this calculation.The remarkable quantitative agreement is only possible, if both lattice calculations (using different actions) give equivalent results and all coefficients of the CEM are sufficiently close to the true QCD values.With all coefficients of the fugacity expansion available, one can study the radius of convergence of CEM, Figure 14.The ratio estimator fails to converge, because of the irregular signs of higher order coefficients (it works for equal or alternating signs).On the other hand, the Mercer-Roberts estimator, converges and extrapolates to a unique radius of convergence, even for coefficients c n pertaining to different observables, as must be the case for a true singularity.A number of improved estimators with faster convergence properties has been proposed in [92].The red points in Figure 14 (right) show the resulting radius of convergence of the CEM for different temperatures.At high temperatures a singularity is predicted at a distance |µ B {T| « π and T ą T pc , which is quantitatively compatible with the first-order Roberge-Weiss transition in the imaginary µ B direction.At lower temperatures there is no Roberge-Weiss transition and the radius of convergence quickly grows.This is an intriguing result, given that no information from imaginary chemical potential has been used as input.It shows that the general approach to detect a non-analytic phase transition is viable, provided that a good estimator for the radius of convergence is used and that sufficiently many coefficients are available (more than 100 in this case!).The correct identification of the Roberge-Weiss transition by the radius of convergence conversely implies that the CEM of lattice QCD does not have another phase transition closer to the origin, i.e., a possible critical endpoint in the real direction must satisfy This is fully consistent with Equation ( 23) but derived by completely different methods.Of course, modelling higher coefficients in terms of lower ones is not unique.The simplest alternative is a description of the baryon number fluctuations up to χ B 8 by a polynomial model, equally without singularity [93].Another one is provided by a rational function model, which can account for singularities and can be applied both to QCD or to a chiral model with a phase transition [94].Equations of state compatible with lattice data and including an Ising critical point in predefined locations have also been constructed [95].While the properties of models are not those of QCD, these analyses altogether do show that there is no sign of criticality in the real µ B direction from the presently available, continuum extrapolated lattice data at zero or imaginary chemical potential, but at best in their higher order completions.It will thus be interesting to test higher order coefficients without modelling by Padé approximants or conformal maps as demonstrated in a Gross-Neveu model [96], or by resummation schemes in pµ B {Tq [97].
Having to go term by term in an expansion can be avoided, if the radius of convergence is instead determined by the Lee-Yang zero [10,11] closest to the origin.Using reweighting to real chemical potential, this was the strategy employed in the first lattice prediction of a critical point on N τ " 4 lattices using unimproved rooted staggered fermions [98].However, it is now understood that the closest Lee-Yang zero was caused by a spectral gap between the unrooted taste quartets of the zeros, after Taylor expanding the reweighting factor, rather than by a phase transition [99].This is related to the general problem of staggered taste quartets splitting up when they cross a branch cut of a rooted determinant [100].It has then been proposed to redefine the rooted staggered determinant by a geometric matching procedure averaging over the quartets, after which it can be represented as an ordinary polynomial [99].A further development [101] concerns the reweighting procedure, which usually is based on sampling with the phase quenched determinant and reweighting in the phase factor.This has a well-known overlap problem between the sampled and reweighted ensembles.To get rid of this, one may neglect the imaginary part of the determinant in the partition function altogether, which is allowed since it will average to zero.One can then sample with the real part of the fermion determinant and reweight in the sign only, which has no overlap problem.Of course, the sign problem remains and one is still faced with a challenging signal to noise ratio.
Application of these new methods using the stout-smeared action on coarse N τ " 4 lattices appears to signal a Lee-Yang zero at µ B " 2.4T [101], which however is is still far from the continuum.Simulation results on N τ " 6 for the renormalised chiral condensate, x ψψy R pT, µq " ´mud f 4 π " x ψψy T,µ ´x ψψy 0,0 ‰ , (31) are shown in Figure 15, both for real and imaginary chemical potential [102].In the left plot, there is no steepening or narrowing of the chiral crossover yet as the real chemical potential is increased.In the right plot the chemical potential is varied for fixed temperature T " 140 MeV.One observes full compatibility of the reweighted real µ B simulations with the analytic continuation from imaginary µ B simulations, but with smaller errors.In this case, also, there is no sign of a non-analyticity, which is again consistent with the results from other methods.

Conclusions
Even though the sign problem of lattice QCD remains unsolved, there has been considerable progress over the last few years towards phenomenologically relevant constraints on the QCD phase diagram.This is due to a number of complementary paths of investigation, each of them refining their methods and having finer lattices as well as different lattice actions at their disposal.
In the Columbia plot at µ B " 0, the region of first-order deconfinement transitions in the heavy mass corner can be directly simulated, and its Zp2q-critical boundary for N f " 2, while not yet continuum extrapolated, is expected to settle in a region around m PS " 4 GeV.The region of first-order chiral phase transitions seen in various lattice discretisations has been observed to shrink continuously with decreasing lattice spacing.There is numerical evidence for several N f ě 3 that their Zp2q-critical boundary extrapolates to the chiral limit at finite lattice spacing, exhibiting tricritical scaling.Such first-order regions are not connected to the continuum and represent lattice artefacts.If this is confirmed by further studies, the Columbia plot will look like in Figure 8.
In another recent development, the critical temperature for the chiral transition in the massless limit has been determined to be in the range T « 126 ´140 MeV, which is significantly lower than the pseudo-critical temperature at the physical point.According to the expected relationship between the massless limit and the physical quark masses, this bounds a possible critical endpoint at finite baryon density to be at µ B ą " 3T.An estimate of the radius of convergence, based on a cluster expansion model consistent with all available zero density baryon number fluctuations for physical QCD, is fully consistent with this.Novel techniques to extract the leading Lee-Yang zero of the partition function, either with improved reweighting methods or by refined series analyses, should be able to provide additional, independent estimates of the proximity of a phase transition to the physical point in the near future.
Funding: Parts of the work reported here were funded by the Deutsche Forschungsgemeinschaft (German Research Foundation), grant number 315477589-TRR 211 (Strong-interaction matter under extreme conditions).

Figure 2 .
Figure 2. (Left): Chiral susceptibility at the physical strange quark mass for a range of decreasing light quark masses on N τ " 8 lattices with HISQ action, from[31].(Right): Magnetic equation of state at the pseudo-critical temperature approaching the chiral limit.Lines represent fits to OpNq-(N " 2, 4 second-order scenario) or Zp2q-scaling with a finite critical quark mass (first-order scenario), from[32].

Figure 5 .
Figure 5. (Left): Phase diagram with a first-order chiral phase transition.For T ă T c the T-axis corresponds to a coexistence line of ˘x ψψy ‰ 0, and T c pm u,d " 0q represents a triple point.(Right): Columbia plot for mass-degenerate quarks.Again every point represents a phase boundary and has an implicitly associated T c pm, N f q.

Figure 6 .
Figure6.The chiral critical surface for unimproved staggered fermions, projected onto the pm{T, N f q plane[53] (left) and onto the pam, N f q plane for N τ " 4[35] (right).Every point of the plane represents a phase boundary with an implicitly tuned βpam, N f , N τ q.Regions above the Zp2q-lines represent crossover, those below first-order transitions.Fits on the left show next-to-leading order tricritical scaling, on the right leading-order tricritical scaling.

Figure 8 .
Figure 8.The Columbia plot in the continuum, as suggested by the tricritical scaling analyses of unimproved staggered and Opaq-improved Wilson fermions.From [53].
e y t k I y w w s T Y d A o 2 h K 9 P 4 f + k 7 b t e 2 f W v K 6 X G x T K O P D g C x + A U e K A K G u A K N E E L E D A E D + A J P D v c e X R e n N d F a 8 5 Z z h y C H 3 D e P g H J 2 Y 4 k < / l a t e x i t > Chiral critical surface goes smoothly from imag. to rea [de Forcrand, O.P. JHEP 07] Chiral+deconfinement transition weaken with real, str Phys.point "deeper" in crossover region than for zero First-order region in RW plane shrinks towards contin t e x i t s h a 1 _ b a s e 6 4 = " Q 5 Q B o J H X a S X 6 K v H m V 3 P N n a 1 a M J c = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m q o M e i F 4 8 V 7 Q e 0 o W y 2 k 3 b p Z h N 2 N 0 I J / Q l e P C j i 1 V / k z X / j t s 1 B W x 8 M P N 6 b Y W Z e k A i u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b /

11 EndFigure 9 .
Figure 9. (Left): Columbia plot with chemical potential as observed on coarse lattices.The bottom plane corresponds to the first center-transition.(Right): The QCD phase diagram at imaginary chemical potential.Vertical lines mark first-order transitions between different center sectors, the dotted lines are the analytic continuation of the thermal transitions at real µ, whose nature depends on the quark masses.

Figure 8 .
Figure 8. Overview of chiral tricritical values of the pion mass in the Roberge-Weiss plane.

Figure 10 .
Figure 10.(Left): Tricritical pseudo-scalar mass values delimiting the first-order chiral region in the RW-plane[71].(Right): Finite size scaling with 3d Zp2q exponents for stout-smeared staggered fermions with m PS « 50 MeV and pm l {m s q phys on N τ " 4 in the RW-plane[73].

Figure 11 .
Figure 11.(Left): Chiral condensate as defined in Equation (7) in the RW plane from HISQ fermions on N τ " 4. (Right): Disconnected chiral susceptibility at fixed quark mass ratio for various lattice sizes.Yellow bands show the location of the RW endpoint transition temperature, extracted from the peak location of χ L .From [74].

Figure 12 .
Figure12.(Left): Relation of the tentative QCD phase diagram with physical light quark masses (back plane) to the chiral limit (front plane) according to[75,76].(Right): If the entire chiral transition line in the massless limit is of second order, the transition at the physical point is crossover everywhere.
sits on the pseudo-critical line of a strengthening crossover.Using the central value from Equation (10) for the chiral critical temperature and imposing the model-independent ordering T cep ă T c " 132 MeV, the chemical potential of a critical point must satisfy µ cep B ą 3.1 T pc p0q « 485 MeV.(23)

Figure 14 .
Figure 14.(Left): Comparison of estimators for the radius of convergence in CEM.(Right): The resulting radius of convergence as a function of T predicts the RW-transition at imaginary µ B .From [90].

2 (Figure 15 .
Figure 15.Renormalised chiral condensate from simulations with stout-smeared staggered fermions at imaginary (blue) and real (black) chemical potential.(Left): Temperature scan at various chemical potentials.(Right): µ B scan at at T " 140 MeV.Coloured bands result from analytic continuation of imaginary chemical potential data.From [102].