Updates on the Columbia plot and its extended/alternative versions

We report on the status of ongoing investigations aiming at locating the deconfinement critical point with standard Wilson fermions and $N_f=2$ flavors towards the continuum limit (standard Columbia plot); locating the tricritical masses at imaginary chemical potential with unimproved staggered fermions at $N_f=2$ (extended Columbia plot); identifying the order of the chiral phase transition at $\mu=0$ for $N_f=2$ via extrapolation from non integer $N_f$ (alternative Columbia plot).


Introduction
Current findings and/or expectations of the lattice QCD community on the order of the thermal phase transition in QCD as function of the two light (assumed degenerate) quark masses m u,d and the strange quark mass m s are wrapped up in the Columbia plot of which we show in Figure 1 two possible versions in agreement with current findings, mostly not yet extrapolated to the continuum limit, findings.
The fact that the location of phase boundaries, is neither qualitatively (for the upper left corner), nor quantitatively established, motivates us to push forward with more investigations on the subject aiming at clarifying, in particular, the picture for N f = 2 degenerate light flavors. In this contribution updates are provided on our attempts to pursue the continuum limit location of the Z 2 critical endpoint in the heavy mass region κ Z 2 heavy at µ = 0 (Sec. 3); use an extrapolation with tricritical exponents from non integer N f for light quarks (Sec. 4) and locate the tricritical heavy and light endpoints m tricr. heavy and m tricr. light at imaginary chemical potential (Sec. 5). Sec. 2 is instead devoted to a description of the strategy adopted in the various projects to locate the relevant phase transitions and establish their order.  In our studies, the chemical potential has been kept fixed at either µ = 0 or at µ i = µ RW i , with the temperature T related to the coupling β according to T = 1/(a(β)N τ ). Studies in Sec. 5 and in Sec. 4 were conducted at a fixed temporal extent of the lattices (no continuum limit is attempted in these cases), while for the study described in Sec. 3 three different values of N τ were considered. The ranges in mass m or hopping parameter κ and gauge coupling constant β were always dictated by our purpose of locating the chiral/deconfinement phase transition with the purpose of mapping down the position and nature of critical boundaries in the phase diagram.
To locate the chiral/deconfinement phase transition and to identify its order, a "common" strategy was adopted, which consisted in a finite size scaling analysis (FSS) of the third and/or fourth standardized moments of the distribution of the (approximate) order parameter. The n th standardized moment, given the distribution of a generic observable O, is expressed as and we will analyze its dependence on some parameter X ∈ {m, κ, β} and on the volume. We will introduce the (approximate) order parameter O for each investigation in the corresponding sections. However, in all cases, in order to extract the order of the transition as a function of the bare quark mass and/or number of flavors, we considered the kurtosis B 4 (O) of the sampled distribution of O.
For the studies described in Sec. 3 and in Sec. 4 it was enough to consider the kurtosis evaluated on the phase boundary i.e. in correspondence to the value of the coupling β c at which the distribution, at various volumes N σ , showed a vanishing skewness B 3 (β = β c ) = 0. For the study described in Sec. 5, the analysis is slightly different due to the more complex phase structure of QCD at imaginary chemical potential, that was first studied by Roberge and Weiss. The order parameter for the Roberge-Weiss phase transition shall be used in this case. Every coupling is, then, critical and the thermal phase transition is located by the position in β of the crossing of the kurtosis datasets at various N σ values.
In the thermodynamic limit N σ → ∞, the universal values taken by the kurtosis B 4 and by the critical exponent ν are well known results. However, the discontinuous step function characterizing the thermodynamic limit is smeared out to a smooth function as soon as a finite volume is considered and a FSS is needed. In all cases we varied the spatial extent of the lattice N σ such that the aspect ratios, governing the size of the box in physical units at finite temperature, was in the range N σ /N τ ∈ [2 − 5]. In the vicinity of a critical point, the kurtosis can be expanded in powers of the scaling variable x = (X − X c )N 1/ν σ and, for large enough volumes, the expansion can be truncated after the linear term, In our case, for the studies described in Sec. 3 and in Sec. 4, the critical value for X c = κ, m corresponds to a second order phase transition in the 3D Ising universality class, so that one can fix B 4 ≈ 1.604 and ν ≈ 0.63 and perform the fit to Eq. (2) with the sole aim of extracting X c (and c).

