Curvature of the pseudocritical line in QCD: Taylor expansion matches analytic continuation

We determine the curvature of the pseudo-critical line of $N_f = 2 + 1$ QCD with physical quark masses via Taylor expansion in the quark chemical potentials. We adopt a discretization based on stout improved staggered fermions and the tree level Symanzik gauge action; the location of the pseudocritical temperature is based on chiral symmetry restoration. Simulations are performed on lattices with different temporal extent ($N_t = 6, 8, 10$), leading to a continuum extrapolated curvature $\kappa = 0.0145(25)$, which is in very good agreement with the continuum extrapolation obtained via analytic continuation and the same discretization, $\kappa = 0.0135(20)$. This result eliminates the possible tension emerging when comparing analytic continuation with earlier results obtained via Taylor expansion.


I. INTRODUCTION
The exploration of the QCD phase diagram is a subject of continuous experimental and theoretical investigation. One of the main issues is represented by the determination of the pseudo-critical line in the T − µ B plane separating the low-T phase, characterized by color confinement and chiral symmetry breaking, from the high-T phase where the so-called Quark-Gluon Plasma (QGP) is thought to be realized. Lattice QCD simulations at non-zero µ B are still hindered by the well known sign problem, however various methods are already effective to circumvent the problem at least for small µ B , where the pseudo-critical line can be approximated at the lowest order in a Taylor expansion in µ 2 The curvature κ of the pseudo-critical line has been determined on the lattice both by analytic continuation [1][2][3][4][5][6][7][8][9][10][11][12], exploiting results obtained at imaginary chemical potentials, and by Taylor expansion [13][14][15][16], i.e. by suitable combinations of expectation values determined at zero chemical potential. The pseudo-critical line for small values of µ B has been investigated also by continuum approaches to the QCD phase diagram (see, e.g., Refs. [17][18][19][20][21] Recently, various lattice investigations have led to a determination of κ by analytic continuation for QCD with physical or almost physical quark masses [22][23][24][25][26]. In particular, Refs. [24] and [25] have provided continuum extrapolated values for κ which are respectively κ = 0.0135 (20) and κ = 0.0149 (21). The two studies adopted a similar discretization (stout improved staggered fermions and the tree level Symanzik gauge action) and slightly different setups for the quark chemical potentials: in Ref. [25], the strangeness neutrality condition reproduced in the heavy-ion experimental environment was enforced explicitly by tuning the strange quark chemical potential µ s appropriately; in Ref. [24], instead, µ s was set to zero while checking at the same time that its influence on κ is negligible (see also Ref. [27]). Results obtained by analytic continuation but adopting a different lattice discretization (HISQ staggered fermions) have led to similar results [22,26].
Such results are typically larger than earlier results obtained via Taylor expansion [28,29], reporting κ ∼ 0.006. In particular, Ref. [29] reported a continuum extrapolated value κ = 0.0066 (20) adopting the same discretization and the same observables (chiral condensate) as in Refs. [23,24], i.e. a value which is more than two standard deviations away from the result from analytic continuation. As discussed in Ref. [23], only a small part of this discrepancy can be accounted for by the different prescriptions used to determine the dependence of T c on µ B , so that a tension remains.
The agreement between results obtained by the two methods is a necessary requirement in order to state that one has a full control over all systematics involved in analytic continuation and in Taylor expansion. Therefore, the importance of clarifying any possible tension cannot be overestimated; indeed, efforts in this direction are already in progress, for instance by adopting the HISQ staggered discretization [27].
In this study, we present a new continuum extrapolation for the curvature obtained via Taylor expansion, considering the same stout staggered discretization adopted in Refs. [29] and [24]. In particular, we consider different prescriptions to determine κ via Taylor expansion and the analysis of the renormalized condensate: for fixed N t = 6, we show that they provide consistent results; then, exploiting simulations on different values of N t (N t = 6, 8, 10), we are able to provide results extrapolated to the continuum limit which is in full agreement with that obtained by analytic continuation.
The paper is organized as follows. In Section II we provide the details regarding the lattice discretization adopted in this study, the various prescriptions to determine κ that we have explored and the observables involved in such prescription. In Section III we illustrate our numerical results and finally, in Section IV we discuss our conclusions.

II. NUMERICAL METHODS
As in our previous studies, Refs. [23,24], we have considered a rooted stout staggered discretization of the N f = 2 + 1 QCD partition function: where S YM is the tree level Symanzik improved gauge action [30,31], written in terms of the original link variables through traces of n × m rectangular loops, W n×m i; µν , while the fermion matrix (M f st ) i, j is built up in terms of the two times stout-smeared [32] links U (2) i; ν , with an isotropic smearing parameter ρ = 0.15; finally, the rooting procedure is used to remove the residual fourth degeneracy of staggered fermions (see Ref. [33] for a discussion of possible related systematics.). Note that the quark chemical potentials are treated as external sources, and are set to zero in the simulations.
The quark mass spectrum has been chosen so as to have two degenerate light quarks, m u = m d ≡ m l . Standard thermal boundary conditions in the temporal direction have been set for bosonic and fermionic degrees of freedom. The temperature of the system, T = 1/(N t a), has been changed, for fixed N t , by changing the lattice spacing a while staying on a line of constant physics [34,35], corresponding to a pseudo-Goldstone pion mass m π 135 MeV and a strange-to-light mass ratio m s /m l = 28.15.

A. Physical observables used to locate Tc and their renormalization
As in Refs. [23,24], the determination of the pseudocritical temperature T c will be based on chiral symmetry restoration, which is the leading phenomenon in the presence of light quark masses. In particular, we will consider the light quark condensate where and V is the spatial volume. The light quark condensate is affected by additive and multiplicative renormalizations and, as in Refs. [23,24], we consider two different renormalization prescriptions. The first one is and has been introduced in Ref. [36]: the leading additive renormalization, which is linear in the quark mass, cancels in the difference with the strange condensate, while the multiplicative renormalization, being independent of T, drops out by normalizing with respect to quantities measured at T = 0 and at the same UV cutoff. The second definition, introduced in Ref. [29], exploits T = 0 quantities to perform the additive renormalization and the value of the bare quark mass to take care of multiplicative ones: The location of T c is usually defined, in terms of the renormalized light condensate, as the point of maximum slope, i.e. the point where ψ ψ r has an inflection point as a function of T and the absolute value of ∂ ψ ψ r /∂T reaches a maximum. Alternatively, one can look at the peak of the chiral susceptibility, i.e. the maximum of χψ ψ ≡ ∂ ψ ψ r /∂m l . Studies exploiting analytic continuation have considered both definitions and then monitored the behavior of T c as a function of the imaginary baryon chemical potential in order to determine κ. In our case, the determination of κ will be based on the matching of derivatives with respect to T and µ B computed at µ B = 0; the main error source will be the statistical one, which is larger and larger as one considers observables representing higher order derivatives. For this reason, we will limit ourselves to the analysis of the renormalized chiral condensate, which is the lowest derivative, while starting from a second order derivative like the chiral susceptibility would be much more difficult.

B. Possible definitions of κ via Taylor expansion
The most natural extension to finite µ B of the prescription to locate T c in terms of ψ ψ r is to still look for an inflection point, i.e. a point where ∂ 2 ψ ψ r /∂T 2 = 0. In order to understand how T c will move, at the lowest order in µ B , following this prescription, we need to consider a Taylor expansion of ψ ψ r (T ) where The prescription is then to require where the quantities A , B , A and B represent second and third order derivatives of A(T ) and B(T ) with respect to T , t ≡ T − T c and we have performed a lowest order Taylor expansion around the pseudo-critical temperature at µ B = 0, T c . Solving Eq. (11) for t one obtains and finally, following the definition of κ in Eq. (1), one obtains From a practical point of view, Eq. (13) means that one needs to evaluate both the renormalized condensate and its µ 2 B -derivative as a function of T around T c , and that must be done with enough precision so that, after a suitable interpolation, one is able to compute numerically their third and second order derivatives with respect to T at T c .
As we shall see, the program above has to face the low statistical accuracy attainable with reasonable statistics, in particular when evaluating the µ 2 B -derivative and especially on the lattices with higher values of N t , which are necessary to take the continuum extrapolation. For this reason, alternative prescriptions for κ have been adopted in the literature. For instance, in Ref. [29], the pseudocritical temperature at finite µ B is defined as the temperature where the renormalized condensate attains the same value as at T c for µ B = 0, i.e.: Then, by definition, the differential d ψ ψ must vanish along the curve T c (µ B ) The advantage of the expression in Eq. (16) with respect to that in Eq. (13) is twofold: one needs to estimate just the first derivative of the renormalized condensate at T c , which is more precise and stable against the choice of the interpolating function than the third one, and one does not need to know the dependence of ∂ ψ ψ r /∂(µ 2 B ) on T , but just its value at T c . However, the prescription is debatable, since there is no strict reason that the condensate should stay constant in value at T c . However, numerical studies at imaginary chemical potential [23] have shown that it gives results for T c which are compatible, within errors, with those obtained by looking at the inflection point.
In the following we will consider both definitions, and refer to them as κ 1 , Eq. (16), and κ 2 , Eq. (13). As we shall see, a detailed comparison between the two definitions will be possible only on N t = 6 lattices, where they will give compatible results, while on lattices with larger N t the statistical errors attained for κ 2 will make it practically useless, so that our present continuum extrapolation will be based on κ 1 alone. Yet, κ 1 is exactly the prescription adopted in Ref. [29], so that a strict comparison will be possible with the results reported there.
Notice that other prescriptions can be found in the literature, which will not be explored in this study. For instance, the determination reported in Ref. [28] (see also Refs. [27,37,38]) assumes a behavior for the pseudocritical temperature which is dictated by the critical scaling around the possible second order point in the O(4) universality class located at m l = 0.
C. Observables needed to determine κ and setup of chemical potentials Apart from the renormalized chiral condensate, which has been already defined above, the other quantity needed for our study is its derivative with respect to µ 2 B . Looking at Eqs. (7) and (8) one realizes that such derivative is trivially obtained combining the derivatives of the finite temperature flavor condensates with respect to µ 2 B , since zero temperature quantities are independent of µ B around µ B = 0. Therefore, we need to compute where the relevant operators entering previous expression are defined as while M f and M f represent first and second derivatives of the fermion matrix, defined in Eq. (4), with respect to µ B , computed at µ B = 0. The way in which such derivatives, M f and M f , are actually taken depends on the quark flavor f and specifies our setup of quark chemical potentials. In particular, as in Refs. [23,24,29], we set µ u = µ d = µ l = µ B /3 and µ s = 0. Therefore, for the strange flavor we have M s = M s = 0, so thatψψ s andψψ s trivially vanish. Instead for f = u, d, considering Eq. (4) and taking into account that ∂/∂µ B = (1/3)∂/∂µ l , we have All traces appearing in Eq. (18) have been computed, as usual, by multiple noisy estimators, paying attention not to combine the same random vectors when estimating product of traces in order to avoid cross-correlations.

III. NUMERICAL RESULTS
We performed numerical simulations on four different lattices, with dimensions 16 3 × 6, 24 3 × 6, 32 3 × 8 and 40 3 × 10, and adopting the bare parameters reported in Table I in order to stay on a line of constant physics. The scale determination is affected by an overall systematic error of the order of 2-3 % [34,35], which however is not relevant to our final results, which are based on the determination of dimensionless ratios of quantities measured at the critical temperature. We have adopted the  [34,35]. The strange-to-light mass ratio has been fixed to ms/m l = 28.15.
A total of ≈ 15-20 K and ≈ 6 K Molecular Dynamics trajectories have been generated, for each value of β, respectively on the two coarsest lattices and on the N t = 8 lattice. On the N t = 10 lattice, the derivative of the chiral condensate has been measured exploiting a dedicated high-statistics simulation performed at the critical temperature and consisting of ≈ 100 K trajectories, while for the values of the chiral condensate around T c we have relied on results obtained in Ref. [24]. Also for the values of the chiral condensate at zero temperature, which are needed to obtain renormalized quantities, we have used results obtained in previous studies [23,24]. Finally, all traces needed in our computations (see Eq. (18)   Eqs. (13) and (16)): In order to estimate the systematic uncertainty related to the choice of the fitting function, we have tried different ansatzes which are summarized in Table II. In particular, for the renormalized condensate we have adopted an arctangent, A = P 1 + P 2 atan(P 3 (T − T c )), an hyperbolic tangent and a cubic polynomial, while for its µ 2 B derivative we have considered a Lorentzian function   Tables III and IV: reported errors include the systematic one related to the choice of the fit range. Results obtained for the pseudocritical temperature T c are instead reported in Table V. In principle, the uncertainty in the determination of T c should contribute to the error given for A , A , B and B computed at T c , however this contribution turns out to be negligible in most cases, apart from the 40 3 × 10 (and in particular for B) where it is marginally appreciable because of the larger uncertainty on T c . For the determinations of T c and A on the N t = 10 lattice we reused the data already obtained in Ref. [24], while B has been obtained from the single dedicated simulation performed at T T c . Results for A and B are reported respectively in Tables VI and VII; the estimate for the error on B(T c ) stemming from the uncertainty on T c has been based on the data available for B as a function of T on the 24 3 × 6 lattice.
As it can be appreciated from Tables III and IV, [24]. Units and conventions are as in Table III. statistical errors. For this reason, in order to obtain our final estimates for κ 1 and κ 2 , which are are based on Eq. (21) and are reported in Table VIII, we have considered the dispersion of values corresponding to all possible combinations of different fitting functions, and added it, when appreciable, to the statistical error. As one can see, present statistics are not enough to reach reliable estimates of B on lattices with N t > 6, where they are affected by errors of the order of 50 %.
A. Discussion of results on Nt = 6 lattices: finite size effects and comparison between κ1 and κ2 Results obtained on N t = 6 for different spatial sizes permit us to make an assessment of the relevance of finite size effects. For a closer comparison, in Fig. 5 we plot together results obtained for T 2 c B on the different volumes. A mild dependence on the spatial volume is visible for some quantities, like for instance the pseudocritical temperature. However results obtained for κ on the 16 3 × 6 and 24 3 × 6 lattices are in perfect agreement within errors, thus confirming the small sensitivity of κ to finite size effects already observed in studies exploiting analytic continuation [23].
Regarding the comparison between κ 1 and κ 2 , one ob-Lattice κ1(ψψr1) κ2(ψψr1) κ1(ψψr2) κ2(ψψr2) 16 3 × 6 0.0122(5) 0.016 (6)    serves a slight tendency for κ 2 to be larger than κ 1 , however this is not significant within errors, which are of the order of 30% for κ 2 ; therefore the two determinations are compatible within our present level of accuracy. Finally, no significant difference is observed between determinations obtained with the two different renormalization prescriptions for the chiral condensate.
B. Continuum extrapolation of κ1 and comparison with analytic continuation A continuum extrapolation is presently possible only for κ 1 , for which results are available for three different lattice spacings, corresponding to N t = 6, 8 and 10, while for κ 2 we can just say that no appreciable variations are observable going from N t = 6 to N t = 8, however errors for N t = 8 are too large to make this statement of any significance. Assuming corrections linear in a 2 , κ 1 (a 2 ) = κ cont 1 + O(a 2 ), and since T = 1 Nta , an extrapolation to the continuum limit for κ 1 can be obtained by a best fit of the function The values of κ 1 obtained on the 24 3 × 6, 32 3 × 8 and 40 3 × 10 lattices are reported in Fig. 6, where are also illustrated the results of the continuum extrapolation, which gives back κ cont 1 (ψψ r1 ) = 0.0147 (22) and κ cont 1 (ψψ r2 ) = 0.0144 (26), with a reduced chi-squared respectively 0.42 and 1.38.
As a final estimate, we quote the average value κ 1 = 0.0145 (25), which is in very good agreement within errors with previous results obtained via analytic continuation, in particular κ = 0.0135 (20) from Ref. [24] where the same discretization and numerical setup of chemical potentials have been adopted. Apart from the continuum limit, a direct comparison with analytic continuation can be made separately for the different lattice sizes and for the different definitions of κ. This is possible in particular exploiting the results reported in Ref. [23], where values of both κ 1 and κ 2 have been reported for 16 3 × 6, 24 3 × 6 and 32 3 × 8 lattices and for both renormalizations of the chiral condensate. An inspection of Table IX, where results from this work and from Ref. [23] are reported together for lattices where both are available, reveals some tension between Taylor expansion and analytic continuation on the coarsest lattices, which tends to disappear on the finer lattices, also because of the larger statistical errors.

IV. DISCUSSION AND CONCLUSIONS
Recently, many numerical investigations have been carried out to determine the curvature κ of the pseudocritical line in the QCD phase diagram departing from the µ B = 0 axis. Estimates obtained by the Taylor expansion technique have been generally lower than those obtained by analytic continuation, however, since the transition is a crossover, care is needed when comparing results obtained by studying different observables.
In this work, the curvature of the pseudocritical line has been studied for N f = 2 + 1 QCD, via Taylor expansion and through numerical simulations performed using the tree-level Symanzik gauge action and the stoutsmeared staggered fermion action. This is the same discretization adopted in Refs. [23,24,29]; moreover, we have adopted the same observables and definitions of κ investigated in those previous studies, in order to make the comparison closer.
In particular, the location of the phase transition has been determined from the inflection point of the chiral condensate and using two renormalization prescriptions,ψψ r1 andψψ r2 , defined respectively in Eqs. (7) and (8). The curvature coefficient has been calculated using two different definitions: the first one, κ 1 , adopted in Ref. [29], assumes that the value of the renormalized condensate stays constant at the critical temperature as the baryon chemical potential is switched on. The second one, κ 2 , which is the same adopted in Refs. [23,24], looks at how the actual inflection point of the condensate moves as a function of µ B : it is preferable because it is does not rely on particular assumptions, however it involves the computation of higher order derivatives, leading to larger numerical uncertainties.
Our main results can be summarized as follows: 1. No statistically significant effect due to the renormalization prescription has been observed and finite size effects have been found to be negligible.  [25]. Bars marked with circles and triangles indicate that estimates have been obtained respectively by Taylor expansion and analytic continuation. More precisely, the red, magenta, blue, cyan and green colors indicate that the estimate has been obtained respectively by Taylor expansion + chiral condensate, Taylor expansion + chiral susceptibility, analytic continuation + chiral condensate, analytic continuation + chiral susceptibility and analytic continuation + combined analysis of various observables.
The values obtained for κ 2 are generally higher than the values obtained for κ 1 , however the difference is well within statistical errors, therefore the two determinations are compatible with our present level of accuracy.
2. The values obtained for κ 2 are compatible with those obtained in Ref. [23] via analytic continuation, however present statistics are not enough to determine a reliable estimate for κ 2 via Taylor expansion on finer lattices and to perform an extrapolation to the continuum. Overall, no discrepancy is observed between the results obtained by Taylor expansion and those obtained by analytic continuation also for κ 1 , though some tension is present on the coarsest lattices, which however tends to disappear on finer lattices.
3. The final continuum extrapolation that we have given for κ 1 is κ cont = 0.0145 (25), which is in agreement with results from analytic continuation [22,23,25]. In Fig. 7 we report a summary of the most recent determinations of κ obtained for QCD at or close to the physical point: the possible tension between analytic continuation and earlier results obtained via Taylor expansion seems to disappear, leaving place to a convergence of the two methods.
Regarding the tension (slightly above 2σ) between our present results and the results reported in Ref. [29], where the same discretization, observables and definition of κ 1 were adopted, a possible explanation could be in the different way adopted to take the continuum limit. In our study, for each N t , we have determined κ 1 at the corresponding pseudo-critical temperature found at µ B = 0, and then we have extrapolated those values to the N t → ∞ limit. In Ref. [29], instead, the definition of κ 1 has been first extended to a wide range of temperatures around T c , still based on monitoring how points (temperatures) where the renormalized condensate assumes a fixed value change as a function of µ B ; then a continuum extrapolation has been performed over the whole range, thereafter taking the value of this extrapolated function at T c .
Two final considerations are in order. First, while a proper extrapolation to the continuum limit for κ 2 is probably out of reach with current computational resources, more statistics would allow to improve the results and make a closer comparison at least on N t = 6 lattices. Second, it must be stressed that the conver-gence towards a common continuum extrapolated value of κ, which is indicated by recent determinations from analytic continuation and Taylor expansion, is still limited to results obtained only from staggered fermion simulations. It would be important, in the future, to have confirmations also from studies adopting different fermion discretizations.