The quantitative data collapse as alternative to the kurtosis fit
The technique of fitting the kurtosis around X c has been employed in many studies in order to locate a phase boundary and/or establish the order of a phase transition. However, especially when it comes to fitting reweighted, rather than just raw, data and due to the necessity of fixing different fit ranges for different volumes, the fit procedure ends up relying on some more or less arbitrary decisions hence being not so solid. Based on these observation an alternative procedure was devised and implemented, which measures in a quantitative and solid way the quality of the collapse among sets of data obtained on different lattice sizes once they are plotted against the scaling variable x, with the critical parameter X c and the exponent ν suitably fixed. More details on this method can be found in [4]. The collapse quality, fixed some criticalX c value and some valueν for the critical exponent and considering all pairs of volumes, is where N V is the number of simulated volumes, the integration is done numerically (after having reweighted with a high enough resolution to safely interpolate) and ∆x = x max − x min is the symmetric interval around x = 0 over which the integration takes place. Once the interpolation and the integration are made, Q(X c , ν) is minimized as a function of its two variables. Since a procedure to get the statistical errors on X c and ν needs to be also set up, what one can do is to use N boot different sets of reweighted kurtosis to minimize Q(X c ,ν), so that as many different estimates of X c and ν are obtained and the corresponding statistical error can be determined. Clearly, a role is also played by the choice of ∆x. However, while using a too wide ∆x can be wrong due to the critical region being smaller, one can extrapolate X c and ν to ∆x → 0. Contrary to the kurtosis fit procedure, in the quantitative data collapse reweighting with a higher resolution can only help getting a better estimate of Q(X c ,ν) 1 and the arbitrariness in the choice of ∆x is removed by the extrapolation. For the above reasons the quantitative data collapse should be used, if not to replace, at least to cross check results from the kurtosis fit.
(a) Shift of the critical m Z 2 masses at N f = 2 towards the continuum limit

Updates on the Columbia plot: Z 2 boundary in in the high masses corner
The entity of the cut-off effects that quantitatively affect our picture of the QCD phase structure have been investigated in previous studies in the upper right corner of the Columbia plot and the Z 2 transitions were already observed to shift to smaller masses for N f = 2, 2 + 1, 3 at µ = 0 [5][6][7][8]. A sketch of this behavior for N f = 2 is given in Figure 2(a). No continuum limit result is known yet. In this case the norm of the Polyakov loop ||L|| was used as approximate order parameter for the deconfinement phase transition. The case of N f = 2 degenerate quarks was addressed and, with a scan in κ, the location of the critical κ Z 2 heavy endpoint on N τ = 6, 8, 10 was investigated with the aim to monitor and possibly model cut-off effects. A previous report on this project is to be found in [9], while here we would like to just focus on the impact of finite size effects in this investigation. It is surely a feature of the heavy mass region that pions, with the adopted discretization, cannot be resolved on our lattices up to N τ ≈ 10. And, while a → 0 with growing N τ , the necessity of keeping the relation 1 N τ N σ satisfied forces us to use larger N σ values to keep the size of the box fixed. This has to be true already for the smallest of the volumes in our FSS analysis.
In view of the increasing cost of simulations some work was invested in devising an alternative and possibly cheaper strategy to locate κ Z 2 heavy . One could, indeed, first identify at any fixed value of the bare mass some minimal physical volume V min characterized by allowing a reliable extraction of κ Z 2 heavy out of a linear fit of the kurtosis. At a different κ or N τ , it should be then enough to e.g. reweight the effective potential V eff at just one fixed V V min as in [7] to locate the phase transition and understand its nature. In practice we start by using a modified fit ansatz for the kurtosis in the vicinity of the which incorporates the finite volume effect for generic observables which are a mixture of energy-like and magnetization-like operator and where the value of the exponent y t − y h is fixed by universality. Then we estimate V min by excluding one-by-one the smallest physical volumes in the fit, until the value of the coefficient of the correction term is compatible with zero. Preliminary results are collected in Table 1 and one example of the performed fits is provided in Figure 2(b).

Updates on an alternative Columbia plot: Z 2 boundary at non-integer N f
In this study, we consider QCD with N f mass-degenerate quarks of mass m at zero density and the partition function reads We formally view this as a partition function of some statistical system characterized by a continuous parameter N f and we try to find out for which (tricritical) value of N f the phase transition displayed by this system changes from first-order to second-order. Of course, the extension of Z N f (m) to noninteger values of N f is not unique, and since we use an interpolation characterized by non-integer powers of the determinant, it does not correspond to any local quantum field theory. However, we are just considering, for any given lattice spacing, a statistical system that represents one particular interpolation between two quantum field theories with integer N f . It is not the specific value of N tric f being relevant, but just its relative location with respect to the integer N f = 2 (see Figure 3). Such result should not depend on the chosen interpolation.
In terms of discretization, we employ staggered fermions, where the first-order region is narrow already on coarse lattices. The RHMC algorithm [11] can be used to simulate any number N f of degenerate flavors of staggered fermions, with N f 4 being the power to which the fermion determinant is raised in the lattice partition function.  We measured the chiral condensate ψ ψ as approximate order parameter at 5 different N f values (N f = 2.8, 2.6, 2.4, 2.2, 2.1). For each parameter set {N f , m, N σ , β}, statistics of about (200k − 400k) trajectories were accumulated over 4 Markov chains, subject to the requirement that the skewness of the chiral condensate distribution is compatible within ∼ 2 − 3 standard deviations among the different chains. The phase transition was located according to the kurtosis fit analysis described in Sec. 2. We observed how, as N f is lowered, m Z 2 decreases towards zero and we used that the critical line, in the vicinity of a tricritical point, is known to display a power law dependence with known critical exponents [13]. The scaling law is The plot of the rescaled mass m 2/5 Z 2 as a function of N f is displayed in Figure 4. The solid line in the mass-rescaled plot shows that the simulated results for N f = 2.2, 2.1 are aligned with the result for N f = 2.0, which we were also able to roughly cross-check via direct simulations, of the extrapolation from imaginary chemical potential from [12]. This preliminarily indicates agreement with the expected scaling relation for N f ≤ 2.2.

Updates on the extended Columbia plot: Roberge-Weiss endpoint
The Columbia plot at µ i = µ RW i displayed in Figure 5(a) looks similar to the one in Figure 1(a), but with the Z 2 lines replaced by tricritical lines, first order triple regions that are wider than at µ = 0 and a second order Z 2 region at intermediate values for the quark masses. In this case the imaginary part of the Polyakov loop L Im was measured as order parameter for the Roberge-Weiss phase transition. Once again, we focused on the case of N f = 2 degenerate unimproved staggered quarks, and tried to locate, with a scan in mass, the tricritical points m tricr. heavy and m tricr. light on N τ = 6 lattices as already done for other discretizations and N τ values [14][15][16]. A previous report on this project is to be found in [17].
For each value of m u,d , simulations were performed at a fixed temporal lattice extent N τ = 6 and at a fixed value of the chemical potential aµ RW i = π/6. The extraction of the critical exponent ν was accomplished both with the kurtosis fit procedure and with the quantitative data collapse described in Sec. 2.1. Results for the critical exponent ν are reported in Figure 5(b). Since results from either kind of analysis happen to agree within a 1σ discrepancy in all (but one) case, they are combined to obtain the final answer on ν. To comment more on our results, it is important to stress that in the FSS, at least three simulated volumes (the largest) should be used. A third volume is essential to cross-check the position of the crossing of the kurtosis of different data sets in correspondence to the critical coupling β c . For each pair of points in    Figure 6. Collection of results on the Z 2 critical line/points in the m π − µ 2 plane. Results from [12,14,15,18,19] are included. Various discretizations and/or results at different N τ can be compared.
In order to decide when to stop accumulating statistics, for large (small) masses the kurtosis of the imaginary part of Polyakov loop was required to be compatible on all the chains within 2 (3) standard deviations. Since this condition can be fulfilled also at a very poor statistics, due to large errors, a further empirical requirement is that values of the kurtosis from different chains must span, errors included, an interval not wider than 0.5. As indicated in Figure 5(b) simulations for a few masses are still running and, at the moment, it is still not possible to give an indication on the position of the two tricritical masses with their statistical error, due to none of the simulated masses falling on the first order triple line. What can be observed, e.g. by the shift in the crossing of subsequent pairs of kurtosis datasets, is that in the heavy (light) mass region larger and larger volumes are needed while the bare mass is increased (decreased) in ranges where the transition becomes tricritical to weak first order. For this reason the preliminary estimate for the location of the two tricritical endpoints in terms of pion masses is unchanged with respect to those indicated in [17].
Here, we can, however, quote a preliminary indication, without errors, in terms of pion masses as m tricr.
π heavy = 2.809 GeV and m tricr. πlight = 350 MeV [4]. One can use this preliminary results to compare with results from other discretizations and/or at other N τ values. The comparison at small masses can be visualized by adding one more point to Fig. 6 in [19], as we do in Figure 6. It is possible, to e.g. consider the Wilson versus staggered discretizations and compare the shift of m tricr.
πlight towards smaller masses going from N τ = 4 to N τ = 6: this is found to amount to 35%(14%) of the value for Wilson (staggered). At large masses the value found for m tricr.
π heavy is, instead, unfortunately still affected by large cut-off effects (am π being still larger than 1).